qpms/beyn.c

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