From 64f77557a7bd68de3d98daef80471213fd42af4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 12 Apr 2020 02:21:56 +0300 Subject: [PATCH 1/4] [WIP;BROKEN] Getting rid of gsl linear algebra in Beyn. Segfaults in LAPACKE_zgesdd. Former-commit-id: 39ce2f70e214bfe2fba6f1ec2960139194dbc761 --- qpms/beyn.c | 405 ++++++++++++++++++++++++++----------------- qpms/beyn.h | 51 +----- tests/CMakeLists.txt | 12 +- tests/tbeyn.c | 1 + 4 files changed, 255 insertions(+), 214 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 36c5605..e69b96d 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -10,30 +10,25 @@ #include #include #include "qpms_error.h" - -// Maybe GSL works? -#include -#include -#include -#include -#include +#include +#include #include "beyn.h" #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 { int M; // dimension of matrices int L; // number of columns of VHat matrix - gsl_vector_complex *eigenvalues, *eigenvalue_errors; - gsl_matrix_complex *eigenvectors; - gsl_matrix_complex *A0, *A1, *A0_coarse, *A1_coarse, *MInvVHat; - gsl_matrix_complex *VHat; - gsl_vector *Sigma, *residuals; + complex double *eigenvalues, *eigenvalue_errors; // [L] + complex double *eigenvectors; // [L][M] (!!!) + complex double *A0, *A1, *A0_coarse, *A1_coarse, *MInvVHat; // [M][L] + complex double *VHat; // [M][L] + double *Sigma, *residuals; // [L] } BeynSolver; // constructor, destructor @@ -253,39 +248,51 @@ beyn_contour_t *beyn_contour_kidney(complex double centre, double rRe, return c; } -void beyn_result_gsl_free(beyn_result_gsl_t *r) { +void beyn_result_free(beyn_result_t *r) { if(r) { - gsl_vector_complex_free(r->eigval); - gsl_vector_complex_free(r->eigval_err); - gsl_vector_free(r->residuals); - gsl_matrix_complex_free(r->eigvec); - gsl_vector_free(r->ranktest_SV); + free(r->eigval); + free(r->eigval_err); + free(r->residuals); + free(r->eigvec); + free(r->ranktest_SV); free(r); } } 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->L = L; QPMS_ENSURE(L <= M, "We expect L <= M, but we got L = %d, M = %d ", L, M); // storage for eigenvalues and eigenvectors - solver->eigenvalues = gsl_vector_complex_calloc(L); - solver->eigenvalue_errors = gsl_vector_complex_calloc(L); - solver->residuals = gsl_vector_calloc(L); - solver->eigenvectors = gsl_matrix_complex_calloc(L, M); + //solver->eigenvalues = gsl_vector_complex_calloc(L); + QPMS_CRASHING_CALLOC(solver->eigenvalues, L, sizeof(complex double)); + //solver->eigenvalue_errors = gsl_vector_complex_calloc(L); + 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 - solver->A0 = gsl_matrix_complex_calloc(M,L); - solver->A1 = gsl_matrix_complex_calloc(M,L); - solver->A0_coarse = gsl_matrix_complex_calloc(M,L); - solver->A1_coarse = gsl_matrix_complex_calloc(M,L); - solver->MInvVHat = gsl_matrix_complex_calloc(M,L); - solver->VHat = gsl_matrix_complex_calloc(M,L); - solver->Sigma = gsl_vector_calloc(L); + //solver->A0 = gsl_matrix_complex_calloc(M,L); + QPMS_CRASHING_CALLOC(solver->A0, M * L, sizeof(complex double)); + //solver->A1 = gsl_matrix_complex_calloc(M,L); + QPMS_CRASHING_CALLOC(solver->A1, M * L, sizeof(complex double)); + //solver->A0_coarse = gsl_matrix_complex_calloc(M,L); + QPMS_CRASHING_CALLOC(solver->A0_coarse, M * L, sizeof(complex double)); + //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 BeynSolver_srandom(solver,(unsigned)time(NULL)); @@ -295,30 +302,30 @@ BeynSolver *BeynSolver_create(int M, int L) void BeynSolver_free(BeynSolver *solver) { - gsl_vector_complex_free(solver->eigenvalues); - gsl_vector_complex_free(solver->eigenvalue_errors); - gsl_matrix_complex_free(solver->eigenvectors); + free(solver->eigenvalues); + free(solver->eigenvalue_errors); + free(solver->eigenvectors); - gsl_matrix_complex_free(solver->A0); - gsl_matrix_complex_free(solver->A1); - gsl_matrix_complex_free(solver->A0_coarse); - gsl_matrix_complex_free(solver->A1_coarse); - gsl_matrix_complex_free(solver->MInvVHat); - gsl_vector_free(solver->Sigma); - gsl_vector_free(solver->residuals); - gsl_matrix_complex_free(solver->VHat); + free(solver->A0); + free(solver->A1); + free(solver->A0_coarse); + free(solver->A1_coarse); + free(solver->MInvVHat); + free(solver->Sigma); + free(solver->residuals); + free(solver->VHat); free(solver); } void BeynSolver_free_partial(BeynSolver *solver) { - gsl_matrix_complex_free(solver->A0); - gsl_matrix_complex_free(solver->A1); - gsl_matrix_complex_free(solver->A0_coarse); - gsl_matrix_complex_free(solver->A1_coarse); - gsl_matrix_complex_free(solver->MInvVHat); - gsl_matrix_complex_free(solver->VHat); + free(solver->A0); + free(solver->A1); + free(solver->A0_coarse); + free(solver->A1_coarse); + free(solver->MInvVHat); + free(solver->VHat); free(solver); } @@ -328,11 +335,8 @@ void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed) if (RandSeed==0) RandSeed=time(0); srandom(RandSeed); - gsl_matrix_complex *VHat=solver->VHat; - for(int nr=0; nrsize1; nr++) - for(int nc=0; ncsize2; nc++) - gsl_matrix_complex_set(VHat,nr,nc,cs2g(zrandN(1, 0))); - + for(size_t i = 0; i < solver->M * solver->L; ++i) + solver->VHat[i] = zrandN(1, 0); } @@ -342,44 +346,47 @@ void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed) * 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, - gsl_matrix_complex *A0, gsl_matrix_complex *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 *A0, complex double *A1, double complex z0, + 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 l = solver->L; - gsl_vector *Sigma = solver->Sigma; + double *Sigma = solver->Sigma; int verbose = 1; // TODO // Beyn Step 3: Compute SVD: A0 = V0_full * Sigma * W0T_full if(verbose) printf(" Beyn: computing SVD...\n"); - gsl_matrix_complex* V0_full = gsl_matrix_complex_alloc(m,l); - gsl_matrix_complex_memcpy(V0_full,A0); - gsl_matrix_complex* W0T_full = gsl_matrix_complex_alloc(l,l); - 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?"); + //gsl_matrix_complex* V0_full = gsl_matrix_complex_alloc(m,l); + //gsl_matrix_complex_memcpy(V0_full,A0); + complex double *V0_full; + QPMS_CRASHING_MALLOCPY(V0_full, A0, m * l * sizeof(complex double)); + //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') '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 */, - V0_full->size2 /* n, number of columns */, - (lapack_complex_double *)(V0_full->data) /*a*/, - V0_full->tda /*lda*/, - Sigma->data /*s*/, + m, // V0_full->size1 /* m, number of rows */, + l, // V0_full->size2 /* n, number of columns */, + V0_full, //(lapack_complex_double *)(V0_full->data) /*a*/, + l, //V0_full->tda /*lda*/, + Sigma, //Sigma->data /*s*/, NULL /*u; not used*/, m /*ldu; should not be really used but lapacke requires it for some obscure reason*/, - (lapack_complex_double *)W0T_full->data /*vt*/, - W0T_full->tda /*ldvt*/ + W0T_full, //(lapack_complex_double *)W0T_full->data /*vt*/, + l //W0T_full->tda /*ldvt*/ )); // Beyn Step 4: Rank test for Sigma // compute effective rank K (number of eigenvalue candidates) int K=0; - for (int k=0; ksize /* this is l, actually */; k++) { - if (verbose) printf("Beyn: SV(%d)=%e\n",k,gsl_vector_get(Sigma, k)); - if (k < rank_sel_min || gsl_vector_get(Sigma, k) > rank_tol) + for (int k=0; k < l/* ksize*/ /* this is l, actually */; k++) { + if (verbose) printf("Beyn: SV(%d)=%e\n",k, Sigma[k] /*gsl_vector_get(Sigma, k)*/); + if (k < rank_sel_min || Sigma[k] /*gsl_vector_get(Sigma, k)*/ > rank_tol) K++; } 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 // set V0, W0T = matrices of first K right, left singular vectors - gsl_matrix_complex *V0 = gsl_matrix_complex_alloc(m,K); - gsl_matrix_complex *W0T= gsl_matrix_complex_alloc(K,l); + //gsl_matrix_complex *V0 = gsl_matrix_complex_alloc(m,K); + //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) { - gsl_vector_complex_view tmp; - tmp = gsl_matrix_complex_column(V0_full, k); - gsl_matrix_complex_set_col(V0, k, &(tmp.vector)); - tmp = gsl_matrix_complex_row(W0T_full, k); - gsl_matrix_complex_set_row(W0T, k, &(tmp.vector)); + //gsl_vector_complex_view tmp; + //tmp = gsl_matrix_complex_column(V0_full, k); + //gsl_matrix_complex_set_col(V0, k, &(tmp.vector)); + for(int i = 0; i < m; ++i) + 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_free(W0T_full); - - gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,l); - gsl_matrix_complex *B = gsl_matrix_complex_calloc(K,K); + //gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,l); + //gsl_matrix_complex *B = gsl_matrix_complex_calloc(K,K); + complex double *TM2, *B; + QPMS_CRASHING_CALLOC(TM2, K * l, sizeof(complex double)); + QPMS_CRASHING_CALLOC(B, K * K, sizeof(complex double)); if(verbose) printf(" Multiplying V0*A1->TM...\n"); - const gsl_complex one = gsl_complex_rect(1,0); - const gsl_complex zero = gsl_complex_rect(0,0); - gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, - V0, A1, zero, TM2); + //const gsl_complex one = gsl_complex_rect(1,0); + //const gsl_complex zero = gsl_complex_rect(0,0); + //gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, + // 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"); - gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, - TM2, W0T, zero, B); + //gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, + // 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(TM2); + //gsl_matrix_complex_free(W0T); + //gsl_matrix_complex_free(TM2); + free(W0T); + free(TM2); 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++) { - 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_matrix_complex_set_col(B,i,tmp); + //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_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; mstride == 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(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); + // dims: B[K,K], S[K,K], Lambda[K] QPMS_ENSURE_SUCCESS(LAPACKE_zgeev( LAPACK_ROW_MAJOR, 'N' /* jobvl; don't compute left eigenvectors */, 'V' /* jobvr; do compute right eigenvectors */, K /* n */, - (lapack_complex_double *)(B->data) /* a */, - B->tda /* lda */, - (lapack_complex_double *) Lambda->data /* w */, + B, //(lapack_complex_double *)(B->data) /* a */, + K, //B->tda /* lda */, + Lambda, //(lapack_complex_double *) Lambda->data /* w */, NULL /* vl */, - m /* ldvl, not used by for some reason needed */, - (lapack_complex_double *)(S->data)/* vr */, - S->tda/* ldvr */ + m /* ldvl, not used but for some reason needed */, + S, //(lapack_complex_double *)(S->data)/* vr */, + K //S->tda/* ldvr */ )); - gsl_matrix_complex_free(B); + //gsl_matrix_complex_free(B); + free(B); // V0S <- V0*S printf("Multiplying V0*S...\n"); - gsl_matrix_complex *V0S = gsl_matrix_complex_alloc(m, K); - QPMS_ENSURE_SUCCESS(gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, - one, V0, S, zero, V0S)); + //gsl_matrix_complex *V0S = gsl_matrix_complex_alloc(m, K); + //QPMS_ENSURE_SUCCESS(gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, + // 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! // (If they differ, there are possibly some spurious eigenvalues.) int KRetained = 0; - gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m, m); - gsl_vector_complex *MVk = gsl_vector_complex_alloc(m); + //gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m, 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) { - 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); - gsl_vector_complex_const_view Vk = gsl_matrix_complex_const_column(V0S, 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 = z0 + Lambda[k]; + //gsl_vector_complex_const_view Vk = gsl_matrix_complex_const_column(V0S, k); double residual = 0; if(res_tol > 0) { - QPMS_ENSURE_SUCCESS(M_function(Mmat, z, Params)); - QPMS_ENSURE_SUCCESS(gsl_blas_zgemv(CblasNoTrans, one, Mmat, &(Vk.vector), zero, MVk)); - residual = gsl_blas_dznrm2(MVk); + QPMS_ENSURE_SUCCESS(M_function(Mmat, m, z, Params)); + //QPMS_ENSURE_SUCCESS(gsl_blas_zgemv(CblasNoTrans, one, Mmat, &(Vk.vector), zero, 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 (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) { - gsl_matrix_complex_set_row(eigenvectors, KRetained, &(Vk.vector)); - gsl_vector_set(solver->residuals, KRetained, residual); + //gsl_matrix_complex_set_row(eigenvectors, KRetained, &(Vk.vector)); + 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; } - gsl_matrix_complex_free(V0S); - gsl_matrix_complex_free(Mmat); - gsl_vector_complex_free(MVk); - gsl_matrix_complex_free(S); - gsl_vector_complex_free(Lambda); + free(V0S); + free(Mmat); + free(MVk); + free(S); + free(Lambda); return KRetained; } -beyn_result_gsl_t *beyn_solve_gsl(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_result_t *beyn_solve(const size_t m, const size_t l, + beyn_function_M_t M_function, beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *params, const beyn_contour_t *contour, double rank_tol, size_t rank_sel_min, double res_tol) { BeynSolver *solver = BeynSolver_create(m, l); - gsl_matrix_complex *A0 = solver->A0; - gsl_matrix_complex *A1 = solver->A1; - gsl_matrix_complex *A0_coarse = solver->A0_coarse; - gsl_matrix_complex *A1_coarse = solver->A1_coarse; - gsl_matrix_complex *MInvVHat = solver->MInvVHat; - gsl_matrix_complex *VHat = solver->VHat; + complex double *A0 = solver->A0; + complex double *A1 = solver->A1; + complex double *A0_coarse = solver->A0_coarse; + complex double *A1_coarse = solver->A1_coarse; + complex double *MInvVHat = solver->MInvVHat; + complex double *VHat = solver->VHat; /***************************************************************/ /* evaluate contour integrals by numerical quadrature to get */ /* A0 and A1 matrices */ /***************************************************************/ - gsl_matrix_complex_set_zero(A0); - gsl_matrix_complex_set_zero(A1); - gsl_matrix_complex_set_zero(A0_coarse); - gsl_matrix_complex_set_zero(A1_coarse); + // TODO zeroing probably redundant (Used calloc...) + memset(A0, 0, m * l * sizeof(complex double)); + memset(A1, 0, m * l * sizeof(complex double)); + memset(A0_coarse, 0, m * l * sizeof(complex double)); + memset(A1_coarse, 0, m * l * sizeof(complex double)); const size_t N = contour->n; if(N & 1) QPMS_WARN("Contour discretisation point number is odd (%zd)," " 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 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) { - 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 { lapack_int *pivot; - gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m,m); - QPMS_ENSURE_SUCCESS(M_function(Mmat, z, params)); + //gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m,m); + 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); +#if 0 QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, m /*m*/, m /*n*/,(lapack_complex_double *) Mmat->data /*a*/, Mmat->tda /*lda*/, pivot /*ipiv*/)); QPMS_ENSURE(MInvVHat->tda == l, "wut?"); 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*/, (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); - gsl_matrix_complex_free(Mmat); + free(Mmat); } - gsl_matrix_complex_scale(MInvVHat, cs2g(dz)); - gsl_matrix_complex_add(A0, MInvVHat); + // gsl_matrix_complex_scale(MInvVHat, cs2g(dz)); + 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) { - 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. - gsl_matrix_complex_add(A1, MInvVHat); + // A_1 scaling as in Beyn's Remark 3.2(b) for numerical stability. + //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) { - gsl_matrix_complex_add(A1_coarse, MInvVHat); - 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); + A1_coarse[i] += 2 * MInvVHat[i]; } } - gsl_vector_complex *eigenvalues = solver->eigenvalues; - gsl_vector_complex *eigenvalue_errors = solver->eigenvalue_errors; - gsl_matrix_complex *eigenvectors = solver->eigenvectors; + complex double *eigenvalues = solver->eigenvalues; + complex double *eigenvalue_errors = solver->eigenvalue_errors; + complex double *eigenvectors = solver->eigenvectors; // 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); @@ -589,10 +668,14 @@ beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l, // 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); - gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors); - beyn_result_gsl_t *result; - QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_gsl_t)); + //gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors); + 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_err = solver->eigenvalue_errors; 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; } +#if 0 // Wrapper of pure C array M-matrix function to GSL. struct beyn_function_M_carr2gsl_param { @@ -664,4 +748,5 @@ void beyn_result_free(beyn_result_t *result) { beyn_result_gsl_free(result->gsl); free(result); } +#endif diff --git a/qpms/beyn.h b/qpms/beyn.h index 86f7702..13a43aa 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -4,20 +4,10 @@ #ifndef BEYN_H #define BEYN_H +#include #include -#include -#include -/// User-supplied function that provides the operator 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. +/// User-supplied function that provides the (row-major) m × m matrix M(z) whose "roots" are to be found. /** Pure C array version */ 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); -/// 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). typedef struct beyn_result_t { size_t neig; ///< Number of eigenvalues found. @@ -88,31 +66,11 @@ typedef struct beyn_result_t { complex double *eigvec; // Rows are the eigenvectors double *ranktest_SV; - /// Private, we wrap it around beyn_result_gsl_t for now. - beyn_result_gsl_t *gsl; - } 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); -beyn_result_gsl_t *beyn_solve_gsl( - 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. - ); - +/// Solve a non-linear eigenproblem using Beyn's algorithm beyn_result_t *beyn_solve( 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). @@ -125,7 +83,4 @@ beyn_result_t *beyn_solve( 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 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index ff416f8..1aacccd 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -72,13 +72,13 @@ target_link_libraries( tbeyn3f qpms gsl lapacke ${QPMS_AMOSLIB} m ) target_include_directories( tbeyn3f PRIVATE .. ) target_compile_definitions( tbeyn3f PRIVATE VARIANTF ) -add_executable( tbeyn_gsl tbeyn_gsl.c ) -target_link_libraries( tbeyn_gsl qpms gsl lapacke ${QPMS_AMOSLIB} m ) -target_include_directories( tbeyn_gsl PRIVATE .. ) +#add_executable( tbeyn_gsl tbeyn_gsl.c ) +#target_link_libraries( tbeyn_gsl qpms gsl lapacke amos m ) +#target_include_directories( tbeyn_gsl PRIVATE .. ) -add_executable( tbeyn_gsl2 tbeyn_gsl2.c ) -target_link_libraries( tbeyn_gsl2 qpms gsl lapacke ${QPMS_AMOSLIB} m ) -target_include_directories( tbeyn_gsl2 PRIVATE .. ) +#add_executable( tbeyn_gsl2 tbeyn_gsl2.c ) +#target_link_libraries( tbeyn_gsl2 qpms gsl lapacke amos m ) +#target_include_directories( tbeyn_gsl2 PRIVATE .. ) add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array tbeyn ) diff --git a/tests/tbeyn.c b/tests/tbeyn.c index 2cef25b..62117de 100644 --- a/tests/tbeyn.c +++ b/tests/tbeyn.c @@ -1,6 +1,7 @@ #include #include #include +#include // Matrix as in Beyn, section 4.11 int M_function(complex double *target, const size_t m, const complex double z, void *no_params) { From 2be35213335824ff7eabcb94f28c9b00ef358ecc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 10 Jun 2021 22:25:56 +0300 Subject: [PATCH 2/4] beyn pure BLAS: fix indices and strides Also fix amos reference in CMakeLists --- qpms/beyn.c | 6 +++--- tests/CMakeLists.txt | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index e69b96d..6d59d1c 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -412,8 +412,8 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_t M_functio 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); + for(int j = 0; j < l; ++j) + MAT(W0T, K, l, k, j) = MAT(W0T_full, l, l, k, j); } //gsl_matrix_complex_free(V0_full); free(V0_full); @@ -535,7 +535,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_t M_functio //QPMS_ENSURE_SUCCESS(gsl_blas_zgemv(CblasNoTrans, one, Mmat, &(Vk.vector), zero, 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 */, + &one, Mmat, m, &MAT(V0S, m, K, 0, k), K /* stride of Vk in V0S */, &zero, MVk, 1); //residual = gsl_blas_dznrm2(MVk); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1aacccd..ccd5fea 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -73,11 +73,11 @@ target_include_directories( tbeyn3f PRIVATE .. ) target_compile_definitions( tbeyn3f PRIVATE VARIANTF ) #add_executable( tbeyn_gsl tbeyn_gsl.c ) -#target_link_libraries( tbeyn_gsl qpms gsl lapacke amos m ) +#target_link_libraries( tbeyn_gsl qpms gsl lapacke ${QPMS_AMOSLIB} m ) #target_include_directories( tbeyn_gsl PRIVATE .. ) #add_executable( tbeyn_gsl2 tbeyn_gsl2.c ) -#target_link_libraries( tbeyn_gsl2 qpms gsl lapacke amos m ) +#target_link_libraries( tbeyn_gsl2 qpms gsl lapacke ${QPMS_AMOSLIB} m ) #target_include_directories( tbeyn_gsl2 PRIVATE .. ) add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array tbeyn ) From f19c590d6ee3974be3497b56c89f82ecaa67c7e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 11 Jun 2021 10:38:49 +0300 Subject: [PATCH 3/4] Pure BLAS Beyn: cleanup commented out GSL code --- qpms/beyn.c | 74 ++--------------------------------------------------- 1 file changed, 2 insertions(+), 72 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 6d59d1c..492bc32 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -269,29 +269,18 @@ BeynSolver *BeynSolver_create(int M, int L) QPMS_ENSURE(L <= M, "We expect L <= M, but we got L = %d, M = %d ", L, M); // storage for eigenvalues and eigenvectors - //solver->eigenvalues = gsl_vector_complex_calloc(L); QPMS_CRASHING_CALLOC(solver->eigenvalues, L, sizeof(complex double)); - //solver->eigenvalue_errors = gsl_vector_complex_calloc(L); 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 - //solver->A0 = gsl_matrix_complex_calloc(M,L); QPMS_CRASHING_CALLOC(solver->A0, M * L, sizeof(complex double)); - //solver->A1 = gsl_matrix_complex_calloc(M,L); QPMS_CRASHING_CALLOC(solver->A1, M * L, sizeof(complex double)); - //solver->A0_coarse = gsl_matrix_complex_calloc(M,L); QPMS_CRASHING_CALLOC(solver->A0_coarse, M * L, sizeof(complex double)); - //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 BeynSolver_srandom(solver,(unsigned)time(NULL)); @@ -359,14 +348,9 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_t M_functio // Beyn Step 3: Compute SVD: A0 = V0_full * Sigma * W0T_full if(verbose) printf(" Beyn: computing SVD...\n"); - //gsl_matrix_complex* V0_full = gsl_matrix_complex_alloc(m,l); - //gsl_matrix_complex_memcpy(V0_full,A0); complex double *V0_full; QPMS_CRASHING_MALLOCPY(V0_full, A0, m * l * sizeof(complex double)); - //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') 'O' /*jobz, 'O' overwrites a with U and saves conjg(V') in vt if m >= n, i.e. if M >= L, which holds*/, m, // V0_full->size1 /* m, number of rows */, @@ -385,8 +369,8 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_t M_functio // compute effective rank K (number of eigenvalue candidates) int K=0; for (int k=0; k < l/* ksize*/ /* this is l, actually */; k++) { - if (verbose) printf("Beyn: SV(%d)=%e\n",k, Sigma[k] /*gsl_vector_get(Sigma, k)*/); - if (k < rank_sel_min || Sigma[k] /*gsl_vector_get(Sigma, k)*/ > rank_tol) + if (verbose) printf("Beyn: SV(%d)=%e\n",k, Sigma[k] ); + if (k < rank_sel_min || Sigma[k] > rank_tol) K++; } if (verbose)printf(" Beyn: %d/%zd relevant singular values\n",K,l); @@ -397,41 +381,26 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_t M_functio // Beyn step 5: B = V0' * A1 * W0 * Sigma^-1 // set V0, W0T = matrices of first K right, left singular vectors - //gsl_matrix_complex *V0 = gsl_matrix_complex_alloc(m,K); - //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) { - //gsl_vector_complex_view tmp; - //tmp = gsl_matrix_complex_column(V0_full, k); - //gsl_matrix_complex_set_col(V0, k, &(tmp.vector)); for(int i = 0; i < m; ++i) 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 < l; ++j) MAT(W0T, K, l, k, j) = MAT(W0T_full, l, l, k, j); } - //gsl_matrix_complex_free(V0_full); free(V0_full); - //gsl_matrix_complex_free(W0T_full); free(W0T_full); - //gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,l); - //gsl_matrix_complex *B = gsl_matrix_complex_calloc(K,K); complex double *TM2, *B; QPMS_CRASHING_CALLOC(TM2, K * l, sizeof(complex double)); QPMS_CRASHING_CALLOC(B, K * K, sizeof(complex double)); if(verbose) printf(" Multiplying V0*A1->TM...\n"); - //const gsl_complex one = gsl_complex_rect(1,0); - //const gsl_complex zero = gsl_complex_rect(0,0); - //gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, - // 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, @@ -440,30 +409,20 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_t M_functio if(verbose) printf(" Multiplying TM*W0T->B...\n"); - //gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, // 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(TM2); free(W0T); free(TM2); if(verbose) printf(" Scaling B <- B*Sigma^{-1}\n"); - //gsl_vector_complex *tmp = gsl_vector_complex_calloc(K); for(int i = 0; i < K; 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_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); - //for(int m=0; m S*Lambda*S' /* According to Beyn's paper (Algorithm 1), one should check conditioning @@ -475,14 +434,10 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_t M_functio */ if(verbose) printf(" Eigensolving (%i,%i)\n",K,K); - //gsl_vector_complex *Lambda = gsl_vector_complex_alloc(K); // eigenvalues - //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(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( LAPACK_ROW_MAJOR, @@ -498,59 +453,44 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_t M_functio K //S->tda/* ldvr */ )); - //gsl_matrix_complex_free(B); free(B); // V0S <- V0*S printf("Multiplying V0*S...\n"); - //gsl_matrix_complex *V0S = gsl_matrix_complex_alloc(m, K); - //QPMS_ENSURE_SUCCESS(gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, - // 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); free(V0); // FIXME!!! make clear relation between KRetained and K in the results! // (If they differ, there are possibly some spurious eigenvalues.) int KRetained = 0; - //gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m, 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) { - //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 = z0 + Lambda[k]; - //gsl_vector_complex_const_view Vk = gsl_matrix_complex_const_column(V0S, k); double residual = 0; if(res_tol > 0) { QPMS_ENSURE_SUCCESS(M_function(Mmat, m, z, Params)); - //QPMS_ENSURE_SUCCESS(gsl_blas_zgemv(CblasNoTrans, one, Mmat, &(Vk.vector), zero, 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), K /* 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 (res_tol > 0 && residual > res_tol) continue; - //gsl_vector_complex_set(eigenvalues, KRetained, zgsl); eigenvalues[KRetained] = z; if(eigenvectors) { - //gsl_matrix_complex_set_row(eigenvectors, KRetained, &(Vk.vector)); 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; @@ -607,7 +547,6 @@ beyn_result_t *beyn_solve(const size_t m, const size_t l, QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, m, l, VHat, z, params)); } else { lapack_int *pivot; - //gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m,m); complex double *Mmat; QPMS_CRASHING_MALLOC(Mmat, m * m * sizeof(complex double)); QPMS_ENSURE_SUCCESS(M_function(Mmat, m, z, params)); @@ -630,30 +569,22 @@ beyn_result_t *beyn_solve(const size_t m, const size_t l, free(Mmat); } - // gsl_matrix_complex_scale(MInvVHat, cs2g(dz)); 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) { - //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]; } // A_1 scaling as in Beyn's Remark 3.2(b) for numerical stability. - //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) { for(size_t i = 0; i < m * l; ++i) - //gsl_matrix_complex_add(A1_coarse, MInvVHat); - //gsl_matrix_complex_add(A1_coarse, MInvVHat); A1_coarse[i] += 2 * MInvVHat[i]; } } @@ -669,7 +600,6 @@ beyn_result_t *beyn_solve(const size_t m, const size_t l, // 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); - //gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors); 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); From 9d556e7b23dedffd073ff163f9e4fc88281c9326 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 11 Jun 2021 10:48:58 +0300 Subject: [PATCH 4/4] Remove all gsl_matrix_... related code in Beyn. Now that pure BLAS implementation is working, this in no longer needed. Getting rid of gsl_matrix and gsl_cblas eliminates header clashes between GSL and other CBLAS implementations. --- qpms/beyn.c | 65 -------------------------------------------- qpms/qpms_cdefs.pxd | 27 ------------------ tests/CMakeLists.txt | 8 ------ tests/tbeyn_gsl.c | 42 ---------------------------- tests/tbeyn_gsl2.c | 37 ------------------------- 5 files changed, 179 deletions(-) delete mode 100644 tests/tbeyn_gsl.c delete mode 100644 tests/tbeyn_gsl2.c diff --git a/qpms/beyn.c b/qpms/beyn.c index 492bc32..8d98ff2 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -1,8 +1,5 @@ #define STATIC_ASSERT(COND,MSG) typedef char static_assertion_##MSG[(COND)?1:-1] -#define cg2s(x) gsl_complex_tostd(x) -#define cs2g(x) gsl_complex_fromstd(x) - #include #include #include @@ -618,65 +615,3 @@ beyn_result_t *beyn_solve(const size_t m, const size_t l, return result; } -#if 0 -// Wrapper of pure C array M-matrix function to GSL. - -struct beyn_function_M_carr2gsl_param { - beyn_function_M_t M_function; - beyn_function_M_inv_Vhat_t M_inv_Vhat_function; - void *param; -}; - -static int beyn_function_M_carr2gsl(gsl_matrix_complex *target_M, complex double z, void *params) -{ - struct beyn_function_M_carr2gsl_param *p = params; - // These could rather be asserts. - QPMS_ENSURE(target_M->size2 == target_M->tda, "Target GSL matrix is not a C-contiguous array. This is a bug, please report!"); - QPMS_ENSURE(target_M->size1 == target_M->size2, "Target is not a square matrix. This is a bug, please report!"); - return p->M_function((complex double *) target_M->data, target_M->size1, z, p->param); -} - -static int beyn_function_M_inv_Vhat_carr2gsl(gsl_matrix_complex *target, - const gsl_matrix_complex *Vhat, complex double z, void *params) -{ - QPMS_UNTESTED; - struct beyn_function_M_carr2gsl_param *p = params; - // These could rather be asserts. - QPMS_ENSURE(target->size2 == target->tda, "Target GSL matrix is not a C-contiguous array. This is a bug, please report!"); - QPMS_ENSURE(Vhat->size2 == Vhat->tda, "Source GSL matrix is not a C-contiguous array. This is a bug, please report!"); - // TODO here I could also check whether Vhat and target have compatible sizes. - return p->M_inv_Vhat_function((complex double *) target->data, target->size1, target->size2, - (complex double *) Vhat->data, z, p->param); -} - -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, - void *params, const beyn_contour_t *contour, double rank_tol, size_t rank_sel_min, double res_tol) { - struct beyn_function_M_carr2gsl_param p = {M, M_inv_Vhat, params}; - return beyn_result_from_beyn_result_gsl( - beyn_solve_gsl(m, l, beyn_function_M_carr2gsl, - (p.M_inv_Vhat_function) ? beyn_function_M_inv_Vhat_carr2gsl : NULL, - (void *) &p, contour, rank_tol, rank_sel_min, res_tol) - ); -} - -beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g) { - struct beyn_result_t *result; - QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_t)); - result->gsl = g; - result->neig = result->gsl->neig; - result->vlen = result->gsl->eigvec->size2; - result->eigval = (complex double *) result->gsl->eigval->data; - result->eigval_err = (complex double *) result->gsl->eigval_err->data; - result->residuals = result->gsl->residuals->data; - result->eigvec = (complex double *) result->gsl->eigvec->data; - result->ranktest_SV = result->gsl->ranktest_SV->data; - return result; -} - -void beyn_result_free(beyn_result_t *result) { - if(result) - beyn_result_gsl_free(result->gsl); - free(result); -} -#endif - diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 2a9f8c9..795323b 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -741,25 +741,10 @@ cdef extern from "ewald.h": double eta, cdouble wavenumber, LatticeDimensionality latdim, PGen *pgen_R, bint pgen_generates_shifted_points, cart3_t k, cart3_t particle_shift) - -cdef extern from "gsl/gsl_complex.h": - ctypedef struct gsl_complex: - double dat[2] - -cdef extern from "gsl/gsl_matrix.h": - ctypedef struct gsl_matrix_complex: - pass - ctypedef struct gsl_vector: - pass - ctypedef struct gsl_vector_complex: - pass - cdef extern from "beyn.h": ctypedef struct beyn_contour_t: bint (*inside_test)(beyn_contour_t *, cdouble z) pass - ctypedef struct beyn_result_gsl_t: - pass ctypedef struct beyn_result_t: size_t neig size_t vlen @@ -768,25 +753,17 @@ cdef extern from "beyn.h": double *residuals cdouble *eigvec double *ranktest_SV - beyn_result_gsl_t *gsl ctypedef enum beyn_contour_halfellipse_orientation: BEYN_CONTOUR_HALFELLIPSE_RE_PLUS BEYN_CONTOUR_HALFELLIPSE_IM_PLUS BEYN_CONTOUR_HALFELLIPSE_RE_MINUS BEYN_CONTOUR_HALFELLIPSE_IM_MINUS - ctypedef int (*beyn_function_M_gsl_t)(gsl_matrix_complex *target_M, cdouble z, void *params) - ctypedef int (*beyn_function_M_inv_Vhat_gsl_t)(gsl_matrix_complex *target, const gsl_matrix_complex *Vhat, cdouble z, void *params) ctypedef int (*beyn_function_M_t)(cdouble *target_M, size_t m, cdouble z, void *params) ctypedef int (*beyn_function_M_inv_Vhat_t)(cdouble *target, size_t m, size_t l, const cdouble *Vhat, cdouble z, void *params) - void beyn_result_gsl_free(beyn_result_gsl_t *result) void beyn_result_free(beyn_result_t *result) - beyn_result_gsl_t *beyn_solve_gsl(size_t m, size_t l, beyn_function_M_gsl_t M, - beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat, void *params, const beyn_contour_t *contour, - double rank_tol, size_t rank_min_sel, double res_tol) - 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, void *params, const beyn_contour_t *contour, double rank_tol, size_t rank_min_sel, double res_tol) @@ -797,7 +774,3 @@ cdef extern from "beyn.h": beyn_contour_t *beyn_contour_kidney(cdouble centre, double halfax_re, double halfax_im, size_t npoints, double rounding, beyn_contour_halfellipse_orientation ori) - - cdouble gsl_comlpex_tostd(gsl_complex z) - gsl_complex gsl_complex_fromstd(cdouble z) - diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index ccd5fea..63df6ae 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -72,14 +72,6 @@ target_link_libraries( tbeyn3f qpms gsl lapacke ${QPMS_AMOSLIB} m ) target_include_directories( tbeyn3f PRIVATE .. ) target_compile_definitions( tbeyn3f PRIVATE VARIANTF ) -#add_executable( tbeyn_gsl tbeyn_gsl.c ) -#target_link_libraries( tbeyn_gsl qpms gsl lapacke ${QPMS_AMOSLIB} m ) -#target_include_directories( tbeyn_gsl PRIVATE .. ) - -#add_executable( tbeyn_gsl2 tbeyn_gsl2.c ) -#target_link_libraries( tbeyn_gsl2 qpms gsl lapacke ${QPMS_AMOSLIB} m ) -#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_test( NAME single_vs_array_translation_coeffs COMMAND test_single_translations_vs_calc ) diff --git a/tests/tbeyn_gsl.c b/tests/tbeyn_gsl.c deleted file mode 100644 index b2dd1d8..0000000 --- a/tests/tbeyn_gsl.c +++ /dev/null @@ -1,42 +0,0 @@ -#include -#include - - -// Matrix as in Beyn, section 4.11 -int M_function(gsl_matrix_complex *target, complex double z, void *no_params) { - int m = target->size1; - - gsl_complex d = gsl_complex_fromstd( 2*m - 4*z / (6*m) ); - gsl_complex od = gsl_complex_fromstd( -(double)m - z / (6*m) ); - - gsl_matrix_complex_set_zero(target); - for (int i = 0; i < m; ++i) { - gsl_matrix_complex_set(target, i, i, d); - if(i > 0) gsl_matrix_complex_set(target, i, i-1, od); - if(i < m - 1) gsl_matrix_complex_set(target, i, i+1, od); - } - gsl_matrix_complex_set(target, m-1, m-1, gsl_complex_fromstd(gsl_complex_tostd(d)/2 + z/(z-1))); - - return 0; -} - -int main() { - complex double z0 = 150+2*I; - double Rx = 148, Ry = 148; - int L = 10, N = 50, dim = 400; - beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); - - beyn_result_gsl_t *result = - beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, - contour, 1e-4, 1, 1e-4); - printf("Found %zd eigenvalues:\n", result->neig); - for (size_t i = 0; i < result->neig; ++i) { - gsl_complex eig = gsl_vector_complex_get(result->eigval, i); - printf("%zd: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); - } - free(contour); - beyn_result_gsl_free(result); - return 0; -} - - diff --git a/tests/tbeyn_gsl2.c b/tests/tbeyn_gsl2.c deleted file mode 100644 index 1da782d..0000000 --- a/tests/tbeyn_gsl2.c +++ /dev/null @@ -1,37 +0,0 @@ -#include -#include - - -// Matrix as in Beyn, section 4.11 -int M_function(gsl_matrix_complex *target, complex double z, void *no_params) { - int m = target->size1; - - gsl_matrix_complex_set_zero(target); - for (int i = 0; i < m; ++i) { - gsl_matrix_complex_set(target, i, i, gsl_complex_rect(i - creal(z), -cimag(z))); - } - - return 0; -} - -int main() { - complex double z0 = 150+2*I; - double Rx = 5, Ry = 5; - int L = 10, N = 600, dim = 200; - beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); - - beyn_result_gsl_t *result = - beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, - contour, 1e-4, 1, 1e-4); - printf("Found %zd eigenvalues:\n", result->neig); - for (size_t i = 0; i < result->neig; ++i) { - gsl_complex eig = gsl_vector_complex_get(result->eigval, i); - printf("%zd: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); - - } - free(contour); - beyn_result_gsl_free(result); - return 0; -} - -