tests/tbeyn3.c

   1 #define _GNU_SOURCE
   2 #include <qpms/beyn.h>
   3 #include <stdio.h>
   4 #include <string.h>
   5 #include <stdlib.h>
   6 #include <math.h>
   7 #include <fenv.h>
   8
   9 static double randU(double a, double b) {return a + (b-a) * random() * (1. / RAND_MAX); }
  10 // Normal distribution via Box-Muller transform
  11 static double randN(double sigma, double mu) {
  12   double u1 = randU(0,1);
  13   double u2 = randU(0,1);
  14   return mu + sigma*sqrt(-2*log(u1))*cos(2.*M_PI*u2);
  15 }
  16
  17 struct param {
  18   double *T0;
  19   double *T1;
  20   double *T2;
  21 };
  22
  23 int M_function(complex double *target, const size_t m, const complex double z, void *params) {
  24   struct param *p = params;
  25
  26   for(size_t i = 0; i < m*m; ++i)
  27     target[i] = p->T0[i] + z*p->T1[i] + cexp(
  28 #ifdef VARIANTB // Also note that this case requires pretty large contour point number (>~ 3000)
  29         (1+3*I)
  30 #else // VARIANTA or VARIANTC
  31         (1+1*I)
  32 #endif
  33         *z*p->T2[i]) +
  34 #ifdef VARIANTC // Essential singularity at zero
  35       cexp(3/z);
  36 #elif defined VARIANTD // Essential singularity outside the contour
  37       cexp(3/(z-1))
  38 #elif defined VARIANTE // High-order pole at zero
  39       3/cpow(z,10);
  40 #elif defined VARIANTF // High-order pole at zero, higher order than dim
  41       .0003/cpow(z,12);
  42 #else // double pole at zero
  43       3/z/z
  44 #endif
  45       ;
  46   return 0;
  47 }
  48
  49 int main(int argc, char **argv) {
  50   feenableexcept(FE_INVALID | FE_OVERFLOW);
  51   complex double z0 = 0+3e-1*I;
  52 #ifdef RXSMALL
  53   double Rx = .1;
  54 #else
  55   double Rx = .3; // Variant B will fail in this case due to large number of eigenvalues (>30)
  56 #endif
  57   double Ry = .25;
  58 #ifdef VARIANTF
  59   int L = 10, N = 150, dim = 10;
  60 #else
  61   int L = 30, N = 150, dim = 60;
  62 #endif
  63   if (argc > 1) N = atoi(argv[1]);
  64   if (argc > 2) L = atoi(argv[2]);
  65 #ifdef IMPLUS
  66   beyn_contour_t *contour = beyn_contour_halfellipse(z0, Rx, Ry, N, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS);
  67 #elif defined IMPLUS_KIDNEY
  68   beyn_contour_t *contour = beyn_contour_kidney(z0, Rx, Ry, 0.3, N, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS);
  69 #else
  70   beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
  71 #endif
  72   struct param p;
  73   p.T0 = malloc(dim*dim*sizeof(double));
  74   p.T1 = malloc(dim*dim*sizeof(double));
  75   p.T2 = malloc(dim*dim*sizeof(double));
  76   for(size_t i = 0; i < dim*dim; ++i) {
  77     p.T0[i] = randN(1,0);
  78     p.T1[i] = randN(1,0);
  79     p.T2[i] = randN(1,0);
  80   }
  81
  82   beyn_result_t *result =
  83     beyn_solve(dim, L, M_function, NULL /*M_inv_Vhat_function*/, &p /*params*/,
  84       contour, 1e-4, 1, 1e-4);
  85   printf("Found %zd eigenvalues:\n", result->neig);
  86   for (size_t i = 0; i < result->neig; ++i) {
  87     complex double eig = result->eigval[i];
  88     printf("%zd: %g%+gj\n", i, creal(eig), cimag(eig));
  89   }
  90   free(contour);
  91   beyn_result_free(result);
  92   free(p.T0);
  93   free(p.T1);
  94   free(p.T2);
  95   return 0;
  96 }
  97
  98