QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

[WIP;BROKEN] Getting rid of gsl linear algebra in Beyn. 

Segfaults in LAPACKE_zgesdd.


Former-commit-id: 39ce2f70e214bfe2fba6f1ec2960139194dbc761
This commit is contained in:
Marek Nečada 2020-04-12 02:21:56 +03:00
parent 7e0d7639e8
commit 64f77557a7
4 changed files with 255 additions and 214 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

405
qpms/beyn.c
View File

@ -10,30 +10,25 @@
#include <math.h> #include <math.h>
#include <time.h> #include <time.h>
#include "qpms_error.h" #include "qpms_error.h"
#include <string.h>
// Maybe GSL works? #include <cblas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_eigen.h>
#include "beyn.h" #include "beyn.h"
#define SQ(x) ((x)*(x)) #define SQ(x) ((x)*(x))
STATIC_ASSERT((sizeof(lapack_complex_double) == sizeof(gsl_complex)), lapack_and_gsl_complex_must_be_consistent); // matrix access
#define MAT(mat_, n_rows_, n_cols_, i_row_, i_col_) (mat_[(n_cols_) * (i_row_) + (i_col_)])
typedef struct BeynSolver typedef struct BeynSolver
{ {
int M; // dimension of matrices int M; // dimension of matrices
int L; // number of columns of VHat matrix int L; // number of columns of VHat matrix
gsl_vector_complex *eigenvalues, *eigenvalue_errors; complex double *eigenvalues, *eigenvalue_errors; // [L]
gsl_matrix_complex *eigenvectors; complex double *eigenvectors; // [L][M] (!!!)
gsl_matrix_complex *A0, *A1, *A0_coarse, *A1_coarse, *MInvVHat; complex double *A0, *A1, *A0_coarse, *A1_coarse, *MInvVHat; // [M][L]
gsl_matrix_complex *VHat; complex double *VHat; // [M][L]
gsl_vector *Sigma, *residuals; double *Sigma, *residuals; // [L]
} BeynSolver; } BeynSolver;
// constructor, destructor // constructor, destructor
@ -253,39 +248,51 @@ beyn_contour_t *beyn_contour_kidney(complex double centre, double rRe,
return c; return c;
} }
void beyn_result_gsl_free(beyn_result_gsl_t *r) { void beyn_result_free(beyn_result_t *r) {
if(r) { if(r) {
gsl_vector_complex_free(r->eigval); free(r->eigval);
gsl_vector_complex_free(r->eigval_err); free(r->eigval_err);
gsl_vector_free(r->residuals); free(r->residuals);
gsl_matrix_complex_free(r->eigvec); free(r->eigvec);
gsl_vector_free(r->ranktest_SV); free(r->ranktest_SV);
free(r); free(r);
} }
} }
BeynSolver *BeynSolver_create(int M, int L) BeynSolver *BeynSolver_create(int M, int L)
{ {
BeynSolver *solver= (BeynSolver *)malloc(sizeof(*solver)); BeynSolver *solver;
QPMS_CRASHING_CALLOC(solver, 1, sizeof(*solver));
solver->M = M; solver->M = M;
solver->L = L; solver->L = L;
QPMS_ENSURE(L <= M, "We expect L <= M, but we got L = %d, M = %d ", L, M); QPMS_ENSURE(L <= M, "We expect L <= M, but we got L = %d, M = %d ", L, M);
// storage for eigenvalues and eigenvectors // storage for eigenvalues and eigenvectors
solver->eigenvalues = gsl_vector_complex_calloc(L); //solver->eigenvalues = gsl_vector_complex_calloc(L);
solver->eigenvalue_errors = gsl_vector_complex_calloc(L); QPMS_CRASHING_CALLOC(solver->eigenvalues, L, sizeof(complex double));
solver->residuals = gsl_vector_calloc(L); //solver->eigenvalue_errors = gsl_vector_complex_calloc(L);
solver->eigenvectors = gsl_matrix_complex_calloc(L, M); QPMS_CRASHING_CALLOC(solver->eigenvalue_errors, L, sizeof(complex double));
//solver->residuals = gsl_vector_calloc(L);
QPMS_CRASHING_CALLOC(solver->residuals, L, sizeof(double));
//solver->eigenvectors = gsl_matrix_complex_calloc(L, M);
QPMS_CRASHING_CALLOC(solver->eigenvectors, L * M, sizeof(complex double));
// storage for singular values, random VHat matrix, etc. used in algorithm // storage for singular values, random VHat matrix, etc. used in algorithm
solver->A0 = gsl_matrix_complex_calloc(M,L); //solver->A0 = gsl_matrix_complex_calloc(M,L);
solver->A1 = gsl_matrix_complex_calloc(M,L); QPMS_CRASHING_CALLOC(solver->A0, M * L, sizeof(complex double));
solver->A0_coarse = gsl_matrix_complex_calloc(M,L); //solver->A1 = gsl_matrix_complex_calloc(M,L);
solver->A1_coarse = gsl_matrix_complex_calloc(M,L); QPMS_CRASHING_CALLOC(solver->A1, M * L, sizeof(complex double));
solver->MInvVHat = gsl_matrix_complex_calloc(M,L); //solver->A0_coarse = gsl_matrix_complex_calloc(M,L);
solver->VHat = gsl_matrix_complex_calloc(M,L); QPMS_CRASHING_CALLOC(solver->A0_coarse, M * L, sizeof(complex double));
solver->Sigma = gsl_vector_calloc(L); //solver->A1_coarse = gsl_matrix_complex_calloc(M,L);
QPMS_CRASHING_CALLOC(solver->A1_coarse, M * L, sizeof(complex double));
//solver->MInvVHat = gsl_matrix_complex_calloc(M,L);
QPMS_CRASHING_CALLOC(solver->MInvVHat, M * L, sizeof(complex double));
//solver->VHat = gsl_matrix_complex_calloc(M,L);
QPMS_CRASHING_CALLOC(solver->VHat, M * L, sizeof(complex double));
//solver->Sigma = gsl_vector_calloc(L);
QPMS_CRASHING_CALLOC(solver->Sigma, L, sizeof(double));
// Beyn Step 1: Generate random matrix VHat // Beyn Step 1: Generate random matrix VHat
BeynSolver_srandom(solver,(unsigned)time(NULL)); BeynSolver_srandom(solver,(unsigned)time(NULL));
@ -295,30 +302,30 @@ BeynSolver *BeynSolver_create(int M, int L)
void BeynSolver_free(BeynSolver *solver) void BeynSolver_free(BeynSolver *solver)
{ {
gsl_vector_complex_free(solver->eigenvalues); free(solver->eigenvalues);
gsl_vector_complex_free(solver->eigenvalue_errors); free(solver->eigenvalue_errors);
gsl_matrix_complex_free(solver->eigenvectors); free(solver->eigenvectors);
gsl_matrix_complex_free(solver->A0); free(solver->A0);
gsl_matrix_complex_free(solver->A1); free(solver->A1);
gsl_matrix_complex_free(solver->A0_coarse); free(solver->A0_coarse);
gsl_matrix_complex_free(solver->A1_coarse); free(solver->A1_coarse);
gsl_matrix_complex_free(solver->MInvVHat); free(solver->MInvVHat);
gsl_vector_free(solver->Sigma); free(solver->Sigma);
gsl_vector_free(solver->residuals); free(solver->residuals);
gsl_matrix_complex_free(solver->VHat); free(solver->VHat);
free(solver); free(solver);
} }
void BeynSolver_free_partial(BeynSolver *solver) void BeynSolver_free_partial(BeynSolver *solver)
{ {
gsl_matrix_complex_free(solver->A0); free(solver->A0);
gsl_matrix_complex_free(solver->A1); free(solver->A1);
gsl_matrix_complex_free(solver->A0_coarse); free(solver->A0_coarse);
gsl_matrix_complex_free(solver->A1_coarse); free(solver->A1_coarse);
gsl_matrix_complex_free(solver->MInvVHat); free(solver->MInvVHat);
gsl_matrix_complex_free(solver->VHat); free(solver->VHat);
free(solver); free(solver);
} }
@ -328,11 +335,8 @@ void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed)
if (RandSeed==0) if (RandSeed==0)
RandSeed=time(0); RandSeed=time(0);
srandom(RandSeed); srandom(RandSeed);
gsl_matrix_complex *VHat=solver->VHat; for(size_t i = 0; i < solver->M * solver->L; ++i)
for(int nr=0; nr<VHat->size1; nr++) solver->VHat[i] = zrandN(1, 0);
for(int nc=0; nc<VHat->size2; nc++)
gsl_matrix_complex_set(VHat,nr,nc,cs2g(zrandN(1, 0)));
} }
@ -342,44 +346,47 @@ void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed)
* and eigenvectors * and eigenvectors
*/ */
static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_function, static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_t M_function,
void *Params, void *Params,
gsl_matrix_complex *A0, gsl_matrix_complex *A1, double complex z0, complex double *A0, complex double *A1, double complex z0,
gsl_vector_complex *eigenvalues, gsl_matrix_complex *eigenvectors, const double rank_tol, size_t rank_sel_min, const double res_tol) complex double *eigenvalues, complex double *eigenvectors, const double rank_tol, size_t rank_sel_min, const double res_tol)
{ {
const size_t m = solver->M; const size_t m = solver->M;
const size_t l = solver->L; const size_t l = solver->L;
gsl_vector *Sigma = solver->Sigma; double *Sigma = solver->Sigma;
int verbose = 1; // TODO int verbose = 1; // TODO
// Beyn Step 3: Compute SVD: A0 = V0_full * Sigma * W0T_full // Beyn Step 3: Compute SVD: A0 = V0_full * Sigma * W0T_full
if(verbose) printf(" Beyn: computing SVD...\n"); if(verbose) printf(" Beyn: computing SVD...\n");
gsl_matrix_complex* V0_full = gsl_matrix_complex_alloc(m,l); //gsl_matrix_complex* V0_full = gsl_matrix_complex_alloc(m,l);
gsl_matrix_complex_memcpy(V0_full,A0); //gsl_matrix_complex_memcpy(V0_full,A0);
gsl_matrix_complex* W0T_full = gsl_matrix_complex_alloc(l,l); complex double *V0_full;
QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride); QPMS_CRASHING_MALLOCPY(V0_full, A0, m * l * sizeof(complex double));
QPMS_ENSURE(V0_full->size1 >= V0_full->size2, "m = %zd, l = %zd, what the hell?"); //gsl_matrix_complex* W0T_full = gsl_matrix_complex_alloc(l,l);
complex double *W0T_full; QPMS_CRASHING_MALLOC(W0T_full, l * l * sizeof(complex double));
//QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride);
//QPMS_ENSURE(V0_full->size1 >= V0_full->size2, "m = %zd, l = %zd, what the hell?");
QPMS_ENSURE_SUCCESS(LAPACKE_zgesdd(LAPACK_ROW_MAJOR, // A = U*Σ*conjg(V') QPMS_ENSURE_SUCCESS(LAPACKE_zgesdd(LAPACK_ROW_MAJOR, // A = U*Σ*conjg(V')
'O' /*jobz, 'O' overwrites a with U and saves conjg(V') in vt if m >= n, i.e. if M >= L, which holds*/, 'O' /*jobz, 'O' overwrites a with U and saves conjg(V') in vt if m >= n, i.e. if M >= L, which holds*/,
V0_full->size1 /* m, number of rows */, m, // V0_full->size1 /* m, number of rows */,
V0_full->size2 /* n, number of columns */, l, // V0_full->size2 /* n, number of columns */,
(lapack_complex_double *)(V0_full->data) /*a*/, V0_full, //(lapack_complex_double *)(V0_full->data) /*a*/,
V0_full->tda /*lda*/, l, //V0_full->tda /*lda*/,
Sigma->data /*s*/, Sigma, //Sigma->data /*s*/,
NULL /*u; not used*/, NULL /*u; not used*/,
m /*ldu; should not be really used but lapacke requires it for some obscure reason*/, m /*ldu; should not be really used but lapacke requires it for some obscure reason*/,
(lapack_complex_double *)W0T_full->data /*vt*/, W0T_full, //(lapack_complex_double *)W0T_full->data /*vt*/,
W0T_full->tda /*ldvt*/ l //W0T_full->tda /*ldvt*/
)); ));
// Beyn Step 4: Rank test for Sigma // Beyn Step 4: Rank test for Sigma
// compute effective rank K (number of eigenvalue candidates) // compute effective rank K (number of eigenvalue candidates)
int K=0; int K=0;
for (int k=0; k<Sigma->size /* this is l, actually */; k++) { for (int k=0; k < l/* k<Sigma->size*/ /* this is l, actually */; k++) {
if (verbose) printf("Beyn: SV(%d)=%e\n",k,gsl_vector_get(Sigma, k)); if (verbose) printf("Beyn: SV(%d)=%e\n",k, Sigma[k] /*gsl_vector_get(Sigma, k)*/);
if (k < rank_sel_min || gsl_vector_get(Sigma, k) > rank_tol) if (k < rank_sel_min || Sigma[k] /*gsl_vector_get(Sigma, k)*/ > rank_tol)
K++; K++;
} }
if (verbose)printf(" Beyn: %d/%zd relevant singular values\n",K,l); if (verbose)printf(" Beyn: %d/%zd relevant singular values\n",K,l);
@ -390,46 +397,70 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
// Beyn step 5: B = V0' * A1 * W0 * Sigma^-1 // Beyn step 5: B = V0' * A1 * W0 * Sigma^-1
// set V0, W0T = matrices of first K right, left singular vectors // set V0, W0T = matrices of first K right, left singular vectors
gsl_matrix_complex *V0 = gsl_matrix_complex_alloc(m,K); //gsl_matrix_complex *V0 = gsl_matrix_complex_alloc(m,K);
gsl_matrix_complex *W0T= gsl_matrix_complex_alloc(K,l); //gsl_matrix_complex *W0T= gsl_matrix_complex_alloc(K,l);
complex double *V0, *W0T;
QPMS_CRASHING_MALLOC(V0, m * K * sizeof(complex double));
QPMS_CRASHING_MALLOC(W0T, K * l * sizeof(complex double));
// TODO this is stupid, some parts could be handled simply by realloc.
for (int k = 0; k < K; ++k) { for (int k = 0; k < K; ++k) {
gsl_vector_complex_view tmp; //gsl_vector_complex_view tmp;
tmp = gsl_matrix_complex_column(V0_full, k); //tmp = gsl_matrix_complex_column(V0_full, k);
gsl_matrix_complex_set_col(V0, k, &(tmp.vector)); //gsl_matrix_complex_set_col(V0, k, &(tmp.vector));
tmp = gsl_matrix_complex_row(W0T_full, k); for(int i = 0; i < m; ++i)
gsl_matrix_complex_set_row(W0T, k, &(tmp.vector)); MAT(V0, m, K, i, k) = MAT(V0_full, m, l, i, k);
//tmp = gsl_matrix_complex_row(W0T_full, k);
//gsl_matrix_complex_set_row(W0T, k, &(tmp.vector));
for(int j = 0; j < K; ++j)
MAT(W0T, K, l, k, j) = MAT(W0T_full, m, l, k, j);
} }
//gsl_matrix_complex_free(V0_full);
free(V0_full);
//gsl_matrix_complex_free(W0T_full);
free(W0T_full);
gsl_matrix_complex_free(V0_full); //gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,l);
gsl_matrix_complex_free(W0T_full); //gsl_matrix_complex *B = gsl_matrix_complex_calloc(K,K);
complex double *TM2, *B;
gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,l); QPMS_CRASHING_CALLOC(TM2, K * l, sizeof(complex double));
gsl_matrix_complex *B = gsl_matrix_complex_calloc(K,K); QPMS_CRASHING_CALLOC(B, K * K, sizeof(complex double));
if(verbose) printf(" Multiplying V0*A1->TM...\n"); if(verbose) printf(" Multiplying V0*A1->TM...\n");
const gsl_complex one = gsl_complex_rect(1,0); //const gsl_complex one = gsl_complex_rect(1,0);
const gsl_complex zero = gsl_complex_rect(0,0); //const gsl_complex zero = gsl_complex_rect(0,0);
gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, //gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one,
V0, A1, zero, TM2); // V0, A1, zero, TM2);
// dims: V0[m,K], A1[m,l], TM2[K,l]
const complex double one = 1, zero = 0;
cblas_zgemm(CblasRowMajor, CblasConjTrans, CblasNoTrans, K, l, m,
&one, V0, K, A1, l, &zero, TM2, l);
if(verbose) printf(" Multiplying TM*W0T->B...\n"); if(verbose) printf(" Multiplying TM*W0T->B...\n");
gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, //gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one,
TM2, W0T, zero, B); // TM2, W0T, zero, B);
// DIMS: TM2[K,l], W0T[K,l], B[K,K]
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans, K, K, l,
&one, TM2, l, W0T, l, &zero, B, K);
gsl_matrix_complex_free(W0T); //gsl_matrix_complex_free(W0T);
gsl_matrix_complex_free(TM2); //gsl_matrix_complex_free(TM2);
free(W0T);
free(TM2);
if(verbose) printf(" Scaling B <- B*Sigma^{-1}\n"); if(verbose) printf(" Scaling B <- B*Sigma^{-1}\n");
gsl_vector_complex *tmp = gsl_vector_complex_calloc(K); //gsl_vector_complex *tmp = gsl_vector_complex_calloc(K);
for(int i = 0; i < K; i++) { for(int i = 0; i < K; i++) {
gsl_matrix_complex_get_col(tmp, B, i); //gsl_matrix_complex_get_col(tmp, B, i);
gsl_vector_complex_scale(tmp, gsl_complex_rect(1.0/gsl_vector_get(Sigma,i), 0)); //gsl_vector_complex_scale(tmp, gsl_complex_rect(1.0/gsl_vector_get(Sigma,i), 0));
gsl_matrix_complex_set_col(B,i,tmp); //gsl_matrix_complex_set_col(B,i,tmp);
for(int j = 0; j < K; j++)
MAT(B, K, K, j, i) /= Sigma[i];
} }
gsl_vector_complex_free(tmp); //gsl_vector_complex_free(tmp);
//for(int m=0; m<K; m++) // B <- B * Sigma^{-1} //for(int m=0; m<K; m++) // B <- B * Sigma^{-1}
@ -444,95 +475,121 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
*/ */
if(verbose) printf(" Eigensolving (%i,%i)\n",K,K); if(verbose) printf(" Eigensolving (%i,%i)\n",K,K);
gsl_vector_complex *Lambda = gsl_vector_complex_alloc(K); // eigenvalues //gsl_vector_complex *Lambda = gsl_vector_complex_alloc(K); // eigenvalues
gsl_matrix_complex *S = gsl_matrix_complex_alloc(K,K); // eigenvectors //gsl_matrix_complex *S = gsl_matrix_complex_alloc(K,K); // eigenvectors
complex double *Lambda /* eigenvalues */ , *S /* eigenvector */;
QPMS_CRASHING_MALLOC(Lambda, K * sizeof(complex double));
QPMS_CRASHING_MALLOC(S, K * K * sizeof(complex double));
QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride); //QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride);
QPMS_ENSURE(Lambda->stride == 1, "Lambda vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride); //QPMS_ENSURE(Lambda->stride == 1, "Lambda vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride);
// dims: B[K,K], S[K,K], Lambda[K]
QPMS_ENSURE_SUCCESS(LAPACKE_zgeev( QPMS_ENSURE_SUCCESS(LAPACKE_zgeev(
LAPACK_ROW_MAJOR, LAPACK_ROW_MAJOR,
'N' /* jobvl; don't compute left eigenvectors */, 'N' /* jobvl; don't compute left eigenvectors */,
'V' /* jobvr; do compute right eigenvectors */, 'V' /* jobvr; do compute right eigenvectors */,
K /* n */, K /* n */,
(lapack_complex_double *)(B->data) /* a */, B, //(lapack_complex_double *)(B->data) /* a */,
B->tda /* lda */, K, //B->tda /* lda */,
(lapack_complex_double *) Lambda->data /* w */, Lambda, //(lapack_complex_double *) Lambda->data /* w */,
NULL /* vl */, NULL /* vl */,
m /* ldvl, not used by for some reason needed */, m /* ldvl, not used but for some reason needed */,
(lapack_complex_double *)(S->data)/* vr */, S, //(lapack_complex_double *)(S->data)/* vr */,
S->tda/* ldvr */ K //S->tda/* ldvr */
)); ));
gsl_matrix_complex_free(B); //gsl_matrix_complex_free(B);
free(B);
// V0S <- V0*S // V0S <- V0*S
printf("Multiplying V0*S...\n"); printf("Multiplying V0*S...\n");
gsl_matrix_complex *V0S = gsl_matrix_complex_alloc(m, K); //gsl_matrix_complex *V0S = gsl_matrix_complex_alloc(m, K);
QPMS_ENSURE_SUCCESS(gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, //QPMS_ENSURE_SUCCESS(gsl_blas_zgemm(CblasNoTrans, CblasNoTrans,
one, V0, S, zero, V0S)); // one, V0, S, zero, V0S));
complex double *V0S;
QPMS_CRASHING_MALLOC(V0S, m * K * sizeof(complex double));
// dims: V0[m,K], S[K,K], V0S[m,K]
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, K, K,
&one, V0, K, S, K, &zero, V0S, K);
gsl_matrix_complex_free(V0); //gsl_matrix_complex_free(V0);
free(V0);
// FIXME!!! make clear relation between KRetained and K in the results! // FIXME!!! make clear relation between KRetained and K in the results!
// (If they differ, there are possibly some spurious eigenvalues.) // (If they differ, there are possibly some spurious eigenvalues.)
int KRetained = 0; int KRetained = 0;
gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m, m); //gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m, m);
gsl_vector_complex *MVk = gsl_vector_complex_alloc(m); //gsl_vector_complex *MVk = gsl_vector_complex_alloc(m);
complex double *Mmat, *MVk;
QPMS_CRASHING_MALLOC(Mmat, m * m * sizeof(complex double));
QPMS_CRASHING_MALLOC(MVk, m * sizeof(complex double));
for (int k = 0; k < K; ++k) { for (int k = 0; k < K; ++k) {
const gsl_complex zgsl = gsl_complex_add(gsl_complex_rect(creal(z0), cimag(z0)), gsl_vector_complex_get(Lambda, k)); //const gsl_complex zgsl = gsl_complex_add(gsl_complex_rect(creal(z0), cimag(z0)), gsl_vector_complex_get(Lambda, k));
const complex double z = GSL_REAL(zgsl) + I*GSL_IMAG(zgsl); //const complex double z = GSL_REAL(zgsl) + I*GSL_IMAG(zgsl);
gsl_vector_complex_const_view Vk = gsl_matrix_complex_const_column(V0S, k); const complex double z = z0 + Lambda[k];
//gsl_vector_complex_const_view Vk = gsl_matrix_complex_const_column(V0S, k);
double residual = 0; double residual = 0;
if(res_tol > 0) { if(res_tol > 0) {
QPMS_ENSURE_SUCCESS(M_function(Mmat, z, Params)); QPMS_ENSURE_SUCCESS(M_function(Mmat, m, z, Params));
QPMS_ENSURE_SUCCESS(gsl_blas_zgemv(CblasNoTrans, one, Mmat, &(Vk.vector), zero, MVk)); //QPMS_ENSURE_SUCCESS(gsl_blas_zgemv(CblasNoTrans, one, Mmat, &(Vk.vector), zero, MVk));
residual = gsl_blas_dznrm2(MVk); // Vk[i] == V0S[i, k]; dims: Mmat[m,m], Vk[m] (V0S[m, K]), MVk[m]
cblas_zgemv(CblasRowMajor, CblasNoTrans, m, m,
&one, Mmat, m, &MAT(V0S, m, K, 0, k), m /* stride of Vk in V0S */,
&zero, MVk, 1);
//residual = gsl_blas_dznrm2(MVk);
residual = cblas_dznrm2(m, MVk, 1);
if (verbose) printf("Beyn: Residual(%i)=%e\n",k,residual); if (verbose) printf("Beyn: Residual(%i)=%e\n",k,residual);
} }
if (res_tol > 0 && residual > res_tol) continue; if (res_tol > 0 && residual > res_tol) continue;
gsl_vector_complex_set(eigenvalues, KRetained, zgsl); //gsl_vector_complex_set(eigenvalues, KRetained, zgsl);
eigenvalues[KRetained] = z;
if(eigenvectors) { if(eigenvectors) {
gsl_matrix_complex_set_row(eigenvectors, KRetained, &(Vk.vector)); //gsl_matrix_complex_set_row(eigenvectors, KRetained, &(Vk.vector));
gsl_vector_set(solver->residuals, KRetained, residual); for(int j = 0; j < m; ++j)
MAT(eigenvectors, l, m, KRetained, j) = MAT(V0S, m, K, j, k);
//gsl_vector_set(solver->residuals, KRetained, residual);
} }
++KRetained; ++KRetained;
} }
gsl_matrix_complex_free(V0S); free(V0S);
gsl_matrix_complex_free(Mmat); free(Mmat);
gsl_vector_complex_free(MVk); free(MVk);
gsl_matrix_complex_free(S); free(S);
gsl_vector_complex_free(Lambda); free(Lambda);
return KRetained; return KRetained;
} }
beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l, beyn_result_t *beyn_solve(const size_t m, const size_t l,
beyn_function_M_gsl_t M_function, beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat_function, beyn_function_M_t M_function, beyn_function_M_inv_Vhat_t M_inv_Vhat_function,
void *params, const beyn_contour_t *contour, void *params, const beyn_contour_t *contour,
double rank_tol, size_t rank_sel_min, double res_tol) double rank_tol, size_t rank_sel_min, double res_tol)
{ {
BeynSolver *solver = BeynSolver_create(m, l); BeynSolver *solver = BeynSolver_create(m, l);
gsl_matrix_complex *A0 = solver->A0; complex double *A0 = solver->A0;
gsl_matrix_complex *A1 = solver->A1; complex double *A1 = solver->A1;
gsl_matrix_complex *A0_coarse = solver->A0_coarse; complex double *A0_coarse = solver->A0_coarse;
gsl_matrix_complex *A1_coarse = solver->A1_coarse; complex double *A1_coarse = solver->A1_coarse;
gsl_matrix_complex *MInvVHat = solver->MInvVHat; complex double *MInvVHat = solver->MInvVHat;
gsl_matrix_complex *VHat = solver->VHat; complex double *VHat = solver->VHat;
/***************************************************************/ /***************************************************************/
/* evaluate contour integrals by numerical quadrature to get */ /* evaluate contour integrals by numerical quadrature to get */
/* A0 and A1 matrices */ /* A0 and A1 matrices */
/***************************************************************/ /***************************************************************/
gsl_matrix_complex_set_zero(A0); // TODO zeroing probably redundant (Used calloc...)
gsl_matrix_complex_set_zero(A1); memset(A0, 0, m * l * sizeof(complex double));
gsl_matrix_complex_set_zero(A0_coarse); memset(A1, 0, m * l * sizeof(complex double));
gsl_matrix_complex_set_zero(A1_coarse); memset(A0_coarse, 0, m * l * sizeof(complex double));
memset(A1_coarse, 0, m * l * sizeof(complex double));
const size_t N = contour->n; const size_t N = contour->n;
if(N & 1) QPMS_WARN("Contour discretisation point number is odd (%zd)," if(N & 1) QPMS_WARN("Contour discretisation point number is odd (%zd),"
" the error estimates might be a bit off.", N); " the error estimates might be a bit off.", N);
@ -544,44 +601,66 @@ beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l,
const complex double z = contour->z_dz[n][0]; const complex double z = contour->z_dz[n][0];
const complex double dz = contour->z_dz[n][1]; const complex double dz = contour->z_dz[n][1];
gsl_matrix_complex_memcpy(MInvVHat, VHat); memcpy(MInvVHat, VHat, m * l * sizeof(complex double));
if(M_inv_Vhat_function) { if(M_inv_Vhat_function) {
QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z, params)); QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, m, l, VHat, z, params));
} else { } else {
lapack_int *pivot; lapack_int *pivot;
gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m,m); //gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m,m);
QPMS_ENSURE_SUCCESS(M_function(Mmat, z, params)); complex double *Mmat;
QPMS_CRASHING_MALLOC(Mmat, m * m * sizeof(complex double));
QPMS_ENSURE_SUCCESS(M_function(Mmat, m, z, params));
QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * m); QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * m);
#if 0
QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR,
m /*m*/, m /*n*/,(lapack_complex_double *) Mmat->data /*a*/, Mmat->tda /*lda*/, pivot /*ipiv*/)); m /*m*/, m /*n*/,(lapack_complex_double *) Mmat->data /*a*/, Mmat->tda /*lda*/, pivot /*ipiv*/));
QPMS_ENSURE(MInvVHat->tda == l, "wut?"); QPMS_ENSURE(MInvVHat->tda == l, "wut?");
QPMS_ENSURE_SUCCESS(LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N' /*trans*/, QPMS_ENSURE_SUCCESS(LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N' /*trans*/,
m /*n*/, l/*nrhs*/, (lapack_complex_double *)(Mmat->data) /*a*/, Mmat->tda /*lda*/, pivot/*ipiv*/, m /*n*/, l/*nrhs*/, (lapack_complex_double *)(Mmat->data) /*a*/, Mmat->tda /*lda*/, pivot/*ipiv*/,
(lapack_complex_double *)(MInvVHat->data) /*b*/, MInvVHat->tda/*ldb*/)); (lapack_complex_double *)(MInvVHat->data) /*b*/, MInvVHat->tda/*ldb*/));
#endif
QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR,
m /*m*/, m /*n*/, Mmat /*a*/, m /*lda*/, pivot /*ipiv*/));
QPMS_ENSURE_SUCCESS(LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N' /*trans*/,
m /*n*/, l/*nrhs*/, Mmat /*a*/, m /*lda*/, pivot/*ipiv*/,
MInvVHat /*b*/, l /*ldb*/));
free(pivot); free(pivot);
gsl_matrix_complex_free(Mmat); free(Mmat);
} }
gsl_matrix_complex_scale(MInvVHat, cs2g(dz)); // gsl_matrix_complex_scale(MInvVHat, cs2g(dz));
gsl_matrix_complex_add(A0, MInvVHat); for(size_t i = 0; i < m * l; ++i)
MInvVHat[i] *= dz;
//gsl_matrix_complex_add(A0, MInvVHat);
for(size_t i = 0; i < m * l; ++i)
A0[i] += MInvVHat[i];
if((n%2)==0) { if((n%2)==0) {
gsl_matrix_complex_add(A0_coarse, MInvVHat); //gsl_matrix_complex_add(A0_coarse, MInvVHat);
gsl_matrix_complex_add(A0_coarse, MInvVHat); //gsl_matrix_complex_add(A0_coarse, MInvVHat);
for(size_t i = 0; i < m * l; ++i)
A0_coarse[i] += 2 * MInvVHat[i];
} }
gsl_matrix_complex_scale(MInvVHat, cs2g(z - z0)); // A_1 scaling as in Beyn's Remark 3.2(b) for numerical stability. // A_1 scaling as in Beyn's Remark 3.2(b) for numerical stability.
gsl_matrix_complex_add(A1, MInvVHat); //gsl_matrix_complex_scale(MInvVHat, cs2g(z - z0));
for(size_t i = 0; i < m * l; ++i)
MInvVHat[i] *= (z - z0);
//gsl_matrix_complex_add(A1, MInvVHat);
for(size_t i = 0; i < m * l; ++i)
A1[i] += MInvVHat[i];
if((n%2)==0) { if((n%2)==0) {
gsl_matrix_complex_add(A1_coarse, MInvVHat); for(size_t i = 0; i < m * l; ++i)
gsl_matrix_complex_add(A1_coarse, MInvVHat); //gsl_matrix_complex_add(A1_coarse, MInvVHat);
//gsl_matrix_complex_add(A1_coarse, MInvVHat);
A1_coarse[i] += 2 * MInvVHat[i];
} }
} }
gsl_vector_complex *eigenvalues = solver->eigenvalues; complex double *eigenvalues = solver->eigenvalues;
gsl_vector_complex *eigenvalue_errors = solver->eigenvalue_errors; complex double *eigenvalue_errors = solver->eigenvalue_errors;
gsl_matrix_complex *eigenvectors = solver->eigenvectors; complex double *eigenvectors = solver->eigenvectors;
// Repeat Steps 3–6 with rougher contour approximation to get an error estimate. // Repeat Steps 3–6 with rougher contour approximation to get an error estimate.
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); 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);
@ -589,10 +668,14 @@ beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l,
// Beyn Steps 3–6 // Beyn Steps 3–6
int K = beyn_process_matrices(solver, M_function, params, A0, A1, z0, eigenvalues, eigenvectors, rank_tol, rank_sel_min, res_tol); int K = beyn_process_matrices(solver, M_function, params, A0, A1, z0, eigenvalues, eigenvectors, rank_tol, rank_sel_min, res_tol);
gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors);
beyn_result_gsl_t *result; //gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors);
QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_gsl_t)); const complex double minusone = -1.;
//TODO maybe change the sizes to correspend to retained eigval count K, not l
cblas_zaxpy(l, &minusone, eigenvalues, 1, eigenvalue_errors, 1);
beyn_result_t *result;
QPMS_CRASHING_MALLOC(result, sizeof(*result));
result->eigval = solver->eigenvalues; result->eigval = solver->eigenvalues;
result->eigval_err = solver->eigenvalue_errors; result->eigval_err = solver->eigenvalue_errors;
result->residuals = solver->residuals; result->residuals = solver->residuals;
@ -605,6 +688,7 @@ beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l,
return result; return result;
} }
#if 0
// Wrapper of pure C array M-matrix function to GSL. // Wrapper of pure C array M-matrix function to GSL.
struct beyn_function_M_carr2gsl_param { struct beyn_function_M_carr2gsl_param {
@ -664,4 +748,5 @@ void beyn_result_free(beyn_result_t *result) {
beyn_result_gsl_free(result->gsl); beyn_result_gsl_free(result->gsl);
free(result); free(result);
} }
#endif

51
qpms/beyn.h
View File

@ -4,20 +4,10 @@
#ifndef BEYN_H #ifndef BEYN_H
#define BEYN_H #define BEYN_H
#include <stddef.h>
#include <complex.h> #include <complex.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex_math.h>
/// User-supplied function that provides the operator M(z) whose "roots" are to be found. /// User-supplied function that provides the (row-major) m × m matrix M(z) whose "roots" are to be found.
/** GSL matrix version */
typedef int (*beyn_function_M_gsl_t)(gsl_matrix_complex *target_M, complex double z, void *params);
/// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$.
/** GSL matrix version */
typedef int (*beyn_function_M_inv_Vhat_gsl_t)(gsl_matrix_complex *target_M_inv_Vhat,
const gsl_matrix_complex *Vhat, complex double z, void *params);
/// User-supplied function that provides the operator M(z) whose "roots" are to be found.
/** Pure C array version */ /** Pure C array version */
typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex double z, void *params); typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex double z, void *params);
@ -66,18 +56,6 @@ beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, dou
size_t n, beyn_contour_halfellipse_orientation or); size_t n, beyn_contour_halfellipse_orientation or);
/// Beyn algorithm result structure (GSL matrix/vector version).
typedef struct beyn_result_gsl_t {
size_t neig; ///< Number of eigenvalues found (a bit redundant?).
gsl_vector_complex *eigval;
gsl_vector_complex *eigval_err;
gsl_vector *residuals;
gsl_matrix_complex *eigvec; // Rows are the eigenvectors
gsl_vector *ranktest_SV;
} beyn_result_gsl_t;
void beyn_result_gsl_free(beyn_result_gsl_t *result);
/// Beyn algorithm result structure (pure C array version). /// Beyn algorithm result structure (pure C array version).
typedef struct beyn_result_t { typedef struct beyn_result_t {
size_t neig; ///< Number of eigenvalues found. size_t neig; ///< Number of eigenvalues found.
@ -88,31 +66,11 @@ typedef struct beyn_result_t {
complex double *eigvec; // Rows are the eigenvectors complex double *eigvec; // Rows are the eigenvectors
double *ranktest_SV; double *ranktest_SV;
/// Private, we wrap it around beyn_result_gsl_t for now.
beyn_result_gsl_t *gsl;
} beyn_result_t; } beyn_result_t;
/// Converts beyn_result_gsl_t from beyn_result_t.
/** After calling this, use beyn_result_free() on the returned pointer;
* do NOT run beyn_result_gsl_free() anymore.
*/
beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g);
void beyn_result_free(beyn_result_t *result); void beyn_result_free(beyn_result_t *result);
beyn_result_gsl_t *beyn_solve_gsl( /// Solve a non-linear eigenproblem using Beyn's algorithm
size_t m, ///< Dimension of the matrix \a M.
size_t l, ///< Number of columns of the random matrix \f$ \hat V \f$ (larger than the expected number of solutions).
beyn_function_M_gsl_t M, ///< Function providing the matrix \f$ M(z) \f$.
beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat, ///< Fuction providing the matrix \f$ M^{-1}(z) \hat V \f$ (optional).
void *params, ///< Parameter pointer passed to M() and M_inv_Vhat().
const beyn_contour_t *contour, ///< Integration contour.
double rank_tol, ///< (default: `1e-4`) TODO DOC.
size_t rank_min_sel, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol.
double res_tol ///< (default: `0.0`) TODO DOC.
);
beyn_result_t *beyn_solve( beyn_result_t *beyn_solve(
size_t m, ///< Dimension of the matrix \a M. size_t m, ///< Dimension of the matrix \a M.
size_t l, ///< Number of columns of the random matrix \f$ \hat V \f$ (larger than the expected number of solutions). size_t l, ///< Number of columns of the random matrix \f$ \hat V \f$ (larger than the expected number of solutions).
@ -125,7 +83,4 @@ beyn_result_t *beyn_solve(
double res_tol ///< (default: `0.0`) TODO DOC. double res_tol ///< (default: `0.0`) TODO DOC.
); );
static inline complex double gsl_complex_tostd(gsl_complex z) { return GSL_REAL(z) + I*GSL_IMAG(z); }
static inline gsl_complex gsl_complex_fromstd(complex double z) { return gsl_complex_rect(creal(z), cimag(z)); }
#endif // BEYN_H #endif // BEYN_H

12
tests/CMakeLists.txt
View File

@ -72,13 +72,13 @@ target_link_libraries( tbeyn3f qpms gsl lapacke ${QPMS_AMOSLIB} m )
target_include_directories( tbeyn3f PRIVATE .. ) target_include_directories( tbeyn3f PRIVATE .. )
target_compile_definitions( tbeyn3f PRIVATE VARIANTF ) target_compile_definitions( tbeyn3f PRIVATE VARIANTF )
add_executable( tbeyn_gsl tbeyn_gsl.c ) #add_executable( tbeyn_gsl tbeyn_gsl.c )
target_link_libraries( tbeyn_gsl qpms gsl lapacke ${QPMS_AMOSLIB} m ) #target_link_libraries( tbeyn_gsl qpms gsl lapacke amos m )
target_include_directories( tbeyn_gsl PRIVATE .. ) #target_include_directories( tbeyn_gsl PRIVATE .. )
add_executable( tbeyn_gsl2 tbeyn_gsl2.c ) #add_executable( tbeyn_gsl2 tbeyn_gsl2.c )
target_link_libraries( tbeyn_gsl2 qpms gsl lapacke ${QPMS_AMOSLIB} m ) #target_link_libraries( tbeyn_gsl2 qpms gsl lapacke amos m )
target_include_directories( tbeyn_gsl2 PRIVATE .. ) #target_include_directories( tbeyn_gsl2 PRIVATE .. )
add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array tbeyn ) add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array tbeyn )

1
tests/tbeyn.c
View File

@ -1,6 +1,7 @@
#include <qpms/beyn.h> #include <qpms/beyn.h>
#include <stdio.h> #include <stdio.h>
#include <string.h> #include <string.h>
#include <stdlib.h>
// Matrix as in Beyn, section 4.11 // Matrix as in Beyn, section 4.11
int M_function(complex double *target, const size_t m, const complex double z, void *no_params) { int M_function(complex double *target, const size_t m, const complex double z, void *no_params) {