diff --git a/qpms/beyn.c b/qpms/beyn.c index 492bc32..8d98ff2 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -1,8 +1,5 @@ #define STATIC_ASSERT(COND,MSG) typedef char static_assertion_##MSG[(COND)?1:-1] -#define cg2s(x) gsl_complex_tostd(x) -#define cs2g(x) gsl_complex_fromstd(x) - #include #include #include @@ -618,65 +615,3 @@ beyn_result_t *beyn_solve(const size_t m, const size_t l, return result; } -#if 0 -// Wrapper of pure C array M-matrix function to GSL. - -struct beyn_function_M_carr2gsl_param { - beyn_function_M_t M_function; - beyn_function_M_inv_Vhat_t M_inv_Vhat_function; - void *param; -}; - -static int beyn_function_M_carr2gsl(gsl_matrix_complex *target_M, complex double z, void *params) -{ - struct beyn_function_M_carr2gsl_param *p = params; - // These could rather be asserts. - QPMS_ENSURE(target_M->size2 == 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); -} - -beyn_result_t *beyn_solve(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, size_t rank_sel_min, double res_tol) { - struct beyn_function_M_carr2gsl_param p = {M, M_inv_Vhat, params}; - return beyn_result_from_beyn_result_gsl( - beyn_solve_gsl(m, l, beyn_function_M_carr2gsl, - (p.M_inv_Vhat_function) ? beyn_function_M_inv_Vhat_carr2gsl : NULL, - (void *) &p, contour, rank_tol, rank_sel_min, res_tol) - ); -} - -beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g) { - struct beyn_result_t *result; - QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_t)); - result->gsl = g; - result->neig = result->gsl->neig; - result->vlen = result->gsl->eigvec->size2; - 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; - result->ranktest_SV = result->gsl->ranktest_SV->data; - return result; -} - -void beyn_result_free(beyn_result_t *result) { - if(result) - beyn_result_gsl_free(result->gsl); - free(result); -} -#endif - diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 2a9f8c9..795323b 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -741,25 +741,10 @@ cdef extern from "ewald.h": double eta, cdouble wavenumber, LatticeDimensionality latdim, PGen *pgen_R, bint pgen_generates_shifted_points, cart3_t k, cart3_t particle_shift) - -cdef extern from "gsl/gsl_complex.h": - ctypedef struct gsl_complex: - double dat[2] - -cdef extern from "gsl/gsl_matrix.h": - ctypedef struct gsl_matrix_complex: - pass - ctypedef struct gsl_vector: - pass - ctypedef struct gsl_vector_complex: - pass - cdef extern from "beyn.h": ctypedef struct beyn_contour_t: bint (*inside_test)(beyn_contour_t *, cdouble z) pass - ctypedef struct beyn_result_gsl_t: - pass ctypedef struct beyn_result_t: size_t neig size_t vlen @@ -768,25 +753,17 @@ cdef extern from "beyn.h": double *residuals cdouble *eigvec double *ranktest_SV - beyn_result_gsl_t *gsl ctypedef enum beyn_contour_halfellipse_orientation: BEYN_CONTOUR_HALFELLIPSE_RE_PLUS BEYN_CONTOUR_HALFELLIPSE_IM_PLUS BEYN_CONTOUR_HALFELLIPSE_RE_MINUS BEYN_CONTOUR_HALFELLIPSE_IM_MINUS - ctypedef int (*beyn_function_M_gsl_t)(gsl_matrix_complex *target_M, cdouble z, void *params) - ctypedef int (*beyn_function_M_inv_Vhat_gsl_t)(gsl_matrix_complex *target, const gsl_matrix_complex *Vhat, cdouble z, void *params) ctypedef int (*beyn_function_M_t)(cdouble *target_M, size_t m, cdouble z, void *params) ctypedef int (*beyn_function_M_inv_Vhat_t)(cdouble *target, size_t m, size_t l, const cdouble *Vhat, cdouble z, void *params) - void beyn_result_gsl_free(beyn_result_gsl_t *result) void beyn_result_free(beyn_result_t *result) - beyn_result_gsl_t *beyn_solve_gsl(size_t m, size_t l, beyn_function_M_gsl_t M, - beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat, void *params, const beyn_contour_t *contour, - double rank_tol, size_t rank_min_sel, double res_tol) - beyn_result_t *beyn_solve(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, size_t rank_min_sel, double res_tol) @@ -797,7 +774,3 @@ cdef extern from "beyn.h": beyn_contour_t *beyn_contour_kidney(cdouble centre, double halfax_re, double halfax_im, size_t npoints, double rounding, beyn_contour_halfellipse_orientation ori) - - cdouble gsl_comlpex_tostd(gsl_complex z) - gsl_complex gsl_complex_fromstd(cdouble z) - diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index ccd5fea..63df6ae 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -72,14 +72,6 @@ target_link_libraries( tbeyn3f qpms gsl lapacke ${QPMS_AMOSLIB} m ) target_include_directories( tbeyn3f PRIVATE .. ) target_compile_definitions( tbeyn3f PRIVATE VARIANTF ) -#add_executable( tbeyn_gsl tbeyn_gsl.c ) -#target_link_libraries( tbeyn_gsl qpms gsl lapacke ${QPMS_AMOSLIB} m ) -#target_include_directories( tbeyn_gsl PRIVATE .. ) - -#add_executable( tbeyn_gsl2 tbeyn_gsl2.c ) -#target_link_libraries( tbeyn_gsl2 qpms gsl lapacke ${QPMS_AMOSLIB} m ) -#target_include_directories( tbeyn_gsl2 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_gsl.c b/tests/tbeyn_gsl.c deleted file mode 100644 index b2dd1d8..0000000 --- a/tests/tbeyn_gsl.c +++ /dev/null @@ -1,42 +0,0 @@ -#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 = - beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, - contour, 1e-4, 1, 1e-4); - printf("Found %zd eigenvalues:\n", result->neig); - for (size_t i = 0; i < result->neig; ++i) { - gsl_complex eig = gsl_vector_complex_get(result->eigval, i); - printf("%zd: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); - } - free(contour); - beyn_result_gsl_free(result); - return 0; -} - - diff --git a/tests/tbeyn_gsl2.c b/tests/tbeyn_gsl2.c deleted file mode 100644 index 1da782d..0000000 --- a/tests/tbeyn_gsl2.c +++ /dev/null @@ -1,37 +0,0 @@ -#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_matrix_complex_set_zero(target); - for (int i = 0; i < m; ++i) { - gsl_matrix_complex_set(target, i, i, gsl_complex_rect(i - creal(z), -cimag(z))); - } - - return 0; -} - -int main() { - complex double z0 = 150+2*I; - double Rx = 5, Ry = 5; - int L = 10, N = 600, dim = 200; - beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); - - beyn_result_gsl_t *result = - beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, - contour, 1e-4, 1, 1e-4); - printf("Found %zd eigenvalues:\n", result->neig); - for (size_t i = 0; i < result->neig; ++i) { - gsl_complex eig = gsl_vector_complex_get(result->eigval, i); - printf("%zd: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); - - } - free(contour); - beyn_result_gsl_free(result); - return 0; -} - -