Implement an API for Beyn's algorithm with standard C arrays.

Former-commit-id: 645490de63f802c8a41f3bad1845cd29e0c3d823

qpms/beyn.c

@@ -18,7 +18,7 @@
 #include <gsl/gsl_blas.h>
 #include <gsl/gsl_eigen.h>
 
-#include <beyn.h>
+#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; n<N; n++) { 
@@ -482,3 +482,54 @@ int beyn_solve_gsl(beyn_result_gsl_t **result, size_t m, size_t l,
   return K;
 }
 
+// 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);
+}
+
+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);
+}
+

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);

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 )

tests/tbeyn.c

@@ -1,21 +1,19 @@
 #include <qpms/beyn.h>
 #include <stdio.h>
-
+#include <string.h>
 
 // Matrix as in Beyn, section 4.11
-int M_function(gsl_matrix_complex *target, complex double z, void *no_params) {
+int M_function(complex double *target, const size_t m, const complex double z, void *no_params) {
-  int m = target->size1;
+  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) );
+  memset(target, 0, m*m*sizeof(complex double));
-  gsl_complex od =  gsl_complex_fromstd( -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);
+    target[i*m + i] = d;
-    if(i > 0) gsl_matrix_complex_set(target, i, i-1, od);
+    if(i > 0) target[i*m + i-1] = od;
-    if(i < m - 1) gsl_matrix_complex_set(target, i, 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;
+  beyn_result_t *result;
-  int K = beyn_solve_gsl(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
+  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);
+    complex double eig = result->eigval[i];
-    printf("%d: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig));
+    printf("%d: %g%+gj\n", i, creal(eig), cimag(eig));
   }
   free(contour);
-  beyn_result_gsl_free(result);
+  beyn_result_free(result);
   return 0;
 }
 

tests/tbeyn_gsl.c

@@ -0,0 +1,42 @@
+#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;
+  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;
+}
+
+