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

Incremental cleanup. 

Former-commit-id: d4d8d41a2edac7cf8ac341cce46e3c976ef68c5e
This commit is contained in:
Marek Nečada 2019-09-02 13:31:22 +03:00
parent 5dc2a44cdd
commit 1585c48071
1 changed files with 42 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

112
qpms/beyn.c
View File

@ -33,39 +33,19 @@ typedef struct BeynSolver
gsl_matrix_complex *A0, *A1, *A0_coarse, *A1_coarse, *MInvVHat; gsl_matrix_complex *A0, *A1, *A0_coarse, *A1_coarse, *MInvVHat;
gsl_matrix_complex *VHat; gsl_matrix_complex *VHat;
gsl_vector *Sigma, *residuals; gsl_vector *Sigma, *residuals;
double complex *Workspace;
} BeynSolver; } BeynSolver;
// constructor, destructor // constructor, destructor
BeynSolver *BeynSolver_create(int M, int L); BeynSolver *BeynSolver_create(int M, int L);
void BeynSolver_free(BeynSolver *solver); void BeynSolver_free(BeynSolver *solver);
// reset the random matrix VHat used in the Beyn algorithm // reset the random matrix VHat used in Beyn's algorithm
//
void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed); void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed);
// for both of the following routines, // Uniformly random number from interval [a, b].
// the return value is the number of eigenvalues found, static double randU(double a, double b) { return a + (b-a) * random() * (1. / RAND_MAX); }
// and the eigenvalues and eigenvectors are stored in the
// Lambda and eigenvectors fields of the BeynSolver structure
// Beyn method for circular contour of radius R,
// centered at z0, using N quadrature points
//int BeynSolve(BeynSolver *solver,
// BeynFunction UserFunction, void *Params,
// double complex z0, double R, int N);
// Beyn method for elliptical contour of horizontal, vertical
// radii Rx, Ry, centered at z0, using N quadrature points
int BeynSolve(BeynSolver *solver,
beyn_function_M_gsl_t M_function, beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat_function, void *params,
double complex z0, double Rx, double Ry, int N);
static double randU(double A, double B) { return A + (B-A) * random() * (1. / RAND_MAX); }
// Random number from normal distribution (via Box-Muller transform, which is enough for our purposes).
static double randN(double Sigma, double Mu) { static double randN(double Sigma, double Mu) {
double u1 = randU(0, 1); double u1 = randU(0, 1);
double u2 = randU(0, 1); double u2 = randU(0, 1);
@ -76,7 +56,6 @@ static complex double zrandN(double sigma, double mu){
return randN(sigma, mu) + I*randN(sigma, mu); return randN(sigma, mu) + I*randN(sigma, mu);
} }
beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double rIm, size_t n) beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double rIm, size_t n)
{ {
beyn_contour_t *c; beyn_contour_t *c;
@ -174,30 +153,29 @@ void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed)
} }
/***************************************************************/ /*
/* perform linear-algebra manipulations on the A0 and A1 */ * linear-algebra manipulations on the A0 and A1 matrices
/* matrices (obtained via numerical quadrature) to extract */ * (obtained via numerical quadrature) to extract eigenvalues
/* eigenvalues 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_gsl_t M_function,
void *Params, void *Params,
gsl_matrix_complex *A0, gsl_matrix_complex *A1, double complex z0, gsl_matrix_complex *A0, gsl_matrix_complex *A1, double complex z0,
gsl_vector_complex *eigenvalues, gsl_matrix_complex *eigenvectors, const double RankTol, const double ResTol) gsl_vector_complex *eigenvalues, gsl_matrix_complex *eigenvectors, const double rank_tol, const double res_tol)
{ {
int M = solver->M; const size_t m = solver->M;
int L = solver->L; const size_t l = solver->L;
gsl_vector *Sigma = solver->Sigma; gsl_vector *Sigma = solver->Sigma;
int verbose = 1; // TODO int verbose = 1; // TODO
// A0 -> V0_full * Sigma * W0T_full' // A0 -> V0_full * Sigma * W0T_full'
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); gsl_matrix_complex* W0T_full = gsl_matrix_complex_alloc(l,l);
//A0->SVD(Sigma, &V0_full, &W0T_full); //A0->SVD(Sigma, &V0_full, &W0T_full);
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);
@ -210,7 +188,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
V0_full->tda /*lda*/, V0_full->tda /*lda*/,
Sigma->data /*s*/, 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*/, (lapack_complex_double *)W0T_full->data /*vt*/,
W0T_full->tda /*ldvt*/ W0T_full->tda /*ldvt*/
)); ));
@ -218,20 +196,20 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
// 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<Sigma->size /* this is l, actually */; k++) {
{ if (verbose) printf("Beyn: SV(%i)=%e\n",k,gsl_vector_get(Sigma, k)); if (verbose) printf("Beyn: SV(%d)=%e\n",k,gsl_vector_get(Sigma, k));
if (gsl_vector_get(Sigma, k) > RankTol ) if (gsl_vector_get(Sigma, k) > rank_tol)
K++; K++;
} }
printf(" Beyn: %i/%i relevant singular values\n",K,L); if (verbose)printf(" Beyn: %d/%zd relevant singular values\n",K,l);
if (K==0) if (K==0) {
{ printf("no singular values found in Beyn eigensolver\n"); QPMS_WARN("no singular values found in Beyn eigensolver\n");
return 0; return 0;
} }
// 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);
for (int k = 0; k < K; ++k) { for (int k = 0; k < K; ++k) {
gsl_vector_complex_view tmp; gsl_vector_complex_view tmp;
@ -245,18 +223,17 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
gsl_matrix_complex_free(W0T_full); gsl_matrix_complex_free(W0T_full);
// B <- V0' * A1 * W0 * Sigma^-1 // B <- V0' * A1 * W0 * Sigma^-1
gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,L); gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,l);
gsl_matrix_complex *B = gsl_matrix_complex_calloc(K,K); gsl_matrix_complex *B = gsl_matrix_complex_calloc(K,K);
printf(" Multiplying V0*A1->TM...\n"); if(verbose) printf(" Multiplying V0*A1->TM...\n");
//V0.Multiply(A1, &TM2, "--transA C"); // TM2 <- V0' * A1 //V0.Multiply(A1, &TM2, "--transA C"); // TM2 <- V0' * A1
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);
printf(" Multiplying TM*W0T->B...\n"); if(verbose) printf(" Multiplying TM*W0T->B...\n");
//TM2.Multiply(&W0T, &B, "--transB C"); // B <- TM2 * W0
gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one,
TM2, W0T, zero, B); TM2, W0T, zero, B);
@ -264,9 +241,9 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
gsl_matrix_complex_free(W0T); gsl_matrix_complex_free(W0T);
gsl_matrix_complex_free(TM2); gsl_matrix_complex_free(TM2);
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);
@ -275,11 +252,8 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
//for(int m=0; m<K; m++) // B <- B * Sigma^{-1} //for(int m=0; m<K; m++) // B <- B * Sigma^{-1}
// for(int n=0; n<K; n++)
// B.ScaleEntry(m,n,1.0/Sigma->GetEntry(n));
// B -> S*Lambda*S' // B -> S*Lambda*S'
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
@ -295,7 +269,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
B->tda /* lda */, B->tda /* lda */,
(lapack_complex_double *) Lambda->data /* w */, (lapack_complex_double *) Lambda->data /* w */,
NULL /* vl */, NULL /* vl */,
M /* ldvl, not used by for some reason needed */, m /* ldvl, not used by for some reason needed */,
(lapack_complex_double *)(S->data)/* vr */, (lapack_complex_double *)(S->data)/* vr */,
S->tda/* ldvr */ S->tda/* ldvr */
)); ));
@ -304,28 +278,28 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
// 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));
gsl_matrix_complex_free(V0); gsl_matrix_complex_free(V0);
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);
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); gsl_vector_complex_const_view Vk = gsl_matrix_complex_const_column(V0S, k);
double residual = 0; double residual = 0;
if(ResTol > 0) { if(res_tol > 0) {
QPMS_ENSURE_SUCCESS(M_function(Mmat, z, Params)); QPMS_ENSURE_SUCCESS(M_function(Mmat, 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); residual = gsl_blas_dznrm2(MVk);
if (verbose) printf("Beyn: Residual(%i)=%e\n",k,residual); if (verbose) printf("Beyn: Residual(%i)=%e\n",k,residual);
} }
if (ResTol > 0 && residual > ResTol) continue; if (res_tol > 0 && residual > res_tol) continue;
gsl_vector_complex_set(eigenvalues, KRetained, zgsl); gsl_vector_complex_set(eigenvalues, KRetained, zgsl);
if(eigenvectors) { if(eigenvectors) {
@ -345,15 +319,13 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
} }
int beyn_solve_gsl(beyn_result_gsl_t **result, size_t m, size_t l, int beyn_solve_gsl(beyn_result_gsl_t **result, 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_gsl_t M_function, beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat_function,
void *params, const beyn_contour_t *contour, void *params, const beyn_contour_t *contour,
double rank_tol, double res_tol) double rank_tol, double res_tol)
{ {
BeynSolver *solver = BeynSolver_create(m, l); BeynSolver *solver = BeynSolver_create(m, l);
const int M = solver->M;
const int L = solver->L;
gsl_matrix_complex *A0 = solver->A0; gsl_matrix_complex *A0 = solver->A0;
gsl_matrix_complex *A1 = solver->A1; gsl_matrix_complex *A1 = solver->A1;
gsl_matrix_complex *A0_coarse = solver->A0_coarse; gsl_matrix_complex *A0_coarse = solver->A0_coarse;
@ -388,14 +360,14 @@ int beyn_solve_gsl(beyn_result_gsl_t **result, size_t m, size_t l,
QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z, params)); QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, 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)); QPMS_ENSURE_SUCCESS(M_function(Mmat, z, params));
QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * M); QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * m);
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*/));
free(pivot); free(pivot);