diff --git a/qpms/beyn.c b/qpms/beyn.c index ed43dd7..9ae71f7 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -18,7 +18,7 @@ #include #include -#include +#include "beyn.h" STATIC_ASSERT((sizeof(lapack_complex_double) == sizeof(gsl_complex)), lapack_and_gsl_complex_must_be_consistent); @@ -412,7 +412,7 @@ int beyn_solve_gsl(beyn_result_gsl_t **result, size_t m, size_t l, gsl_matrix_complex_set_zero(A0Coarse); gsl_matrix_complex_set_zero(A1Coarse); size_t N = contour->n; - double DeltaTheta = 2.0*M_PI / ((double)N); + //double DeltaTheta = 2.0*M_PI / ((double)N); printf(" Evaluating contour integral (%zd points)...\n",N); const complex double z0 = contour->centre; for(int n=0; nsize2 == target_M->tda, "Target GSL matrix is not a C-contiguous array. This is a bug, please report!"); + QPMS_ENSURE(target_M->size1 == target_M->size2, "Target is not a square matrix. This is a bug, please report!"); + return p->M_function((complex double *) target_M->data, target_M->size1, z, p->param); +} + +static int beyn_function_M_inv_Vhat_carr2gsl(gsl_matrix_complex *target, + const gsl_matrix_complex *Vhat, complex double z, void *params) +{ + QPMS_UNTESTED; + struct beyn_function_M_carr2gsl_param *p = params; + // These could rather be asserts. + QPMS_ENSURE(target->size2 == target->tda, "Target GSL matrix is not a C-contiguous array. This is a bug, please report!"); + QPMS_ENSURE(Vhat->size2 == Vhat->tda, "Source GSL matrix is not a C-contiguous array. This is a bug, please report!"); + // TODO here I could also check whether Vhat and target have compatible sizes. + return p->M_inv_Vhat_function((complex double *) target->data, target->size1, target->size2, + (complex double *) Vhat->data, z, p->param); +} + +int beyn_solve(beyn_result_t **result, size_t m, size_t l, beyn_function_M_t M, beyn_function_M_inv_Vhat_t M_inv_Vhat, + void *params, const beyn_contour_t *contour, double rank_tol, double res_tol) { + struct beyn_function_M_carr2gsl_param p = {M, M_inv_Vhat, params}; + QPMS_CRASHING_MALLOC(*result, sizeof(beyn_result_t)); + int retval = beyn_solve_gsl(&((*result)->gsl), m, l, beyn_function_M_carr2gsl, + (p.M_inv_Vhat_function) ? beyn_function_M_inv_Vhat_carr2gsl : NULL, + (void *) &p, contour, rank_tol, res_tol); + (*result)->neig = (*result)->gsl->neig; + (*result)->eigval = (complex double *) (*result)->gsl->eigval->data; + (*result)->eigval_err = (complex double *) (*result)->gsl->eigval_err->data; + (*result)->residuals = (*result)->gsl->residuals->data; + (*result)->eigvec = (complex double *) (*result)->gsl->eigvec->data; + return retval; +} + +void beyn_result_free(beyn_result_t *result) { + if(result) + beyn_result_gsl_free(result->gsl); + free(result); +} + diff --git a/qpms/beyn.h b/qpms/beyn.h index bfa6c0a..df605a7 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -24,7 +24,7 @@ typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex dou /// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$. /** Pure C array version */ typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, size_t m, size_t l, - const gsl_matrix_complex *Vhat, complex double z, void *params); + const complex double *Vhat, complex double z, void *params); /// Complex plane integration contour structure. typedef struct beyn_contour_t { @@ -55,6 +55,10 @@ typedef struct beyn_result_t { complex double *eigval_err; double *residuals; complex double *eigvec; + + /// Private, we wrap it around beyn_result_gsl_t for now. + beyn_result_gsl_t *gsl; + } beyn_result_t; void beyn_result_free(beyn_result_t *result); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index eb0b17e..5413990 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -19,6 +19,10 @@ add_executable( tbeyn tbeyn.c ) target_link_libraries( tbeyn qpms gsl lapacke amos m ) target_include_directories( tbeyn PRIVATE .. ) +add_executable( tbeyn_gsl tbeyn_gsl.c ) +target_link_libraries( tbeyn_gsl qpms gsl lapacke amos m ) +target_include_directories( tbeyn_gsl PRIVATE .. ) + add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array tbeyn ) add_test( NAME single_vs_array_translation_coeffs COMMAND test_single_translations_vs_calc ) diff --git a/tests/tbeyn.c b/tests/tbeyn.c index bbf3f98..1e0322c 100644 --- a/tests/tbeyn.c +++ b/tests/tbeyn.c @@ -1,21 +1,19 @@ #include #include - +#include // Matrix as in Beyn, section 4.11 -int M_function(gsl_matrix_complex *target, complex double z, void *no_params) { - int m = target->size1; +int M_function(complex double *target, const size_t m, const complex double z, void *no_params) { + complex double d = 2*m - 4*z / (6*m); + complex double od = -((double)m) - z / (6*m); - gsl_complex d = gsl_complex_fromstd( 2*m - 4*z / (6*m) ); - gsl_complex od = gsl_complex_fromstd( -m - z / (6*m) ); - - gsl_matrix_complex_set_zero(target); + memset(target, 0, m*m*sizeof(complex double)); for (int i = 0; i < m; ++i) { - gsl_matrix_complex_set(target, i, i, d); - if(i > 0) gsl_matrix_complex_set(target, i, i-1, od); - if(i < m - 1) gsl_matrix_complex_set(target, i, i+1, od); + target[i*m + i] = d; + if(i > 0) target[i*m + i-1] = od; + if(i < m - 1) target[i*m + i+1] = od; } - gsl_matrix_complex_set(target, m-1, m-1, gsl_complex_fromstd(gsl_complex_tostd(d)/2 + z/(z-1))); + target[m*(m-1) + m-1] = d/2 + z/(z-1); return 0; } @@ -26,16 +24,16 @@ int main() { int L = 10, N = 50, dim = 400; beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); - beyn_result_gsl_t *result; - int K = beyn_solve_gsl(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, + beyn_result_t *result; + int K = beyn_solve(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, contour, 1e-4, 0); printf("Found %d eigenvalues:\n", K); for (int i = 0; i < K; ++i) { - gsl_complex eig = gsl_vector_complex_get(result->eigval, i); - printf("%d: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); + complex double eig = result->eigval[i]; + printf("%d: %g%+gj\n", i, creal(eig), cimag(eig)); } free(contour); - beyn_result_gsl_free(result); + beyn_result_free(result); return 0; } diff --git a/tests/tbeyn_gsl.c b/tests/tbeyn_gsl.c new file mode 100644 index 0000000..8d492d3 --- /dev/null +++ b/tests/tbeyn_gsl.c @@ -0,0 +1,42 @@ +#include +#include + + +// Matrix as in Beyn, section 4.11 +int M_function(gsl_matrix_complex *target, complex double z, void *no_params) { + int m = target->size1; + + gsl_complex d = gsl_complex_fromstd( 2*m - 4*z / (6*m) ); + gsl_complex od = gsl_complex_fromstd( -(double)m - z / (6*m) ); + + gsl_matrix_complex_set_zero(target); + for (int i = 0; i < m; ++i) { + gsl_matrix_complex_set(target, i, i, d); + if(i > 0) gsl_matrix_complex_set(target, i, i-1, od); + if(i < m - 1) gsl_matrix_complex_set(target, i, i+1, od); + } + gsl_matrix_complex_set(target, m-1, m-1, gsl_complex_fromstd(gsl_complex_tostd(d)/2 + z/(z-1))); + + return 0; +} + +int main() { + complex double z0 = 150+2*I; + double Rx = 148, Ry = 148; + int L = 10, N = 50, dim = 400; + beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); + + beyn_result_gsl_t *result; + int K = beyn_solve_gsl(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, + contour, 1e-4, 0); + printf("Found %d eigenvalues:\n", K); + for (int i = 0; i < K; ++i) { + gsl_complex eig = gsl_vector_complex_get(result->eigval, i); + printf("%d: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); + } + free(contour); + beyn_result_gsl_free(result); + return 0; +} + +