From 3fb9f23af5348a9b55c6856bb9bbf71a4152056c Mon Sep 17 00:00:00 2001 From: Kristian Arjas Date: Wed, 21 Aug 2019 11:55:05 +0300 Subject: [PATCH 01/34] Add Beyn eigenmode solver (cherry picked from commit 9ed6ee319a1d0b295ccbd690e7c0c0261111b456 [formerly 896e706a2da075c178b7e736761bdc01e4ea9800]) Former-commit-id: 6de5c14d48258b6105b76eb4fab5b555284caaa4 --- qpms/libBeyn.c | 431 +++++++++++++++++++++++++++++++++++++++++++++++++ qpms/libBeyn.h | 95 +++++++++++ 2 files changed, 526 insertions(+) create mode 100644 qpms/libBeyn.c create mode 100644 qpms/libBeyn.h diff --git a/qpms/libBeyn.c b/qpms/libBeyn.c new file mode 100644 index 0000000..0ad4698 --- /dev/null +++ b/qpms/libBeyn.c @@ -0,0 +1,431 @@ +/* Copyright (C) 2005-2011 M. T. Homer Reid + * + * This file is part of SCUFF-EM. + * + * SCUFF-EM is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * SCUFF-EM is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +/* + * libBeyn.cc -- implementation of Beyn's method for + * -- nonlinear eigenvalue problems + * + * Homer Reid -- 6/2016 + * + */ +#include +#include +#include +#include + +#include + + +//#include +//#include + +// Maybe GSL works? +#include +#include +#include +#include +#include + +#include + +//#define II cdouble(0.0,1.0) + +// Uniformly random number between -2 and 2 +double complex zrandN(){ + double a = (rand()*4.0/RAND_MAX) - 2; + double b = (rand()*4.0/RAND_MAX) - 2; + return a + b*I; +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +BeynSolver *CreateBeynSolver(int M, int L) +{ + BeynSolver *Solver= (BeynSolver *)malloc(sizeof(*Solver)); + + Solver->M = M; + Solver->L = L; + + int MLMax = (M>L) ? M : L; + int MLMin = (MEigenvalues = gsl_vector_complex_calloc(L); + Solver->EVErrors = gsl_vector_complex_calloc(L); + Solver->Residuals = gsl_vector_complex_calloc(L); + Solver->Eigenvectors = gsl_matrix_complex_calloc(M, L); + + // 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->A0Coarse = gsl_matrix_complex_calloc(M,L); + Solver->A1Coarse = 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_complex_calloc(MLMin); + ReRandomize(Solver,(unsigned)time(NULL)); + + // internal workspace: need storage for 2 MxL matrices + // plus 3 LxL matrices + #define MLBUFFERS 2 + #define LLBUFFERS 3 + size_t ML = MLMax*L, LL = L*L; + size_t WorkspaceSize = (MLBUFFERS*ML + LLBUFFERS*LL)*sizeof(double complex); + + Solver->Workspace = (double complex*)malloc( WorkspaceSize ); + + return Solver; + +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +void DestroyBeynSolver(BeynSolver *Solver) +{ + gsl_vector_complex_free(Solver->Eigenvalues); + gsl_vector_complex_free(Solver->EVErrors); + gsl_matrix_complex_free(Solver->Eigenvectors); + + gsl_matrix_complex_free(Solver->A0); + gsl_matrix_complex_free(Solver->A1); + gsl_matrix_complex_free(Solver->A0Coarse); + gsl_matrix_complex_free(Solver->A1Coarse); + gsl_matrix_complex_free(Solver->MInvVHat); + gsl_vector_complex_free(Solver->Sigma); + gsl_matrix_complex_free(Solver->VHat); + + free(Solver->Workspace); + + free(Solver); +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ + +void ReRandomize(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_set(VHat,nr,nc,zrandN()); + +} + + +/***************************************************************/ +/* perform linear-algebra manipulations on the A0 and A1 */ +/* matrices (obtained via numerical quadrature) to extract */ +/* eigenvalues and eigenvectors */ +/***************************************************************/ + +int ProcessAMatrices(BeynSolver *Solver, BeynFunction UserFunc, + void *Params, + gsl_matrix_complex *A0, gsl_matrix_complex *A1, double complex z0, + gsl_vector_complex *Eigenvalues, gsl_matrix_complex *Eigenvectors) +{ + int M = Solver->M; + int L = Solver->L; + gsl_vector_complex *Sigma = Solver->Sigma; + + + int Verbose = 0;//CheckEnv("SCUFF_BEYN_VERBOSE"); + double RankTol=1.0e-4; //CheckEnv("SCUFF_BEYN_RANK_TOL",&RankTol); + double ResTol=0.0; // CheckEnv("SCUFF_BEYN_RES_TOL",&ResTol); + + // A0 -> V0Full * Sigma * W0TFull' + printf(" Beyn: computing SVD...\n"); + gsl_matrix_complex* V0Full = gsl_matrix_complex_calloc(M,L); + gsl_matrix_memcpy(V0Full,A0); + + gsl_matrix_complex* W0TFull = gsl_matrix_complex_calloc(L,L); + //A0->SVD(Sigma, &V0Full, &W0TFull); + gsl_vector_complex *work = gsl_vector_complex_alloc(M); + + gsl_linalg_SV_decomp(V0Full, W0TFull, Sigma, work); + + + // compute effective rank K (number of eigenvalue candidates) + int K=0; + for(int k=0; ksize; k++) + { if (Verbose) printf("Beyn: SV(%i)=%e",k,gsl_vector_get(Sigma, k)); + if (gsl_vector_get(Sigma, k) > RankTol ) + K++; + } + printf(" Beyn: %i/%i relevant singular values\n",K,L); + if (K==0) + { printf("no singular values found in Beyn eigensolver\n"); + return 0; + } + + + // 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_vector_complex *TempM = gsl_vector_complex_calloc(M); + gsl_vector_complex *TempL = gsl_vector_complex_calloc(L); + for(int k=0; kTM...\n"); + //V0.Multiply(A1, &TM2, "--transA C"); // TM2 <- V0' * A1 + gsl_complex one = gsl_complex_rect(1,0); + gsl_complex zero = gsl_complex_rect(0,0); + gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, + V0, A1, zero, TM2); + + printf(" Multiplying TM*W0T->B...\n"); + //TM2.Multiply(&W0T, &B, "--transB C"); // B <- TM2 * W0 + + gsl_blas_zgemm(CblasNoTrans, CblasTrans, one, + TM2, W0T, zero, B); + + gsl_matrix_free(W0T); + gsl_matrix_free(TM2); + + 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_get_col(tmp, B, i); + gsl_vector_scale(tmp, 1.0/gsl_vector_get(Sigma,i)); + gsl_matrix_set_col(B,i,tmp); + } + gsl_vector_complex_free(tmp); + + //for(int m=0; mGetEntry(n)); + + + // B -> S*Lambda*S' + printf(" Eigensolving (%i,%i)\n",K,K); + + gsl_vector_complex *Lambda = gsl_vector_complex_calloc(K); // Eigenvalues + gsl_matrix_complex *S = gsl_matrix_complex_calloc(K,K); // Eigenvectors + gsl_matrix_complex *Eye = gsl_matrix_complex_alloc(K,K); + + gsl_vector_complex *alph = gsl_vector_complex_calloc(K); + gsl_vector_complex *beta = gsl_vector_complex_calloc(K); + + gsl_matrix_set_identity(Eye); + + gsl_eigen_genv_workspace * W = gsl_eigen_genv_alloc(K); + gsl_eigen_genv(B, Eye, alph, beta, S,W); + gsl_eigen_genv_free(W); + + gsl_complex tmpa; + gsl_complex tmpb; + for(int i = 0; i < K; i++){ + tmpb = gsl_vector_complex_get(beta,i); + tmpa = gsl_vector_complex_get(alph,i); + if(gsl_complex_abs(tmpb)){ + gsl_vector_complex_set(Lambda, i, gsl_complex_div(tmpa,tmpb)); + printf("Eigenvalue %e + %e i found\n",GSL_REAL(gsl_complex_div(tmpa,tmpb)), GSL_IMAG(gsl_complex_div(tmpa,tmpb))); + } else{ + printf("Beta %d is zero \n",i); + } + if(!gsl_complex_abs(tmpa)){ + printf("Alpha %d is zero \n",i); + } + } + + gsl_vector_complex_free(alph); + gsl_vector_complex_free(beta); + gsl_matrix_free(Eye); + + //B.NSEig(&Lambda, &S); + + + // V0S <- V0*S + printf("Multiplying V0*S...\n"); + + gsl_vector_complex *V = gsl_vector_complex_alloc(K); + gsl_vector_complex *s = gsl_vector_complex_alloc(K); + + printf("Evaluating retained values \n"); + int KRetained=0; + gsl_vector_complex * om = gsl_vector_alloc(1); + for(int k=0; k0.0) + { /*gsl_matrix_complex Vk(M,1,V); + gsl_matrix_complex MVk(M,1,MLBuffers[0]); + UserFunc(z, Params, &Vk, &MVk); + Residual=VecNorm(MVk.ZM, M); + */ + gsl_vector_complex_set(om,1,tmp_c); + Residual = min_sing(om,Params); + if (1) printf("Beyn: Residual(%i)=%e\n",k,Residual); + } + if (ResTol>0.0 && Residual>ResTol) continue; + + //Eigenvalues->SetEntry(KRetained, z); + gsl_vector_complex_set(Eigenvalues, KRetained, tmp_c); + gsl_matrix_complex_set_col(Eigenvectors, KRetained, V); + /*if (Eigenvectors) + { + //Eigenvectors->SetEntries(":", KRetained, V); + //Solver->Residuals->SetEntry(KRetained,Residual); + } + */ + KRetained++; + } + } + + printf("%d eigenvalues found \n",KRetained); + + gsl_matrix_free(S); + gsl_matrix_free(V0); + gsl_vector_free(Lambda); + gsl_vector_complex_free(om); + return KRetained; +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ + +int BeynSolve(BeynSolver *Solver, + BeynFunction UserFunc, void *Params, + double complex z0, double Rx, double Ry, int N) +{ + /***************************************************************/ + /* force N to be even so we can simultaneously evaluate */ + /* the integral with N/2 quadrature points */ + /***************************************************************/ + + if ( (N%2)==1 ) N++; + + /*if (Rx==Ry) + printf("Applying Beyn method with z0=%s,R=%e,N=%i...\n",z2s(z0),Rx,N); + else + printf("Applying Beyn method with z0=%s,Rx=%e,Ry=%e,N=%i...\n",z2s(z0),Rx,Ry,N); +*/ + int M = Solver->M; + int L = Solver->L; + gsl_matrix_complex *A0 = Solver->A0; + gsl_matrix_complex *A1 = Solver->A1; + gsl_matrix_complex *A0Coarse = Solver->A0Coarse; + gsl_matrix_complex *A1Coarse = Solver->A1Coarse; + gsl_matrix_complex *MInvVHat = Solver->MInvVHat; + gsl_matrix_complex *VHat = Solver->VHat; + + /***************************************************************/ + /* evaluate contour integrals by numerical quadrature to get */ + /* A0 and A1 matrices */ + /***************************************************************/ + + gsl_matrix_set_zero(A0); + gsl_matrix_set_zero(A1); + gsl_matrix_set_zero(A0Coarse); + gsl_matrix_set_zero(A1Coarse); + double DeltaTheta = 2.0*M_PI / ((double)N); + printf(" Evaluating contour integral (%i points)...\n",N); + for(int n=0; nCopy(VHat); + gsl_matrix_memcpy(MInvVHat, VHat); + + // Tän pitäis kutsua eval_WT + // Output ilmeisesti tallentuun neljänteen parametriin + + UserFunc(z0+zz, Params, VHat, MInvVHat); + + gsl_matrix_complex_scale(MInvVHat, dz); + gsl_matrix_complex_add(A0, MInvVHat); + if((n%2)==0) { + gsl_matrix_complex_add(A0Coarse, MInvVHat); + gsl_matrix_complex_add(A0Coarse, MInvVHat); + } + + gsl_matrix_complex_scale(MInvVHat, z1); + gsl_matrix_complex_add(A1, MInvVHat); + if((n%2)==0) { + gsl_matrix_complex_add(A1Coarse, MInvVHat); + gsl_matrix_complex_add(A1Coarse, MInvVHat); + } + } + + gsl_vector_complex *Eigenvalues = Solver->Eigenvalues; + gsl_vector_complex *EVErrors = Solver->EVErrors; + gsl_matrix_complex *Eigenvectors = Solver->Eigenvectors; + + int K = ProcessAMatrices(Solver, UserFunc, Params, A0, A1, z0, Eigenvalues, Eigenvectors); + //int KCoarse = ProcessAMatrices(Solver, UserFunc, Params, A0Coarse, A1Coarse, z0, EVErrors, Eigenvectors); +// Log("{K,KCoarse}={%i,%i}",K,KCoarse); + /* + for(int k=0; kN && kN; k++) + { EVErrors->ZV[k] -= Eigenvalues->ZV[k]; + EVErrors->ZV[k] = cdouble( fabs(real(EVErrors->ZV[k])), + fabs(imag(EVErrors->ZV[k])) + ); + } + +*/ + return 0; +} + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +/* +int BeynSolve(BeynSolver *Solver, + BeynFunction UserFunction, void *Params, + cdouble z0, double R, int N) +{ return BeynSolve(Solver, UserFunction, Params, z0, R, R, N); } +*/ diff --git a/qpms/libBeyn.h b/qpms/libBeyn.h new file mode 100644 index 0000000..5230bef --- /dev/null +++ b/qpms/libBeyn.h @@ -0,0 +1,95 @@ +/* Copyright (C) 2005-2011 M. T. Homer Reid + * + * This file is part of SCUFF-EM. + * + * SCUFF-EM is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * SCUFF-EM is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +/* + * libBeyn.h -- header file for libBeyn, a simple implementation of + * -- Beyn's algorithm for nonlinear eigenproblems + * -- + * -- This is packaged together with SCUFF-EM, but + * -- it is really a standalone independent entity + * -- for general-purpose use in solving nonlinear + * -- eigenproblems. + */ + + +#ifndef LIBBEYN_H +#define LIBBEYN_H + +#include +#include +#include +#include +#include + +#include + +// Needs to be changed to gsl or something +//#include + +#include +/***************************************************************/ +/* prototype for user-supplied function passed to BeynMethod. */ +/* The user's function should replace VHat with */ +/* Inverse[ M(z) ] * VHat. */ +/***************************************************************/ +typedef void (*BeynFunction)(double complex z, void *Params, gsl_matrix_complex *VHat, gsl_matrix_complex *MVHat); + +/***************************************************************/ +/***************************************************************/ +/***************************************************************/ +typedef struct BeynSolver +{ + int M; // dimension of matrices + int L; // number of columns of VHat matrix + + gsl_vector_complex *Eigenvalues, *EVErrors, *Residuals; + gsl_matrix_complex *Eigenvectors; + gsl_matrix_complex *A0, *A1, *A0Coarse, *A1Coarse, *MInvVHat; + gsl_matrix_complex *VHat; + gsl_vector_complex *Sigma; + double complex *Workspace; + + } BeynSolver; + +// constructor, destructor +BeynSolver *CreateBeynSolver(int M, int L); +void DestroyBeynSolver(BeynSolver *Solver); + +// reset the random matrix VHat used in the Beyn algorithm +// +void ReRandomize(BeynSolver *Solver, unsigned int RandSeed); + +// for both of the following routines, +// the return value is the number of eigenvalues found, +// 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, + BeynFunction UserFunction, void *Params, + double complex z0, double Rx, double Ry, int N); + +#endif From b52942a5e596ecfa730b679829960d7bc645200f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 21 Aug 2019 13:05:17 +0300 Subject: [PATCH 02/34] WIP real->complex cleanup in beyn algorithm. Former-commit-id: 2cfe5b0c87924ab04255d36e46ee3c603f25ba6d --- qpms/CMakeLists.txt | 4 +- qpms/{libBeyn.c => beyn.c} | 88 ++++++++++++++++++-------------------- qpms/{libBeyn.h => beyn.h} | 8 ++-- 3 files changed, 49 insertions(+), 51 deletions(-) rename qpms/{libBeyn.c => beyn.c} (88%) rename qpms/{libBeyn.h => beyn.h} (97%) diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index 3dffad0..3ec1a24 100644 --- a/qpms/CMakeLists.txt +++ b/qpms/CMakeLists.txt @@ -13,7 +13,9 @@ include_directories(${DIRS}) add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c ewald.c ewaldsf.c pointgroups.c latticegens.c lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c - bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c) + bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c + beyn.c + ) use_c99() set(LIBS ${LIBS} ${GSL_LIBRARIES} ${GSLCBLAS_LIBRARIES}) diff --git a/qpms/libBeyn.c b/qpms/beyn.c similarity index 88% rename from qpms/libBeyn.c rename to qpms/beyn.c index 0ad4698..6b4ef58 100644 --- a/qpms/libBeyn.c +++ b/qpms/beyn.c @@ -31,10 +31,6 @@ #include - -//#include -//#include - // Maybe GSL works? #include #include @@ -42,15 +38,13 @@ #include #include -#include - -//#define II cdouble(0.0,1.0) +#include // Uniformly random number between -2 and 2 -double complex zrandN(){ +gsl_complex zrandN(){ double a = (rand()*4.0/RAND_MAX) - 2; double b = (rand()*4.0/RAND_MAX) - 2; - return a + b*I; + return gsl_complex_rect(a, b); } /***************************************************************/ @@ -79,7 +73,7 @@ BeynSolver *CreateBeynSolver(int M, int L) Solver->A1Coarse = 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_complex_calloc(MLMin); + Solver->Sigma = gsl_vector_calloc(MLMin); ReRandomize(Solver,(unsigned)time(NULL)); // internal workspace: need storage for 2 MxL matrices @@ -109,7 +103,7 @@ void DestroyBeynSolver(BeynSolver *Solver) gsl_matrix_complex_free(Solver->A0Coarse); gsl_matrix_complex_free(Solver->A1Coarse); gsl_matrix_complex_free(Solver->MInvVHat); - gsl_vector_complex_free(Solver->Sigma); + gsl_vector_free(Solver->Sigma); gsl_matrix_complex_free(Solver->VHat); free(Solver->Workspace); @@ -129,7 +123,7 @@ void ReRandomize(BeynSolver *Solver, unsigned int RandSeed) gsl_matrix_complex *VHat=Solver->VHat; for(int nr=0; nrsize1; nr++) for(int nc=0; ncsize2; nc++) - gsl_matrix_set(VHat,nr,nc,zrandN()); + gsl_matrix_complex_set(VHat,nr,nc,zrandN()); } @@ -147,7 +141,7 @@ int ProcessAMatrices(BeynSolver *Solver, BeynFunction UserFunc, { int M = Solver->M; int L = Solver->L; - gsl_vector_complex *Sigma = Solver->Sigma; + gsl_vector *Sigma = Solver->Sigma; int Verbose = 0;//CheckEnv("SCUFF_BEYN_VERBOSE"); @@ -157,13 +151,14 @@ int ProcessAMatrices(BeynSolver *Solver, BeynFunction UserFunc, // A0 -> V0Full * Sigma * W0TFull' printf(" Beyn: computing SVD...\n"); gsl_matrix_complex* V0Full = gsl_matrix_complex_calloc(M,L); - gsl_matrix_memcpy(V0Full,A0); + gsl_matrix_complex_memcpy(V0Full,A0); gsl_matrix_complex* W0TFull = gsl_matrix_complex_calloc(L,L); //A0->SVD(Sigma, &V0Full, &W0TFull); gsl_vector_complex *work = gsl_vector_complex_alloc(M); - - gsl_linalg_SV_decomp(V0Full, W0TFull, Sigma, work); + + // FIXME not supported by GSL; use LAPACKE_zgesdd + gsl_linalg_complex_SV_decomp(V0Full, W0TFull, Sigma, work); // compute effective rank K (number of eigenvalue candidates) @@ -188,17 +183,17 @@ int ProcessAMatrices(BeynSolver *Solver, BeynFunction UserFunc, gsl_vector_complex *TempL = gsl_vector_complex_calloc(L); for(int k=0; k0.0 && Residual>ResTol) continue; @@ -322,9 +318,9 @@ int ProcessAMatrices(BeynSolver *Solver, BeynFunction UserFunc, printf("%d eigenvalues found \n",KRetained); - gsl_matrix_free(S); - gsl_matrix_free(V0); - gsl_vector_free(Lambda); + gsl_matrix_complex_free(S); + gsl_matrix_complex_free(V0); + gsl_vector_complex_free(Lambda); gsl_vector_complex_free(om); return KRetained; } @@ -349,8 +345,8 @@ int BeynSolve(BeynSolver *Solver, else printf("Applying Beyn method with z0=%s,Rx=%e,Ry=%e,N=%i...\n",z2s(z0),Rx,Ry,N); */ - int M = Solver->M; - int L = Solver->L; + //int M = Solver->M; + //int L = Solver->L; gsl_matrix_complex *A0 = Solver->A0; gsl_matrix_complex *A1 = Solver->A1; gsl_matrix_complex *A0Coarse = Solver->A0Coarse; @@ -363,10 +359,10 @@ int BeynSolve(BeynSolver *Solver, /* A0 and A1 matrices */ /***************************************************************/ - gsl_matrix_set_zero(A0); - gsl_matrix_set_zero(A1); - gsl_matrix_set_zero(A0Coarse); - gsl_matrix_set_zero(A1Coarse); + gsl_matrix_complex_set_zero(A0); + gsl_matrix_complex_set_zero(A1); + gsl_matrix_complex_set_zero(A0Coarse); + gsl_matrix_complex_set_zero(A1Coarse); double DeltaTheta = 2.0*M_PI / ((double)N); printf(" Evaluating contour integral (%i points)...\n",N); for(int n=0; nCopy(VHat); - gsl_matrix_memcpy(MInvVHat, VHat); + gsl_matrix_complex_memcpy(MInvVHat, VHat); // Tän pitäis kutsua eval_WT // Output ilmeisesti tallentuun neljänteen parametriin @@ -402,7 +398,7 @@ int BeynSolve(BeynSolver *Solver, } gsl_vector_complex *Eigenvalues = Solver->Eigenvalues; - gsl_vector_complex *EVErrors = Solver->EVErrors; + //gsl_vector_complex *EVErrors = Solver->EVErrors; gsl_matrix_complex *Eigenvectors = Solver->Eigenvectors; int K = ProcessAMatrices(Solver, UserFunc, Params, A0, A1, z0, Eigenvalues, Eigenvectors); diff --git a/qpms/libBeyn.h b/qpms/beyn.h similarity index 97% rename from qpms/libBeyn.h rename to qpms/beyn.h index 5230bef..98404fe 100644 --- a/qpms/libBeyn.h +++ b/qpms/beyn.h @@ -28,8 +28,8 @@ */ -#ifndef LIBBEYN_H -#define LIBBEYN_H +#ifndef BEYN_H +#define BEYN_H #include #include @@ -62,7 +62,7 @@ typedef struct BeynSolver gsl_matrix_complex *Eigenvectors; gsl_matrix_complex *A0, *A1, *A0Coarse, *A1Coarse, *MInvVHat; gsl_matrix_complex *VHat; - gsl_vector_complex *Sigma; + gsl_vector *Sigma; double complex *Workspace; } BeynSolver; @@ -92,4 +92,4 @@ int BeynSolve(BeynSolver *Solver, BeynFunction UserFunction, void *Params, double complex z0, double Rx, double Ry, int N); -#endif +#endif // BEYN_H From 6e875d65aefe1413252bae1e87229b1c1ba9fdd5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 22 Aug 2019 17:02:53 +0300 Subject: [PATCH 03/34] WIP Dealing with the Beyn clusterfuck (compiles now). TODO inverse M -matrix Former-commit-id: eb2a37128c04c10406dc65eca7d47152b4d93db9 --- qpms/beyn.c | 195 +++++++++++++++++++++++++--------------------------- qpms/beyn.h | 29 +++----- 2 files changed, 104 insertions(+), 120 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 6b4ef58..14b05eb 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -1,6 +1,6 @@ -/* Copyright (C) 2005-2011 M. T. Homer Reid - * - * This file is part of SCUFF-EM. +/* + * This file was originally part of SCUFF-EM by M. T. Homer Reid. + * Modified by Kristian Arjas and Marek Nečada to work without libhmat and libhrutil. * * SCUFF-EM is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -17,19 +17,15 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ -/* - * libBeyn.cc -- implementation of Beyn's method for - * -- nonlinear eigenvalue problems - * - * Homer Reid -- 6/2016 - * - */ +#define STATIC_ASSERT(COND,MSG) typedef char static_assertion_##MSG[(COND)?1:-1] + +#include +#include #include #include #include #include - -#include +#include "qpms_error.h" // Maybe GSL works? #include @@ -40,6 +36,8 @@ #include +STATIC_ASSERT((sizeof(lapack_complex_double) == sizeof(gsl_complex)), lapack_and_gsl_complex_must_be_consistent); + // Uniformly random number between -2 and 2 gsl_complex zrandN(){ double a = (rand()*4.0/RAND_MAX) - 2; @@ -56,14 +54,17 @@ BeynSolver *CreateBeynSolver(int M, int L) Solver->M = M; Solver->L = L; + QPMS_ENSURE(L <= M, "We expect L <= M, but we got L = %d, M = %d ", L, M); +#if 0 int MLMax = (M>L) ? M : L; +#endif int MLMin = (MEigenvalues = gsl_vector_complex_calloc(L); Solver->EVErrors = gsl_vector_complex_calloc(L); - Solver->Residuals = gsl_vector_complex_calloc(L); + Solver->Residuals = gsl_vector_calloc(L); Solver->Eigenvectors = gsl_matrix_complex_calloc(M, L); // storage for singular values, random VHat matrix, etc. used in algorithm @@ -76,15 +77,14 @@ BeynSolver *CreateBeynSolver(int M, int L) Solver->Sigma = gsl_vector_calloc(MLMin); ReRandomize(Solver,(unsigned)time(NULL)); +#if 0 // internal workspace: need storage for 2 MxL matrices // plus 3 LxL matrices #define MLBUFFERS 2 #define LLBUFFERS 3 size_t ML = MLMax*L, LL = L*L; - size_t WorkspaceSize = (MLBUFFERS*ML + LLBUFFERS*LL)*sizeof(double complex); +#endif - Solver->Workspace = (double complex*)malloc( WorkspaceSize ); - return Solver; } @@ -104,10 +104,9 @@ void DestroyBeynSolver(BeynSolver *Solver) gsl_matrix_complex_free(Solver->A1Coarse); gsl_matrix_complex_free(Solver->MInvVHat); gsl_vector_free(Solver->Sigma); + gsl_vector_free(Solver->Residuals); gsl_matrix_complex_free(Solver->VHat); - free(Solver->Workspace); - free(Solver); } @@ -134,7 +133,7 @@ void ReRandomize(BeynSolver *Solver, unsigned int RandSeed) /* eigenvalues and eigenvectors */ /***************************************************************/ -int ProcessAMatrices(BeynSolver *Solver, BeynFunction UserFunc, +int ProcessAMatrices(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) @@ -155,15 +154,27 @@ int ProcessAMatrices(BeynSolver *Solver, BeynFunction UserFunc, gsl_matrix_complex* W0TFull = gsl_matrix_complex_calloc(L,L); //A0->SVD(Sigma, &V0Full, &W0TFull); - gsl_vector_complex *work = gsl_vector_complex_alloc(M); + QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride); // FIXME not supported by GSL; use LAPACKE_zgesdd - gsl_linalg_complex_SV_decomp(V0Full, W0TFull, Sigma, work); + //gsl_linalg_complex_SV_decomp(V0Full, W0TFull, Sigma, work); + 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*/, + V0Full->size1 /* m, number of rows */, + V0Full->size2 /* n, number of columns */, + (lapack_complex_double *)(V0Full->data) /*a*/, + V0Full->tda /*lda*/, + Sigma->data /*s*/, + NULL /*u; not used*/, + 0 /*ldu; not used*/, + (lapack_complex_double *)W0TFull->data /*vt*/, + W0TFull->tda /*ldvt*/ + )); // compute effective rank K (number of eigenvalue candidates) int K=0; - for(int k=0; ksize; k++) + for(int k=0; ksize /* this is L, actually */; k++) { if (Verbose) printf("Beyn: SV(%i)=%e",k,gsl_vector_get(Sigma, k)); if (gsl_vector_get(Sigma, k) > RankTol ) K++; @@ -191,7 +202,6 @@ int ProcessAMatrices(BeynSolver *Solver, BeynFunction UserFunc, gsl_matrix_complex_free(V0Full); gsl_matrix_complex_free(W0TFull); - gsl_vector_complex_free(work); gsl_vector_complex_free(TempM); gsl_vector_complex_free(TempL); @@ -202,8 +212,8 @@ int ProcessAMatrices(BeynSolver *Solver, BeynFunction UserFunc, printf(" Multiplying V0*A1->TM...\n"); //V0.Multiply(A1, &TM2, "--transA C"); // TM2 <- V0' * A1 - gsl_complex one = gsl_complex_rect(1,0); - gsl_complex zero = gsl_complex_rect(0,0); + 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); @@ -234,94 +244,72 @@ int ProcessAMatrices(BeynSolver *Solver, BeynFunction UserFunc, // B -> S*Lambda*S' printf(" Eigensolving (%i,%i)\n",K,K); - gsl_vector_complex *Lambda = gsl_vector_complex_calloc(K); // Eigenvalues - gsl_matrix_complex *S = gsl_matrix_complex_calloc(K,K); // Eigenvectors - gsl_matrix_complex *Eye = gsl_matrix_complex_alloc(K,K); - - gsl_vector_complex *alph = gsl_vector_complex_calloc(K); - gsl_vector_complex *beta = gsl_vector_complex_calloc(K); - - gsl_matrix_complex_set_identity(Eye); + gsl_vector_complex *Lambda = gsl_vector_complex_alloc(K); // Eigenvalues + gsl_matrix_complex *S = gsl_matrix_complex_alloc(K,K); // Eigenvectors // FIXME general complex eigensystems not supported by GSL (use LAPACKE_zgee?) - gsl_eigen_genv_workspace * W = gsl_eigen_genv_alloc(K); - gsl_eigen_genv(B, Eye, alph, beta, S,W); - gsl_eigen_genv_free(W); + //gsl_eigen_genv_workspace * W = gsl_eigen_genv_alloc(K); + //gsl_eigen_genv(B, Eye, alph, beta, S,W); + //gsl_eigen_genv_free(W); + 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_SUCCESS(LAPACKE_zgeev( + LAPACK_ROW_MAJOR, + 'N' /* jobvl; don't compute left eigenvectors */, + 'V' /* jobvr; do compute right eigenvectors */, + K /* n */, + (lapack_complex_double *)(S->data) /* a */, + S->tda /* lda */, + (lapack_complex_double *) Lambda->data /* w */, + NULL /* vl */, + 0 /* ldvl */, + (lapack_complex_double *)(B->data)/* vr */, + B->tda/* ldvr */ + )); - gsl_complex tmpa; - gsl_complex tmpb; - for(int i = 0; i < K; i++){ - tmpb = gsl_vector_complex_get(beta,i); - tmpa = gsl_vector_complex_get(alph,i); - if(gsl_complex_abs(tmpb)){ - gsl_vector_complex_set(Lambda, i, gsl_complex_div(tmpa,tmpb)); - printf("Eigenvalue %e + %e i found\n",GSL_REAL(gsl_complex_div(tmpa,tmpb)), GSL_IMAG(gsl_complex_div(tmpa,tmpb))); - } else{ - printf("Beta %d is zero \n",i); - } - if(!gsl_complex_abs(tmpa)){ - printf("Alpha %d is zero \n",i); - } - } - gsl_vector_complex_free(alph); - gsl_vector_complex_free(beta); - gsl_matrix_complex_free(Eye); - //B.NSEig(&Lambda, &S); + //B.NSEig(&Lambda, &S); - // 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_vector_complex *V = gsl_vector_complex_alloc(K); - gsl_vector_complex *s = gsl_vector_complex_alloc(K); - printf("Evaluating retained values \n"); - int KRetained=0; - gsl_vector_complex * om = gsl_vector_complex_alloc(1); - for(int k=0; k0.0) - { /*gsl_matrix_complex Vk(M,1,V); - gsl_matrix_complex MVk(M,1,MLBuffers[0]); - UserFunc(z, Params, &Vk, &MVk); - Residual=VecNorm(MVk.ZM, M); - */ - gsl_vector_complex_set(om,1,tmp_c); - Residual = min_sing(om,Params); // in unitcell.c - if (1) printf("Beyn: Residual(%i)=%e\n",k,Residual); - } - if (ResTol>0.0 && Residual>ResTol) continue; - - //Eigenvalues->SetEntry(KRetained, z); - gsl_vector_complex_set(Eigenvalues, KRetained, tmp_c); - gsl_matrix_complex_set_col(Eigenvectors, KRetained, V); - /*if (Eigenvectors) - { - //Eigenvectors->SetEntries(":", KRetained, V); - //Solver->Residuals->SetEntry(KRetained,Residual); - } - */ - KRetained++; + int KRetained = 0; + gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(M,L); + gsl_vector_complex *MVk = gsl_vector_complex_alloc(M); + 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); + + double residual = 0; + if(ResTol > 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); + if (Verbose) printf("Beyn: Residual(%i)=%e\n",k,residual); } - } + if (ResTol > 0 && residual > ResTol) continue; - printf("%d eigenvalues found \n",KRetained); + gsl_vector_complex_set(Eigenvalues, KRetained, zgsl); + if(Eigenvectors) { + gsl_matrix_complex_set_col(Eigenvectors, KRetained, &(Vk.vector)); + 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_matrix_complex_free(V0); gsl_vector_complex_free(Lambda); - gsl_vector_complex_free(om); + return KRetained; } @@ -329,8 +317,8 @@ int ProcessAMatrices(BeynSolver *Solver, BeynFunction UserFunc, /***************************************************************/ /***************************************************************/ -int BeynSolve(BeynSolver *Solver, - BeynFunction UserFunc, void *Params, +int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, + beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *Params, double complex z0, double Rx, double Ry, int N) { /***************************************************************/ @@ -380,7 +368,12 @@ int BeynSolve(BeynSolver *Solver, // Tän pitäis kutsua eval_WT // Output ilmeisesti tallentuun neljänteen parametriin - UserFunc(z0+zz, Params, VHat, MInvVHat); + if(M_inv_Vhat_function) { + QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z0+zz, Params)); + } else { + QPMS_NOT_IMPLEMENTED("TODO"); + } + //UserFunc(z0+zz, Params, VHat, MInvVHat); gsl_matrix_complex_scale(MInvVHat, dz); gsl_matrix_complex_add(A0, MInvVHat); @@ -401,7 +394,7 @@ int BeynSolve(BeynSolver *Solver, //gsl_vector_complex *EVErrors = Solver->EVErrors; gsl_matrix_complex *Eigenvectors = Solver->Eigenvectors; - int K = ProcessAMatrices(Solver, UserFunc, Params, A0, A1, z0, Eigenvalues, Eigenvectors); + int K = ProcessAMatrices(Solver, M_function, Params, A0, A1, z0, Eigenvalues, Eigenvectors); //int KCoarse = ProcessAMatrices(Solver, UserFunc, Params, A0Coarse, A1Coarse, z0, EVErrors, Eigenvectors); // Log("{K,KCoarse}={%i,%i}",K,KCoarse); /* diff --git a/qpms/beyn.h b/qpms/beyn.h index 98404fe..61ec3b4 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -31,24 +31,15 @@ #ifndef BEYN_H #define BEYN_H -#include -#include -#include -#include -#include - #include - -// Needs to be changed to gsl or something -//#include - #include -/***************************************************************/ -/* prototype for user-supplied function passed to BeynMethod. */ -/* The user's function should replace VHat with */ -/* Inverse[ M(z) ] * VHat. */ -/***************************************************************/ -typedef void (*BeynFunction)(double complex z, void *Params, gsl_matrix_complex *VHat, gsl_matrix_complex *MVHat); + +/// User-supplied function that provides the operator M(z) whose "roots" are to be found. +typedef int (*beyn_function_M_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$. +typedef int (*beyn_function_M_inv_Vhat_t)(gsl_matrix_complex *target_M_inv_Vhat, + const gsl_matrix_complex *Vhat, complex double z, void *params); /***************************************************************/ /***************************************************************/ @@ -58,11 +49,11 @@ typedef struct BeynSolver int M; // dimension of matrices int L; // number of columns of VHat matrix - gsl_vector_complex *Eigenvalues, *EVErrors, *Residuals; + gsl_vector_complex *Eigenvalues, *EVErrors; gsl_matrix_complex *Eigenvectors; gsl_matrix_complex *A0, *A1, *A0Coarse, *A1Coarse, *MInvVHat; gsl_matrix_complex *VHat; - gsl_vector *Sigma; + gsl_vector *Sigma, *Residuals; double complex *Workspace; } BeynSolver; @@ -89,7 +80,7 @@ void ReRandomize(BeynSolver *Solver, unsigned int RandSeed); // Beyn method for elliptical contour of horizontal, vertical // radii Rx, Ry, centered at z0, using N quadrature points int BeynSolve(BeynSolver *Solver, - BeynFunction UserFunction, void *Params, + beyn_function_M_t M_function, beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *params, double complex z0, double Rx, double Ry, int N); #endif // BEYN_H From cf9892279dcf5d37ae2b230aab070c4256499022 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 23 Aug 2019 09:34:19 +0300 Subject: [PATCH 04/34] Beyn ready for testing Former-commit-id: b058f06d1f4d80e366c2375764b0c7a0f252568b --- qpms/beyn.c | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 14b05eb..fbe9763 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -280,7 +280,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, int KRetained = 0; - gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(M,L); + gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(M,M); gsl_vector_complex *MVk = gsl_vector_complex_alloc(M); 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)); @@ -333,8 +333,8 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, else printf("Applying Beyn method with z0=%s,Rx=%e,Ry=%e,N=%i...\n",z2s(z0),Rx,Ry,N); */ - //int M = Solver->M; - //int L = Solver->L; + const int M = Solver->M; + const int L = Solver->L; gsl_matrix_complex *A0 = Solver->A0; gsl_matrix_complex *A1 = Solver->A1; gsl_matrix_complex *A0Coarse = Solver->A0Coarse; @@ -363,6 +363,7 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, double complex zz = Rx*CT + Ry*ST*I; //MInvVHat->Copy(VHat); + // Mitä varten tämä kopiointi on? gsl_matrix_complex_memcpy(MInvVHat, VHat); // Tän pitäis kutsua eval_WT @@ -371,7 +372,17 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, if(M_inv_Vhat_function) { QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z0+zz, Params)); } else { - QPMS_NOT_IMPLEMENTED("TODO"); + lapack_int *pivot; + gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(M,M); + QPMS_ENSURE_SUCCESS(M_function(Mmat, z0+zz, Params)); + QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * M); + 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_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 *)VHat->data /*b*/, VHat->tda/*ldb*/)); + free(pivot); + gsl_matrix_complex_free(Mmat); } //UserFunc(z0+zz, Params, VHat, MInvVHat); From 142d97d484b724c0d8beef6672363b9da8292f64 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 23 Aug 2019 11:16:46 +0300 Subject: [PATCH 05/34] Beyn's method test (fails) Former-commit-id: 535a94df45d1ddd6b23b6c94a8167ca969389099 --- qpms/beyn.c | 2 +- qpms/beyn.h | 5 +++++ tests/CMakeLists.txt | 6 +++++- tests/tbeyn.c | 39 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 50 insertions(+), 2 deletions(-) create mode 100644 tests/tbeyn.c diff --git a/qpms/beyn.c b/qpms/beyn.c index fbe9763..3dce3a8 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -417,7 +417,7 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, } */ - return 0; + return K; } /***************************************************************/ diff --git a/qpms/beyn.h b/qpms/beyn.h index 61ec3b4..6b9f69f 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -33,6 +33,7 @@ #include #include +#include /// User-supplied function that provides the operator M(z) whose "roots" are to be found. typedef int (*beyn_function_M_t)(gsl_matrix_complex *target_M, complex double z, void *params); @@ -83,4 +84,8 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *params, double complex z0, double Rx, double Ry, int N); + +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 6e6797b..eb0b17e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -15,7 +15,11 @@ add_executable( test_scalar_ewald32 test_scalar_ewald32.c ) target_link_libraries( test_scalar_ewald32 qpms gsl lapacke amos m ) target_include_directories( test_scalar_ewald32 PRIVATE .. ) -add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array ) +add_executable( tbeyn tbeyn.c ) +target_link_libraries( tbeyn qpms gsl lapacke amos m ) +target_include_directories( tbeyn 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.c b/tests/tbeyn.c new file mode 100644 index 0000000..a78113f --- /dev/null +++ b/tests/tbeyn.c @@ -0,0 +1,39 @@ +#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( -m - z / (6*m) ); + + gsl_matrix_complex_set_zero(target); + for (int i = 0; i < m; ++i) { + gsl_matrix_complex_set(target, i, i, (gsl_complex) d); + if(i > 0) gsl_matrix_complex_set(target, i, i-1, (gsl_complex) od); + if(i < m - 1) gsl_matrix_complex_set(target, i, i+1, (gsl_complex) od); + } + + return 0; +} + +int main() { + complex double z0 = 150+2*I; + double Rx = 148, Ry = 148; + int L = 10, N = 50, dim = 400; + BeynSolver * solver = CreateBeynSolver(dim, L); + + int K = BeynSolve(solver, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, + z0, Rx, Ry, N); + printf("Found %d eigenvalues:\n", K); + for (int i = 0; i < K; ++i) { + gsl_complex eig = gsl_vector_complex_get(solver->Eigenvalues, i); + printf("%d: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); + } + DestroyBeynSolver(solver); + return 0; +} + + From e898fd5ad89cf1005bb1e7ef3f2d0939b951d5d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 23 Aug 2019 12:20:51 +0300 Subject: [PATCH 06/34] tbeyn does not crash anymore but the result is wrong Former-commit-id: 3288074674f1a3f61b43d905a96c865fa32faccd --- qpms/beyn.c | 28 ++++++++++++---------------- 1 file changed, 12 insertions(+), 16 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 3dce3a8..f889109 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -156,8 +156,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, //A0->SVD(Sigma, &V0Full, &W0TFull); QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride); - // FIXME not supported by GSL; use LAPACKE_zgesdd - //gsl_linalg_complex_SV_decomp(V0Full, W0TFull, Sigma, work); + QPMS_ENSURE(V0Full->size1 >= V0Full->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*/, V0Full->size1 /* m, number of rows */, @@ -166,7 +165,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, V0Full->tda /*lda*/, Sigma->data /*s*/, NULL /*u; not used*/, - 0 /*ldu; not used*/, + M /*ldu; should not be really used but lapacke requires it for some obscure reason*/, (lapack_complex_double *)W0TFull->data /*vt*/, W0TFull->tda /*ldvt*/ )); @@ -186,24 +185,21 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, } + // 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_vector_complex *TempM = gsl_vector_complex_calloc(M); - gsl_vector_complex *TempL = gsl_vector_complex_calloc(L); - for(int k=0; ktda /* lda */, (lapack_complex_double *) Lambda->data /* w */, NULL /* vl */, - 0 /* ldvl */, + M /* ldvl, not used by for some reason needed */, (lapack_complex_double *)(B->data)/* vr */, B->tda/* ldvr */ )); From 8666657fcc30d704f25b3ab855fd24d2143ad6a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 23 Aug 2019 13:48:28 +0300 Subject: [PATCH 07/34] Beyn fix memory leaks and swapped arguments in zgeev. Former-commit-id: a63ba1bbd2ce50f664a419e4742fcead2fb355e1 --- qpms/beyn.c | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index f889109..93023fe 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -212,6 +212,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, const gsl_complex zero = gsl_complex_rect(0,0); gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, V0, A1, zero, TM2); + printf(" Multiplying TM*W0T->B...\n"); //TM2.Multiply(&W0T, &B, "--transB C"); // B <- TM2 * W0 @@ -254,16 +255,16 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, 'N' /* jobvl; don't compute left eigenvectors */, 'V' /* jobvr; do compute right eigenvectors */, K /* n */, - (lapack_complex_double *)(S->data) /* a */, - S->tda /* lda */, + (lapack_complex_double *)(B->data) /* a */, + B->tda /* lda */, (lapack_complex_double *) Lambda->data /* w */, NULL /* vl */, M /* ldvl, not used by for some reason needed */, - (lapack_complex_double *)(B->data)/* vr */, - B->tda/* ldvr */ + (lapack_complex_double *)(S->data)/* vr */, + S->tda/* ldvr */ )); - + gsl_matrix_complex_free(B); //B.NSEig(&Lambda, &S); @@ -273,6 +274,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, QPMS_ENSURE_SUCCESS(gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, V0, S, zero, V0S)); + gsl_matrix_complex_free(V0); int KRetained = 0; From 93d34c9830b22ef8625552dd959ef00c91781b28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 23 Aug 2019 14:30:39 +0300 Subject: [PATCH 08/34] Probably wrong target to zgetrs Former-commit-id: 4d22626121778d77d61f4fa9f600fc87b91c1bec --- qpms/beyn.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 93023fe..a6b9f48 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -143,7 +143,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, gsl_vector *Sigma = Solver->Sigma; - int Verbose = 0;//CheckEnv("SCUFF_BEYN_VERBOSE"); + int Verbose = 1;//CheckEnv("SCUFF_BEYN_VERBOSE"); double RankTol=1.0e-4; //CheckEnv("SCUFF_BEYN_RANK_TOL",&RankTol); double ResTol=0.0; // CheckEnv("SCUFF_BEYN_RES_TOL",&ResTol); @@ -378,7 +378,7 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, M /*m*/, M /*n*/,(lapack_complex_double *) Mmat->data /*a*/, Mmat->tda /*lda*/, pivot /*ipiv*/)); 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 *)VHat->data /*b*/, VHat->tda/*ldb*/)); + (lapack_complex_double *)MInvVHat->data /*b*/, MInvVHat->tda/*ldb*/)); free(pivot); gsl_matrix_complex_free(Mmat); } From 51fa6f1dd6a7bc23e6a42bad2ee60c8cf30fe36b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 23 Aug 2019 15:05:23 +0300 Subject: [PATCH 09/34] Try to reproduce Reid's RNG Former-commit-id: b607cb40b9ac77970e0c4250cf033d22554de815 --- qpms/beyn.c | 33 ++++++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 9 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index a6b9f48..ebf6f05 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -19,6 +19,9 @@ #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 @@ -38,12 +41,26 @@ STATIC_ASSERT((sizeof(lapack_complex_double) == sizeof(gsl_complex)), lapack_and_gsl_complex_must_be_consistent); +double randU(double A, double B) { return A + (B-A) * random() * (1. / RAND_MAX); } +double randN(double Sigma, double Mu) { + double u1 = randU(0, 1); + double u2 = randU(0, 1); + return Mu + Sigma*sqrt(-2*log(u1))*cos(2.*M_PI*u2); +} + +#if 0 // Uniformly random number between -2 and 2 gsl_complex zrandN(){ double a = (rand()*4.0/RAND_MAX) - 2; double b = (rand()*4.0/RAND_MAX) - 2; return gsl_complex_rect(a, b); } +#endif + +complex double zrandN(double sigma, double mu){ + return randN(sigma, mu) + I*randN(sigma, mu); +} + /***************************************************************/ /***************************************************************/ @@ -122,7 +139,7 @@ void ReRandomize(BeynSolver *Solver, unsigned int 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,zrandN()); + gsl_matrix_complex_set(VHat,nr,nc,cs2g(zrandN(1, 0))); } @@ -355,10 +372,8 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, { double Theta = ((double)n)*DeltaTheta; double CT = cos(Theta), ST=sin(Theta); - gsl_complex z1 = gsl_complex_rect(Rx*CT, Ry*ST); - gsl_complex dz = gsl_complex_rect(Ry*CT/((double)N),(Rx*ST/((double)N))); - - double complex zz = Rx*CT + Ry*ST*I; + complex double z1 = Rx*CT + I*Ry*ST; + complex double dz = (I*Rx*ST + Ry*CT) / N; //MInvVHat->Copy(VHat); // Mitä varten tämä kopiointi on? @@ -368,11 +383,11 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, // Output ilmeisesti tallentuun neljänteen parametriin if(M_inv_Vhat_function) { - QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z0+zz, Params)); + QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z0+z1, Params)); } else { lapack_int *pivot; gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(M,M); - QPMS_ENSURE_SUCCESS(M_function(Mmat, z0+zz, Params)); + QPMS_ENSURE_SUCCESS(M_function(Mmat, z0+z1, Params)); QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * M); QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, M /*m*/, M /*n*/,(lapack_complex_double *) Mmat->data /*a*/, Mmat->tda /*lda*/, pivot /*ipiv*/)); @@ -384,14 +399,14 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, } //UserFunc(z0+zz, Params, VHat, MInvVHat); - gsl_matrix_complex_scale(MInvVHat, dz); + gsl_matrix_complex_scale(MInvVHat, cs2g(dz)); gsl_matrix_complex_add(A0, MInvVHat); if((n%2)==0) { gsl_matrix_complex_add(A0Coarse, MInvVHat); gsl_matrix_complex_add(A0Coarse, MInvVHat); } - gsl_matrix_complex_scale(MInvVHat, z1); + gsl_matrix_complex_scale(MInvVHat, cs2g(z1)); gsl_matrix_complex_add(A1, MInvVHat); if((n%2)==0) { gsl_matrix_complex_add(A1Coarse, MInvVHat); From c0c44c11b6fbe0e0e02344f9f64b867e3922d76b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 23 Aug 2019 20:09:21 +0300 Subject: [PATCH 10/34] Add some excessive checks for matrix inversion results Former-commit-id: 229cdff48f8103cbadf00bbc92c11fa19e51a4b0 --- qpms/beyn.c | 38 +++++++++++++++++++++++++++++++++----- 1 file changed, 33 insertions(+), 5 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index ebf6f05..357fe71 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -166,10 +166,10 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, // A0 -> V0Full * Sigma * W0TFull' printf(" Beyn: computing SVD...\n"); - gsl_matrix_complex* V0Full = gsl_matrix_complex_calloc(M,L); + gsl_matrix_complex* V0Full = gsl_matrix_complex_alloc(M,L); gsl_matrix_complex_memcpy(V0Full,A0); - gsl_matrix_complex* W0TFull = gsl_matrix_complex_calloc(L,L); + gsl_matrix_complex* W0TFull = gsl_matrix_complex_alloc(L,L); //A0->SVD(Sigma, &V0Full, &W0TFull); QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride); @@ -368,8 +368,7 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, gsl_matrix_complex_set_zero(A1Coarse); double DeltaTheta = 2.0*M_PI / ((double)N); printf(" Evaluating contour integral (%i points)...\n",N); - for(int n=0; ndata /*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*/)); + +#ifndef NDEBUG + // Check the inversion result + complex double minvhat_check[M][L]; + for(int i = 0; i < M; ++i) for(int j = 0; j < L; ++j) + minvhat_check[i][j] = cg2s(gsl_matrix_complex_get(MInvVHat, i, j)); + complex double vhat_recon[M][L]; + + cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, + M, L, M, &one, Mmat_check[0], M, minvhat_check[0], L, &zero, vhat_recon[0], L); + + for(int i = 0; i < M; ++i) for(int j = 0; j < L; ++j) + if (cabs(vhat_recon[i][j] - vhat_copy[i][j]) > 1e-8*(1+cabs(vhat_recon[i][j] + vhat_copy[i][j]))) + QPMS_WTF; +#endif + + free(pivot); gsl_matrix_complex_free(Mmat); } @@ -412,7 +440,7 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, gsl_matrix_complex_add(A1Coarse, MInvVHat); gsl_matrix_complex_add(A1Coarse, MInvVHat); } - } + } gsl_vector_complex *Eigenvalues = Solver->Eigenvalues; //gsl_vector_complex *EVErrors = Solver->EVErrors; From 00e3a9ce09b00dffd3b9f83151afb6ba9dc5ba7c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 27 Aug 2019 19:06:37 +0300 Subject: [PATCH 11/34] Fix test matrix definition Former-commit-id: f5874c5cf632b9fe5c561154370a3460fc1c7b46 --- tests/tbeyn.c | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/tbeyn.c b/tests/tbeyn.c index a78113f..cfd1703 100644 --- a/tests/tbeyn.c +++ b/tests/tbeyn.c @@ -11,10 +11,11 @@ int M_function(gsl_matrix_complex *target, complex double z, void *no_params) { gsl_matrix_complex_set_zero(target); for (int i = 0; i < m; ++i) { - gsl_matrix_complex_set(target, i, i, (gsl_complex) d); - if(i > 0) gsl_matrix_complex_set(target, i, i-1, (gsl_complex) od); - if(i < m - 1) gsl_matrix_complex_set(target, i, i+1, (gsl_complex) od); + 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; } @@ -24,6 +25,7 @@ int main() { double Rx = 148, Ry = 148; int L = 10, N = 50, dim = 400; BeynSolver * solver = CreateBeynSolver(dim, L); + ReRandomize(solver, 666); int K = BeynSolve(solver, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, z0, Rx, Ry, N); From 5471367aad5fce073b847975684039443a2258eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 27 Aug 2019 20:10:28 +0300 Subject: [PATCH 12/34] Finally reproduce Reid's implementation. Former-commit-id: 215088edddd93c79b4a4ad3ee9836595bcc69167 --- qpms/beyn.c | 34 ++++------------------------------ 1 file changed, 4 insertions(+), 30 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 357fe71..9024a21 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -191,7 +191,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, // 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(%i)=%e",k,gsl_vector_get(Sigma, k)); + { if (Verbose) printf("Beyn: SV(%i)=%e\n",k,gsl_vector_get(Sigma, k)); if (gsl_vector_get(Sigma, k) > RankTol ) K++; } @@ -234,7 +234,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, printf(" Multiplying TM*W0T->B...\n"); //TM2.Multiply(&W0T, &B, "--transB C"); // B <- TM2 * W0 - gsl_blas_zgemm(CblasNoTrans, CblasTrans, one, + gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, TM2, W0T, zero, B); gsl_matrix_complex_free(W0T); @@ -380,46 +380,20 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, // Tän pitäis kutsua eval_WT // Output ilmeisesti tallentuun neljänteen parametriin -#ifndef NDEBUG - complex double vhat_copy[M][L]; //dbg - for(int i = 0; i < M; ++i) for(int j = 0; j < L; ++j) //dbg - vhat_copy[i][j] = cg2s(gsl_matrix_complex_get(VHat, i, j)); //dbg -#endif if(M_inv_Vhat_function) { QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z0+z1, Params)); } else { - const complex double zero = 0, one = 1; lapack_int *pivot; gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(M,M); QPMS_ENSURE_SUCCESS(M_function(Mmat, z0+z1, Params)); -#ifndef NDEBUG - complex double Mmat_check[M][M]; //dbg - for(int i = 0; i < M; ++i) for(int j = 0; j < M; ++j) //dbg - Mmat_check[i][j] = cg2s(gsl_matrix_complex_get(Mmat, i, j)); //dbg -#endif QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * M); 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*/)); - -#ifndef NDEBUG - // Check the inversion result - complex double minvhat_check[M][L]; - for(int i = 0; i < M; ++i) for(int j = 0; j < L; ++j) - minvhat_check[i][j] = cg2s(gsl_matrix_complex_get(MInvVHat, i, j)); - complex double vhat_recon[M][L]; - - cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, - M, L, M, &one, Mmat_check[0], M, minvhat_check[0], L, &zero, vhat_recon[0], L); - - for(int i = 0; i < M; ++i) for(int j = 0; j < L; ++j) - if (cabs(vhat_recon[i][j] - vhat_copy[i][j]) > 1e-8*(1+cabs(vhat_recon[i][j] + vhat_copy[i][j]))) - QPMS_WTF; -#endif + M /*n*/, L/*nrhs*/, (lapack_complex_double *)(Mmat->data) /*a*/, Mmat->tda /*lda*/, pivot/*ipiv*/, + (lapack_complex_double *)(MInvVHat->data) /*b*/, MInvVHat->tda/*ldb*/)); free(pivot); From 1aa989015546a0ad38476f8deab7b0f54bed2ab2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 27 Aug 2019 20:41:30 +0300 Subject: [PATCH 13/34] Reindend, add also the 'coarse values' calculation Former-commit-id: 2dc73a2875823cae187585787bd4d344dea232f9 --- qpms/beyn.c | 201 ++++++++++++++++++++++++++-------------------------- 1 file changed, 101 insertions(+), 100 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 9024a21..fc3d797 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -51,9 +51,9 @@ double randN(double Sigma, double Mu) { #if 0 // Uniformly random number between -2 and 2 gsl_complex zrandN(){ - double a = (rand()*4.0/RAND_MAX) - 2; - double b = (rand()*4.0/RAND_MAX) - 2; - return gsl_complex_rect(a, b); + double a = (rand()*4.0/RAND_MAX) - 2; + double b = (rand()*4.0/RAND_MAX) - 2; + return gsl_complex_rect(a, b); } #endif @@ -97,13 +97,13 @@ BeynSolver *CreateBeynSolver(int M, int L) #if 0 // internal workspace: need storage for 2 MxL matrices // plus 3 LxL matrices - #define MLBUFFERS 2 - #define LLBUFFERS 3 +#define MLBUFFERS 2 +#define LLBUFFERS 3 size_t ML = MLMax*L, LL = L*L; #endif - + return Solver; - + } /***************************************************************/ @@ -134,12 +134,12 @@ void DestroyBeynSolver(BeynSolver *Solver) void ReRandomize(BeynSolver *Solver, unsigned int RandSeed) { if (RandSeed==0) - RandSeed=time(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(int nc=0; ncsize2; nc++) + gsl_matrix_complex_set(VHat,nr,nc,cs2g(zrandN(1, 0))); } @@ -163,7 +163,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, int Verbose = 1;//CheckEnv("SCUFF_BEYN_VERBOSE"); double RankTol=1.0e-4; //CheckEnv("SCUFF_BEYN_RANK_TOL",&RankTol); double ResTol=0.0; // CheckEnv("SCUFF_BEYN_RES_TOL",&ResTol); - + // A0 -> V0Full * Sigma * W0TFull' printf(" Beyn: computing SVD...\n"); gsl_matrix_complex* V0Full = gsl_matrix_complex_alloc(M,L); @@ -171,37 +171,37 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, gsl_matrix_complex* W0TFull = gsl_matrix_complex_alloc(L,L); //A0->SVD(Sigma, &V0Full, &W0TFull); - + QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride); QPMS_ENSURE(V0Full->size1 >= V0Full->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*/, - V0Full->size1 /* m, number of rows */, - V0Full->size2 /* n, number of columns */, - (lapack_complex_double *)(V0Full->data) /*a*/, - V0Full->tda /*lda*/, - 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 *)W0TFull->data /*vt*/, - W0TFull->tda /*ldvt*/ - )); + 'O' /*jobz, 'O' overwrites a with U and saves conjg(V') in vt if m >= n, i.e. if M >= L, which holds*/, + V0Full->size1 /* m, number of rows */, + V0Full->size2 /* n, number of columns */, + (lapack_complex_double *)(V0Full->data) /*a*/, + V0Full->tda /*lda*/, + 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 *)W0TFull->data /*vt*/, + W0TFull->tda /*ldvt*/ + )); + - // 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(%i)=%e\n",k,gsl_vector_get(Sigma, k)); - if (gsl_vector_get(Sigma, k) > RankTol ) + { if (Verbose) printf("Beyn: SV(%i)=%e\n",k,gsl_vector_get(Sigma, k)); + if (gsl_vector_get(Sigma, k) > RankTol ) K++; - } + } printf(" Beyn: %i/%i relevant singular values\n",K,L); if (K==0) - { printf("no singular values found in Beyn eigensolver\n"); - return 0; - } + { printf("no singular values found in Beyn eigensolver\n"); + return 0; + } + - // set V0, W0T = matrices of first K right, left singular vectors gsl_matrix_complex *V0 = gsl_matrix_complex_alloc(M,K); @@ -217,7 +217,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, gsl_matrix_complex_free(V0Full); gsl_matrix_complex_free(W0TFull); - + // B <- V0' * A1 * W0 * Sigma^-1 gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,L); @@ -228,14 +228,14 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, 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); + V0, A1, zero, TM2); + - printf(" Multiplying TM*W0T->B...\n"); //TM2.Multiply(&W0T, &B, "--transB C"); // B <- TM2 * W0 - + gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, - TM2, W0T, zero, B); + TM2, W0T, zero, B); gsl_matrix_complex_free(W0T); gsl_matrix_complex_free(TM2); @@ -243,14 +243,14 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, 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); + 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_vector_complex_free(tmp); //for(int m=0; mGetEntry(n)); @@ -282,7 +282,7 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, )); gsl_matrix_complex_free(B); - + //B.NSEig(&Lambda, &S); // V0S <- V0*S @@ -333,21 +333,21 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, /***************************************************************/ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, - beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *Params, - double complex z0, double Rx, double Ry, int N) + beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *Params, + double complex z0, double Rx, double Ry, int N) { /***************************************************************/ /* force N to be even so we can simultaneously evaluate */ /* the integral with N/2 quadrature points */ /***************************************************************/ - + if ( (N%2)==1 ) N++; /*if (Rx==Ry) - printf("Applying Beyn method with z0=%s,R=%e,N=%i...\n",z2s(z0),Rx,N); - else - printf("Applying Beyn method with z0=%s,Rx=%e,Ry=%e,N=%i...\n",z2s(z0),Rx,Ry,N); -*/ + printf("Applying Beyn method with z0=%s,R=%e,N=%i...\n",z2s(z0),Rx,N); + else + printf("Applying Beyn method with z0=%s,Rx=%e,Ry=%e,N=%i...\n",z2s(z0),Rx,Ry,N); + */ const int M = Solver->M; const int L = Solver->L; gsl_matrix_complex *A0 = Solver->A0; @@ -361,7 +361,7 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, /* 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(A0Coarse); @@ -369,69 +369,70 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, double DeltaTheta = 2.0*M_PI / ((double)N); printf(" Evaluating contour integral (%i points)...\n",N); for(int n=0; nCopy(VHat); - // Mitä varten tämä kopiointi on? - gsl_matrix_complex_memcpy(MInvVHat, VHat); + //MInvVHat->Copy(VHat); + // Mitä varten tämä kopiointi on? + gsl_matrix_complex_memcpy(MInvVHat, VHat); - // Tän pitäis kutsua eval_WT - // Output ilmeisesti tallentuun neljänteen parametriin - - if(M_inv_Vhat_function) { - QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z0+z1, Params)); - } else { - lapack_int *pivot; - gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(M,M); - QPMS_ENSURE_SUCCESS(M_function(Mmat, z0+z1, Params)); - QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * M); - 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*/)); + // Tän pitäis kutsua eval_WT + // Output ilmeisesti tallentuun neljänteen parametriin + + if(M_inv_Vhat_function) { + QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z0+z1, Params)); + } else { + lapack_int *pivot; + gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(M,M); + QPMS_ENSURE_SUCCESS(M_function(Mmat, z0+z1, Params)); + QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * M); + 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*/)); - free(pivot); - gsl_matrix_complex_free(Mmat); - } - //UserFunc(z0+zz, Params, VHat, MInvVHat); + free(pivot); + gsl_matrix_complex_free(Mmat); + } + //UserFunc(z0+zz, Params, VHat, MInvVHat); - gsl_matrix_complex_scale(MInvVHat, cs2g(dz)); - gsl_matrix_complex_add(A0, MInvVHat); - if((n%2)==0) { - gsl_matrix_complex_add(A0Coarse, MInvVHat); - gsl_matrix_complex_add(A0Coarse, MInvVHat); - } - - gsl_matrix_complex_scale(MInvVHat, cs2g(z1)); - gsl_matrix_complex_add(A1, MInvVHat); - if((n%2)==0) { - gsl_matrix_complex_add(A1Coarse, MInvVHat); - gsl_matrix_complex_add(A1Coarse, MInvVHat); - } + gsl_matrix_complex_scale(MInvVHat, cs2g(dz)); + gsl_matrix_complex_add(A0, MInvVHat); + if((n%2)==0) { + gsl_matrix_complex_add(A0Coarse, MInvVHat); + gsl_matrix_complex_add(A0Coarse, MInvVHat); + } + + gsl_matrix_complex_scale(MInvVHat, cs2g(z1)); + gsl_matrix_complex_add(A1, MInvVHat); + if((n%2)==0) { + gsl_matrix_complex_add(A1Coarse, MInvVHat); + gsl_matrix_complex_add(A1Coarse, MInvVHat); + } } gsl_vector_complex *Eigenvalues = Solver->Eigenvalues; - //gsl_vector_complex *EVErrors = Solver->EVErrors; + gsl_vector_complex *EVErrors = Solver->EVErrors; gsl_matrix_complex *Eigenvectors = Solver->Eigenvectors; int K = ProcessAMatrices(Solver, M_function, Params, A0, A1, z0, Eigenvalues, Eigenvectors); - //int KCoarse = ProcessAMatrices(Solver, UserFunc, Params, A0Coarse, A1Coarse, z0, EVErrors, Eigenvectors); -// Log("{K,KCoarse}={%i,%i}",K,KCoarse); - /* - for(int k=0; kN && kN; k++) - { EVErrors->ZV[k] -= Eigenvalues->ZV[k]; - EVErrors->ZV[k] = cdouble( fabs(real(EVErrors->ZV[k])), - fabs(imag(EVErrors->ZV[k])) - ); - } + int KCoarse = ProcessAMatrices(Solver, M_function, Params, A0Coarse, A1Coarse, z0, EVErrors, Eigenvectors); + // Log("{K,KCoarse}={%i,%i}",K,KCoarse); + gsl_blas_zaxpy(gsl_complex_rect(-1,0), Eigenvalues, EVErrors); +#if 0 + for(size_t k = 0; k < EVErrors->size && k < Eigenvalues->size; ++k) { -*/ + EVErrors->ZV[k] -= Eigenvalues->ZV[k]; + EVErrors->ZV[k] = cdouble( fabs(real(EVErrors->ZV[k])), + fabs(imag(EVErrors->ZV[k])) + ); + } +#endif return K; } From 4c21fde628cc20eb0ad19179dfc47c5c14b613b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 28 Aug 2019 15:12:42 +0300 Subject: [PATCH 14/34] WIP rewriting Beyn API Former-commit-id: 0579928c1aac32ec82db0364080d46fcce7721a6 --- qpms/beyn.c | 149 ++++++++++++++++++++++++++++++++++------------------ qpms/beyn.h | 97 ++++++++++------------------------ 2 files changed, 125 insertions(+), 121 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index fc3d797..0725722 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -1,22 +1,3 @@ -/* - * This file was originally part of SCUFF-EM by M. T. Homer Reid. - * Modified by Kristian Arjas and Marek Nečada to work without libhmat and libhrutil. - * - * SCUFF-EM is free software; you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation; either version 2 of the License, or - * (at your option) any later version. - * - * SCUFF-EM is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - */ - #define STATIC_ASSERT(COND,MSG) typedef char static_assertion_##MSG[(COND)?1:-1] #define cg2s(x) gsl_complex_tostd(x) @@ -41,27 +22,85 @@ STATIC_ASSERT((sizeof(lapack_complex_double) == sizeof(gsl_complex)), lapack_and_gsl_complex_must_be_consistent); -double randU(double A, double B) { return A + (B-A) * random() * (1. / RAND_MAX); } -double randN(double Sigma, double Mu) { + +typedef struct BeynSolver +{ + int M; // dimension of matrices + int L; // number of columns of VHat matrix + + gsl_vector_complex *Eigenvalues, *EVErrors; + gsl_matrix_complex *Eigenvectors; + gsl_matrix_complex *A0, *A1, *A0Coarse, *A1Coarse, *MInvVHat; + gsl_matrix_complex *VHat; + gsl_vector *Sigma, *Residuals; + double complex *Workspace; +} BeynSolver; + +// constructor, destructor +BeynSolver *CreateBeynSolver(int M, int L); +void DestroyBeynSolver(BeynSolver *Solver); + +// reset the random matrix VHat used in the Beyn algorithm +// +void ReRandomize(BeynSolver *Solver, unsigned int RandSeed); + +// for both of the following routines, +// the return value is the number of eigenvalues found, +// 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_t M_function, beyn_function_M_inv_Vhat_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); } + +static double randN(double Sigma, double Mu) { double u1 = randU(0, 1); double u2 = randU(0, 1); return Mu + Sigma*sqrt(-2*log(u1))*cos(2.*M_PI*u2); } -#if 0 -// Uniformly random number between -2 and 2 -gsl_complex zrandN(){ - double a = (rand()*4.0/RAND_MAX) - 2; - double b = (rand()*4.0/RAND_MAX) - 2; - return gsl_complex_rect(a, b); -} -#endif - -complex double zrandN(double sigma, double mu){ +static complex double zrandN(double sigma, double 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 *c; + QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + n*sizeof(c->z_dz[0])); + c->n = n; + for(size_t i = 0; i < n; ++i) { + double t = i*2*M_PI/n; + double st = sin(t), ct = cos(t); + c->z_zd[i][0] = centre + ct*rRe + st*rIm; + c->z_zd[i][1] = (I*rRe*st + rIm*ct) / n; + } + return c; +} + +void beyn_result_gsl_free(beyn_result_gsl_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); + free(r); + } +} + /***************************************************************/ /***************************************************************/ /***************************************************************/ @@ -73,11 +112,6 @@ BeynSolver *CreateBeynSolver(int M, int L) Solver->L = L; QPMS_ENSURE(L <= M, "We expect L <= M, but we got L = %d, M = %d ", L, M); -#if 0 - int MLMax = (M>L) ? M : L; -#endif - int MLMin = (MEigenvalues = gsl_vector_complex_calloc(L); Solver->EVErrors = gsl_vector_complex_calloc(L); @@ -91,7 +125,7 @@ BeynSolver *CreateBeynSolver(int M, int L) Solver->A1Coarse = 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(MLMin); + Solver->Sigma = gsl_vector_calloc(L); ReRandomize(Solver,(unsigned)time(NULL)); #if 0 @@ -127,6 +161,20 @@ void DestroyBeynSolver(BeynSolver *Solver) free(Solver); } +void DestroyBeynSolverPartial(BeynSolver *solver) +{ + gsl_matrix_complex_free(Solver->A0); + gsl_matrix_complex_free(Solver->A1); + gsl_matrix_complex_free(Solver->A0Coarse); + gsl_matrix_complex_free(Solver->A1Coarse); + gsl_matrix_complex_free(Solver->MInvVHat); + gsl_vector_free(Solver->Sigma); + gsl_matrix_complex_free(Solver->VHat); + + free(Solver); +} + + /***************************************************************/ /***************************************************************/ /***************************************************************/ @@ -261,10 +309,6 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, gsl_vector_complex *Lambda = gsl_vector_complex_alloc(K); // Eigenvalues gsl_matrix_complex *S = gsl_matrix_complex_alloc(K,K); // Eigenvectors - // FIXME general complex eigensystems not supported by GSL (use LAPACKE_zgee?) - //gsl_eigen_genv_workspace * W = gsl_eigen_genv_alloc(K); - //gsl_eigen_genv(B, Eye, alph, beta, S,W); - //gsl_eigen_genv_free(W); 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_SUCCESS(LAPACKE_zgeev( @@ -283,8 +327,6 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, gsl_matrix_complex_free(B); - //B.NSEig(&Lambda, &S); - // V0S <- V0*S printf("Multiplying V0*S...\n"); gsl_matrix_complex *V0S = gsl_matrix_complex_alloc(M, K); @@ -436,12 +478,15 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, return K; } -/***************************************************************/ -/***************************************************************/ -/***************************************************************/ -/* -int BeynSolve(BeynSolver *Solver, - BeynFunction UserFunction, void *Params, - cdouble z0, double R, int N) -{ return BeynSolve(Solver, UserFunction, Params, z0, R, R, N); } -*/ +// This is currently just a wrapper over the old mess that is to be rewritten gradually. +int beyn_solve_gsl(beyn_result_gsl_t **result, 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, double res_tol) +{ + + + return 0; +} + + + diff --git a/qpms/beyn.h b/qpms/beyn.h index 6b9f69f..aab2763 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -1,33 +1,3 @@ -/* Copyright (C) 2005-2011 M. T. Homer Reid - * - * This file is part of SCUFF-EM. - * - * SCUFF-EM is free software; you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation; either version 2 of the License, or - * (at your option) any later version. - * - * SCUFF-EM is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - */ - -/* - * libBeyn.h -- header file for libBeyn, a simple implementation of - * -- Beyn's algorithm for nonlinear eigenproblems - * -- - * -- This is packaged together with SCUFF-EM, but - * -- it is really a standalone independent entity - * -- for general-purpose use in solving nonlinear - * -- eigenproblems. - */ - - #ifndef BEYN_H #define BEYN_H @@ -36,54 +6,43 @@ #include /// User-supplied function that provides the operator M(z) whose "roots" are to be found. -typedef int (*beyn_function_M_t)(gsl_matrix_complex *target_M, complex double z, void *params); +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$. -typedef int (*beyn_function_M_inv_Vhat_t)(gsl_matrix_complex *target_M_inv_Vhat, +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); -/***************************************************************/ -/***************************************************************/ -/***************************************************************/ -typedef struct BeynSolver -{ - int M; // dimension of matrices - int L; // number of columns of VHat matrix - gsl_vector_complex *Eigenvalues, *EVErrors; - gsl_matrix_complex *Eigenvectors; - gsl_matrix_complex *A0, *A1, *A0Coarse, *A1Coarse, *MInvVHat; - gsl_matrix_complex *VHat; - gsl_vector *Sigma, *Residuals; - double complex *Workspace; +/// Complex plane integration contour structure. +typedef struct beyn_contour_t { + size_t n; ///< Number of discretisation points. + complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points. +} beyn_contour_t; - } BeynSolver; +/// Complex plane elliptic integration contour with axes parallel to the real, imaginary axes. +/** Free using free(). */ +beyn_contour_t *beyn_contour_ellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints); -// constructor, destructor -BeynSolver *CreateBeynSolver(int M, int L); -void DestroyBeynSolver(BeynSolver *Solver); +typedef struct beyn_result_gsl_t { + size_t neig; ///< Number of eigenvalues found + gsl_vector_complex *eigval; + gsl_vector_complex *eigval_err; + gsl_vector *residuals; + gsl_matrix_complex *eigvec; +} beyn_result_gsl_t; -// reset the random matrix VHat used in the Beyn algorithm -// -void ReRandomize(BeynSolver *Solver, unsigned int RandSeed); - -// for both of the following routines, -// the return value is the number of eigenvalues found, -// 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_t M_function, beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *params, - double complex z0, double Rx, double Ry, int N); +void beyn_result_gsl_free(beyn_result_gsl_t *result); +int beyn_solve_gsl(beyn_result_gsl_t **result, + 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. + 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)); } From 537c4b2d377313c68929f799c886d47430757655 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 29 Aug 2019 17:18:59 +0300 Subject: [PATCH 15/34] Beyn solver new API Former-commit-id: cb242415d66acebd42aa6c12ef74696f5dea85de --- qpms/beyn.c | 88 +++++++++++++++++++++++---------------------------- qpms/beyn.h | 1 + tests/tbeyn.c | 13 ++++---- 3 files changed, 48 insertions(+), 54 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 0725722..ed43dd7 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -58,7 +58,7 @@ void ReRandomize(BeynSolver *Solver, unsigned int RandSeed); // 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_t M_function, beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *params, + 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); @@ -81,12 +81,13 @@ beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double r { beyn_contour_t *c; QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + n*sizeof(c->z_dz[0])); + c->centre = centre; c->n = n; for(size_t i = 0; i < n; ++i) { double t = i*2*M_PI/n; double st = sin(t), ct = cos(t); - c->z_zd[i][0] = centre + ct*rRe + st*rIm; - c->z_zd[i][1] = (I*rRe*st + rIm*ct) / n; + c->z_dz[i][0] = centre + ct*rRe + I*st*rIm; + c->z_dz[i][1] = (I*rRe*st + rIm*ct) / n; } return c; } @@ -161,7 +162,7 @@ void DestroyBeynSolver(BeynSolver *Solver) free(Solver); } -void DestroyBeynSolverPartial(BeynSolver *solver) +void DestroyBeynSolverPartial(BeynSolver *Solver) { gsl_matrix_complex_free(Solver->A0); gsl_matrix_complex_free(Solver->A1); @@ -198,10 +199,10 @@ void ReRandomize(BeynSolver *Solver, unsigned int RandSeed) /* eigenvalues and eigenvectors */ /***************************************************************/ -int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, +int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_gsl_t M_function, void *Params, gsl_matrix_complex *A0, gsl_matrix_complex *A1, double complex z0, - gsl_vector_complex *Eigenvalues, gsl_matrix_complex *Eigenvectors) + gsl_vector_complex *Eigenvalues, gsl_matrix_complex *Eigenvectors, const double RankTol, const double ResTol) { int M = Solver->M; int L = Solver->L; @@ -209,8 +210,8 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, int Verbose = 1;//CheckEnv("SCUFF_BEYN_VERBOSE"); - double RankTol=1.0e-4; //CheckEnv("SCUFF_BEYN_RANK_TOL",&RankTol); - double ResTol=0.0; // CheckEnv("SCUFF_BEYN_RES_TOL",&ResTol); + //double RankTol=1.0e-4; //CheckEnv("SCUFF_BEYN_RANK_TOL",&RankTol); + //double ResTol=0.0; // CheckEnv("SCUFF_BEYN_RES_TOL",&ResTol); // A0 -> V0Full * Sigma * W0TFull' printf(" Beyn: computing SVD...\n"); @@ -249,8 +250,6 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, return 0; } - - // 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); @@ -266,7 +265,6 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, gsl_matrix_complex_free(V0Full); gsl_matrix_complex_free(W0TFull); - // B <- V0' * A1 * W0 * Sigma^-1 gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,L); gsl_matrix_complex *B = gsl_matrix_complex_calloc(K,K); @@ -278,7 +276,6 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, V0, A1, zero, TM2); - printf(" Multiplying TM*W0T->B...\n"); //TM2.Multiply(&W0T, &B, "--transB C"); // B <- TM2 * W0 @@ -302,7 +299,6 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, // for(int n=0; nGetEntry(n)); - // B -> S*Lambda*S' printf(" Eigensolving (%i,%i)\n",K,K); @@ -335,7 +331,6 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, gsl_matrix_complex_free(V0); - int KRetained = 0; gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(M,M); gsl_vector_complex *MVk = gsl_vector_complex_alloc(M); @@ -373,17 +368,25 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_t M_function, /***************************************************************/ /***************************************************************/ /***************************************************************/ - +#if 0 int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *Params, double complex z0, double Rx, double Ry, int N) +#endif + +int beyn_solve_gsl(beyn_result_gsl_t **result, size_t m, size_t l, + 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, + double rank_tol, double res_tol) { + BeynSolver *Solver = CreateBeynSolver(m, l); + /***************************************************************/ /* force N to be even so we can simultaneously evaluate */ /* the integral with N/2 quadrature points */ /***************************************************************/ - if ( (N%2)==1 ) N++; + // if ( (N%2)==1 ) N++; /*if (Rx==Ry) printf("Applying Beyn method with z0=%s,R=%e,N=%i...\n",z2s(z0),Rx,N); @@ -408,13 +411,13 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, gsl_matrix_complex_set_zero(A1); gsl_matrix_complex_set_zero(A0Coarse); gsl_matrix_complex_set_zero(A1Coarse); + size_t N = contour->n; double DeltaTheta = 2.0*M_PI / ((double)N); - printf(" Evaluating contour integral (%i points)...\n",N); + printf(" Evaluating contour integral (%zd points)...\n",N); + const complex double z0 = contour->centre; for(int n=0; nz_dz[n][0]; + const complex double dz = contour->z_dz[n][1]; //MInvVHat->Copy(VHat); // Mitä varten tämä kopiointi on? @@ -424,11 +427,11 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, // Output ilmeisesti tallentuun neljänteen parametriin if(M_inv_Vhat_function) { - QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z0+z1, Params)); + QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z, params)); } else { lapack_int *pivot; gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(M,M); - QPMS_ENSURE_SUCCESS(M_function(Mmat, z0+z1, Params)); + QPMS_ENSURE_SUCCESS(M_function(Mmat, z, params)); QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * M); QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, M /*m*/, M /*n*/,(lapack_complex_double *) Mmat->data /*a*/, Mmat->tda /*lda*/, pivot /*ipiv*/)); @@ -437,11 +440,10 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, M /*n*/, L/*nrhs*/, (lapack_complex_double *)(Mmat->data) /*a*/, Mmat->tda /*lda*/, pivot/*ipiv*/, (lapack_complex_double *)(MInvVHat->data) /*b*/, MInvVHat->tda/*ldb*/)); - free(pivot); gsl_matrix_complex_free(Mmat); } - //UserFunc(z0+zz, Params, VHat, MInvVHat); + //UserFunc(z0+zz, params, VHat, MInvVHat); gsl_matrix_complex_scale(MInvVHat, cs2g(dz)); gsl_matrix_complex_add(A0, MInvVHat); @@ -450,7 +452,7 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, gsl_matrix_complex_add(A0Coarse, MInvVHat); } - gsl_matrix_complex_scale(MInvVHat, cs2g(z1)); + gsl_matrix_complex_scale(MInvVHat, cs2g(z - z0)); gsl_matrix_complex_add(A1, MInvVHat); if((n%2)==0) { gsl_matrix_complex_add(A1Coarse, MInvVHat); @@ -462,31 +464,21 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, gsl_vector_complex *EVErrors = Solver->EVErrors; gsl_matrix_complex *Eigenvectors = Solver->Eigenvectors; - int K = ProcessAMatrices(Solver, M_function, Params, A0, A1, z0, Eigenvalues, Eigenvectors); - int KCoarse = ProcessAMatrices(Solver, M_function, Params, A0Coarse, A1Coarse, z0, EVErrors, Eigenvectors); + int K = ProcessAMatrices(Solver, M_function, params, A0, A1, z0, Eigenvalues, Eigenvectors, rank_tol, res_tol); + int KCoarse = ProcessAMatrices(Solver, M_function, params, A0Coarse, A1Coarse, z0, EVErrors, Eigenvectors, rank_tol, res_tol); // Log("{K,KCoarse}={%i,%i}",K,KCoarse); + gsl_blas_zaxpy(gsl_complex_rect(-1,0), Eigenvalues, EVErrors); -#if 0 - for(size_t k = 0; k < EVErrors->size && k < Eigenvalues->size; ++k) { + // TODO Original did also fabs on the complex and real parts ^^^. + + QPMS_CRASHING_MALLOC(*result, sizeof(beyn_result_gsl_t)); + (*result)->eigval = Solver->Eigenvalues; + (*result)->eigval_err = Solver->EVErrors; + (*result)->residuals = Solver->Residuals; + (*result)->eigvec = Solver->Eigenvectors; + + DestroyBeynSolverPartial(Solver); - EVErrors->ZV[k] -= Eigenvalues->ZV[k]; - EVErrors->ZV[k] = cdouble( fabs(real(EVErrors->ZV[k])), - fabs(imag(EVErrors->ZV[k])) - ); - } -#endif return K; } -// This is currently just a wrapper over the old mess that is to be rewritten gradually. -int beyn_solve_gsl(beyn_result_gsl_t **result, 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, double res_tol) -{ - - - return 0; -} - - - diff --git a/qpms/beyn.h b/qpms/beyn.h index aab2763..0d4b45e 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -16,6 +16,7 @@ typedef int (*beyn_function_M_inv_Vhat_gsl_t)(gsl_matrix_complex *target_M_inv_V /// Complex plane integration contour structure. typedef struct beyn_contour_t { size_t n; ///< Number of discretisation points. + complex double centre; ///< TODO what is the exact purpose of this? complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points. } beyn_contour_t; diff --git a/tests/tbeyn.c b/tests/tbeyn.c index cfd1703..bbf3f98 100644 --- a/tests/tbeyn.c +++ b/tests/tbeyn.c @@ -24,17 +24,18 @@ int main() { complex double z0 = 150+2*I; double Rx = 148, Ry = 148; int L = 10, N = 50, dim = 400; - BeynSolver * solver = CreateBeynSolver(dim, L); - ReRandomize(solver, 666); + beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); - int K = BeynSolve(solver, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, - z0, Rx, Ry, N); + beyn_result_gsl_t *result; + int K = beyn_solve_gsl(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, + contour, 1e-4, 0); printf("Found %d eigenvalues:\n", K); for (int i = 0; i < K; ++i) { - gsl_complex eig = gsl_vector_complex_get(solver->Eigenvalues, i); + gsl_complex eig = gsl_vector_complex_get(result->eigval, i); printf("%d: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); } - DestroyBeynSolver(solver); + free(contour); + beyn_result_gsl_free(result); return 0; } From bc703cbcca84f58c19e6ed26aff7b77b8c2799ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 2 Sep 2019 10:33:07 +0300 Subject: [PATCH 16/34] Beyn pure C api headers Former-commit-id: 0e1da18ccf552d496eed4a073611c57d5e9619ba --- qpms/beyn.h | 38 +++++++++++++++++++++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/qpms/beyn.h b/qpms/beyn.h index 0d4b45e..bfa6c0a 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -1,3 +1,6 @@ +/** \file beyn.h + * \brief Beyn's algorithm for nonlinear eigenvalue problems. + */ #ifndef BEYN_H #define BEYN_H @@ -6,12 +9,22 @@ #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. +/** Pure C array version */ +typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex double z, void *params); + +/// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$. +/** Pure C array version */ +typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, size_t m, size_t l, + const gsl_matrix_complex *Vhat, complex double z, void *params); /// Complex plane integration contour structure. typedef struct beyn_contour_t { @@ -24,8 +37,9 @@ typedef struct beyn_contour_t { /** Free using free(). */ beyn_contour_t *beyn_contour_ellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints); +/// Beyn algorithm result structure (GSL matrix/vector version). typedef struct beyn_result_gsl_t { - size_t neig; ///< Number of eigenvalues found + size_t neig; ///< Number of eigenvalues found (a bit redundant?). gsl_vector_complex *eigval; gsl_vector_complex *eigval_err; gsl_vector *residuals; @@ -34,6 +48,17 @@ typedef struct 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. + complex double *eigval; + complex double *eigval_err; + double *residuals; + complex double *eigvec; +} beyn_result_t; + +void beyn_result_free(beyn_result_t *result); + int beyn_solve_gsl(beyn_result_gsl_t **result, 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). @@ -45,6 +70,17 @@ int beyn_solve_gsl(beyn_result_gsl_t **result, double res_tol ///< (default: `0.0`) TODO DOC. ); +int beyn_solve(beyn_result_t **result, + 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_t M, ///< Function providing the matrix \f$ M(z) \f$. + beyn_function_M_inv_Vhat_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. + 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)); } From cabe7646407e2951c928845afec4fb551dd3d47b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 2 Sep 2019 11:57:25 +0300 Subject: [PATCH 17/34] Implement an API for Beyn's algorithm with standard C arrays. Former-commit-id: 645490de63f802c8a41f3bad1845cd29e0c3d823 --- qpms/beyn.c | 55 ++++++++++++++++++++++++++++++++++++++++++-- qpms/beyn.h | 6 ++++- tests/CMakeLists.txt | 4 ++++ tests/tbeyn.c | 30 +++++++++++------------- tests/tbeyn_gsl.c | 42 +++++++++++++++++++++++++++++++++ 5 files changed, 118 insertions(+), 19 deletions(-) create mode 100644 tests/tbeyn_gsl.c diff --git a/qpms/beyn.c b/qpms/beyn.c index ed43dd7..9ae71f7 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -18,7 +18,7 @@ #include #include -#include +#include "beyn.h" STATIC_ASSERT((sizeof(lapack_complex_double) == sizeof(gsl_complex)), lapack_and_gsl_complex_must_be_consistent); @@ -412,7 +412,7 @@ int beyn_solve_gsl(beyn_result_gsl_t **result, size_t m, size_t l, gsl_matrix_complex_set_zero(A0Coarse); gsl_matrix_complex_set_zero(A1Coarse); size_t N = contour->n; - double DeltaTheta = 2.0*M_PI / ((double)N); + //double DeltaTheta = 2.0*M_PI / ((double)N); printf(" Evaluating contour integral (%zd points)...\n",N); const complex double z0 = contour->centre; for(int n=0; nsize2 == 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); +} + +int beyn_solve(beyn_result_t **result, 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, double res_tol) { + struct beyn_function_M_carr2gsl_param p = {M, M_inv_Vhat, params}; + QPMS_CRASHING_MALLOC(*result, sizeof(beyn_result_t)); + int retval = beyn_solve_gsl(&((*result)->gsl), m, l, beyn_function_M_carr2gsl, + (p.M_inv_Vhat_function) ? beyn_function_M_inv_Vhat_carr2gsl : NULL, + (void *) &p, contour, rank_tol, res_tol); + (*result)->neig = (*result)->gsl->neig; + (*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; + return retval; +} + +void beyn_result_free(beyn_result_t *result) { + if(result) + beyn_result_gsl_free(result->gsl); + free(result); +} + diff --git a/qpms/beyn.h b/qpms/beyn.h index bfa6c0a..df605a7 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -24,7 +24,7 @@ typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex dou /// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$. /** Pure C array version */ typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, size_t m, size_t l, - const gsl_matrix_complex *Vhat, complex double z, void *params); + const complex double *Vhat, complex double z, void *params); /// Complex plane integration contour structure. typedef struct beyn_contour_t { @@ -55,6 +55,10 @@ typedef struct beyn_result_t { complex double *eigval_err; double *residuals; complex double *eigvec; + + /// Private, we wrap it around beyn_result_gsl_t for now. + beyn_result_gsl_t *gsl; + } beyn_result_t; void beyn_result_free(beyn_result_t *result); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index eb0b17e..5413990 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -19,6 +19,10 @@ add_executable( tbeyn tbeyn.c ) target_link_libraries( tbeyn qpms gsl lapacke amos m ) target_include_directories( tbeyn 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_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.c b/tests/tbeyn.c index bbf3f98..1e0322c 100644 --- a/tests/tbeyn.c +++ b/tests/tbeyn.c @@ -1,21 +1,19 @@ #include #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; +int M_function(complex double *target, const size_t m, const complex double z, void *no_params) { + complex double d = 2*m - 4*z / (6*m); + complex double od = -((double)m) - z / (6*m); - gsl_complex d = gsl_complex_fromstd( 2*m - 4*z / (6*m) ); - gsl_complex od = gsl_complex_fromstd( -m - z / (6*m) ); - - gsl_matrix_complex_set_zero(target); + memset(target, 0, m*m*sizeof(complex double)); 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); + target[i*m + i] = d; + if(i > 0) target[i*m + i-1] = od; + if(i < m - 1) target[i*m + i+1] = od; } - gsl_matrix_complex_set(target, m-1, m-1, gsl_complex_fromstd(gsl_complex_tostd(d)/2 + z/(z-1))); + target[m*(m-1) + m-1] = d/2 + z/(z-1); return 0; } @@ -26,16 +24,16 @@ int main() { int L = 10, N = 50, dim = 400; beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); - beyn_result_gsl_t *result; - int K = beyn_solve_gsl(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, + beyn_result_t *result; + int K = beyn_solve(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, contour, 1e-4, 0); printf("Found %d eigenvalues:\n", K); for (int i = 0; i < K; ++i) { - gsl_complex eig = gsl_vector_complex_get(result->eigval, i); - printf("%d: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); + complex double eig = result->eigval[i]; + printf("%d: %g%+gj\n", i, creal(eig), cimag(eig)); } free(contour); - beyn_result_gsl_free(result); + beyn_result_free(result); return 0; } diff --git a/tests/tbeyn_gsl.c b/tests/tbeyn_gsl.c new file mode 100644 index 0000000..8d492d3 --- /dev/null +++ b/tests/tbeyn_gsl.c @@ -0,0 +1,42 @@ +#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; + int K = beyn_solve_gsl(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, + contour, 1e-4, 0); + printf("Found %d eigenvalues:\n", K); + for (int i = 0; i < K; ++i) { + gsl_complex eig = gsl_vector_complex_get(result->eigval, i); + printf("%d: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); + } + free(contour); + beyn_result_gsl_free(result); + return 0; +} + + From 5dc2a44cdd9c714e7656cabd93c5ddd347d30306 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 2 Sep 2019 12:25:26 +0300 Subject: [PATCH 18/34] Incremental cleanup and style assimilation Former-commit-id: eb31cf03313edbad256331592afa4c6d212feab8 --- qpms/beyn.c | 284 ++++++++++++++++++++++------------------------------ 1 file changed, 120 insertions(+), 164 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 9ae71f7..c50b625 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -25,41 +25,41 @@ STATIC_ASSERT((sizeof(lapack_complex_double) == sizeof(gsl_complex)), lapack_and typedef struct BeynSolver { - int M; // dimension of matrices - int L; // number of columns of VHat matrix + int M; // dimension of matrices + int L; // number of columns of VHat matrix - gsl_vector_complex *Eigenvalues, *EVErrors; - gsl_matrix_complex *Eigenvectors; - gsl_matrix_complex *A0, *A1, *A0Coarse, *A1Coarse, *MInvVHat; - gsl_matrix_complex *VHat; - gsl_vector *Sigma, *Residuals; - double complex *Workspace; + 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; + double complex *Workspace; } BeynSolver; // constructor, destructor -BeynSolver *CreateBeynSolver(int M, int L); -void DestroyBeynSolver(BeynSolver *Solver); +BeynSolver *BeynSolver_create(int M, int L); +void BeynSolver_free(BeynSolver *solver); // reset the random matrix VHat used in the Beyn algorithm // -void ReRandomize(BeynSolver *Solver, unsigned int RandSeed); +void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed); // for both of the following routines, // the return value is the number of eigenvalues found, // and the eigenvalues and eigenvectors are stored in the -// Lambda and Eigenvectors fields of the BeynSolver structure +// 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, +//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); +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); @@ -102,90 +102,71 @@ void beyn_result_gsl_free(beyn_result_gsl_t *r) { } } -/***************************************************************/ -/***************************************************************/ -/***************************************************************/ -BeynSolver *CreateBeynSolver(int M, int L) +BeynSolver *BeynSolver_create(int M, int L) { - BeynSolver *Solver= (BeynSolver *)malloc(sizeof(*Solver)); + BeynSolver *solver= (BeynSolver *)malloc(sizeof(*solver)); - Solver->M = M; - Solver->L = L; + 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->EVErrors = gsl_vector_complex_calloc(L); - Solver->Residuals = gsl_vector_calloc(L); - Solver->Eigenvectors = gsl_matrix_complex_calloc(M, L); + 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(M, L); // 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->A0Coarse = gsl_matrix_complex_calloc(M,L); - Solver->A1Coarse = 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); - ReRandomize(Solver,(unsigned)time(NULL)); + 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); + BeynSolver_srandom(solver,(unsigned)time(NULL)); -#if 0 - // internal workspace: need storage for 2 MxL matrices - // plus 3 LxL matrices -#define MLBUFFERS 2 -#define LLBUFFERS 3 - size_t ML = MLMax*L, LL = L*L; -#endif - - return Solver; + return solver; } -/***************************************************************/ -/***************************************************************/ -/***************************************************************/ -void DestroyBeynSolver(BeynSolver *Solver) +void BeynSolver_free(BeynSolver *solver) { - gsl_vector_complex_free(Solver->Eigenvalues); - gsl_vector_complex_free(Solver->EVErrors); - gsl_matrix_complex_free(Solver->Eigenvectors); + gsl_vector_complex_free(solver->eigenvalues); + gsl_vector_complex_free(solver->eigenvalue_errors); + gsl_matrix_complex_free(solver->eigenvectors); - gsl_matrix_complex_free(Solver->A0); - gsl_matrix_complex_free(Solver->A1); - gsl_matrix_complex_free(Solver->A0Coarse); - gsl_matrix_complex_free(Solver->A1Coarse); - gsl_matrix_complex_free(Solver->MInvVHat); - gsl_vector_free(Solver->Sigma); - gsl_vector_free(Solver->Residuals); - gsl_matrix_complex_free(Solver->VHat); + 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); + free(solver); } -void DestroyBeynSolverPartial(BeynSolver *Solver) +void BeynSolver_free_partial(BeynSolver *solver) { - gsl_matrix_complex_free(Solver->A0); - gsl_matrix_complex_free(Solver->A1); - gsl_matrix_complex_free(Solver->A0Coarse); - gsl_matrix_complex_free(Solver->A1Coarse); - gsl_matrix_complex_free(Solver->MInvVHat); - gsl_vector_free(Solver->Sigma); - gsl_matrix_complex_free(Solver->VHat); + 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_matrix_complex_free(solver->VHat); - free(Solver); + free(solver); } - -/***************************************************************/ -/***************************************************************/ -/***************************************************************/ - -void ReRandomize(BeynSolver *Solver, unsigned int RandSeed) +void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed) { if (RandSeed==0) RandSeed=time(0); srandom(RandSeed); - gsl_matrix_complex *VHat=Solver->VHat; + 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))); @@ -199,48 +180,46 @@ void ReRandomize(BeynSolver *Solver, unsigned int RandSeed) /* eigenvalues and eigenvectors */ /***************************************************************/ -int ProcessAMatrices(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, 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 RankTol, const double ResTol) { - int M = Solver->M; - int L = Solver->L; - gsl_vector *Sigma = Solver->Sigma; + int M = solver->M; + int L = solver->L; + gsl_vector *Sigma = solver->Sigma; - int Verbose = 1;//CheckEnv("SCUFF_BEYN_VERBOSE"); - //double RankTol=1.0e-4; //CheckEnv("SCUFF_BEYN_RANK_TOL",&RankTol); - //double ResTol=0.0; // CheckEnv("SCUFF_BEYN_RES_TOL",&ResTol); + int verbose = 1; // TODO - // A0 -> V0Full * Sigma * W0TFull' + // A0 -> V0_full * Sigma * W0T_full' printf(" Beyn: computing SVD...\n"); - gsl_matrix_complex* V0Full = gsl_matrix_complex_alloc(M,L); - gsl_matrix_complex_memcpy(V0Full,A0); + gsl_matrix_complex* V0_full = gsl_matrix_complex_alloc(M,L); + gsl_matrix_complex_memcpy(V0_full,A0); - gsl_matrix_complex* W0TFull = gsl_matrix_complex_alloc(L,L); - //A0->SVD(Sigma, &V0Full, &W0TFull); + gsl_matrix_complex* W0T_full = gsl_matrix_complex_alloc(L,L); + //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(V0Full->size1 >= V0Full->size2, "m = %zd, l = %zd, what the hell?"); + 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*/, - V0Full->size1 /* m, number of rows */, - V0Full->size2 /* n, number of columns */, - (lapack_complex_double *)(V0Full->data) /*a*/, - V0Full->tda /*lda*/, + 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*/, NULL /*u; not used*/, M /*ldu; should not be really used but lapacke requires it for some obscure reason*/, - (lapack_complex_double *)W0TFull->data /*vt*/, - W0TFull->tda /*ldvt*/ + (lapack_complex_double *)W0T_full->data /*vt*/, + W0T_full->tda /*ldvt*/ )); // 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(%i)=%e\n",k,gsl_vector_get(Sigma, k)); + { if (verbose) printf("Beyn: SV(%i)=%e\n",k,gsl_vector_get(Sigma, k)); if (gsl_vector_get(Sigma, k) > RankTol ) K++; } @@ -256,14 +235,14 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_gsl_t M_function, for (int k = 0; k < K; ++k) { gsl_vector_complex_view tmp; - tmp = gsl_matrix_complex_column(V0Full, k); + tmp = gsl_matrix_complex_column(V0_full, k); gsl_matrix_complex_set_col(V0, k, &(tmp.vector)); - tmp = gsl_matrix_complex_row(W0TFull, k); + tmp = gsl_matrix_complex_row(W0T_full, k); gsl_matrix_complex_set_row(W0T, k, &(tmp.vector)); } - gsl_matrix_complex_free(V0Full); - gsl_matrix_complex_free(W0TFull); + gsl_matrix_complex_free(V0_full); + gsl_matrix_complex_free(W0T_full); // B <- V0' * A1 * W0 * Sigma^-1 gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,L); @@ -302,8 +281,8 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_gsl_t M_function, // B -> S*Lambda*S' 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 + gsl_vector_complex *Lambda = gsl_vector_complex_alloc(K); // eigenvalues + gsl_matrix_complex *S = gsl_matrix_complex_alloc(K,K); // eigenvectors 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); @@ -344,14 +323,14 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_gsl_t M_function, 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); - 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; - gsl_vector_complex_set(Eigenvalues, KRetained, zgsl); - if(Eigenvectors) { - gsl_matrix_complex_set_col(Eigenvectors, KRetained, &(Vk.vector)); - gsl_vector_set(Solver->Residuals, KRetained, residual); + gsl_vector_complex_set(eigenvalues, KRetained, zgsl); + if(eigenvectors) { + gsl_matrix_complex_set_col(eigenvectors, KRetained, &(Vk.vector)); + gsl_vector_set(solver->residuals, KRetained, residual); } ++KRetained; } @@ -365,42 +344,22 @@ int ProcessAMatrices(BeynSolver *Solver, beyn_function_M_gsl_t M_function, return KRetained; } -/***************************************************************/ -/***************************************************************/ -/***************************************************************/ -#if 0 -int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function, - beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *Params, - double complex z0, double Rx, double Ry, int N) -#endif int beyn_solve_gsl(beyn_result_gsl_t **result, size_t m, size_t l, 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, double rank_tol, double res_tol) { - BeynSolver *Solver = CreateBeynSolver(m, l); + BeynSolver *solver = BeynSolver_create(m, l); - /***************************************************************/ - /* force N to be even so we can simultaneously evaluate */ - /* the integral with N/2 quadrature points */ - /***************************************************************/ - - // if ( (N%2)==1 ) N++; - - /*if (Rx==Ry) - printf("Applying Beyn method with z0=%s,R=%e,N=%i...\n",z2s(z0),Rx,N); - else - printf("Applying Beyn method with z0=%s,Rx=%e,Ry=%e,N=%i...\n",z2s(z0),Rx,Ry,N); - */ - const int M = Solver->M; - const int L = Solver->L; - gsl_matrix_complex *A0 = Solver->A0; - gsl_matrix_complex *A1 = Solver->A1; - gsl_matrix_complex *A0Coarse = Solver->A0Coarse; - gsl_matrix_complex *A1Coarse = Solver->A1Coarse; - gsl_matrix_complex *MInvVHat = Solver->MInvVHat; - gsl_matrix_complex *VHat = Solver->VHat; + const int M = solver->M; + const int L = solver->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; /***************************************************************/ /* evaluate contour integrals by numerical quadrature to get */ @@ -409,18 +368,17 @@ int beyn_solve_gsl(beyn_result_gsl_t **result, size_t m, size_t l, gsl_matrix_complex_set_zero(A0); gsl_matrix_complex_set_zero(A1); - gsl_matrix_complex_set_zero(A0Coarse); - gsl_matrix_complex_set_zero(A1Coarse); - size_t N = contour->n; - //double DeltaTheta = 2.0*M_PI / ((double)N); - printf(" Evaluating contour integral (%zd points)...\n",N); + gsl_matrix_complex_set_zero(A0_coarse); + gsl_matrix_complex_set_zero(A1_coarse); + 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); + const complex double z0 = contour->centre; for(int n=0; nz_dz[n][0]; const complex double dz = contour->z_dz[n][1]; - //MInvVHat->Copy(VHat); - // Mitä varten tämä kopiointi on? gsl_matrix_complex_memcpy(MInvVHat, VHat); // Tän pitäis kutsua eval_WT @@ -443,41 +401,39 @@ int beyn_solve_gsl(beyn_result_gsl_t **result, size_t m, size_t l, free(pivot); gsl_matrix_complex_free(Mmat); } - //UserFunc(z0+zz, params, VHat, MInvVHat); gsl_matrix_complex_scale(MInvVHat, cs2g(dz)); gsl_matrix_complex_add(A0, MInvVHat); if((n%2)==0) { - gsl_matrix_complex_add(A0Coarse, MInvVHat); - gsl_matrix_complex_add(A0Coarse, MInvVHat); + gsl_matrix_complex_add(A0_coarse, MInvVHat); + gsl_matrix_complex_add(A0_coarse, MInvVHat); } gsl_matrix_complex_scale(MInvVHat, cs2g(z - z0)); gsl_matrix_complex_add(A1, MInvVHat); if((n%2)==0) { - gsl_matrix_complex_add(A1Coarse, MInvVHat); - gsl_matrix_complex_add(A1Coarse, MInvVHat); + gsl_matrix_complex_add(A1_coarse, MInvVHat); + gsl_matrix_complex_add(A1_coarse, MInvVHat); } } - gsl_vector_complex *Eigenvalues = Solver->Eigenvalues; - gsl_vector_complex *EVErrors = Solver->EVErrors; - gsl_matrix_complex *Eigenvectors = Solver->Eigenvectors; - - int K = ProcessAMatrices(Solver, M_function, params, A0, A1, z0, Eigenvalues, Eigenvectors, rank_tol, res_tol); - int KCoarse = ProcessAMatrices(Solver, M_function, params, A0Coarse, A1Coarse, z0, EVErrors, Eigenvectors, rank_tol, res_tol); - // Log("{K,KCoarse}={%i,%i}",K,KCoarse); - - gsl_blas_zaxpy(gsl_complex_rect(-1,0), Eigenvalues, EVErrors); + gsl_vector_complex *eigenvalues = solver->eigenvalues; + gsl_vector_complex *eigenvalue_errors = solver->eigenvalue_errors; + gsl_matrix_complex *eigenvectors = solver->eigenvectors; + + int K = beyn_process_matrices(solver, M_function, params, A0, A1, z0, eigenvalues, eigenvectors, rank_tol, res_tol); + int K_coarse = beyn_process_matrices(solver, M_function, params, A0_coarse, A1_coarse, z0, eigenvalue_errors, eigenvectors, rank_tol, res_tol); + + gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors); // TODO Original did also fabs on the complex and real parts ^^^. QPMS_CRASHING_MALLOC(*result, sizeof(beyn_result_gsl_t)); - (*result)->eigval = Solver->Eigenvalues; - (*result)->eigval_err = Solver->EVErrors; - (*result)->residuals = Solver->Residuals; - (*result)->eigvec = Solver->Eigenvectors; - - DestroyBeynSolverPartial(Solver); + (*result)->eigval = solver->eigenvalues; + (*result)->eigval_err = solver->eigenvalue_errors; + (*result)->residuals = solver->residuals; + (*result)->eigvec = solver->eigenvectors; + + BeynSolver_free_partial(solver); return K; } From 1585c48071a68bfed7d2f7f560fe4e7beb45f900 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 2 Sep 2019 13:31:22 +0300 Subject: [PATCH 19/34] Incremental cleanup. Former-commit-id: d4d8d41a2edac7cf8ac341cce46e3c976ef68c5e --- qpms/beyn.c | 112 ++++++++++++++++++++-------------------------------- 1 file changed, 42 insertions(+), 70 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index c50b625..9b05de7 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -33,39 +33,19 @@ typedef struct BeynSolver gsl_matrix_complex *A0, *A1, *A0_coarse, *A1_coarse, *MInvVHat; gsl_matrix_complex *VHat; gsl_vector *Sigma, *residuals; - double complex *Workspace; } BeynSolver; // constructor, destructor BeynSolver *BeynSolver_create(int M, int L); 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); -// for both of the following routines, -// the return value is the number of eigenvalues found, -// 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); } +// Uniformly random number from interval [a, b]. +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) { double u1 = 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); } - beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double rIm, size_t n) { 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 */ -/* matrices (obtained via numerical quadrature) to extract */ -/* eigenvalues and eigenvectors */ -/***************************************************************/ +/* + * linear-algebra manipulations on the A0 and A1 matrices + * (obtained via numerical quadrature) to extract eigenvalues + * and eigenvectors + */ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_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 RankTol, const double ResTol) + gsl_vector_complex *eigenvalues, gsl_matrix_complex *eigenvectors, const double rank_tol, const double res_tol) { - int M = solver->M; - int L = solver->L; + const size_t m = solver->M; + const size_t l = solver->L; gsl_vector *Sigma = solver->Sigma; - int verbose = 1; // TODO // A0 -> V0_full * Sigma * W0T_full' - printf(" Beyn: computing SVD...\n"); - gsl_matrix_complex* V0_full = gsl_matrix_complex_alloc(M,L); + 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); + gsl_matrix_complex* W0T_full = gsl_matrix_complex_alloc(l,l); //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); @@ -210,7 +188,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun V0_full->tda /*lda*/, Sigma->data /*s*/, 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->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) int K=0; - for(int k=0; ksize /* this is L, actually */; k++) - { if (verbose) printf("Beyn: SV(%i)=%e\n",k,gsl_vector_get(Sigma, k)); - if (gsl_vector_get(Sigma, k) > RankTol ) + 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 (gsl_vector_get(Sigma, k) > rank_tol) K++; } - printf(" Beyn: %i/%i relevant singular values\n",K,L); - if (K==0) - { printf("no singular values found in Beyn eigensolver\n"); + if (verbose)printf(" Beyn: %d/%zd relevant singular values\n",K,l); + if (K==0) { + QPMS_WARN("no singular values found in Beyn eigensolver\n"); return 0; } // 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); for (int k = 0; k < K; ++k) { 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); // 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); - printf(" Multiplying V0*A1->TM...\n"); + if(verbose) printf(" Multiplying V0*A1->TM...\n"); //V0.Multiply(A1, &TM2, "--transA C"); // TM2 <- V0' * A1 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); - printf(" Multiplying TM*W0T->B...\n"); - //TM2.Multiply(&W0T, &B, "--transB C"); // B <- TM2 * W0 + if(verbose) printf(" Multiplying TM*W0T->B...\n"); gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, 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(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); - for(int i = 0; i < K; i++){ + 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); @@ -275,11 +252,8 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun //for(int m=0; mGetEntry(n)); - // 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_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 */, (lapack_complex_double *) Lambda->data /* w */, 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 */, S->tda/* ldvr */ )); @@ -304,28 +278,28 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun // V0S <- V0*S 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, one, V0, S, zero, V0S)); gsl_matrix_complex_free(V0); 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); 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); double residual = 0; - if(ResTol > 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); 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); 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, void *params, const beyn_contour_t *contour, double rank_tol, double res_tol) { 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 *A1 = solver->A1; 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)); } else { 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_CRASHING_MALLOC(pivot, sizeof(lapack_int) * M); + QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * m); 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?"); + 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*/, + m /*n*/, l/*nrhs*/, (lapack_complex_double *)(Mmat->data) /*a*/, Mmat->tda /*lda*/, pivot/*ipiv*/, (lapack_complex_double *)(MInvVHat->data) /*b*/, MInvVHat->tda/*ldb*/)); free(pivot); From cbf039710f2afc06b40e609c9866c4acf0bb6bd0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 9 Sep 2019 20:50:31 +0300 Subject: [PATCH 20/34] cdefs for beyn.h Former-commit-id: 846f26b8a6e1a2829c379b1bcc9edae1a0965b26 --- qpms/qpms_cdefs.pxd | 44 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 2efe72c..941e51c 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -535,3 +535,47 @@ cdef extern from "ewald.h": int cx_gamma_inc_CF_e(double a, cdouble x, qpms_csf_result *result) int complex_gamma_inc_e(double a, cdouble x, qpms_csf_result *result) +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: + pass + ctypedef struct beyn_result_gsl_t: + pass + ctypedef struct beyn_result_t: + size_t neig + cdouble *eigval + cdouble *eigval_err + double *residuals + cdouble *eigvec + beyn_result_gsl_t *gsl + + 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) + + int beyn_solve_gsl(beyn_result_gsl_t **result, 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, double res_tol) + + int beyn_solve(beyn_result_t **result, 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, double res_tol) + + cdouble gsl_comlpex_tostd(gsl_complex z) + gsl_complex gsl_complex_fromstd(cdouble z) + From 3bc59096bcaa74eefe4c0e8f279c584756619633 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 10 Sep 2019 14:04:56 +0300 Subject: [PATCH 21/34] Public beyn_result_gsl_t -> beyn_result_gsl_t convertor. Former-commit-id: 56fb70384c20b93dbc8b119c6b5b4128ce4fcaf9 --- qpms/beyn.c | 21 +++++++++++++++------ qpms/beyn.h | 6 ++++++ 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 9b05de7..3c3b1ed 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -444,17 +444,26 @@ int beyn_solve(beyn_result_t **result, size_t m, size_t l, beyn_function_M_t M, void *params, const beyn_contour_t *contour, double rank_tol, double res_tol) { struct beyn_function_M_carr2gsl_param p = {M, M_inv_Vhat, params}; QPMS_CRASHING_MALLOC(*result, sizeof(beyn_result_t)); - int retval = beyn_solve_gsl(&((*result)->gsl), m, l, beyn_function_M_carr2gsl, + struct beyn_result_gsl_t *result_gsl; + int retval = beyn_solve_gsl(&result_gsl, m, l, beyn_function_M_carr2gsl, (p.M_inv_Vhat_function) ? beyn_function_M_inv_Vhat_carr2gsl : NULL, (void *) &p, contour, rank_tol, res_tol); - (*result)->neig = (*result)->gsl->neig; - (*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 = beyn_result_from_beyn_result_gsl(result_gsl); return retval; } +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->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; + return result; +} + void beyn_result_free(beyn_result_t *result) { if(result) beyn_result_gsl_free(result->gsl); diff --git a/qpms/beyn.h b/qpms/beyn.h index df605a7..b32025e 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -61,6 +61,12 @@ typedef struct 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); int beyn_solve_gsl(beyn_result_gsl_t **result, From 11aa48d2da0b90822adc5d86a4f0e322bba78b45 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 10 Sep 2019 14:54:08 +0300 Subject: [PATCH 22/34] Beyn minor fixes, cleaner API. Former-commit-id: 44d7147a14194a7b8e5a66552dd3855029b9e370 --- qpms/beyn.c | 30 ++++++++++++++++-------------- qpms/beyn.h | 5 +++-- qpms/qpms_cdefs.pxd | 5 +++-- tests/tbeyn.c | 10 +++++----- tests/tbeyn_gsl.c | 10 +++++----- 5 files changed, 32 insertions(+), 28 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 3c3b1ed..d2495af 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -284,6 +284,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun gsl_matrix_complex_free(V0); + // FIXME!!! make clear relation between KRetained and K in the results! int KRetained = 0; gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m, m); gsl_vector_complex *MVk = gsl_vector_complex_alloc(m); @@ -319,7 +320,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun } -int beyn_solve_gsl(beyn_result_gsl_t **result, const size_t m, const size_t l, +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, void *params, const beyn_contour_t *contour, double rank_tol, double res_tol) @@ -399,15 +400,17 @@ int beyn_solve_gsl(beyn_result_gsl_t **result, const size_t m, const size_t l, gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors); // TODO Original did also fabs on the complex and real parts ^^^. - QPMS_CRASHING_MALLOC(*result, sizeof(beyn_result_gsl_t)); - (*result)->eigval = solver->eigenvalues; - (*result)->eigval_err = solver->eigenvalue_errors; - (*result)->residuals = solver->residuals; - (*result)->eigvec = solver->eigenvectors; + beyn_result_gsl_t *result; + QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_gsl_t)); + result->eigval = solver->eigenvalues; + result->eigval_err = solver->eigenvalue_errors; + result->residuals = solver->residuals; + result->eigvec = solver->eigenvectors; + result->neig = K; BeynSolver_free_partial(solver); - return K; + return result; } // Wrapper of pure C array M-matrix function to GSL. @@ -440,16 +443,14 @@ static int beyn_function_M_inv_Vhat_carr2gsl(gsl_matrix_complex *target, (complex double *) Vhat->data, z, p->param); } -int beyn_solve(beyn_result_t **result, size_t m, size_t l, beyn_function_M_t M, beyn_function_M_inv_Vhat_t M_inv_Vhat, +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, double res_tol) { struct beyn_function_M_carr2gsl_param p = {M, M_inv_Vhat, params}; - QPMS_CRASHING_MALLOC(*result, sizeof(beyn_result_t)); - struct beyn_result_gsl_t *result_gsl; - int retval = beyn_solve_gsl(&result_gsl, m, l, beyn_function_M_carr2gsl, + 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, res_tol); - *result = beyn_result_from_beyn_result_gsl(result_gsl); - return retval; + (void *) &p, contour, rank_tol, res_tol) + ); } beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g) { @@ -457,6 +458,7 @@ beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g) { 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; diff --git a/qpms/beyn.h b/qpms/beyn.h index b32025e..cbcb64d 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -51,6 +51,7 @@ 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. + size_t vlen; ///< Vector space dimension (also the leading dimension of eigvec). complex double *eigval; complex double *eigval_err; double *residuals; @@ -69,7 +70,7 @@ beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g); void beyn_result_free(beyn_result_t *result); -int beyn_solve_gsl(beyn_result_gsl_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$. @@ -80,7 +81,7 @@ int beyn_solve_gsl(beyn_result_gsl_t **result, double res_tol ///< (default: `0.0`) TODO DOC. ); -int beyn_solve(beyn_result_t **result, +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). beyn_function_M_t M, ///< Function providing the matrix \f$ M(z) \f$. diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 941e51c..046c343 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -554,6 +554,7 @@ cdef extern from "beyn.h": pass ctypedef struct beyn_result_t: size_t neig + size_t vlen cdouble *eigval cdouble *eigval_err double *residuals @@ -568,11 +569,11 @@ cdef extern from "beyn.h": void beyn_result_gsl_free(beyn_result_gsl_t *result) void beyn_result_free(beyn_result_t *result) - int beyn_solve_gsl(beyn_result_gsl_t **result, size_t m, size_t l, beyn_function_M_gsl_t M, + 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, double res_tol) - int beyn_solve(beyn_result_t **result, size_t m, size_t l, beyn_function_M_t M, + 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, double res_tol) diff --git a/tests/tbeyn.c b/tests/tbeyn.c index 1e0322c..8cdb4ea 100644 --- a/tests/tbeyn.c +++ b/tests/tbeyn.c @@ -24,13 +24,13 @@ int main() { int L = 10, N = 50, dim = 400; beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); - beyn_result_t *result; - int K = beyn_solve(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, + beyn_result_t *result = + beyn_solve(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, contour, 1e-4, 0); - printf("Found %d eigenvalues:\n", K); - for (int i = 0; i < K; ++i) { + printf("Found %zd eigenvalues:\n", result->neig); + for (size_t i = 0; i < result->neig; ++i) { complex double eig = result->eigval[i]; - printf("%d: %g%+gj\n", i, creal(eig), cimag(eig)); + printf("%zd: %g%+gj\n", i, creal(eig), cimag(eig)); } free(contour); beyn_result_free(result); diff --git a/tests/tbeyn_gsl.c b/tests/tbeyn_gsl.c index 8d492d3..36c46a7 100644 --- a/tests/tbeyn_gsl.c +++ b/tests/tbeyn_gsl.c @@ -26,13 +26,13 @@ int main() { int L = 10, N = 50, dim = 400; beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); - beyn_result_gsl_t *result; - int K = beyn_solve_gsl(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, + beyn_result_gsl_t *result = + beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, contour, 1e-4, 0); - printf("Found %d eigenvalues:\n", K); - for (int i = 0; i < K; ++i) { + 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("%d: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); + printf("%zd: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig)); } free(contour); beyn_result_gsl_free(result); From 2762ada49c97961e5837dd42c407d8879cf5d519 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 16 Sep 2019 08:02:31 +0300 Subject: [PATCH 23/34] Beyn: keep singular values from rank test Former-commit-id: eb82025e460575335f9d6e3d9acb3cfdfe8867f0 --- qpms/beyn.c | 4 +++- qpms/beyn.h | 2 ++ qpms/qpms_cdefs.pxd | 1 + 3 files changed, 6 insertions(+), 1 deletion(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index d2495af..9cdf408 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -77,6 +77,7 @@ void beyn_result_gsl_free(beyn_result_gsl_t *r) { 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); } } @@ -134,7 +135,6 @@ void BeynSolver_free_partial(BeynSolver *solver) 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_matrix_complex_free(solver->VHat); free(solver); @@ -406,6 +406,7 @@ beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l, result->eigval_err = solver->eigenvalue_errors; result->residuals = solver->residuals; result->eigvec = solver->eigenvectors; + result->ranktest_SV = solver->Sigma; result->neig = K; BeynSolver_free_partial(solver); @@ -463,6 +464,7 @@ beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g) { 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; } diff --git a/qpms/beyn.h b/qpms/beyn.h index cbcb64d..bd17928 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -44,6 +44,7 @@ typedef struct beyn_result_gsl_t { gsl_vector_complex *eigval_err; gsl_vector *residuals; gsl_matrix_complex *eigvec; + gsl_vector *ranktest_SV; } beyn_result_gsl_t; void beyn_result_gsl_free(beyn_result_gsl_t *result); @@ -56,6 +57,7 @@ typedef struct beyn_result_t { complex double *eigval_err; double *residuals; complex double *eigvec; + double *ranktest_SV; /// Private, we wrap it around beyn_result_gsl_t for now. beyn_result_gsl_t *gsl; diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 046c343..7ea32e9 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -559,6 +559,7 @@ cdef extern from "beyn.h": cdouble *eigval_err double *residuals cdouble *eigvec + double *ranktest_SV beyn_result_gsl_t *gsl ctypedef int (*beyn_function_M_gsl_t)(gsl_matrix_complex *target_M, cdouble z, void *params) From 9d0f910bba0752374684599567cfa3e2fc4dd006 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 17 Sep 2019 07:17:04 +0300 Subject: [PATCH 24/34] Some comments to beyn.c Former-commit-id: cc756b3343a9290f3b3b8c79cec75559b11a952f --- qpms/beyn.c | 35 ++++++++++++++++++++++------------- qpms/beyn.h | 9 ++++++++- 2 files changed, 30 insertions(+), 14 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 9cdf408..1623666 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -104,6 +104,7 @@ BeynSolver *BeynSolver_create(int M, int L) solver->MInvVHat = gsl_matrix_complex_calloc(M,L); solver->VHat = gsl_matrix_complex_calloc(M,L); solver->Sigma = gsl_vector_calloc(L); + // Beyn Step 1: Generate random matrix VHat BeynSolver_srandom(solver,(unsigned)time(NULL)); return solver; @@ -170,14 +171,11 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun int verbose = 1; // TODO - // A0 -> V0_full * Sigma * W0T_full' + // 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); - //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(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') @@ -194,6 +192,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun )); + // 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++) { @@ -207,6 +206,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun return 0; } + // 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); @@ -222,12 +222,11 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun gsl_matrix_complex_free(V0_full); gsl_matrix_complex_free(W0T_full); - // B <- V0' * A1 * W0 * Sigma^-1 gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,l); gsl_matrix_complex *B = gsl_matrix_complex_calloc(K,K); if(verbose) printf(" Multiplying V0*A1->TM...\n"); - //V0.Multiply(A1, &TM2, "--transA C"); // TM2 <- V0' * A1 + const gsl_complex one = gsl_complex_rect(1,0); const gsl_complex zero = gsl_complex_rect(0,0); gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, @@ -251,8 +250,16 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun gsl_vector_complex_free(tmp); //for(int m=0; m S*Lambda*S' + + // Beyn step 6. + // Eigenvalue decomposition B -> S*Lambda*S' + /* According to Beyn's paper (Algorithm 1), one should check conditioning + * of the eigenvalues; if they are ill-conditioned, one should perform + * a procedure involving Schur decomposition and reordering. + * + * Beyn refers to MATLAB routines eig, condeig, schur, ordschur. + * I am not sure about the equivalents in LAPACK, TODO check zgeevx, zgeesx. + */ if(verbose) printf(" Eigensolving (%i,%i)\n",K,K); gsl_vector_complex *Lambda = gsl_vector_complex_alloc(K); // eigenvalues @@ -285,6 +292,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun gsl_matrix_complex_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); @@ -347,6 +355,8 @@ beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l, if(N & 1) QPMS_WARN("Contour discretisation point number is odd (%zd)," " the error estimates might be a bit off.", N); + + // Beyn Step 2: Computa A0, A1 const complex double z0 = contour->centre; for(int n=0; nz_dz[n][0]; @@ -354,9 +364,6 @@ beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l, gsl_matrix_complex_memcpy(MInvVHat, VHat); - // Tän pitäis kutsua eval_WT - // Output ilmeisesti tallentuun neljänteen parametriin - if(M_inv_Vhat_function) { QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z, params)); } else { @@ -382,7 +389,7 @@ beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l, gsl_matrix_complex_add(A0_coarse, MInvVHat); } - gsl_matrix_complex_scale(MInvVHat, cs2g(z - z0)); + 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); if((n%2)==0) { gsl_matrix_complex_add(A1_coarse, MInvVHat); @@ -394,11 +401,13 @@ beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l, gsl_vector_complex *eigenvalue_errors = solver->eigenvalue_errors; gsl_matrix_complex *eigenvectors = solver->eigenvectors; + // Beyn Steps 3–6 int K = beyn_process_matrices(solver, M_function, params, A0, A1, z0, eigenvalues, eigenvectors, rank_tol, res_tol); + // 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, rank_tol, res_tol); gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors); - // TODO Original did also fabs on the complex and real parts ^^^. + // Reid did also fabs on the complex and real parts ^^^. beyn_result_gsl_t *result; QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_gsl_t)); diff --git a/qpms/beyn.h b/qpms/beyn.h index bd17928..4cf9213 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -29,7 +29,14 @@ typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, siz /// Complex plane integration contour structure. typedef struct beyn_contour_t { size_t n; ///< Number of discretisation points. - complex double centre; ///< TODO what is the exact purpose of this? + /// "Centre" of the contour. + /** + * This point is used in the rescaling of the \f$ A_1 \f$ matrix as in + * Beyn's Remark 3.2 (b) in order to improve the numerical stability. + * It does not have to be a centre in some strictly defined sense, + * but it should be "somewhere around" where the contour is. + */ + complex double centre; complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points. } beyn_contour_t; From bd0b8a83e96d24919780e223666e1e0547a147d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 17 Sep 2019 07:58:59 +0300 Subject: [PATCH 25/34] Beyn: functionality for testing if points are inside contour. Former-commit-id: 5c88b6e9aa4308201871a9b19c8d20c04b93a718 --- qpms/beyn.c | 15 ++++++++++++++- qpms/beyn.h | 2 ++ qpms/qpms_cdefs.pxd | 1 + 3 files changed, 17 insertions(+), 1 deletion(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 1623666..b33291c 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -56,10 +56,19 @@ static complex double zrandN(double sigma, double mu){ return randN(sigma, mu) + I*randN(sigma, mu); } +static inline double dsq(double a) { return a * a; } + +static _Bool beyn_contour_ellipse_inside_test(struct beyn_contour_t *c, complex double z) { + double rRe = c->z_dz[c->n][0]; + double rIm = c->z_dz[c->n][1]; + complex double zn = z - c->centre; + return dsq(creal(zn)/rRe) + dsq(cimag(zn)/rIm) < 1; +} + beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double rIm, size_t n) { beyn_contour_t *c; - QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + n*sizeof(c->z_dz[0])); + QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + (n+1)*sizeof(c->z_dz[0])); c->centre = centre; c->n = n; for(size_t i = 0; i < n; ++i) { @@ -68,6 +77,10 @@ beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double r c->z_dz[i][0] = centre + ct*rRe + I*st*rIm; c->z_dz[i][1] = (I*rRe*st + rIm*ct) / n; } + // We hide the half-axes metadata after the discretisation points. + c->z_dz[n][0] = rRe; + c->z_dz[n][1] = rIm; + c->inside_test = beyn_contour_ellipse_inside_test; return c; } diff --git a/qpms/beyn.h b/qpms/beyn.h index 4cf9213..5eeccf9 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -37,6 +37,8 @@ typedef struct beyn_contour_t { * but it should be "somewhere around" where the contour is. */ complex double centre; + /// Function testing that a point \a z lies inside the contour (optional). + _Bool (*inside_test)(struct beyn_contour_t *, complex double z); complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points. } beyn_contour_t; diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 7ea32e9..1e7712e 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -549,6 +549,7 @@ cdef extern from "gsl/gsl_matrix.h": 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 From 42931eb0a0f321ed222b2e63222ca28e6c59b09c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 17 Sep 2019 14:11:23 +0300 Subject: [PATCH 26/34] Additional test to Beyn solver (Example 4.9). Former-commit-id: 63a78fa09fd3875e5d10b6b4fb6eea72010cfcee --- tests/CMakeLists.txt | 4 +++ tests/tbeyn.c | 2 +- tests/tbeyn2.c | 62 ++++++++++++++++++++++++++++++++++++++++++++ tests/tbeyn_gsl.c | 2 +- 4 files changed, 68 insertions(+), 2 deletions(-) create mode 100644 tests/tbeyn2.c diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5413990..ac7a18e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -19,6 +19,10 @@ add_executable( tbeyn tbeyn.c ) target_link_libraries( tbeyn qpms gsl lapacke amos m ) target_include_directories( tbeyn PRIVATE .. ) +add_executable( tbeyn2 tbeyn2.c ) +target_link_libraries( tbeyn2 qpms gsl lapacke amos m ) +target_include_directories( tbeyn2 PRIVATE .. ) + add_executable( tbeyn_gsl tbeyn_gsl.c ) target_link_libraries( tbeyn_gsl qpms gsl lapacke amos m ) target_include_directories( tbeyn_gsl PRIVATE .. ) diff --git a/tests/tbeyn.c b/tests/tbeyn.c index 8cdb4ea..e2f52f4 100644 --- a/tests/tbeyn.c +++ b/tests/tbeyn.c @@ -26,7 +26,7 @@ int main() { beyn_result_t *result = beyn_solve(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, - contour, 1e-4, 0); + contour, 1e-4, 1e-4); printf("Found %zd eigenvalues:\n", result->neig); for (size_t i = 0; i < result->neig; ++i) { complex double eig = result->eigval[i]; diff --git a/tests/tbeyn2.c b/tests/tbeyn2.c new file mode 100644 index 0000000..43b37d5 --- /dev/null +++ b/tests/tbeyn2.c @@ -0,0 +1,62 @@ +#include +#include +#include +#include +#include + +static double randU(double a, double b) {return a + (b-a) * random() * (1. / RAND_MAX); } +// Normal distribution via Box-Muller transform +static double randN(double sigma, double mu) { + double u1 = randU(0,1); + double u2 = randU(0,1); + return mu + sigma*sqrt(-2*log(u1))*cos(2.*M_PI*u2); +} + +struct param49 { + double *T0; + double *T1; + double *T2; +}; + +// Matrix as in Beyn, example 4.9 +int M_function(complex double *target, const size_t m, const complex double z, void *params) { + struct param49 *p = params; + + for(size_t i = 0; i < m*m; ++i) + target[i] = p->T0[i] + z*p->T1[i] + z*z*p->T2[i]; + return 0; +} + +int main(int argc, char **argv) { + complex double z0 = 0; + double Rx = .3, Ry = .3; + int L = 30, N = 150, dim = 60; + if (argc > 1) N = atoi(argv[1]); + beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); + struct param49 p; + p.T0 = malloc(dim*dim*sizeof(double)); + p.T1 = malloc(dim*dim*sizeof(double)); + p.T2 = malloc(dim*dim*sizeof(double)); + for(size_t i = 0; i < dim*dim; ++i) { + p.T0[i] = randN(1,0); + p.T1[i] = randN(1,0); + p.T2[i] = randN(1,0); + } + + beyn_result_t *result = + beyn_solve(dim, L, M_function, NULL /*M_inv_Vhat_function*/, &p /*params*/, + contour, 1e-4, 1e-4); + printf("Found %zd eigenvalues:\n", result->neig); + for (size_t i = 0; i < result->neig; ++i) { + complex double eig = result->eigval[i]; + printf("%zd: %g%+gj\n", i, creal(eig), cimag(eig)); + } + free(contour); + beyn_result_free(result); + free(p.T0); + free(p.T1); + free(p.T2); + return 0; +} + + diff --git a/tests/tbeyn_gsl.c b/tests/tbeyn_gsl.c index 36c46a7..9f73b0c 100644 --- a/tests/tbeyn_gsl.c +++ b/tests/tbeyn_gsl.c @@ -28,7 +28,7 @@ int main() { beyn_result_gsl_t *result = beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, - contour, 1e-4, 0); + contour, 1e-4, 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); From 6445c3523ecfea3c5e8a362e1bdbfaecff6b19f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 17 Sep 2019 15:14:09 +0300 Subject: [PATCH 27/34] Experiments with Beyn's algorithm and singularities. It seems that everything is fine in the end except essential singularities inside the contour. Former-commit-id: 9fe7135cb30d1fff1dc3ff9bf8a9152c6b0ef9b4 --- tests/CMakeLists.txt | 35 +++++++++++++++++ tests/tbeyn3.c | 89 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 124 insertions(+) create mode 100644 tests/tbeyn3.c diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index ac7a18e..8185d79 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -23,6 +23,41 @@ add_executable( tbeyn2 tbeyn2.c ) target_link_libraries( tbeyn2 qpms gsl lapacke amos m ) target_include_directories( tbeyn2 PRIVATE .. ) +add_executable( tbeyn3a tbeyn3.c ) +target_link_libraries( tbeyn3a qpms gsl lapacke amos m ) +target_include_directories( tbeyn3a PRIVATE .. ) +target_compile_definitions( tbeyn3a PRIVATE VARIANTA ) + +add_executable( tbeyn3b tbeyn3.c ) +target_link_libraries( tbeyn3b qpms gsl lapacke amos m ) +target_include_directories( tbeyn3b PRIVATE .. ) +target_compile_definitions( tbeyn3b PRIVATE VARIANTB RXSMALL ) + +add_executable( tbeyn3bfail tbeyn3.c ) +target_link_libraries( tbeyn3bfail qpms gsl lapacke amos m ) +target_include_directories( tbeyn3bfail PRIVATE .. ) +target_compile_definitions( tbeyn3bfail PRIVATE VARIANTB ) + +add_executable( tbeyn3c tbeyn3.c ) +target_link_libraries( tbeyn3c qpms gsl lapacke amos m ) +target_include_directories( tbeyn3c PRIVATE .. ) +target_compile_definitions( tbeyn3c PRIVATE VARIANTC RXSMALL ) + +add_executable( tbeyn3d tbeyn3.c ) +target_link_libraries( tbeyn3d qpms gsl lapacke amos m ) +target_include_directories( tbeyn3d PRIVATE .. ) +target_compile_definitions( tbeyn3d PRIVATE VARIANTD RXSMALL ) + +add_executable( tbeyn3e tbeyn3.c ) +target_link_libraries( tbeyn3e qpms gsl lapacke amos m ) +target_include_directories( tbeyn3e PRIVATE .. ) +target_compile_definitions( tbeyn3e PRIVATE VARIANTE RXSMALL ) + +add_executable( tbeyn3f tbeyn3.c ) +target_link_libraries( tbeyn3f qpms gsl lapacke amos 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 amos m ) target_include_directories( tbeyn_gsl PRIVATE .. ) diff --git a/tests/tbeyn3.c b/tests/tbeyn3.c new file mode 100644 index 0000000..0c983b6 --- /dev/null +++ b/tests/tbeyn3.c @@ -0,0 +1,89 @@ +#include +#include +#include +#include +#include + +static double randU(double a, double b) {return a + (b-a) * random() * (1. / RAND_MAX); } +// Normal distribution via Box-Muller transform +static double randN(double sigma, double mu) { + double u1 = randU(0,1); + double u2 = randU(0,1); + return mu + sigma*sqrt(-2*log(u1))*cos(2.*M_PI*u2); +} + +struct param { + double *T0; + double *T1; + double *T2; +}; + +int M_function(complex double *target, const size_t m, const complex double z, void *params) { + struct param *p = params; + + for(size_t i = 0; i < m*m; ++i) + target[i] = p->T0[i] + z*p->T1[i] + cexp( +#ifdef VARIANTB // Also note that this case requires pretty large contour point number (>~ 3000) + (1+3*I) +#else // VARIANTA or VARIANTC + (1+1*I) +#endif + *z*p->T2[i]) + +#ifdef VARIANTC // Essential singularity at zero + cexp(3/z); +#elif defined VARIANTD // Essential singularity outside the contour + cexp(3/(z-1)) +#elif defined VARIANTE // High-order pole at zero + 3/cpow(z,10); +#elif defined VARIANTF // High-order pole at zero, higher order than dim + .0003/cpow(z,12); +#else // double pole at zero + 3/z/z +#endif + ; + return 0; +} + +int main(int argc, char **argv) { + complex double z0 = 0; +#ifdef RXSMALL + double Rx = .1; +#else + double Rx = .3; // Variant B will fail in this case due to large number of eigenvalues (>30) +#endif + double Ry = .3; +#ifdef VARIANTF + int L = 10, N = 150, dim = 10; +#else + int L = 30, N = 150, dim = 60; +#endif + if (argc > 1) N = atoi(argv[1]); + if (argc > 2) L = atoi(argv[2]); + beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); + struct param p; + p.T0 = malloc(dim*dim*sizeof(double)); + p.T1 = malloc(dim*dim*sizeof(double)); + p.T2 = malloc(dim*dim*sizeof(double)); + for(size_t i = 0; i < dim*dim; ++i) { + p.T0[i] = randN(1,0); + p.T1[i] = randN(1,0); + p.T2[i] = randN(1,0); + } + + beyn_result_t *result = + beyn_solve(dim, L, M_function, NULL /*M_inv_Vhat_function*/, &p /*params*/, + contour, 1e-4, 1e-4); + printf("Found %zd eigenvalues:\n", result->neig); + for (size_t i = 0; i < result->neig; ++i) { + complex double eig = result->eigval[i]; + printf("%zd: %g%+gj\n", i, creal(eig), cimag(eig)); + } + free(contour); + beyn_result_free(result); + free(p.T0); + free(p.T1); + free(p.T2); + return 0; +} + + From 2fcc6282faa8b9aeb94116745d2d728593a107ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 26 Sep 2019 07:30:27 +0300 Subject: [PATCH 28/34] "Half-ellipse" integration contour. Former-commit-id: 800079fba837d65f84af62d5adfad177a6aa3cdf --- qpms/beyn.c | 73 +++++++++++++++++++++++++++++++++++++++++++++ qpms/beyn.h | 14 +++++++++ qpms/qpms_cdefs.pxd | 5 ++++ 3 files changed, 92 insertions(+) diff --git a/qpms/beyn.c b/qpms/beyn.c index b33291c..a24a59e 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -84,6 +84,79 @@ beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double r return c; } + +beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe, + double rIm, size_t n, beyn_contour_halfellipse_orientation or) +{ + beyn_contour_t *c; + QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + (n+1)*sizeof(c->z_dz[0]) + + sizeof(beyn_contour_halfellipse_orientation)); + c->centre = centre; + c->n = n; + const size_t nline = n/2; + const size_t narc = n - nline; + complex double faktor; + const double positive_zero = copysign(0., +1.); + const double negative_zero = copysign(0., -1.); + double l = rRe, h = rIm; + switch(or) { + case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS: + faktor = -I; + l = rIm; h = rRe; + break; + case BEYN_CONTOUR_HALFELLIPSE_RE_MINUS: + faktor = I; + l = rIm; h = rRe; + break; + case BEYN_CONTOUR_HALFELLIPSE_IM_PLUS: + faktor = 1; + break; + case BEYN_CONTOUR_HALFELLIPSE_IM_MINUS: + faktor = -1; + break; + default: QPMS_WTF; + } + + for(size_t i = 0; i < narc; ++i) { + double t = (i+0.5)*M_PI/narc; + double st = sin(t), ct = cos(t); + c->z_dz[i][0] = centre + faktor*(ct*l + I*st*h); + c->z_dz[i][1] = faktor * (I*l*st + h*ct) / narc; + } + for(size_t i = 0; i < nline; ++i) { + double t = 0.5 * (1 - (double) nline) + i; + complex double z = centre + faktor * t * l; + switch(or) { // ensure correct zero signs; CHECKME!!! + case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS: + if(creal(z) == 0 && signbit(creal(z))) + z = positive_zero + I * cimag(z); + break; + case BEYN_CONTOUR_HALFELLIPSE_RE_MINUS: + if(creal(z) == 0 && !signbit(creal(z))) + z = negative_zero + I * cimag(z); + break; + case BEYN_CONTOUR_HALFELLIPSE_IM_PLUS: + if(cimag(z) == 0 && signbit(cimag(z))) + z = creal(z) + I * positive_zero; + break; + case BEYN_CONTOUR_HALFELLIPSE_IM_MINUS: + if(cimag(z) == 0 && !signbit(cimag(z))) + z = creal(z) + I * negative_zero; + break; + default: QPMS_WTF; + } + c->z_dz[narc + i][0] = z; + c->z_dz[narc + i][1] = faktor * l / narc; + } + // We hide the half-axes metadata after the discretisation points. + c->z_dz[n][0] = rRe; + c->z_dz[n][1] = rIm; + // ugly... + *((beyn_contour_halfellipse_orientation *) &c->z_dz[n+1][0]) = or; + c->inside_test = NULL; // TODO beyn_contour_halfellipse_inside_test; + return c; +} + void beyn_result_gsl_free(beyn_result_gsl_t *r) { if(r) { gsl_vector_complex_free(r->eigval); diff --git a/qpms/beyn.h b/qpms/beyn.h index 5eeccf9..c822ecd 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -46,6 +46,20 @@ typedef struct beyn_contour_t { /** Free using free(). */ beyn_contour_t *beyn_contour_ellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints); + +typedef enum { + BEYN_CONTOUR_HALFELLIPSE_RE_PLUS = 3, + BEYN_CONTOUR_HALFELLIPSE_RE_MINUS = 1, + BEYN_CONTOUR_HALFELLIPSE_IM_PLUS = 0, + BEYN_CONTOUR_HALFELLIPSE_IM_MINUS = 2, +} beyn_contour_halfellipse_orientation; + + +/// Complex plane "half-elliptic" integration contour with axes parallel to the real, imaginary axes. +/** Free using free(). */ +beyn_contour_t *beyn_contour_halfellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints, + 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?). diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 1e7712e..4e4c634 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -562,6 +562,11 @@ cdef extern from "beyn.h": 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) From b3c3eeb2a28a2011c047551904566eeacb5df020 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 26 Sep 2019 10:06:31 +0300 Subject: [PATCH 29/34] Fix derivatives in integration contours. Former-commit-id: 7575a69a82eb19126aaac9aede9f170815b6015a --- qpms/beyn.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index a24a59e..bbc9f2b 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -75,7 +75,7 @@ beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double r double t = i*2*M_PI/n; double st = sin(t), ct = cos(t); c->z_dz[i][0] = centre + ct*rRe + I*st*rIm; - c->z_dz[i][1] = (I*rRe*st + rIm*ct) / n; + c->z_dz[i][1] = (-rRe*st + I*rIm*ct) * (2*M_PI / n); } // We hide the half-axes metadata after the discretisation points. c->z_dz[n][0] = rRe; @@ -121,7 +121,7 @@ beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe, double t = (i+0.5)*M_PI/narc; double st = sin(t), ct = cos(t); c->z_dz[i][0] = centre + faktor*(ct*l + I*st*h); - c->z_dz[i][1] = faktor * (I*l*st + h*ct) / narc; + c->z_dz[i][1] = faktor * (-l*st + I*h*ct) * (M_PI / narc); } for(size_t i = 0; i < nline; ++i) { double t = 0.5 * (1 - (double) nline) + i; @@ -146,7 +146,7 @@ beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe, default: QPMS_WTF; } c->z_dz[narc + i][0] = z; - c->z_dz[narc + i][1] = faktor * l / narc; + c->z_dz[narc + i][1] = faktor * l / nline; } // We hide the half-axes metadata after the discretisation points. c->z_dz[n][0] = rRe; From 2e7b925be2fc53ae7debee2f6a6fd825ca549389 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 26 Sep 2019 11:36:41 +0300 Subject: [PATCH 30/34] WIP trying to fix the half-ellipse contour. Former-commit-id: cbba6fffacbb38322a590478899f38c7d7dafe8e --- qpms/beyn.c | 4 ++-- tests/CMakeLists.txt | 5 +++++ tests/tbeyn3.c | 11 +++++++++-- 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index bbc9f2b..9bb0974 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -125,7 +125,7 @@ beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe, } for(size_t i = 0; i < nline; ++i) { double t = 0.5 * (1 - (double) nline) + i; - complex double z = centre + faktor * t * l; + complex double z = centre + faktor * t * 2 * l / nline; switch(or) { // ensure correct zero signs; CHECKME!!! case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS: if(creal(z) == 0 && signbit(creal(z))) @@ -146,7 +146,7 @@ beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe, default: QPMS_WTF; } c->z_dz[narc + i][0] = z; - c->z_dz[narc + i][1] = faktor * l / nline; + c->z_dz[narc + i][1] = faktor * 2 * l / nline; } // We hide the half-axes metadata after the discretisation points. c->z_dz[n][0] = rRe; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8185d79..5d39526 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -28,6 +28,11 @@ target_link_libraries( tbeyn3a qpms gsl lapacke amos m ) target_include_directories( tbeyn3a PRIVATE .. ) target_compile_definitions( tbeyn3a PRIVATE VARIANTA ) +add_executable( tbeyn3a_implus tbeyn3.c ) +target_link_libraries( tbeyn3a_implus qpms gsl lapacke amos m ) +target_include_directories( tbeyn3a_implus PRIVATE .. ) +target_compile_definitions( tbeyn3a_implus PRIVATE VARIANTA IMPLUS ) + add_executable( tbeyn3b tbeyn3.c ) target_link_libraries( tbeyn3b qpms gsl lapacke amos m ) target_include_directories( tbeyn3b PRIVATE .. ) diff --git a/tests/tbeyn3.c b/tests/tbeyn3.c index 0c983b6..fd0260e 100644 --- a/tests/tbeyn3.c +++ b/tests/tbeyn3.c @@ -1,8 +1,10 @@ +#define _GNU_SOURCE #include #include #include #include #include +#include static double randU(double a, double b) {return a + (b-a) * random() * (1. / RAND_MAX); } // Normal distribution via Box-Muller transform @@ -45,13 +47,14 @@ int M_function(complex double *target, const size_t m, const complex double z, v } int main(int argc, char **argv) { - complex double z0 = 0; + feenableexcept(FE_INVALID | FE_OVERFLOW); + complex double z0 = 0+3e-1*I; #ifdef RXSMALL double Rx = .1; #else double Rx = .3; // Variant B will fail in this case due to large number of eigenvalues (>30) #endif - double Ry = .3; + double Ry = .25; #ifdef VARIANTF int L = 10, N = 150, dim = 10; #else @@ -59,7 +62,11 @@ int main(int argc, char **argv) { #endif if (argc > 1) N = atoi(argv[1]); if (argc > 2) L = atoi(argv[2]); +#ifdef IMPLUS + beyn_contour_t *contour = beyn_contour_halfellipse(z0, Rx, Ry, N, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS); +#else beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); +#endif struct param p; p.T0 = malloc(dim*dim*sizeof(double)); p.T1 = malloc(dim*dim*sizeof(double)); From d6084f3264cfb149fc9f41c3e699b6c9e6bfced5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 29 Sep 2019 22:32:53 +0300 Subject: [PATCH 31/34] WIP Alternative integration contour (untested) Former-commit-id: c82da58ac33364797d7175a69feed475184937d9 --- qpms/beyn.c | 142 +++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 119 insertions(+), 23 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 9bb0974..9f6b081 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -19,6 +19,7 @@ #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); @@ -85,6 +86,36 @@ beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double r } +// Sets correct sign to zero for a given branch cut orientation +static inline complex double +align_zero(complex double z, beyn_contour_halfellipse_orientation or) +{ + // Maybe redundant, TODO check the standard. + const double positive_zero = copysign(0., +1.); + const double negative_zero = copysign(0., -1.); + switch(or) { // ensure correct zero signs; CHECKME!!! + case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS: + if(creal(z) == 0 && signbit(creal(z))) + z = positive_zero + I * cimag(z); + break; + case BEYN_CONTOUR_HALFELLIPSE_RE_MINUS: + if(creal(z) == 0 && !signbit(creal(z))) + z = negative_zero + I * cimag(z); + break; + case BEYN_CONTOUR_HALFELLIPSE_IM_PLUS: + if(cimag(z) == 0 && signbit(cimag(z))) + z = creal(z) + I * positive_zero; + break; + case BEYN_CONTOUR_HALFELLIPSE_IM_MINUS: + if(cimag(z) == 0 && !signbit(cimag(z))) + z = creal(z) + I * negative_zero; + break; + default: QPMS_WTF; + } + return z; +} + + beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe, double rIm, size_t n, beyn_contour_halfellipse_orientation or) { @@ -96,8 +127,6 @@ beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe, const size_t nline = n/2; const size_t narc = n - nline; complex double faktor; - const double positive_zero = copysign(0., +1.); - const double negative_zero = copysign(0., -1.); double l = rRe, h = rIm; switch(or) { case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS: @@ -125,27 +154,7 @@ beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe, } for(size_t i = 0; i < nline; ++i) { double t = 0.5 * (1 - (double) nline) + i; - complex double z = centre + faktor * t * 2 * l / nline; - switch(or) { // ensure correct zero signs; CHECKME!!! - case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS: - if(creal(z) == 0 && signbit(creal(z))) - z = positive_zero + I * cimag(z); - break; - case BEYN_CONTOUR_HALFELLIPSE_RE_MINUS: - if(creal(z) == 0 && !signbit(creal(z))) - z = negative_zero + I * cimag(z); - break; - case BEYN_CONTOUR_HALFELLIPSE_IM_PLUS: - if(cimag(z) == 0 && signbit(cimag(z))) - z = creal(z) + I * positive_zero; - break; - case BEYN_CONTOUR_HALFELLIPSE_IM_MINUS: - if(cimag(z) == 0 && !signbit(cimag(z))) - z = creal(z) + I * negative_zero; - break; - default: QPMS_WTF; - } - c->z_dz[narc + i][0] = z; + c->z_dz[narc + i][0] = align_zero(centre + faktor * t * 2 * l / nline, or); c->z_dz[narc + i][1] = faktor * 2 * l / nline; } // We hide the half-axes metadata after the discretisation points. @@ -157,6 +166,93 @@ beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe, return c; } +beyn_contour_t *beyn_contour_kidney(complex double centre, double rRe, + double rIm, const double rounding, const size_t n, beyn_contour_halfellipse_orientation or) +{ + beyn_contour_t *c; + QPMS_ENSURE(rounding >= 0 && rounding < .5, "rounding must lie in the interval [0, 0.5)"); + QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + (n+1)*sizeof(c->z_dz[0]) + + sizeof(beyn_contour_halfellipse_orientation)); + c->centre = centre; + c->n = n; + complex double faktor; + double l = rRe, h = rIm; + switch(or) { + case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS: + faktor = -I; + l = rIm; h = rRe; + break; + case BEYN_CONTOUR_HALFELLIPSE_RE_MINUS: + faktor = I; + l = rIm; h = rRe; + break; + case BEYN_CONTOUR_HALFELLIPSE_IM_PLUS: + faktor = 1; + break; + case BEYN_CONTOUR_HALFELLIPSE_IM_MINUS: + faktor = -1; + break; + default: QPMS_WTF; + } + + // Small circle centre coordinates. + const double y = rounding * h; // distance from the cut / straight line + const double x = sqrt(SQ(h - y) - SQ(y)); + + const double alpha = asin(y/(h-y)); + const double ar = l/h; // aspect ratio + + // Parameter range (equal to the contour length if ar == 1) + const double tmax = 2 * (x + (M_PI_2 + alpha) * y + (M_PI_2 - alpha) * h); + const double dt = tmax / n; + + size_t i = 0; + double t; + // Straight line, first part + double t_lo = 0, t_hi = x; + for(; t = i * dt, t <= t_hi; ++i) { + c->z_dz[i][0] = align_zero(centre + (t - t_lo) * ar * faktor, or); + c->z_dz[i][1] = dt * ar * faktor; + } + // First small arc + t_lo = t_hi; t_hi = t_lo + (M_PI_2 + alpha) * y; + for(; t = i * dt, dt < t_hi; ++i) { + double phi = (t - t_lo) / y - M_PI_2; + c->z_dz[i][0] = centre + (ar * (x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor; + c->z_dz[i][1] = dt * y * (- ar * sin(phi) + cos(phi) * I) * faktor; + } + // Big arc + t_lo = t_hi; t_hi = t_lo + (M_PI + 2 * alpha) * h; + for(; t = i * dt, dt < t_hi; ++i) { + double phi = (t - t_lo) / h + M_PI_2 + alpha; + c->z_dz[i][0] = centre + (ar * (h * cos(phi)) + h * (1 + sin(phi)) * I) * faktor; + c->z_dz[i][1] = dt * h * (- ar * sin(phi) + cos(phi) * I) * faktor; + } + // Second small arc + t_lo = t_hi; t_hi = t_lo + (M_PI_2 + alpha) * y; + for(; t = i * dt, dt < t_hi; ++i) { + double phi = (t - t_lo) / y + M_PI - alpha; + c->z_dz[i][0] = centre + (ar * (x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor; + c->z_dz[i][1] = dt * y * (- ar * sin(phi) + cos(phi) * I) * faktor; + } + // Straight line, second part + t_lo = t_hi; t_hi = tmax; + for(; t = i * dt, i < n; ++i) { + c->z_dz[i][0] = align_zero(centre + (t - t_lo - x) * ar * faktor, or); + c->z_dz[i][1] = dt * ar * faktor; + } + +#if 0 // TODO later + // We hide the half-axes metadata after the discretisation points. + c->z_dz[n][0] = rRe; + c->z_dz[n][1] = rIm; + // ugly... + *((beyn_contour_halfellipse_orientation *) &c->z_dz[n+1][0]) = or; +#endif + c->inside_test = NULL; // TODO beyn_contour_halfellipse_inside_test; + return c; +} + void beyn_result_gsl_free(beyn_result_gsl_t *r) { if(r) { gsl_vector_complex_free(r->eigval); From f940e62e5272a1df5dcb39a79181146fdcc9b621 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 30 Sep 2019 08:53:57 +0300 Subject: [PATCH 32/34] Test and fix the rounded half-ellipse contour.. Former-commit-id: 448c30d375b3f9d0abab9442ae00bfcd123e5cb9 --- qpms/beyn.c | 20 ++++++++++---------- qpms/beyn.h | 6 ++++++ tests/CMakeLists.txt | 4 ++++ tests/kidneycontour.c | 19 +++++++++++++++++++ 4 files changed, 39 insertions(+), 10 deletions(-) create mode 100644 tests/kidneycontour.c diff --git a/qpms/beyn.c b/qpms/beyn.c index 9f6b081..28911be 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -216,24 +216,24 @@ beyn_contour_t *beyn_contour_kidney(complex double centre, double rRe, } // First small arc t_lo = t_hi; t_hi = t_lo + (M_PI_2 + alpha) * y; - for(; t = i * dt, dt < t_hi; ++i) { + for(; t = i * dt, t < t_hi; ++i) { double phi = (t - t_lo) / y - M_PI_2; c->z_dz[i][0] = centre + (ar * (x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor; - c->z_dz[i][1] = dt * y * (- ar * sin(phi) + cos(phi) * I) * faktor; + c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor; } // Big arc - t_lo = t_hi; t_hi = t_lo + (M_PI + 2 * alpha) * h; - for(; t = i * dt, dt < t_hi; ++i) { - double phi = (t - t_lo) / h + M_PI_2 + alpha; - c->z_dz[i][0] = centre + (ar * (h * cos(phi)) + h * (1 + sin(phi)) * I) * faktor; - c->z_dz[i][1] = dt * h * (- ar * sin(phi) + cos(phi) * I) * faktor; + t_lo = t_hi; t_hi = t_lo + (M_PI - 2 * alpha) * h; + for(; t = i * dt, t < t_hi; ++i) { + double phi = (t - t_lo) / h + alpha; + c->z_dz[i][0] = centre + (ar * (h * cos(phi)) + h * sin(phi) * I) * faktor; + c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor; } // Second small arc t_lo = t_hi; t_hi = t_lo + (M_PI_2 + alpha) * y; - for(; t = i * dt, dt < t_hi; ++i) { + for(; t = i * dt, t < t_hi; ++i) { double phi = (t - t_lo) / y + M_PI - alpha; - c->z_dz[i][0] = centre + (ar * (x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor; - c->z_dz[i][1] = dt * y * (- ar * sin(phi) + cos(phi) * I) * faktor; + c->z_dz[i][0] = centre + (ar * (- x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor; + c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor; } // Straight line, second part t_lo = t_hi; t_hi = tmax; diff --git a/qpms/beyn.h b/qpms/beyn.h index c822ecd..25bc1c7 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -60,6 +60,12 @@ typedef enum { beyn_contour_t *beyn_contour_halfellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints, beyn_contour_halfellipse_orientation or); +/// Similar to halfellipse but with rounded corners. +beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, double halfax_im, + double rounding, ///< Must be in interval [0, 0.5) + 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?). diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5d39526..a8513dc 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -15,6 +15,10 @@ add_executable( test_scalar_ewald32 test_scalar_ewald32.c ) target_link_libraries( test_scalar_ewald32 qpms gsl lapacke amos m ) target_include_directories( test_scalar_ewald32 PRIVATE .. ) +add_executable( kidneycontour kidneycontour.c ) +target_link_libraries( kidneycontour qpms gsl lapacke amos m ) +target_include_directories( kidneycontour PRIVATE .. ) + add_executable( tbeyn tbeyn.c ) target_link_libraries( tbeyn qpms gsl lapacke amos m ) target_include_directories( tbeyn PRIVATE .. ) diff --git a/tests/kidneycontour.c b/tests/kidneycontour.c new file mode 100644 index 0000000..8f4e635 --- /dev/null +++ b/tests/kidneycontour.c @@ -0,0 +1,19 @@ +#include +#include + +#define CPAIR(x) creal(x), cimag(x) + +int main(int argc, char **argv) +{ + double rRe = 2e3, rIm = 1.5e3, rounding = 0.2; + complex double centre = 1e3 * I; + size_t n = 100; + + beyn_contour_t *c = beyn_contour_kidney(centre, rRe, rIm, rounding, n, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS); + + for(size_t i = 0; i < n; ++i) + printf("%g\t%g\t%g\t%g\n", CPAIR(c->z_dz[i][0]), CPAIR(c->z_dz[i][1])); + + free(c); + return 0; +} From dad170b8febee65c389b6292b41f375e65f22c37 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 30 Sep 2019 10:54:10 +0300 Subject: [PATCH 33/34] Kidney contour basic test Former-commit-id: 861cf5f9a8409db360e48db74274f20069b94391 --- tests/CMakeLists.txt | 5 +++++ tests/tbeyn3.c | 2 ++ 2 files changed, 7 insertions(+) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index a8513dc..ad15e0f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -37,6 +37,11 @@ target_link_libraries( tbeyn3a_implus qpms gsl lapacke amos m ) target_include_directories( tbeyn3a_implus PRIVATE .. ) target_compile_definitions( tbeyn3a_implus PRIVATE VARIANTA IMPLUS ) +add_executable( tbeyn3a_kidney tbeyn3.c ) +target_link_libraries( tbeyn3a_kidney qpms gsl lapacke amos m ) +target_include_directories( tbeyn3a_kidney PRIVATE .. ) +target_compile_definitions( tbeyn3a_kidney PRIVATE VARIANTA IMPLUS_KIDNEY ) + add_executable( tbeyn3b tbeyn3.c ) target_link_libraries( tbeyn3b qpms gsl lapacke amos m ) target_include_directories( tbeyn3b PRIVATE .. ) diff --git a/tests/tbeyn3.c b/tests/tbeyn3.c index fd0260e..52ede9a 100644 --- a/tests/tbeyn3.c +++ b/tests/tbeyn3.c @@ -64,6 +64,8 @@ int main(int argc, char **argv) { if (argc > 2) L = atoi(argv[2]); #ifdef IMPLUS beyn_contour_t *contour = beyn_contour_halfellipse(z0, Rx, Ry, N, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS); +#elif defined IMPLUS_KIDNEY + beyn_contour_t *contour = beyn_contour_kidney(z0, Rx, Ry, 0.3, N, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS); #else beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); #endif From 79e6c47f94c13e5b39e541ec741e13bb365e476d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 1 Oct 2019 11:52:29 +0300 Subject: [PATCH 34/34] Beyn integration contours to pxd Former-commit-id: d33de9a694ea7d5600a627062c7119e8483d753e --- qpms/qpms_cdefs.pxd | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 4e4c634..75770ba 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -584,6 +584,13 @@ cdef extern from "beyn.h": beyn_function_M_inv_Vhat_t M_inv_Vhat, void *params, const beyn_contour_t *contour, double rank_tol, double res_tol) + beyn_contour_t *beyn_contour_ellipse(cdouble centre, double halfax_re, double halfax_im, size_t npoints) + beyn_contour_t *beyn_contour_halfellipse(cdouble centre, double halfax_re, double halfax_im, size_t npoints, + beyn_contour_halfellipse_orientation ori) + 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)