tests/tbeyn_gsl.c

   1 #include <qpms/beyn.h>
   2 #include <stdio.h>
   3
   4
   5 // Matrix as in Beyn, section 4.11
   6 int M_function(gsl_matrix_complex *target, complex double z, void *no_params) {
   7   int m = target->size1;
   8
   9   gsl_complex d = gsl_complex_fromstd(  2*m - 4*z / (6*m) );
  10   gsl_complex od =  gsl_complex_fromstd( -(double)m - z / (6*m) );
  11
  12   gsl_matrix_complex_set_zero(target);
  13   for (int i = 0; i < m; ++i) {
  14     gsl_matrix_complex_set(target, i, i, d);
  15     if(i > 0) gsl_matrix_complex_set(target, i, i-1, od);
  16     if(i < m - 1) gsl_matrix_complex_set(target, i, i+1, od);
  17   }
  18   gsl_matrix_complex_set(target, m-1, m-1, gsl_complex_fromstd(gsl_complex_tostd(d)/2 + z/(z-1)));
  19
  20   return 0;
  21 }
  22
  23 int main() {
  24   complex double z0 = 150+2*I;
  25   double Rx = 148, Ry = 148;
  26   int L = 10, N = 50, dim = 400;
  27   beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
  28
  29   beyn_result_gsl_t *result =
  30     beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
  31       contour, 1e-4, 1, 1e-4);
  32   printf("Found %zd eigenvalues:\n", result->neig);
  33   for (size_t i = 0; i < result->neig; ++i) {
  34     gsl_complex eig = gsl_vector_complex_get(result->eigval, i);
  35     printf("%zd: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig));
  36   }
  37   free(contour);
  38   beyn_result_gsl_free(result);
  39   return 0;
  40 }
  41
  42