tests/tbeyn2.c

   1 #include <qpms/beyn.h>
   2 #include <stdio.h>
   3 #include <string.h>
   4 #include <stdlib.h>
   5 #include <math.h>
   6
   7 static double randU(double a, double b) {return a + (b-a) * random() * (1. / RAND_MAX); }
   8 // Normal distribution via Box-Muller transform
   9 static double randN(double sigma, double mu) {
  10   double u1 = randU(0,1);
  11   double u2 = randU(0,1);
  12   return mu + sigma*sqrt(-2*log(u1))*cos(2.*M_PI*u2);
  13 }
  14
  15 struct param49 {
  16   double *T0;
  17   double *T1;
  18   double *T2;
  19 };
  20
  21 // Matrix as in Beyn, example 4.9
  22 int M_function(complex double *target, const size_t m, const complex double z, void *params) {
  23   struct param49 *p = params;
  24
  25   for(size_t i = 0; i < m*m; ++i)
  26     target[i] = p->T0[i] + z*p->T1[i] + z*z*p->T2[i];
  27   return 0;
  28 }
  29
  30 int main(int argc, char **argv) {
  31   complex double z0 = 0;
  32   double Rx = .3, Ry = .3;
  33   int L = 30, N = 150, dim = 60;
  34   if (argc > 1) N = atoi(argv[1]);
  35   beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
  36   struct param49 p;
  37   p.T0 = malloc(dim*dim*sizeof(double));
  38   p.T1 = malloc(dim*dim*sizeof(double));
  39   p.T2 = malloc(dim*dim*sizeof(double));
  40   for(size_t i = 0; i < dim*dim; ++i) {
  41     p.T0[i] = randN(1,0);
  42     p.T1[i] = randN(1,0);
  43     p.T2[i] = randN(1,0);
  44   }
  45
  46   beyn_result_t *result =
  47     beyn_solve(dim, L, M_function, NULL /*M_inv_Vhat_function*/, &p /*params*/,
  48       contour, 1e-4, 1, 1e-4);
  49   printf("Found %zd eigenvalues:\n", result->neig);
  50   for (size_t i = 0; i < result->neig; ++i) {
  51     complex double eig = result->eigval[i];
  52     printf("%zd: %g%+gj\n", i, creal(eig), cimag(eig));
  53   }
  54   free(contour);
  55   beyn_result_free(result);
  56   free(p.T0);
  57   free(p.T1);
  58   free(p.T2);
  59   return 0;
  60 }
  61
  62