From fe406e230d6f723df5f4b06151842e08841db8b4 Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Mon, 18 Dec 2017 11:26:42 -0500 Subject: [PATCH] block with a single vector case is working with IRBL --- .../ImplicitlyRestartedBlockLanczos.h | 811 +++++++++++------ .../ImplicitlyRestartedBlockLanczos.h.bak | 835 ++++++++++++++++++ .../iterative/ImplicitlyRestartedLanczos.h | 67 +- .../ImplicitlyRestartedLanczos.h.bak | 625 +++++++++++++ tests/lanczos/Test_dwf_block_lanczos.cc | 18 +- tests/lanczos/Test_dwf_lanczos.cc | 10 +- 6 files changed, 2085 insertions(+), 281 deletions(-) create mode 100644 lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h.bak create mode 100644 lib/algorithms/iterative/ImplicitlyRestartedLanczos.h.bak diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index 92acb67e..8d7a41c4 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -78,6 +78,7 @@ public: : _Linop(Linop), _poly(poly), Nstop(_Nstop), Nu(_Nu), Nk(_Nk), Nm(_Nm), Nblock_m(_Nm/_Nu), Nblock_k(_Nk/_Nu), + //eresid(_eresid), MaxIter(10), eresid(_eresid), MaxIter(_MaxIter), diagonalisation(_diagonalisation) { assert( (Nk%Nu==0) && (Nm%Nu==0) ); }; @@ -139,12 +140,10 @@ until convergence clog <<" -- total Nm = "<< Nm <<" vectors"<< std::endl; clog <<" -- size of eval = "<< eval.size() << std::endl; clog <<" -- size of evec = "<< evec.size() << std::endl; - if ( diagonalisation == IRLdiagonaliseWithDSTEGR ) { - clog << "Diagonalisation is DSTEGR "<< std::endl; - } else if ( diagonalisation == IRLdiagonaliseWithQR ) { - clog << "Diagonalisation is QR "<< std::endl; - } else if ( diagonalisation == IRLdiagonaliseWithEigen ) { + if ( diagonalisation == IRLdiagonaliseWithEigen ) { clog << "Diagonalisation is Eigen "<< std::endl; + } else { + abort(); } clog << std::string(74,'*') << std::endl; @@ -157,6 +156,7 @@ until convergence std::vector eval2(Nm); Eigen::MatrixXcd Qt = Eigen::MatrixXcd::Zero(Nm,Nm); + Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm); std::vector Iconv(Nm); std::vector B(Nm,grid); // waste of space replicating @@ -165,14 +165,11 @@ until convergence std::vector f_copy(Nu,grid); Field v(grid); - int k1 = 1; - int k2 = Nk; - Nconv = 0; RealD beta_k; - // Set initial vector + // set initial vector for (int i=0; i=Nstop ){ goto converged; } + } // end of iter loop clog <<"**************************************************************************"<< std::endl; @@ -308,7 +396,6 @@ until convergence clog << " -- beta(k) = "<< beta_k << "\n"; clog << " -- Nconv = "<< Nconv << "\n"; clog <<"**************************************************************************"<< std::endl; -#endif } private: @@ -341,31 +428,33 @@ private: Real beta; - clog << "A: b = " << b << std::endl; // 3. wk:=Avk−βkv_{k−1} for (int k=L, u=0; k0) { - clog << "B: b = " << b << std::endl; for (int u=0; u>& lmd, + std::vector>& lme, + int Nu, int Nb, int Nk, int Nm, + Eigen::MatrixXcd& M) + { + //clog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; + assert( Nk%Nu == 0 && Nm%Nu == 0 ); + assert( Nk <= Nm ); + + // rearrange + for ( int u=0; u QRD(Mtmp); + Q = QRD.householderQ(); + for (int j=0; j QRD(Mtmp); + Q = QRD.householderQ(); + R = QRD.matrixQR(); // upper triangular part is the R matrix. + // lower triangular part used to represent series + // of Q sequence. + + // equivalent operation of Qprod *= Q + //M = Eigen::MatrixXcd::Zero(Nm,Nm); + + //for (int i=0; i Nm) kmax = Nm; + // for (int k=i; ki) M(i,j) = conj(M(j,i)); + if (i-j > Nu || j-i > Nu) M(i,j) = 0.; + } + } + + //clog << "shiftedQRDecompEigen() end" << endl; + } +#endif +#if 0 + void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nu, int Nm, + RealD Dsh, + Eigen::MatrixXcd& Qprod) + { + //clog << "shiftedQRDecompEigen() begin" << '\n'; + Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm); + Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm); + + Mtmp = M; + for (int i=0; i QRD(Mtmp); + Q = QRD.householderQ(); + + M = Q.adjoint()*(M*Q); + for (int i=0; i QRD2(Mtmp); + Qprod = QRD2.householderQ(); + + Mtmp -= Qprod; + clog << "Frobenius norm ||Qprod(after) - Qprod|| = " << Mtmp.norm() << std::endl; +#endif + //clog << "shiftedQRDecompEigen() end" << endl; + } +#endif +#if 0 + void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nm, + RealD Dsh, + Eigen::MatrixXcd& Qprod) + { + //clog << "shiftedQRDecompEigen() begin" << '\n'; + Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm); + //Eigen::MatrixXcd Qtmp = Eigen::MatrixXcd::Zero(Nm,Nm); + + Mtmp = Qprod.adjoint()*(M*Qprod); + for (int i=0; i QRD(Mtmp); + //Qtmp = Qprod*QRD.householderQ(); + + //Eigen::HouseholderQR QRD2(Qtmp); + //Qprod = QRD2.householderQ(); + + Qprod *= QRD.householderQ(); + //ComplexD p; + //RealD r; + + //r = 0.; + //for (int k=0; k //memset + +#define clog std::cout << GridLogMessage + +namespace Grid { + +///////////////////////////////////////////////////////////// +// Implicitly restarted block lanczos +///////////////////////////////////////////////////////////// +template +class ImplicitlyRestartedBlockLanczos { + +private: + + std::string cname = std::string("ImplicitlyRestartedBlockLanczos"); + int MaxIter; // Max iterations + int Nstop; // Number of evecs checked for convergence + int Nu; // Numbeer of vecs in the unit block + int Nk; // Number of converged sought + int Nm; // total number of vectors + int Nblock_k; // Nk/Nu + int Nblock_m; // Nm/Nu + RealD eresid; + IRLdiagonalisation diagonalisation; + //////////////////////////////////// + // Embedded objects + //////////////////////////////////// + SortEigen _sort; + LinearOperatorBase &_Linop; + OperatorFunction &_poly; + + ///////////////////////// + // Constructor + ///////////////////////// +public: + ImplicitlyRestartedBlockLanczos(LinearOperatorBase &Linop, // op + OperatorFunction & poly, // polynomial + int _Nstop, // really sought vecs + int _Nu, // vecs in the unit block + int _Nk, // sought vecs + int _Nm, // total vecs + RealD _eresid, // resid in lmd deficit + int _MaxIter, // Max iterations + IRLdiagonalisation _diagonalisation = IRLdiagonaliseWithEigen) + : _Linop(Linop), _poly(poly), + Nstop(_Nstop), Nu(_Nu), Nk(_Nk), Nm(_Nm), + Nblock_m(_Nm/_Nu), Nblock_k(_Nk/_Nu), + //eresid(_eresid), MaxIter(10), + eresid(_eresid), MaxIter(_MaxIter), + diagonalisation(_diagonalisation) + { assert( (Nk%Nu==0) && (Nm%Nu==0) ); }; + + //////////////////////////////// + // Helpers + //////////////////////////////// + static RealD normalize(Field& v) + { + RealD nn = norm2(v); + nn = sqrt(nn); + v = v * (1.0/nn); + return nn; + } + + void orthogonalize(Field& w, std::vector& evec, int k) + { + typedef typename Field::scalar_type MyComplex; + MyComplex ip; + + for(int j=0; j K P = M − K † +Compute the factorization AVM = VM HM + fM eM +repeat + Q=I + for i = 1,...,P do + QiRi =HM −θiI Q = QQi + H M = Q †i H M Q i + end for + βK =HM(K+1,K) σK =Q(M,K) + r=vK+1βK +rσK + VK =VM(1:M)Q(1:M,1:K) + HK =HM(1:K,1:K) + →AVK =VKHK +fKe†K † Extend to an M = K + P step factorization AVM = VMHM + fMeM +until convergence +*/ + void calc(std::vector& eval, + std::vector& evec, + const std::vector& src, int& Nconv) + { + std::string fname = std::string(cname+"::calc()"); + GridBase *grid = evec[0]._grid; + assert(grid == src[0]._grid); + assert( Nu = src.size() ); + + clog << std::string(74,'*') << std::endl; + clog << fname + " starting iteration 0 / "<< MaxIter<< std::endl; + clog << std::string(74,'*') << std::endl; + clog <<" -- seek Nk = "<< Nk <<" vectors"<< std::endl; + clog <<" -- accept Nstop = "<< Nstop <<" vectors"<< std::endl; + clog <<" -- total Nm = "<< Nm <<" vectors"<< std::endl; + clog <<" -- size of eval = "<< eval.size() << std::endl; + clog <<" -- size of evec = "<< evec.size() << std::endl; + if ( diagonalisation == IRLdiagonaliseWithEigen ) { + clog << "Diagonalisation is Eigen "<< std::endl; + } else { + abort(); + } + clog << std::string(74,'*') << std::endl; + + assert(Nm == evec.size() && Nm == eval.size()); + + std::vector> lmd(Nu,std::vector(Nm,0.0)); + std::vector> lme(Nu,std::vector(Nm,0.0)); + std::vector> lmd2(Nu,std::vector(Nm,0.0)); + std::vector> lme2(Nu,std::vector(Nm,0.0)); + std::vector eval2(Nm); + + Eigen::MatrixXcd Qt = Eigen::MatrixXcd::Zero(Nm,Nm); + Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm); + + std::vector Iconv(Nm); + std::vector B(Nm,grid); // waste of space replicating + + std::vector f(Nu,grid); + std::vector f_copy(Nu,grid); + Field v(grid); + + Nconv = 0; + + RealD beta_k; + + // set initial vector + for (int i=0; i=Nstop ){ + goto converged; + } + + } // end of iter loop + + clog <<"**************************************************************************"<< std::endl; + std::cout<< GridLogError << fname + " NOT converged."; + clog <<"**************************************************************************"<< std::endl; + abort(); + + converged: + // Sorting + eval.resize(Nconv); + evec.resize(Nconv,grid); + for(int i=0; i>& lmd, + std::vector>& lme, + std::vector& evec, + std::vector& w, + std::vector& w_copy, + int b) + { + const RealD tiny = 1.0e-20; + + int Nu = w.size(); + int Nm = evec.size(); + assert( b < Nm/Nu ); + + // converts block index to full indicies for an interval [L,R) + int L = Nu*b; + int R = Nu*(b+1); + + Real beta; + + // 3. wk:=Avk−βkv_{k−1} + for (int k=L, u=0; k0) { + for (int u=0; u>& lmd, + std::vector>& lme, + int Nu, int Nb, int Nk, int Nm, + Eigen::MatrixXcd& M) + { + //clog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; + assert( Nk%Nu == 0 && Nm%Nu == 0 ); + assert( Nk <= Nm ); + + // rearrange + for ( int u=0; u QRD(Mtmp); +// Q = QRD.householderQ(); +// +// M = Q.adjoint()*(M*Q); +//#if 0 +// Qprod *= Q; +//#else +// Mtmp = Qprod*Q; +// +// Eigen::HouseholderQR QRD2(Mtmp); +// Qprod = QRD2.householderQ(); +// +// Mtmp -= Qprod; +// clog << "Frobenius norm ||Qprod(after) - Qprod|| = " << Mtmp.norm() << std::endl; +//#endif +// //clog << "shiftedQRDecompEigen() end" << endl; +// } + void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nm, + RealD Dsh, + Eigen::MatrixXcd& Qprod) + { + //clog << "shiftedQRDecompEigen() begin" << '\n'; + Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm); + //Eigen::MatrixXcd Qtmp = Eigen::MatrixXcd::Zero(Nm,Nm); + + Mtmp = Qprod.adjoint()*(M*Qprod); + for (int i=0; i QRD(Mtmp); + //Qtmp = Qprod*QRD.householderQ(); + + //Eigen::HouseholderQR QRD2(Qtmp); + //Qprod = QRD2.householderQ(); + + Qprod *= QRD.householderQ(); + //ComplexD p; + //RealD r; + + //r = 0.; + //for (int k=0; k0) w -= lme[k-1] * evec[k-1]; + if(k>0) { + w -= lme[k-1] * evec[k-1]; + //clog << "ckpt A (k= " << k << ")" << '\n'; + //clog << "lme = " << lme[k-1] << '\n'; + //clog << "norm(w) = " << norm2(w) << std::endl; + } ComplexD zalph = innerProduct(evec[k],w); // 4. αk:=(wk,vk) RealD alph = real(zalph); + //clog << "ckpt B (k= " << k << ")" << '\n'; + //clog << "lmd = " << alph << std::endl; w = w - alph * evec[k];// 5. wk:=wk−αkvk @@ -622,4 +680,5 @@ void diagonalize_lapack(std::vector& lmd, }; } +#undef clog #endif diff --git a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h.bak b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h.bak new file mode 100644 index 00000000..a8723f32 --- /dev/null +++ b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h.bak @@ -0,0 +1,625 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/algorithms/iterative/ImplicitlyRestartedLanczos.h + + Copyright (C) 2015 + +Author: Peter Boyle +Author: Chulwoo Jung +Author: Guido Cossu + + This program 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. + + This program 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., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#ifndef GRID_IRL_H +#define GRID_IRL_H + +#include //memset + +namespace Grid { + + enum IRLdiagonalisation { + IRLdiagonaliseWithDSTEGR, + IRLdiagonaliseWithQR, + IRLdiagonaliseWithEigen + }; + +//////////////////////////////////////////////////////////////////////////////// +// Helper class for sorting the evalues AND evectors by Field +// Use pointer swizzle on vectors +//////////////////////////////////////////////////////////////////////////////// +template +class SortEigen { + private: + static bool less_lmd(RealD left,RealD right){ + return left > right; + } + static bool less_pair(std::pair& left, + std::pair& right){ + return left.first > (right.first); + } + + public: + void push(std::vector& lmd,std::vector& evec,int N) { + + //////////////////////////////////////////////////////////////////////// + // PAB: FIXME: VERY VERY VERY wasteful: takes a copy of the entire vector set. + // : The vector reorder should be done by pointer swizzle somehow + //////////////////////////////////////////////////////////////////////// + std::vector cpy(lmd.size(),evec[0]._grid); + for(int i=0;i > emod(lmd.size()); + + for(int i=0;i(lmd[i],&cpy[i]); + + partial_sort(emod.begin(),emod.begin()+N,emod.end(),less_pair); + + typename std::vector >::iterator it = emod.begin(); + for(int i=0;ifirst; + evec[i]=*(it->second); + ++it; + } + } + void push(std::vector& lmd,int N) { + std::partial_sort(lmd.begin(),lmd.begin()+N,lmd.end(),less_lmd); + } + bool saturated(RealD lmd, RealD thrs) { + return fabs(lmd) > fabs(thrs); + } +}; + +///////////////////////////////////////////////////////////// +// Implicitly restarted lanczos +///////////////////////////////////////////////////////////// +template +class ImplicitlyRestartedLanczos { + +private: + + int MaxIter; // Max iterations + int Nstop; // Number of evecs checked for convergence + int Nk; // Number of converged sought + int Nm; // Nm -- total number of vectors + RealD eresid; + IRLdiagonalisation diagonalisation; + //////////////////////////////////// + // Embedded objects + //////////////////////////////////// + SortEigen _sort; + LinearOperatorBase &_Linop; + OperatorFunction &_poly; + + ///////////////////////// + // Constructor + ///////////////////////// +public: + ImplicitlyRestartedLanczos(LinearOperatorBase &Linop, // op + OperatorFunction & poly, // polynomial + int _Nstop, // really sought vecs + int _Nk, // sought vecs + int _Nm, // total vecs + RealD _eresid, // resid in lmd deficit + int _MaxIter, // Max iterations + IRLdiagonalisation _diagonalisation= IRLdiagonaliseWithEigen ) : + _Linop(Linop), _poly(poly), + Nstop(_Nstop), Nk(_Nk), Nm(_Nm), + eresid(_eresid), MaxIter(_MaxIter), + diagonalisation(_diagonalisation) + { }; + + //////////////////////////////// + // Helpers + //////////////////////////////// + static RealD normalise(Field& v) + { + RealD nn = norm2(v); + nn = sqrt(nn); + v = v * (1.0/nn); + return nn; + } + + void orthogonalize(Field& w, std::vector& evec, int k) + { + typedef typename Field::scalar_type MyComplex; + MyComplex ip; + + for(int j=0; j K P = M − K † +Compute the factorization AVM = VM HM + fM eM +repeat + Q=I + for i = 1,...,P do + QiRi =HM −θiI Q = QQi + H M = Q †i H M Q i + end for + βK =HM(K+1,K) σK =Q(M,K) + r=vK+1βK +rσK + VK =VM(1:M)Q(1:M,1:K) + HK =HM(1:K,1:K) + →AVK =VKHK +fKe†K † Extend to an M = K + P step factorization AVM = VMHM + fMeM +until convergence +*/ + void calc(std::vector& eval, std::vector& evec, const Field& src, int& Nconv) + { + + GridBase *grid = evec[0]._grid; + assert(grid == src._grid); + + std::cout << GridLogMessage <<"**************************************************************************"<< std::endl; + std::cout << GridLogMessage <<" ImplicitlyRestartedLanczos::calc() starting iteration 0 / "<< MaxIter<< std::endl; + std::cout << GridLogMessage <<"**************************************************************************"<< std::endl; + std::cout << GridLogMessage <<" -- seek Nk = " << Nk <<" vectors"<< std::endl; + std::cout << GridLogMessage <<" -- accept Nstop = " << Nstop <<" vectors"<< std::endl; + std::cout << GridLogMessage <<" -- total Nm = " << Nm <<" vectors"<< std::endl; + std::cout << GridLogMessage <<" -- size of eval = " << eval.size() << std::endl; + std::cout << GridLogMessage <<" -- size of evec = " << evec.size() << std::endl; + if ( diagonalisation == IRLdiagonaliseWithDSTEGR ) { + std::cout << GridLogMessage << "Diagonalisation is DSTEGR "< lme(Nm); + std::vector lme2(Nm); + std::vector eval2(Nm); + + Eigen::MatrixXd Qt = Eigen::MatrixXd::Zero(Nm,Nm); + + std::vector Iconv(Nm); + std::vector B(Nm,grid); // waste of space replicating + + Field f(grid); + Field v(grid); + + int k1 = 1; + int k2 = Nk; + + Nconv = 0; + + RealD beta_k; + + // Set initial vector + evec[0] = src; + std::cout << GridLogMessage <<"norm2(src)= " << norm2(src)<=Nstop ){ + goto converged; + } + } // end of iter loop + + std::cout << GridLogMessage <<"**************************************************************************"<< std::endl; + std::cout<< GridLogError <<" ImplicitlyRestartedLanczos::calc() NOT converged."; + std::cout << GridLogMessage <<"**************************************************************************"<< std::endl; + abort(); + + converged: + // Sorting + eval.resize(Nconv); + evec.resize(Nconv,grid); + for(int i=0; i& lmd, + std::vector& lme, + std::vector& evec, + Field& w,int Nm,int k) + { + const RealD tiny = 1.0e-20; + assert( k< Nm ); + + _poly(_Linop,evec[k],w); // 3. wk:=Avk−βkv_{k−1} + + if(k>0) w -= lme[k-1] * evec[k-1]; + + ComplexD zalph = innerProduct(evec[k],w); // 4. αk:=(wk,vk) + RealD alph = real(zalph); + + w = w - alph * evec[k];// 5. wk:=wk−αkvk + + RealD beta = normalise(w); // 6. βk+1 := ∥wk∥2. If βk+1 = 0 then Stop + // 7. vk+1 := wk/βk+1 + + lmd[k] = alph; + lme[k] = beta; + + if ( k > 0 ) orthogonalize(w,evec,k); // orthonormalise + if ( k < Nm-1) evec[k+1] = w; + + if ( beta < tiny ) std::cout << GridLogMessage << " beta is tiny "<& lmd, std::vector& lme, + int Nk, int Nm, + Eigen::MatrixXd & Qt, // Nm x Nm + GridBase *grid) + { + Eigen::MatrixXd TriDiag = Eigen::MatrixXd::Zero(Nk,Nk); + + for(int i=0;i eigensolver(TriDiag); + + for (int i = 0; i < Nk; i++) { + lmd[Nk-1-i] = eigensolver.eigenvalues()(i); + } + for (int i = 0; i < Nk; i++) { + for (int j = 0; j < Nk; j++) { + Qt(Nk-1-i,j) = eigensolver.eigenvectors()(j,i); + } + } + } + /////////////////////////////////////////////////////////////////////////// + // File could end here if settle on Eigen ??? + /////////////////////////////////////////////////////////////////////////// + + void qr_decomp(std::vector& lmd, // Nm + std::vector& lme, // Nm + int Nk, int Nm, // Nk, Nm + Eigen::MatrixXd& Qt, // Nm x Nm matrix + RealD Dsh, int kmin, int kmax) + { + int k = kmin-1; + RealD x; + + RealD Fden = 1.0/hypot(lmd[k]-Dsh,lme[k]); + RealD c = ( lmd[k] -Dsh) *Fden; + RealD s = -lme[k] *Fden; + + RealD tmpa1 = lmd[k]; + RealD tmpa2 = lmd[k+1]; + RealD tmpb = lme[k]; + + lmd[k] = c*c*tmpa1 +s*s*tmpa2 -2.0*c*s*tmpb; + lmd[k+1] = s*s*tmpa1 +c*c*tmpa2 +2.0*c*s*tmpb; + lme[k] = c*s*(tmpa1-tmpa2) +(c*c-s*s)*tmpb; + x =-s*lme[k+1]; + lme[k+1] = c*lme[k+1]; + + for(int i=0; i& lmd, std::vector& lme, + int Nk, int Nm, + Eigen::MatrixXd & Qt, + GridBase *grid) + { + Qt = Eigen::MatrixXd::Identity(Nm,Nm); + if ( diagonalisation == IRLdiagonaliseWithDSTEGR ) { + diagonalize_lapack(lmd,lme,Nk,Nm,Qt,grid); + } else if ( diagonalisation == IRLdiagonaliseWithQR ) { + diagonalize_QR(lmd,lme,Nk,Nm,Qt,grid); + } else if ( diagonalisation == IRLdiagonaliseWithEigen ) { + diagonalize_Eigen(lmd,lme,Nk,Nm,Qt,grid); + } else { + assert(0); + } + } + +#ifdef USE_LAPACK +void LAPACK_dstegr(char *jobz, char *range, int *n, double *d, double *e, + double *vl, double *vu, int *il, int *iu, double *abstol, + int *m, double *w, double *z, int *ldz, int *isuppz, + double *work, int *lwork, int *iwork, int *liwork, + int *info); +#endif + +void diagonalize_lapack(std::vector& lmd, + std::vector& lme, + int Nk, int Nm, + Eigen::MatrixXd& Qt, + GridBase *grid) +{ +#ifdef USE_LAPACK + const int size = Nm; + int NN = Nk; + double evals_tmp[NN]; + double evec_tmp[NN][NN]; + memset(evec_tmp[0],0,sizeof(double)*NN*NN); + double DD[NN]; + double EE[NN]; + for (int i = 0; i< NN; i++) { + for (int j = i - 1; j <= i + 1; j++) { + if ( j < NN && j >= 0 ) { + if (i==j) DD[i] = lmd[i]; + if (i==j) evals_tmp[i] = lmd[i]; + if (j==(i-1)) EE[j] = lme[j]; + } + } + } + int evals_found; + int lwork = ( (18*NN) > (1+4*NN+NN*NN)? (18*NN):(1+4*NN+NN*NN)) ; + int liwork = 3+NN*10 ; + int iwork[liwork]; + double work[lwork]; + int isuppz[2*NN]; + char jobz = 'V'; // calculate evals & evecs + char range = 'I'; // calculate all evals + // char range = 'A'; // calculate all evals + char uplo = 'U'; // refer to upper half of original matrix + char compz = 'I'; // Compute eigenvectors of tridiagonal matrix + int ifail[NN]; + int info; + int total = grid->_Nprocessors; + int node = grid->_processor; + int interval = (NN/total)+1; + double vl = 0.0, vu = 0.0; + int il = interval*node+1 , iu = interval*(node+1); + if (iu > NN) iu=NN; + double tol = 0.0; + if (1) { + memset(evals_tmp,0,sizeof(double)*NN); + if ( il <= NN){ + LAPACK_dstegr(&jobz, &range, &NN, + (double*)DD, (double*)EE, + &vl, &vu, &il, &iu, // these four are ignored if second parameteris 'A' + &tol, // tolerance + &evals_found, evals_tmp, (double*)evec_tmp, &NN, + isuppz, + work, &lwork, iwork, &liwork, + &info); + for (int i = iu-1; i>= il-1; i--){ + evals_tmp[i] = evals_tmp[i - (il-1)]; + if (il>1) evals_tmp[i-(il-1)]=0.; + for (int j = 0; j< NN; j++){ + evec_tmp[i][j] = evec_tmp[i - (il-1)][j]; + if (il>1) evec_tmp[i-(il-1)][j]=0.; + } + } + } + { + grid->GlobalSumVector(evals_tmp,NN); + grid->GlobalSumVector((double*)evec_tmp,NN*NN); + } + } + // Safer to sort instead of just reversing it, + // but the document of the routine says evals are sorted in increasing order. + // qr gives evals in decreasing order. + for(int i=0;i& lmd, std::vector& lme, + int Nk, int Nm, + Eigen::MatrixXd & Qt, + GridBase *grid) + { + int Niter = 100*Nm; + int kmin = 1; + int kmax = Nk; + + // (this should be more sophisticated) + for(int iter=0; iter= kmin; --j){ + RealD dds = fabs(lmd[j-1])+fabs(lmd[j]); + if(fabs(lme[j-1])+dds > dds){ + kmax = j+1; + goto continued; + } + } + Niter = iter; + return; + + continued: + for(int j=0; j dds){ + kmin = j+1; + break; + } + } + } + std::cout << GridLogError << "[QL method] Error - Too many iteration: "< seeds5({5,6,7,8}); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); - //GridParallelRNG RNG5rb(FrbGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG5rb(FrbGrid); RNG5.SeedFixedIntegers(seeds5); // ypj [note] why seed RNG5 again? bug? In this case, run with a default seed(). - GridParallelRNG RNG5rb(FrbGrid); //RNG5rb.SeedFixedIntegers(seeds5); + //GridParallelRNG RNG5rb(FrbGrid); //RNG5rb.SeedFixedIntegers(seeds5); LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4, Umu); @@ -77,19 +77,19 @@ int main (int argc, char ** argv) SchurDiagTwoOperator HermOp(Ddwf); // SchurDiagMooeeOperator HermOp(Ddwf); - const int Nstop = 30; - const int Nu = 4; - const int Nk = 60; - const int Np = 60; + const int Nstop = 50; + const int Nu = 1; + const int Nk = 200; + const int Np = 200; const int Nm = Nk+Np; - const int MaxIt= 10000; + const int MaxIt= 10; RealD resid = 1.0e-8; //std::vector Coeffs { 0.,-1.}; // ypj [note] this may not be supported by some compilers - std::vector Coeffs({ 0.,1.}); + std::vector Coeffs({ 0.,-1.}); Polynomial PolyX(Coeffs); - Chebyshev Cheb(0.2,5.,11); + Chebyshev Cheb(0.2,5.5,11); // ChebyshevLanczos Cheb(9.,1.,0.,20); // Cheb.csv(std::cout); // exit(-24); diff --git a/tests/lanczos/Test_dwf_lanczos.cc b/tests/lanczos/Test_dwf_lanczos.cc index 1dd5dae3..6c65608d 100644 --- a/tests/lanczos/Test_dwf_lanczos.cc +++ b/tests/lanczos/Test_dwf_lanczos.cc @@ -75,16 +75,16 @@ int main (int argc, char ** argv) SchurDiagTwoOperator HermOp(Ddwf); // SchurDiagMooeeOperator HermOp(Ddwf); - const int Nstop = 30; - const int Nk = 40; - const int Np = 40; + const int Nstop = 50; + const int Nk = 200; + const int Np = 200; const int Nm = Nk+Np; - const int MaxIt= 10000; + const int MaxIt= 100; RealD resid = 1.0e-8; std::vector Coeffs { 0.,-1.}; Polynomial PolyX(Coeffs); - Chebyshev Cheb(0.2,5.,11); + Chebyshev Cheb(0.2,5.5,11); // ChebyshevLanczos Cheb(9.,1.,0.,20); // Cheb.csv(std::cout); // exit(-24);