qpms/beyn.c

   1 #define STATIC_ASSERT(COND,MSG) typedef char static_assertion_##MSG[(COND)?1:-1]
   2
   3 #define cg2s(x) gsl_complex_tostd(x)
   4 #define cs2g(x) gsl_complex_fromstd(x)
   5
   6 #include <complex.h>
   7 #include <lapacke.h>
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <math.h>
  11 #include <time.h>
  12 #include "qpms_error.h"
  13
  14 // Maybe GSL works?
  15 #include <gsl/gsl_matrix.h>
  16 #include <gsl/gsl_complex_math.h>
  17 #include <gsl/gsl_linalg.h>
  18 #include <gsl/gsl_blas.h>
  19 #include <gsl/gsl_eigen.h>
  20
  21 #include "beyn.h"
  22 #define SQ(x) ((x)*(x))
  23
  24 STATIC_ASSERT((sizeof(lapack_complex_double) == sizeof(gsl_complex)), lapack_and_gsl_complex_must_be_consistent);
  25
  26
  27 typedef struct BeynSolver
  28 {
  29   int M;   // dimension of matrices
  30   int L;   // number of columns of VHat matrix
  31
  32   gsl_vector_complex *eigenvalues, *eigenvalue_errors;
  33   gsl_matrix_complex *eigenvectors;
  34   gsl_matrix_complex *A0, *A1, *A0_coarse, *A1_coarse, *MInvVHat;
  35   gsl_matrix_complex *VHat;
  36   gsl_vector *Sigma, *residuals;
  37 } BeynSolver;
  38
  39 // constructor, destructor
  40 BeynSolver *BeynSolver_create(int M, int L);
  41 void BeynSolver_free(BeynSolver *solver);
  42
  43 // reset the random matrix VHat used in Beyn's algorithm
  44 void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed);
  45
  46 // Uniformly random number from interval [a, b].
  47 static double randU(double a, double b) { return a + (b-a) * random() * (1. / RAND_MAX); }
  48
  49 // Random number from normal distribution (via Box-Muller transform, which is enough for our purposes).
  50 static double randN(double Sigma, double Mu) {
  51   double u1 = randU(0, 1);
  52   double u2 = randU(0, 1);
  53   return Mu + Sigma*sqrt(-2*log(u1))*cos(2.*M_PI*u2);
  54 }
  55
  56 static complex double zrandN(double sigma, double mu){
  57   return randN(sigma, mu) + I*randN(sigma, mu);
  58 }
  59
  60 static inline double dsq(double a) { return a * a; }
  61
  62 static _Bool beyn_contour_ellipse_inside_test(struct beyn_contour_t *c, complex double z) {
  63   double rRe = c->z_dz[c->n][0];
  64   double rIm = c->z_dz[c->n][1];
  65   complex double zn = z - c->centre;
  66   return dsq(creal(zn)/rRe) + dsq(cimag(zn)/rIm) < 1;
  67 }
  68
  69 beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double rIm, size_t n)
  70 {
  71   beyn_contour_t *c;
  72   QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + (n+1)*sizeof(c->z_dz[0]));
  73   c->centre = centre;
  74   c->n = n;
  75   for(size_t i = 0; i < n; ++i) {
  76     double t = i*2*M_PI/n;
  77     double st = sin(t), ct = cos(t);
  78     c->z_dz[i][0] = centre + ct*rRe + I*st*rIm;
  79     c->z_dz[i][1] = (-rRe*st + I*rIm*ct) * (2*M_PI / n);
  80   }
  81   // We hide the half-axes metadata after the discretisation points.
  82   c->z_dz[n][0] = rRe;
  83   c->z_dz[n][1] = rIm;
  84   c->inside_test = beyn_contour_ellipse_inside_test;
  85   return c;
  86 }
  87
  88
  89 // Sets correct sign to zero for a given branch cut orientation
  90 static inline complex double
  91 align_zero(complex double z, beyn_contour_halfellipse_orientation or)
  92 {
  93   // Maybe redundant, TODO check the standard.
  94   const double positive_zero = copysign(0., +1.);
  95   const double negative_zero = copysign(0., -1.);
  96   switch(or) { // ensure correct zero signs; CHECKME!!!
  97     case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS:
  98       if(creal(z) == 0 && signbit(creal(z)))
  99         z = positive_zero + I * cimag(z);
 100       break;
 101     case BEYN_CONTOUR_HALFELLIPSE_RE_MINUS:
 102       if(creal(z) == 0 && !signbit(creal(z)))
 103         z = negative_zero + I * cimag(z);
 104       break;
 105     case BEYN_CONTOUR_HALFELLIPSE_IM_PLUS:
 106       if(cimag(z) == 0 && signbit(cimag(z)))
 107         z = creal(z) + I * positive_zero;
 108       break;
 109     case BEYN_CONTOUR_HALFELLIPSE_IM_MINUS:
 110       if(cimag(z) == 0 && !signbit(cimag(z)))
 111         z = creal(z) + I * negative_zero;
 112       break;
 113     default: QPMS_WTF;
 114   }
 115   return z;
 116 }
 117
 118
 119 beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe,
 120     double rIm, size_t n, beyn_contour_halfellipse_orientation or)
 121 {
 122   beyn_contour_t *c;
 123   QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + (n+1)*sizeof(c->z_dz[0])
 124       + sizeof(beyn_contour_halfellipse_orientation));
 125   c->centre = centre;
 126   c->n = n;
 127   const size_t nline = n/2;
 128   const size_t narc = n - nline;
 129   complex double faktor;
 130   double l = rRe, h = rIm;
 131   switch(or) {
 132     case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS:
 133       faktor = -I;
 134       l = rIm; h = rRe;
 135       break;
 136     case BEYN_CONTOUR_HALFELLIPSE_RE_MINUS:
 137       faktor = I;
 138       l = rIm; h = rRe;
 139       break;
 140     case BEYN_CONTOUR_HALFELLIPSE_IM_PLUS:
 141       faktor = 1;
 142       break;
 143     case BEYN_CONTOUR_HALFELLIPSE_IM_MINUS:
 144       faktor = -1;
 145       break;
 146     default: QPMS_WTF;
 147   }
 148
 149   for(size_t i = 0; i < narc; ++i) {
 150     double t = (i+0.5)*M_PI/narc;
 151     double st = sin(t), ct = cos(t);
 152     c->z_dz[i][0] = centre + faktor*(ct*l + I*st*h);
 153     c->z_dz[i][1] = faktor * (-l*st + I*h*ct) * (M_PI / narc);
 154   }
 155   for(size_t i = 0; i < nline; ++i) {
 156     double t = 0.5 * (1 - (double) nline) + i;
 157     c->z_dz[narc + i][0] = align_zero(centre + faktor * t * 2 * l / nline, or);
 158     c->z_dz[narc + i][1] = faktor * 2 * l / nline;
 159   }
 160   // We hide the half-axes metadata after the discretisation points.
 161   c->z_dz[n][0] = rRe;
 162   c->z_dz[n][1] = rIm;
 163   // ugly...
 164   *((beyn_contour_halfellipse_orientation *) &c->z_dz[n+1][0]) = or;
 165   c->inside_test = NULL; // TODO beyn_contour_halfellipse_inside_test;
 166   return c;
 167 }
 168
 169 beyn_contour_t *beyn_contour_kidney(complex double centre, double rRe,
 170     double rIm, const double rounding, const size_t n, beyn_contour_halfellipse_orientation or)
 171 {
 172   beyn_contour_t *c;
 173   QPMS_ENSURE(rounding >= 0 && rounding < .5, "rounding must lie in the interval [0, 0.5)");
 174   QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + (n+1)*sizeof(c->z_dz[0])
 175       + sizeof(beyn_contour_halfellipse_orientation));
 176   c->centre = centre;
 177   c->n = n;
 178   complex double faktor;
 179   double l = rRe, h = rIm;
 180   switch(or) {
 181     case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS:
 182       faktor = -I;
 183       l = rIm; h = rRe;
 184       break;
 185     case BEYN_CONTOUR_HALFELLIPSE_RE_MINUS:
 186       faktor = I;
 187       l = rIm; h = rRe;
 188       break;
 189     case BEYN_CONTOUR_HALFELLIPSE_IM_PLUS:
 190       faktor = 1;
 191       break;
 192     case BEYN_CONTOUR_HALFELLIPSE_IM_MINUS:
 193       faktor = -1;
 194       break;
 195     default: QPMS_WTF;
 196   }
 197
 198   // Small circle centre coordinates.
 199   const double y = rounding * h; // distance from the cut / straight line
 200   const double x = sqrt(SQ(h - y) - SQ(y));
 201
 202   const double alpha = asin(y/(h-y));
 203   const double ar = l/h; // aspect ratio
 204
 205   // Parameter range (equal to the contour length if ar == 1)
 206   const double tmax = 2 * (x + (M_PI_2 + alpha) * y + (M_PI_2 - alpha) * h);
 207   const double dt = tmax / n;
 208
 209   size_t i = 0;
 210   double t;
 211   // Straight line, first part
 212   double t_lo = 0, t_hi = x;
 213   for(; t = i * dt, t <= t_hi; ++i) {
 214     c->z_dz[i][0] = align_zero(centre + (t - t_lo) * ar * faktor, or);
 215     c->z_dz[i][1] = dt * ar * faktor;
 216   }
 217   // First small arc
 218   t_lo = t_hi; t_hi = t_lo + (M_PI_2 + alpha) * y;
 219   for(; t = i * dt, t < t_hi; ++i) {
 220     double phi = (t - t_lo) / y - M_PI_2;
 221     c->z_dz[i][0] = centre + (ar * (x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor;
 222     c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor;
 223   }
 224   // Big arc
 225   t_lo = t_hi; t_hi = t_lo + (M_PI - 2 * alpha) * h;
 226   for(; t = i * dt, t < t_hi; ++i) {
 227     double phi = (t - t_lo) / h + alpha;
 228     c->z_dz[i][0] = centre + (ar * (h * cos(phi)) + h * sin(phi) * I) * faktor;
 229     c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor;
 230   }
 231   // Second small arc
 232   t_lo = t_hi; t_hi = t_lo + (M_PI_2 + alpha) * y;
 233   for(; t = i * dt, t < t_hi; ++i) {
 234     double phi = (t - t_lo) / y + M_PI - alpha;
 235     c->z_dz[i][0] = centre + (ar * (- x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor;
 236     c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor;
 237   }
 238   // Straight line, second part
 239   t_lo = t_hi; t_hi = tmax;
 240   for(; t = i * dt, i < n; ++i) {
 241     c->z_dz[i][0] = align_zero(centre + (t - t_lo - x) * ar * faktor, or);
 242     c->z_dz[i][1] = dt * ar * faktor;
 243   }
 244
 245 #if 0 // TODO later
 246   // We hide the half-axes metadata after the discretisation points.
 247   c->z_dz[n][0] = rRe;
 248   c->z_dz[n][1] = rIm;
 249   // ugly...
 250   *((beyn_contour_halfellipse_orientation *) &c->z_dz[n+1][0]) = or;
 251 #endif
 252   c->inside_test = NULL; // TODO beyn_contour_halfellipse_inside_test;
 253   return c;
 254 }
 255
 256 void beyn_result_gsl_free(beyn_result_gsl_t *r) {
 257   if(r) {
 258     gsl_vector_complex_free(r->eigval);
 259     gsl_vector_complex_free(r->eigval_err);
 260     gsl_vector_free(r->residuals);
 261     gsl_matrix_complex_free(r->eigvec);
 262     gsl_vector_free(r->ranktest_SV);
 263     free(r);
 264   }
 265 }
 266
 267 BeynSolver *BeynSolver_create(int M, int L)
 268 {
 269   BeynSolver *solver= (BeynSolver *)malloc(sizeof(*solver));
 270
 271   solver->M = M;
 272   solver->L = L;
 273   QPMS_ENSURE(L <= M, "We expect L <= M, but we got L = %d, M = %d ", L, M);
 274
 275   // storage for eigenvalues and eigenvectors
 276   solver->eigenvalues  = gsl_vector_complex_calloc(L);
 277   solver->eigenvalue_errors     = gsl_vector_complex_calloc(L);
 278   solver->residuals    = gsl_vector_calloc(L);
 279   solver->eigenvectors = gsl_matrix_complex_calloc(L, M);
 280
 281   // storage for singular values, random VHat matrix, etc. used in algorithm
 282   solver->A0           = gsl_matrix_complex_calloc(M,L);
 283   solver->A1           = gsl_matrix_complex_calloc(M,L);
 284   solver->A0_coarse     = gsl_matrix_complex_calloc(M,L);
 285   solver->A1_coarse     = gsl_matrix_complex_calloc(M,L);
 286   solver->MInvVHat     = gsl_matrix_complex_calloc(M,L);
 287   solver->VHat         = gsl_matrix_complex_calloc(M,L);
 288   solver->Sigma        = gsl_vector_calloc(L);
 289   // Beyn Step 1: Generate random matrix VHat
 290   BeynSolver_srandom(solver,(unsigned)time(NULL));
 291
 292   return solver;
 293
 294 }
 295
 296 void BeynSolver_free(BeynSolver *solver)
 297 {
 298   gsl_vector_complex_free(solver->eigenvalues);
 299   gsl_vector_complex_free(solver->eigenvalue_errors);
 300   gsl_matrix_complex_free(solver->eigenvectors);
 301
 302   gsl_matrix_complex_free(solver->A0);
 303   gsl_matrix_complex_free(solver->A1);
 304   gsl_matrix_complex_free(solver->A0_coarse);
 305   gsl_matrix_complex_free(solver->A1_coarse);
 306   gsl_matrix_complex_free(solver->MInvVHat);
 307   gsl_vector_free(solver->Sigma);
 308   gsl_vector_free(solver->residuals);
 309   gsl_matrix_complex_free(solver->VHat);
 310
 311   free(solver);
 312 }
 313
 314 void BeynSolver_free_partial(BeynSolver *solver)
 315 {
 316   gsl_matrix_complex_free(solver->A0);
 317   gsl_matrix_complex_free(solver->A1);
 318   gsl_matrix_complex_free(solver->A0_coarse);
 319   gsl_matrix_complex_free(solver->A1_coarse);
 320   gsl_matrix_complex_free(solver->MInvVHat);
 321   gsl_matrix_complex_free(solver->VHat);
 322
 323   free(solver);
 324 }
 325
 326 void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed)
 327 {
 328   if (RandSeed==0)
 329     RandSeed=time(0);
 330   srandom(RandSeed);
 331   gsl_matrix_complex *VHat=solver->VHat;
 332   for(int nr=0; nr<VHat->size1; nr++)
 333     for(int nc=0; nc<VHat->size2; nc++)
 334       gsl_matrix_complex_set(VHat,nr,nc,cs2g(zrandN(1, 0)));
 335
 336 }
 337
 338
 339 /*
 340  * linear-algebra manipulations on the A0 and A1 matrices
 341  * (obtained via numerical quadrature) to extract eigenvalues
 342  * and eigenvectors
 343  */
 344
 345 static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_function,
 346     void *Params,
 347     gsl_matrix_complex *A0, gsl_matrix_complex *A1, double complex z0,
 348     gsl_vector_complex *eigenvalues, gsl_matrix_complex *eigenvectors, const double rank_tol, size_t rank_sel_min, const double res_tol)
 349 {
 350   const size_t m = solver->M;
 351   const size_t l = solver->L;
 352   gsl_vector *Sigma = solver->Sigma;
 353
 354   int verbose = 1; // TODO
 355
 356   // Beyn Step 3: Compute SVD: A0 = V0_full * Sigma * W0T_full
 357   if(verbose) printf(" Beyn: computing SVD...\n");
 358   gsl_matrix_complex* V0_full = gsl_matrix_complex_alloc(m,l);
 359   gsl_matrix_complex_memcpy(V0_full,A0);
 360   gsl_matrix_complex* W0T_full = gsl_matrix_complex_alloc(l,l);
 361   QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride);
 362   QPMS_ENSURE(V0_full->size1 >= V0_full->size2, "m = %zd, l = %zd, what the hell?");
 363   QPMS_ENSURE_SUCCESS(LAPACKE_zgesdd(LAPACK_ROW_MAJOR, // A = U*Σ*conjg(V')
 364         'O' /*jobz, 'O' overwrites a with U and saves conjg(V') in vt if m >= n, i.e. if M >= L, which holds*/,
 365         V0_full->size1 /* m, number of rows */,
 366         V0_full->size2 /* n, number of columns */,
 367         (lapack_complex_double *)(V0_full->data) /*a*/,
 368         V0_full->tda /*lda*/,
 369         Sigma->data /*s*/,
 370         NULL /*u; not used*/,
 371         m /*ldu; should not be really used but lapacke requires it for some obscure reason*/,
 372         (lapack_complex_double *)W0T_full->data /*vt*/,
 373         W0T_full->tda /*ldvt*/
 374         ));
 375
 376
 377   // Beyn Step 4: Rank test for Sigma
 378   // compute effective rank K (number of eigenvalue candidates)
 379   int K=0;
 380   for (int k=0; k<Sigma->size /* this is l, actually */; k++) {
 381     if (verbose) printf("Beyn: SV(%d)=%e\n",k,gsl_vector_get(Sigma, k));
 382     if (k < rank_sel_min || gsl_vector_get(Sigma, k) > rank_tol)
 383       K++;
 384   }
 385   if (verbose)printf(" Beyn: %d/%zd relevant singular values\n",K,l);
 386   if (K==0) {
 387     QPMS_WARN("no singular values found in Beyn eigensolver\n");
 388     return 0;
 389   }
 390
 391   // Beyn step 5:  B = V0' * A1 * W0 * Sigma^-1
 392   // set V0, W0T = matrices of first K right, left singular vectors
 393   gsl_matrix_complex *V0 = gsl_matrix_complex_alloc(m,K);
 394   gsl_matrix_complex *W0T= gsl_matrix_complex_alloc(K,l);
 395
 396   for (int k = 0; k < K; ++k) {
 397     gsl_vector_complex_view tmp;
 398     tmp = gsl_matrix_complex_column(V0_full, k);
 399     gsl_matrix_complex_set_col(V0, k, &(tmp.vector));
 400     tmp = gsl_matrix_complex_row(W0T_full, k);
 401     gsl_matrix_complex_set_row(W0T, k, &(tmp.vector));
 402   }
 403
 404   gsl_matrix_complex_free(V0_full);
 405   gsl_matrix_complex_free(W0T_full);
 406
 407   gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,l);
 408   gsl_matrix_complex *B = gsl_matrix_complex_calloc(K,K);
 409
 410   if(verbose) printf(" Multiplying V0*A1->TM...\n");
 411
 412   const gsl_complex one = gsl_complex_rect(1,0);
 413   const gsl_complex zero = gsl_complex_rect(0,0);
 414   gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one,
 415       V0, A1, zero, TM2);
 416
 417   if(verbose) printf(" Multiplying TM*W0T->B...\n");
 418
 419   gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one,
 420       TM2, W0T, zero, B);
 421
 422   gsl_matrix_complex_free(W0T);
 423   gsl_matrix_complex_free(TM2);
 424
 425   if(verbose) printf(" Scaling B <- B*Sigma^{-1}\n");
 426   gsl_vector_complex *tmp = gsl_vector_complex_calloc(K);
 427   for(int i = 0; i < K; i++) {
 428     gsl_matrix_complex_get_col(tmp, B, i);
 429     gsl_vector_complex_scale(tmp, gsl_complex_rect(1.0/gsl_vector_get(Sigma,i), 0));
 430     gsl_matrix_complex_set_col(B,i,tmp);
 431   }
 432   gsl_vector_complex_free(tmp);
 433
 434   //for(int m=0; m<K; m++)                  //  B <- B * Sigma^{-1}
 435
 436   // Beyn step 6.
 437   // Eigenvalue decomposition B -> S*Lambda*S'
 438   /* According to Beyn's paper (Algorithm 1), one should check conditioning
 439    * of the eigenvalues; if they are ill-conditioned, one should perform
 440    * a procedure involving Schur decomposition and reordering.
 441    *
 442    * Beyn refers to MATLAB routines eig, condeig, schur, ordschur.
 443    * I am not sure about the equivalents in LAPACK, TODO check zgeevx, zgeesx.
 444    */
 445   if(verbose) printf(" Eigensolving (%i,%i)\n",K,K);
 446
 447   gsl_vector_complex *Lambda = gsl_vector_complex_alloc(K); // eigenvalues
 448   gsl_matrix_complex *S = gsl_matrix_complex_alloc(K,K); // eigenvectors
 449
 450   QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride);
 451   QPMS_ENSURE(Lambda->stride == 1, "Lambda vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride);
 452   QPMS_ENSURE_SUCCESS(LAPACKE_zgeev(
 453         LAPACK_ROW_MAJOR,
 454         'N' /* jobvl; don't compute left eigenvectors */,
 455         'V' /* jobvr; do compute right eigenvectors */,
 456         K   /* n */,
 457         (lapack_complex_double *)(B->data) /* a */,
 458         B->tda  /* lda */,
 459         (lapack_complex_double *) Lambda->data /* w */,
 460         NULL /* vl */,
 461         m /* ldvl, not used by for some reason needed */,
 462         (lapack_complex_double *)(S->data)/* vr */,
 463         S->tda/* ldvr */
 464         ));
 465
 466   gsl_matrix_complex_free(B);
 467
 468   // V0S <- V0*S
 469   printf("Multiplying V0*S...\n");
 470   gsl_matrix_complex *V0S = gsl_matrix_complex_alloc(m, K);
 471   QPMS_ENSURE_SUCCESS(gsl_blas_zgemm(CblasNoTrans, CblasNoTrans,
 472         one, V0, S, zero, V0S));
 473
 474   gsl_matrix_complex_free(V0);
 475
 476   // FIXME!!! make clear relation between KRetained and K in the results!
 477   // (If they differ, there are possibly some spurious eigenvalues.)
 478   int KRetained = 0;
 479   gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m, m);
 480   gsl_vector_complex *MVk = gsl_vector_complex_alloc(m);
 481   for (int k = 0; k < K; ++k) {
 482     const gsl_complex zgsl = gsl_complex_add(gsl_complex_rect(creal(z0), cimag(z0)), gsl_vector_complex_get(Lambda, k));
 483     const complex double z = GSL_REAL(zgsl) + I*GSL_IMAG(zgsl);
 484     gsl_vector_complex_const_view Vk = gsl_matrix_complex_const_column(V0S, k);
 485
 486     double residual = 0;
 487     if(res_tol > 0) {
 488       QPMS_ENSURE_SUCCESS(M_function(Mmat, z, Params));
 489       QPMS_ENSURE_SUCCESS(gsl_blas_zgemv(CblasNoTrans, one, Mmat, &(Vk.vector), zero, MVk));
 490       residual = gsl_blas_dznrm2(MVk);
 491       if (verbose) printf("Beyn: Residual(%i)=%e\n",k,residual);
 492     }
 493     if (res_tol > 0 && residual > res_tol) continue;
 494
 495     gsl_vector_complex_set(eigenvalues, KRetained, zgsl);
 496     if(eigenvectors) {
 497       gsl_matrix_complex_set_row(eigenvectors, KRetained, &(Vk.vector));
 498       gsl_vector_set(solver->residuals, KRetained, residual);
 499     }
 500     ++KRetained;
 501   }
 502
 503   gsl_matrix_complex_free(V0S);
 504   gsl_matrix_complex_free(Mmat);
 505   gsl_vector_complex_free(MVk);
 506   gsl_matrix_complex_free(S);
 507   gsl_vector_complex_free(Lambda);
 508
 509   return KRetained;
 510 }
 511
 512
 513 beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l,
 514     beyn_function_M_gsl_t M_function, beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat_function,
 515     void *params, const beyn_contour_t *contour,
 516     double rank_tol, size_t rank_sel_min, double res_tol)
 517 {
 518   BeynSolver *solver = BeynSolver_create(m, l);
 519
 520   gsl_matrix_complex *A0           = solver->A0;
 521   gsl_matrix_complex *A1           = solver->A1;
 522   gsl_matrix_complex *A0_coarse     = solver->A0_coarse;
 523   gsl_matrix_complex *A1_coarse     = solver->A1_coarse;
 524   gsl_matrix_complex *MInvVHat     = solver->MInvVHat;
 525   gsl_matrix_complex *VHat         = solver->VHat;
 526
 527   /***************************************************************/
 528   /* evaluate contour integrals by numerical quadrature to get   */
 529   /* A0 and A1 matrices                                          */
 530   /***************************************************************/
 531
 532   gsl_matrix_complex_set_zero(A0);
 533   gsl_matrix_complex_set_zero(A1);
 534   gsl_matrix_complex_set_zero(A0_coarse);
 535   gsl_matrix_complex_set_zero(A1_coarse);
 536   const size_t N = contour->n;
 537   if(N & 1) QPMS_WARN("Contour discretisation point number is odd (%zd),"
 538      " the error estimates might be a bit off.", N);
 539
 540
 541   // Beyn Step 2: Computa A0, A1
 542   const complex double z0 = contour->centre;
 543   for(int n=0; n<N; n++) {
 544     const complex double z = contour->z_dz[n][0];
 545     const complex double dz = contour->z_dz[n][1];
 546
 547     gsl_matrix_complex_memcpy(MInvVHat, VHat);
 548
 549     if(M_inv_Vhat_function) {
 550       QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z, params));
 551     } else {
 552       lapack_int *pivot;
 553       gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m,m);
 554       QPMS_ENSURE_SUCCESS(M_function(Mmat, z, params));
 555       QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * m);
 556       QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR,
 557             m /*m*/, m /*n*/,(lapack_complex_double *) Mmat->data /*a*/, Mmat->tda /*lda*/, pivot /*ipiv*/));
 558       QPMS_ENSURE(MInvVHat->tda == l, "wut?");
 559       QPMS_ENSURE_SUCCESS(LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N' /*trans*/,
 560             m /*n*/, l/*nrhs*/, (lapack_complex_double *)(Mmat->data) /*a*/, Mmat->tda /*lda*/, pivot/*ipiv*/,
 561             (lapack_complex_double *)(MInvVHat->data) /*b*/, MInvVHat->tda/*ldb*/));
 562
 563       free(pivot);
 564       gsl_matrix_complex_free(Mmat);
 565     }
 566
 567     gsl_matrix_complex_scale(MInvVHat, cs2g(dz));
 568     gsl_matrix_complex_add(A0, MInvVHat);
 569     if((n%2)==0) {
 570       gsl_matrix_complex_add(A0_coarse, MInvVHat);
 571       gsl_matrix_complex_add(A0_coarse, MInvVHat);
 572     }
 573
 574     gsl_matrix_complex_scale(MInvVHat, cs2g(z - z0)); // A_1 scaling as in Beyn's Remark 3.2(b) for numerical stability.
 575     gsl_matrix_complex_add(A1, MInvVHat);
 576     if((n%2)==0) {
 577       gsl_matrix_complex_add(A1_coarse, MInvVHat);
 578       gsl_matrix_complex_add(A1_coarse, MInvVHat);
 579     }
 580   }
 581
 582   gsl_vector_complex *eigenvalues  = solver->eigenvalues;
 583   gsl_vector_complex *eigenvalue_errors     = solver->eigenvalue_errors;
 584   gsl_matrix_complex *eigenvectors = solver->eigenvectors;
 585
 586   // Repeat Steps 3–6 with rougher contour approximation to get an error estimate.
 587   int K_coarse = beyn_process_matrices(solver, M_function, params, A0_coarse, A1_coarse, z0, eigenvalue_errors, /*eigenvectors_coarse*/ NULL, rank_tol, rank_sel_min, res_tol);
 588   // Reid did also fabs on the complex and real parts ^^^.
 589
 590   // Beyn Steps 3–6
 591   int K = beyn_process_matrices(solver, M_function, params, A0, A1, z0, eigenvalues, eigenvectors, rank_tol, rank_sel_min, res_tol);
 592   gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors);
 593
 594   beyn_result_gsl_t *result;
 595   QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_gsl_t));
 596   result->eigval = solver->eigenvalues;
 597   result->eigval_err = solver->eigenvalue_errors;
 598   result->residuals = solver->residuals;
 599   result->eigvec = solver->eigenvectors;
 600   result->ranktest_SV = solver->Sigma;
 601   result->neig = K;
 602
 603   BeynSolver_free_partial(solver);
 604
 605   return result;
 606 }
 607
 608 // Wrapper of pure C array M-matrix function to GSL.
 609
 610 struct beyn_function_M_carr2gsl_param {
 611   beyn_function_M_t M_function;
 612   beyn_function_M_inv_Vhat_t M_inv_Vhat_function;
 613   void *param;
 614 };
 615
 616 static int beyn_function_M_carr2gsl(gsl_matrix_complex *target_M, complex double z, void *params)
 617 {
 618   struct beyn_function_M_carr2gsl_param *p = params;
 619   // These could rather be asserts.
 620   QPMS_ENSURE(target_M->size2 == target_M->tda, "Target GSL matrix is not a C-contiguous array. This is a bug, please report!");
 621   QPMS_ENSURE(target_M->size1 == target_M->size2, "Target is not a square matrix. This is a bug, please report!");
 622   return p->M_function((complex double *) target_M->data, target_M->size1, z, p->param);
 623 }
 624
 625 static int beyn_function_M_inv_Vhat_carr2gsl(gsl_matrix_complex *target,
 626     const gsl_matrix_complex *Vhat, complex double z, void *params)
 627 {
 628   QPMS_UNTESTED;
 629   struct beyn_function_M_carr2gsl_param *p = params;
 630   // These could rather be asserts.
 631   QPMS_ENSURE(target->size2 == target->tda, "Target GSL matrix is not a C-contiguous array. This is a bug, please report!");
 632   QPMS_ENSURE(Vhat->size2 == Vhat->tda, "Source GSL matrix is not a C-contiguous array. This is a bug, please report!");
 633   // TODO here I could also check whether Vhat and target have compatible sizes.
 634   return p->M_inv_Vhat_function((complex double *) target->data, target->size1, target->size2,
 635       (complex double *) Vhat->data, z, p->param);
 636 }
 637
 638 beyn_result_t *beyn_solve(size_t m, size_t l, beyn_function_M_t M, beyn_function_M_inv_Vhat_t M_inv_Vhat,
 639     void *params, const beyn_contour_t *contour, double rank_tol, size_t rank_sel_min, double res_tol) {
 640   struct beyn_function_M_carr2gsl_param p = {M, M_inv_Vhat, params};
 641   return beyn_result_from_beyn_result_gsl(
 642     beyn_solve_gsl(m, l, beyn_function_M_carr2gsl,
 643       (p.M_inv_Vhat_function) ? beyn_function_M_inv_Vhat_carr2gsl : NULL,
 644       (void *) &p, contour, rank_tol, rank_sel_min, res_tol)
 645     );
 646 }
 647
 648 beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g) {
 649   struct beyn_result_t *result;
 650   QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_t));
 651   result->gsl = g;
 652   result->neig = result->gsl->neig;
 653   result->vlen = result->gsl->eigvec->size2;
 654   result->eigval = (complex double *) result->gsl->eigval->data;
 655   result->eigval_err = (complex double *) result->gsl->eigval_err->data;
 656   result->residuals = result->gsl->residuals->data;
 657   result->eigvec = (complex double *) result->gsl->eigvec->data;
 658   result->ranktest_SV = result->gsl->ranktest_SV->data;
 659   return result;
 660 }
 661
 662 void beyn_result_free(beyn_result_t *result) {
 663   if(result)
 664     beyn_result_gsl_free(result->gsl);
 665   free(result);
 666 }
 667