tests/tbeyn.c

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