QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Remove all gsl_matrix_... related code in Beyn. 

Now that pure BLAS implementation is working, this in no longer needed.

Getting rid of gsl_matrix and gsl_cblas eliminates header clashes between
GSL and other CBLAS implementations.
This commit is contained in:
Marek Nečada 2021-06-11 10:48:58 +03:00
parent f19c590d6e
commit 9d556e7b23
5 changed files with 0 additions and 179 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
qpms/beyn.c
View File

@ -1,8 +1,5 @@
#define STATIC_ASSERT(COND,MSG) typedef char static_assertion_##MSG[(COND)?1:-1] #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 <complex.h> #include <complex.h>
#include <lapacke.h> #include <lapacke.h>
#include <stdio.h> #include <stdio.h>
@ -618,65 +615,3 @@ beyn_result_t *beyn_solve(const size_t m, const size_t l,
return result; 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

27
qpms/qpms_cdefs.pxd
View File

@ -741,25 +741,10 @@ cdef extern from "ewald.h":
double eta, cdouble wavenumber, LatticeDimensionality latdim, PGen *pgen_R, double eta, cdouble wavenumber, LatticeDimensionality latdim, PGen *pgen_R,
bint pgen_generates_shifted_points, cart3_t k, cart3_t particle_shift) 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": cdef extern from "beyn.h":
ctypedef struct beyn_contour_t: ctypedef struct beyn_contour_t:
bint (*inside_test)(beyn_contour_t *, cdouble z) bint (*inside_test)(beyn_contour_t *, cdouble z)
pass pass
ctypedef struct beyn_result_gsl_t:
pass
ctypedef struct beyn_result_t: ctypedef struct beyn_result_t:
size_t neig size_t neig
size_t vlen size_t vlen
@ -768,25 +753,17 @@ cdef extern from "beyn.h":
double *residuals double *residuals
cdouble *eigvec cdouble *eigvec
double *ranktest_SV double *ranktest_SV
beyn_result_gsl_t *gsl
ctypedef enum beyn_contour_halfellipse_orientation: ctypedef enum beyn_contour_halfellipse_orientation:
BEYN_CONTOUR_HALFELLIPSE_RE_PLUS BEYN_CONTOUR_HALFELLIPSE_RE_PLUS
BEYN_CONTOUR_HALFELLIPSE_IM_PLUS BEYN_CONTOUR_HALFELLIPSE_IM_PLUS
BEYN_CONTOUR_HALFELLIPSE_RE_MINUS BEYN_CONTOUR_HALFELLIPSE_RE_MINUS
BEYN_CONTOUR_HALFELLIPSE_IM_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_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) 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) 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_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, 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) 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, beyn_contour_t *beyn_contour_kidney(cdouble centre, double halfax_re, double halfax_im, size_t npoints,
double rounding, beyn_contour_halfellipse_orientation ori) double rounding, beyn_contour_halfellipse_orientation ori)
cdouble gsl_comlpex_tostd(gsl_complex z)
gsl_complex gsl_complex_fromstd(cdouble z)

8
tests/CMakeLists.txt
View File

@ -72,14 +72,6 @@ target_link_libraries( tbeyn3f qpms gsl lapacke ${QPMS_AMOSLIB} m )
target_include_directories( tbeyn3f PRIVATE .. ) target_include_directories( tbeyn3f PRIVATE .. )
target_compile_definitions( tbeyn3f PRIVATE VARIANTF ) 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_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 ) add_test( NAME single_vs_array_translation_coeffs COMMAND test_single_translations_vs_calc )

42
tests/tbeyn_gsl.c
View File

@ -1,42 +0,0 @@
#include <qpms/beyn.h>
#include <stdio.h>
// 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;
}

37
tests/tbeyn_gsl2.c
View File

@ -1,37 +0,0 @@
#include <qpms/beyn.h>
#include <stdio.h>
// 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;
}