From 91cc33e90710ade34ed8d794c8d7a0bdcbb844ab Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Sun, 5 Nov 2017 23:46:03 -0500 Subject: [PATCH 01/15] IRBL development test bed is added; bootstrap.sh is updated to avoid repeated download of the EIGEN package after the first build. --- .gitignore | 6 ++ bootstrap.sh | 15 ++-- tests/lanczos/Test_dwf_block_lanczos.cc | 107 ++++++++++++++++++++++++ tests/lanczos/Test_wilson_lanczos.cc | 2 +- 4 files changed, 124 insertions(+), 6 deletions(-) create mode 100644 tests/lanczos/Test_dwf_block_lanczos.cc diff --git a/.gitignore b/.gitignore index d743ee06..44e063b9 100644 --- a/.gitignore +++ b/.gitignore @@ -88,11 +88,17 @@ Thumbs.db ################### build*/* +# bootstrap # +############# +*.tar.bz2* + # IDE related files # ##################### *.xcodeproj/* build.sh .vscode +.ctags +tags # Eigen source # ################ diff --git a/bootstrap.sh b/bootstrap.sh index ac27ef25..1510e3a1 100755 --- a/bootstrap.sh +++ b/bootstrap.sh @@ -1,11 +1,16 @@ #!/usr/bin/env bash -EIGEN_URL='http://bitbucket.org/eigen/eigen/get/3.3.3.tar.bz2' +EIGEN_SRC='3.3.3.tar.bz2' +EIGEN_URL="http://bitbucket.org/eigen/eigen/get/${EIGEN_SRC}" -echo "-- deploying Eigen source..." -wget ${EIGEN_URL} --no-check-certificate -./scripts/update_eigen.sh `basename ${EIGEN_URL}` -#rm `basename ${EIGEN_URL}` +if [ -f ${EIGEN_SRC} ]; then + echo "-- skip deploying Eigen source..." +else + echo "-- deploying Eigen source..." + wget ${EIGEN_URL} --no-check-certificate + ./scripts/update_eigen.sh `basename ${EIGEN_URL}` + #rm `basename ${EIGEN_URL}` +fi echo '-- generating Make.inc files...' ./scripts/filelist diff --git a/tests/lanczos/Test_dwf_block_lanczos.cc b/tests/lanczos/Test_dwf_block_lanczos.cc new file mode 100644 index 00000000..da721efa --- /dev/null +++ b/tests/lanczos/Test_dwf_block_lanczos.cc @@ -0,0 +1,107 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_lanczos.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 */ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +typedef typename GparityDomainWallFermionR::FermionField FermionField; + +RealD AllZero(RealD x){ return 0.;} + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls=8; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n",UGrid,UrbGrid,FGrid,FrbGrid); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG RNG5rb(FrbGrid); RNG5.SeedFixedIntegers(seeds5); + + LatticeGaugeField Umu(UGrid); + SU3::HotConfiguration(RNG4, Umu); + + std::vector U(4,UGrid); + for(int mu=0;mu(Umu,mu); + } + + RealD mass=0.01; + RealD M5=1.8; + RealD mob_b=1.5; +// DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + GparityMobiusFermionD ::ImplParams params; + std::vector twists({1,1,1,0}); + params.twists = twists; + GparityMobiusFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params); + +// MdagMLinearOperator HermOp(Ddwf); +// SchurDiagTwoOperator HermOp(Ddwf); + SchurDiagTwoOperator HermOp(Ddwf); +// SchurDiagMooeeOperator HermOp(Ddwf); + + const int Nstop = 30; + const int Nk = 40; + const int Np = 40; + const int Nm = Nk+Np; + const int MaxIt= 10000; + RealD resid = 1.0e-8; + + std::vector Coeffs { 0.,-1.}; + Polynomial PolyX(Coeffs); + Chebyshev Cheb(0.2,5.,11); +// ChebyshevLanczos Cheb(9.,1.,0.,20); +// Cheb.csv(std::cout); +// exit(-24); + ImplicitlyRestartedLanczos IRL(HermOp,Cheb,Nstop,Nk,Nm,resid,MaxIt); + + + std::vector eval(Nm); + FermionField src(FrbGrid); + gaussian(RNG5rb,src); + std::vector evec(Nm,FrbGrid); + for(int i=0;i<1;i++){ + std::cout << GridLogMessage < Date: Mon, 27 Nov 2017 13:39:52 -0500 Subject: [PATCH 02/15] fix vector assign bug --- lib/algorithms/approx/Chebyshev.h | 8 ++++++-- tests/lanczos/Test_dwf_block_lanczos.cc | 14 +++++++++----- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/lib/algorithms/approx/Chebyshev.h b/lib/algorithms/approx/Chebyshev.h index f8c21a05..e8b3e113 100644 --- a/lib/algorithms/approx/Chebyshev.h +++ b/lib/algorithms/approx/Chebyshev.h @@ -84,7 +84,8 @@ namespace Grid { public: void csv(std::ostream &out){ RealD diff = hi-lo; - for (RealD x=lo-0.2*diff; x 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); LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4, Umu); @@ -76,13 +78,15 @@ int main (int argc, char ** argv) // SchurDiagMooeeOperator HermOp(Ddwf); const int Nstop = 30; - const int Nk = 40; - const int Np = 40; + const int Nk = 60; + const int Np = 60; const int Nm = Nk+Np; const int MaxIt= 10000; RealD resid = 1.0e-8; - std::vector Coeffs { 0.,-1.}; + //std::vector Coeffs { 0.,-1.}; + // ypj [note] this may not be supported by some compilers + std::vector Coeffs({ 0.,1.}); Polynomial PolyX(Coeffs); Chebyshev Cheb(0.2,5.,11); // ChebyshevLanczos Cheb(9.,1.,0.,20); @@ -100,7 +104,7 @@ int main (int argc, char ** argv) }; int Nconv; - //IRL.calc(eval,evec,src,Nconv); + IRL.calc(eval,evec,src,Nconv); Grid_finalize(); From 5139eaf4917334ba88c6d0c9f292129ed321af90 Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Sun, 3 Dec 2017 23:55:22 -0500 Subject: [PATCH 03/15] block Lanczos construction is added. --- lib/algorithms/Algorithms.h | 1 + .../ImplicitlyRestartedBlockLanczos.h | 700 ++++++++++++++++++ tests/lanczos/Test_dwf_block_lanczos.cc | 17 +- 3 files changed, 711 insertions(+), 7 deletions(-) create mode 100644 lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h diff --git a/lib/algorithms/Algorithms.h b/lib/algorithms/Algorithms.h index e6a3f530..28bdfe24 100644 --- a/lib/algorithms/Algorithms.h +++ b/lib/algorithms/Algorithms.h @@ -53,6 +53,7 @@ Author: Peter Boyle #include #include #include +#include #include #include #include diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h new file mode 100644 index 00000000..92acb67e --- /dev/null +++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -0,0 +1,700 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h + + Copyright (C) 2015 + +Author: Peter Boyle +Author: Chulwoo Jung +Author: Yong-Chull Jang +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_IRBL_H +#define GRID_IRBL_H + +#include //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(_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 == IRLdiagonaliseWithDSTEGR ) { + clog << "Diagonalisation is DSTEGR "<< std::endl; + } else if ( diagonalisation == IRLdiagonaliseWithQR ) { + clog << "Diagonalisation is QR "<< std::endl; + } else if ( diagonalisation == IRLdiagonaliseWithEigen ) { + clog << "Diagonalisation is Eigen "<< std::endl; + } + 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); + + 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); + + int k1 = 1; + int k2 = Nk; + + 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; + + 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 HermOp(Ddwf); const int Nstop = 30; + const int Nu = 4; const int Nk = 60; const int Np = 60; const int Nm = Nk+Np; @@ -92,19 +93,21 @@ int main (int argc, char ** argv) // ChebyshevLanczos Cheb(9.,1.,0.,20); // Cheb.csv(std::cout); // exit(-24); - ImplicitlyRestartedLanczos IRL(HermOp,Cheb,Nstop,Nk,Nm,resid,MaxIt); + ImplicitlyRestartedBlockLanczos IRBL(HermOp,Cheb,Nstop,Nu,Nk,Nm,resid,MaxIt); std::vector eval(Nm); - FermionField src(FrbGrid); - gaussian(RNG5rb,src); + + std::vector src(Nu,FrbGrid); + for ( int i=0; i evec(Nm,FrbGrid); - for(int i=0;i<1;i++){ - std::cout << GridLogMessage < Date: Mon, 18 Dec 2017 11:26:42 -0500 Subject: [PATCH 04/15] 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); From 89c4e9b16863fcd684404cc405d26676aefc5331 Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Thu, 21 Dec 2017 23:13:39 -0500 Subject: [PATCH 05/15] first complete version of IRBL; requires practical test and clean up --- .../ImplicitlyRestartedBlockLanczos.h | 31 ++++++++++--------- tests/lanczos/Test_dwf_block_lanczos.cc | 6 ++-- tests/lanczos/Test_dwf_lanczos.cc | 8 ++--- 3 files changed, 24 insertions(+), 21 deletions(-) diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index 8d7a41c4..694f14fc 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -198,18 +198,7 @@ until convergence // clog << "ckpt A2: lmd[" << k << "] = " << lmd[0][k] << '\n'; //} - // residual vector -#if 1 // ypj[fixme] temporary to check a case when block has one vector - for ( int i=0; i HermOp(Ddwf); // SchurDiagMooeeOperator HermOp(Ddwf); - const int Nstop = 50; + const int Nstop = 120; const int Nu = 1; - const int Nk = 200; - const int Np = 200; + const int Nk = 240; + const int Np = 240; const int Nm = Nk+Np; const int MaxIt= 10; RealD resid = 1.0e-8; diff --git a/tests/lanczos/Test_dwf_lanczos.cc b/tests/lanczos/Test_dwf_lanczos.cc index 6c65608d..58219007 100644 --- a/tests/lanczos/Test_dwf_lanczos.cc +++ b/tests/lanczos/Test_dwf_lanczos.cc @@ -75,11 +75,11 @@ int main (int argc, char ** argv) SchurDiagTwoOperator HermOp(Ddwf); // SchurDiagMooeeOperator HermOp(Ddwf); - const int Nstop = 50; - const int Nk = 200; - const int Np = 200; + const int Nstop = 120; + const int Nk = 240; + const int Np = 240; const int Nm = Nk+Np; - const int MaxIt= 100; + const int MaxIt= 10; RealD resid = 1.0e-8; std::vector Coeffs { 0.,-1.}; From 3cb8cb72824edbf084785c785278d6861eb87387 Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Sat, 23 Dec 2017 14:54:07 -0500 Subject: [PATCH 06/15] 'typename' is added to compile with AVX512 using GCC7.2.0; a semicolon was missing in Grid_avx512.h and the bug is fixed. Option SKL is added to configure script for skylake processor specific AVX512 operations. Code can be compiled with --enable-simd=SKL using GCC 7.2.0, but Test_simd fails. AVX512 support for complex double type with non-intel compilers makes this error; it needs a review. --- lib/algorithms/iterative/ImplicitlyRestartedLanczos.h | 2 +- lib/qcd/action/fermion/CayleyFermion5Dvec.cc | 6 ++++-- lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc | 3 ++- lib/qcd/action/fermion/MobiusEOFAFermionvec.cc | 3 ++- lib/simd/Grid_avx512.h | 2 +- 5 files changed, 10 insertions(+), 6 deletions(-) diff --git a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h index b20afcb9..7541a544 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h @@ -267,7 +267,7 @@ until convergence Qt = Eigen::MatrixXd::Identity(Nm,Nm); for(int ip=k2; ip(Umu,mu); + } +#endif + + //RealD mass=0.01; + RealD mass=0.00107; RealD M5=1.8; RealD mob_b=1.5; // DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); @@ -77,10 +92,10 @@ int main (int argc, char ** argv) SchurDiagTwoOperator HermOp(Ddwf); // SchurDiagMooeeOperator HermOp(Ddwf); - const int Nstop = 120; - const int Nu = 1; - const int Nk = 240; - const int Np = 240; + const int Nstop = 100; + const int Nu = 10; + const int Nk = 200; + const int Np = 200; const int Nm = Nk+Np; const int MaxIt= 10; RealD resid = 1.0e-8; @@ -89,7 +104,8 @@ int main (int argc, char ** argv) // ypj [note] this may not be supported by some compilers std::vector Coeffs({ 0.,-1.}); Polynomial PolyX(Coeffs); - Chebyshev Cheb(0.2,5.5,11); + //Chebyshev Cheb(0.2,5.5,11); + Chebyshev Cheb(0.03,5.5,21); // ChebyshevLanczos Cheb(9.,1.,0.,20); // Cheb.csv(std::cout); // exit(-24); From bfc0306a439725c5884324c5292f1dbbf73db27a Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Fri, 29 Dec 2017 00:58:24 -0500 Subject: [PATCH 08/15] job parameters can be provided by cmdarg. --- lib/util/Init.cc | 8 ++ lib/util/Init.h | 3 + tests/lanczos/Test_dwf_block_lanczos.cc | 170 +++++++++++++++++++----- 3 files changed, 145 insertions(+), 36 deletions(-) diff --git a/lib/util/Init.cc b/lib/util/Init.cc index 1266d34d..3bdd52ad 100644 --- a/lib/util/Init.cc +++ b/lib/util/Init.cc @@ -157,6 +157,14 @@ void GridCmdOptionInt(std::string &str,int & val) return; } +// ypj [add] +void GridCmdOptionFloat(std::string &str,double & val) +{ + std::stringstream ss(str); + ss>>val; + return; +} + void GridParseLayout(char **argv,int argc, std::vector &latt, diff --git a/lib/util/Init.h b/lib/util/Init.h index 3da00742..e71c95c6 100644 --- a/lib/util/Init.h +++ b/lib/util/Init.h @@ -54,6 +54,9 @@ namespace Grid { std::string GridCmdVectorIntToString(const std::vector & vec); void GridCmdOptionCSL(std::string str,std::vector & vec); void GridCmdOptionIntVector(std::string &str,std::vector & vec); + // ypj [add] + void GridCmdOptionInt(std::string &str,int & val); + void GridCmdOptionFloat(std::string &str,double & val); void GridParseLayout(char **argv,int argc, diff --git a/tests/lanczos/Test_dwf_block_lanczos.cc b/tests/lanczos/Test_dwf_block_lanczos.cc index 64b72c98..71879d25 100644 --- a/tests/lanczos/Test_dwf_block_lanczos.cc +++ b/tests/lanczos/Test_dwf_block_lanczos.cc @@ -35,17 +35,125 @@ typedef typename GparityDomainWallFermionR::FermionField FermionField; RealD AllZero(RealD x){ return 0.;} +class CmdJobParams +{ + public: + std::string gaugefile; + + int Ls; + double mass; + double M5; + double mob_b; + + int Nu; + int Nk; + int Np; + int Nstop; + int MaxIter; + double resid; + + double low; + double high; + int order; + + CmdJobParams() + : gaugefile("Hot"), + Ls(8), mass(0.01), M5(1.8), mob_b(1.5), + Nu(4), Nk(200), Np(200), Nstop(100), MaxIter(10), resid(1.0e-8), + low(0.2), high(5.5), order(11) + {}; + + void Parse(char **argv, int argc); +}; + + +void CmdJobParams::Parse(char **argv,int argc) +{ + std::string arg; + std::vector vi; + + if( GridCmdOptionExists(argv,argv+argc,"--gconf") ){ + gaugefile = GridCmdOptionPayload(argv,argv+argc,"--gconf"); + } + + if( GridCmdOptionExists(argv,argv+argc,"--Ls") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--Ls"); + GridCmdOptionInt(arg,Ls); + } + + if( GridCmdOptionExists(argv,argv+argc,"--mass") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--mass"); + GridCmdOptionFloat(arg,mass); + } + + if( GridCmdOptionExists(argv,argv+argc,"--M5") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--M5"); + GridCmdOptionFloat(arg,M5); + } + + if( GridCmdOptionExists(argv,argv+argc,"--mob_b") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--mob_b"); + GridCmdOptionFloat(arg,mob_b); + } + + if( GridCmdOptionExists(argv,argv+argc,"--irbl") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--irbl"); + GridCmdOptionIntVector(arg,vi); + Nu = vi[0]; + Nk = vi[1]; + Np = vi[2]; + Nstop = vi[3]; + MaxIter = vi[4]; + } + + if( GridCmdOptionExists(argv,argv+argc,"--resid") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--resid"); + GridCmdOptionFloat(arg,resid); + } + + if( GridCmdOptionExists(argv,argv+argc,"--cheby_l") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_l"); + GridCmdOptionFloat(arg,low); + } + + if( GridCmdOptionExists(argv,argv+argc,"--cheby_u") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_u"); + GridCmdOptionFloat(arg,high); + } + + if( GridCmdOptionExists(argv,argv+argc,"--cheby_n") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_n"); + GridCmdOptionInt(arg,order); + } + + if ( CartesianCommunicator::RankWorld() == 0 ) { + clog <<" Gauge Configuration "<< gaugefile << '\n'; + clog <<" Ls "<< Ls << '\n'; + clog <<" mass "<< mass << '\n'; + clog <<" M5 "<< M5 << '\n'; + clog <<" mob_b "<< mob_b << '\n'; + clog <<" Nu "<< Nu << '\n'; + clog <<" Nk "<< Nk << '\n'; + clog <<" Np "<< Np << '\n'; + clog <<" Nstop "<< Nstop << '\n'; + clog <<" MaxIter "<< MaxIter << '\n'; + clog <<" resid "<< resid << '\n'; + clog <<" Cheby Poly "<< low << "," << high << "," << order << std::endl; + } +} + + int main (int argc, char ** argv) { Grid_init(&argc,&argv); - - const int Ls=8; - //const int Ls=16; + + CmdJobParams JP; + JP.Parse(argv,argc); GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); - GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); - GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(JP.Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(JP.Ls,UGrid); printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n",UGrid,UrbGrid,FGrid,FrbGrid); std::vector seeds4({1,2,3,4}); @@ -56,31 +164,23 @@ int main (int argc, char ** argv) // ypj [note] why seed RNG5 again? bug? In this case, run with a default seed(). //GridParallelRNG RNG5rb(FrbGrid); //RNG5rb.SeedFixedIntegers(seeds5); -#if 0 // ypj [note] generate a random configuration - LatticeGaugeField Umu(UGrid); - SU3::HotConfiguration(RNG4, Umu); - - std::vector U(4,UGrid); - for(int mu=0;mu(Umu,mu); - } -#else // read configuration from a file LatticeGaugeField Umu(UGrid); std::vector U(4,UGrid); - FieldMetaData header; - std::string file("./ckpoint_lat"); - NerscIO::readConfiguration(Umu,header,file); - + if ( JP.gaugefile.compare("Hot") == 0 ) { + SU3::HotConfiguration(RNG4, Umu); + } else { + FieldMetaData header; + NerscIO::readConfiguration(Umu,header,JP.gaugefile); + } + for(int mu=0;mu(Umu,mu); } -#endif - //RealD mass=0.01; - RealD mass=0.00107; - RealD M5=1.8; - RealD mob_b=1.5; + RealD mass = JP.mass; + RealD M5 = JP.M5; + RealD mob_b = JP.mob_b; // DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); GparityMobiusFermionD ::ImplParams params; std::vector twists({1,1,1,0}); @@ -92,27 +192,25 @@ int main (int argc, char ** argv) SchurDiagTwoOperator HermOp(Ddwf); // SchurDiagMooeeOperator HermOp(Ddwf); - const int Nstop = 100; - const int Nu = 10; - const int Nk = 200; - const int Np = 200; - const int Nm = Nk+Np; - const int MaxIt= 10; - RealD resid = 1.0e-8; + int Nu = JP.Nu; + int Nk = JP.Nk; + int Nm = Nk+JP.Np; //std::vector Coeffs { 0.,-1.}; // ypj [note] this may not be supported by some compilers std::vector Coeffs({ 0.,-1.}); Polynomial PolyX(Coeffs); //Chebyshev Cheb(0.2,5.5,11); - Chebyshev Cheb(0.03,5.5,21); -// ChebyshevLanczos Cheb(9.,1.,0.,20); + Chebyshev Cheb(JP.low,JP.high,JP.order); // Cheb.csv(std::cout); -// exit(-24); - ImplicitlyRestartedBlockLanczos IRBL(HermOp,Cheb,Nstop,Nu,Nk,Nm,resid,MaxIt); - + ImplicitlyRestartedBlockLanczos IRBL(HermOp, + Cheb, + JP.Nstop, + Nu,Nk,Nm, + JP.resid, + JP.MaxIter); - std::vector eval(Nm); + std::vector eval(Nm); std::vector src(Nu,FrbGrid); for ( int i=0; i Date: Fri, 29 Dec 2017 23:26:17 -0500 Subject: [PATCH 09/15] force hermicity to the block alpha and force diagonal of the block beta to be real --- .../ImplicitlyRestartedBlockLanczos.h | 61 +++++++++++++------ 1 file changed, 42 insertions(+), 19 deletions(-) diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index 694f14fc..ff8050b7 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -105,6 +105,18 @@ public: } normalize(w); } + + void orthogonalize_blockhead(Field& w, std::vector& evec, int k, int Nu) + { + typedef typename Field::scalar_type MyComplex; + MyComplex ip; + + for(int j=0; j0) { for (int u=0; u=Nstop ){ goto converged; } + + if ( iter < MaxIter-1 ) { + if ( Nu == 1 ) { + // reconstruct initial vector for additional pole space + blockwiseStep(lmd,lme,evec,f,f_copy,Nblock_k-1); + } + else { + // update the first block + for (int i=0; i QRD(Mtmp); - Q = QRD.householderQ(); - for (int j=0; j Nm) kmax = Nm; - // for (int k=i; k 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.; - } + Mtmp(i,i) = real(Mtmp(i,i)) + Dsh; } - //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 = Mtmp; - 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; ki) M(i,j) = conj(M(j,i)); + // if (i-j > Nu || j-i > Nu) M(i,j) = 0.; // } - // r = 0.; - // for (int k=0; k +Author: Chulwoo Jung +Author: Yong-Chull Jang +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_IRBL_H +#define GRID_IRBL_H + +#include //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& evec, int k, int Nu) + { + 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(); + 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 vi; + double re,im; + int expect, idx; + std::string vstr; + std::ifstream pfile; if( GridCmdOptionExists(argv,argv+argc,"--gconf") ){ gaugefile = GridCmdOptionPayload(argv,argv+argc,"--gconf"); } - if( GridCmdOptionExists(argv,argv+argc,"--Ls") ){ - arg = GridCmdOptionPayload(argv,argv+argc,"--Ls"); - GridCmdOptionInt(arg,Ls); + if( GridCmdOptionExists(argv,argv+argc,"--phase") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--phase"); + pfile.open(arg); + assert(pfile); + expect = 0; + while( pfile >> vstr ) { + if ( vstr.compare("boundary_phase") == 0 ) { + pfile >> vstr; + GridCmdOptionInt(vstr,idx); + assert(expect==idx); + pfile >> vstr; + GridCmdOptionFloat(vstr,re); + pfile >> vstr; + GridCmdOptionFloat(vstr,im); + boundary_phase.push_back({re,im}); + expect++; + } + } + pfile.close(); + } else { + for (int i=0; i<4; ++i) boundary_phase.push_back({1.,0.}); + } + + if( GridCmdOptionExists(argv,argv+argc,"--omega") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--omega"); + pfile.open(arg); + assert(pfile); + Ls = 0; + while( pfile >> vstr ) { + if ( vstr.compare("omega") == 0 ) { + pfile >> vstr; + GridCmdOptionInt(vstr,idx); + assert(Ls==idx); + pfile >> vstr; + GridCmdOptionFloat(vstr,re); + pfile >> vstr; + GridCmdOptionFloat(vstr,im); + omega.push_back({re,im}); + Ls++; + } + } + pfile.close(); + } else { + if( GridCmdOptionExists(argv,argv+argc,"--Ls") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--Ls"); + GridCmdOptionInt(arg,Ls); + } } if( GridCmdOptionExists(argv,argv+argc,"--mass") ){ @@ -127,18 +178,25 @@ void CmdJobParams::Parse(char **argv,int argc) } if ( CartesianCommunicator::RankWorld() == 0 ) { - clog <<" Gauge Configuration "<< gaugefile << '\n'; - clog <<" Ls "<< Ls << '\n'; - clog <<" mass "<< mass << '\n'; - clog <<" M5 "<< M5 << '\n'; - clog <<" mob_b "<< mob_b << '\n'; - clog <<" Nu "<< Nu << '\n'; - clog <<" Nk "<< Nk << '\n'; - clog <<" Np "<< Np << '\n'; - clog <<" Nstop "<< Nstop << '\n'; - clog <<" MaxIter "<< MaxIter << '\n'; - clog <<" resid "<< resid << '\n'; - clog <<" Cheby Poly "<< low << "," << high << "," << order << std::endl; + std::streamsize ss = std::cout.precision(); + std::cout << GridLogMessage <<" Gauge Configuration "<< gaugefile << '\n'; + std::cout.precision(15); + for ( int i=0; i<4; ++i ) std::cout << GridLogMessage <<" boundary_phase["<< i << "] = " << boundary_phase[i] << '\n'; + std::cout.precision(ss); + std::cout << GridLogMessage <<" Ls "<< Ls << '\n'; + std::cout << GridLogMessage <<" mass "<< mass << '\n'; + std::cout << GridLogMessage <<" M5 "<< M5 << '\n'; + std::cout << GridLogMessage <<" mob_b "<< mob_b << '\n'; + std::cout.precision(15); + for ( int i=0; i twists({1,1,1,0}); - params.twists = twists; - GparityMobiusFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params); +// RealD mob_b = JP.mob_b; // Gparity +// std::vector omega; // ZMobius + +// GparityMobiusFermionD ::ImplParams params; +// std::vector twists({1,1,1,0}); +// params.twists = twists; +// GparityMobiusFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params); +// SchurDiagTwoOperator HermOp(Ddwf); -// MdagMLinearOperator HermOp(Ddwf); -// SchurDiagTwoOperator HermOp(Ddwf); - SchurDiagTwoOperator HermOp(Ddwf); -// SchurDiagMooeeOperator HermOp(Ddwf); + //WilsonFermionR::ImplParams params; + ZMobiusFermionR::ImplParams params; + params.overlapCommsCompute = true; + params.boundary_phases = JP.boundary_phase; + ZMobiusFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,JP.omega,1.,0.,params); + SchurDiagTwoOperator HermOp(Ddwf); int Nu = JP.Nu; int Nk = JP.Nk; @@ -217,7 +279,7 @@ int main (int argc, char ** argv) std::vector evec(Nm,FrbGrid); for(int i=0;i<1;++i){ - clog << i <<" / "<< Nm <<" grid pointer "<< evec[i]._grid << std::endl; + std::cout << GridLogMessage << i <<" / "<< Nm <<" grid pointer "<< evec[i]._grid << std::endl; }; int Nconv; From fbe1209f7e94a8a4fb0615820c2f176838ad9c20 Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Wed, 31 Jan 2018 12:10:24 -0500 Subject: [PATCH 11/15] count converged eigenvalues not assuming candidates are sorted --- .../iterative/ImplicitlyRestartedBlockLanczos.h | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index 933f2882..70c9afeb 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -166,6 +166,7 @@ until convergence std::vector> lmd2(Nu,std::vector(Nm,0.0)); std::vector> lme2(Nu,std::vector(Nm,0.0)); std::vector eval2(Nm); + std::vector resid(Nk); Eigen::MatrixXcd Qt = Eigen::MatrixXcd::Zero(Nm,Nm); Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm); @@ -267,6 +268,7 @@ until convergence eval2[i] = vnum/vden; v -= eval2[i]*B[i]; RealD vv = norm2(v); + resid[i] = vv; std::cout.precision(13); clog << "[" << std::setw(3)<< std::setiosflags(std::ios_base::right) < Date: Wed, 31 Jan 2018 12:40:33 -0500 Subject: [PATCH 12/15] relocate the printing block of converged eigenvalues --- .../iterative/ImplicitlyRestartedBlockLanczos.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index 70c9afeb..bf849a0f 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -285,10 +285,6 @@ until convergence } // i-loop end clog <<" #modes converged: "<=Nstop ){ - goto converged; - } for(int i=0; i=Nstop ){ + goto converged; + } if ( iter < MaxIter-1 ) { if ( Nu == 1 ) { From 3caf0e8b0969ef337312c8d9b17493a27752b00a Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Fri, 2 Feb 2018 18:43:37 -0500 Subject: [PATCH 13/15] print Ritz values before and after the shift; rollback to simple last lanczos vector rebuild; print formatting is changed --- .../ImplicitlyRestartedBlockLanczos.h | 183 ++++++++---------- 1 file changed, 76 insertions(+), 107 deletions(-) diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index bf849a0f..71aa2891 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -118,23 +118,6 @@ public: normalize(w); } -/* Rudy Arthur's thesis pp.137 ------------------------- -Require: M > 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) @@ -197,10 +180,7 @@ until convergence int iter; for(iter = 0; iter>& lmd, std::vector>& lme, std::vector& evec, From 11219a8f7a200e644160d16fd51cbf3b33c9763d Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Mon, 5 Mar 2018 18:12:42 -0500 Subject: [PATCH 14/15] add explicit restart method rbl --- .../ImplicitlyRestartedBlockLanczos.h | 340 ++++++++++++++---- tests/lanczos/Test_dwf_block_lanczos.cc | 57 ++- 2 files changed, 309 insertions(+), 88 deletions(-) diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index 71aa2891..2de35244 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -33,10 +33,12 @@ Author: Guido Cossu #include //memset -#define clog std::cout << GridLogMessage +#define Glog std::cout << GridLogMessage namespace Grid { +enum class LanczosType { irbl, rbl }; + ///////////////////////////////////////////////////////////// // Implicitly restarted block lanczos ///////////////////////////////////////////////////////////// @@ -53,6 +55,7 @@ private: int Nm; // total number of vectors int Nblock_k; // Nk/Nu int Nblock_m; // Nm/Nu + int Nconv_test_interval; // Number of skipped vectors when checking a convergence RealD eresid; IRLdiagonalisation diagonalisation; //////////////////////////////////// @@ -69,6 +72,7 @@ public: ImplicitlyRestartedBlockLanczos(LinearOperatorBase &Linop, // op OperatorFunction & poly, // polynomial int _Nstop, // really sought vecs + int _Nconv_test_interval, // conv check interval int _Nu, // vecs in the unit block int _Nk, // sought vecs int _Nm, // total vecs @@ -76,7 +80,8 @@ public: int _MaxIter, // Max iterations IRLdiagonalisation _diagonalisation = IRLdiagonaliseWithEigen) : _Linop(Linop), _poly(poly), - Nstop(_Nstop), Nu(_Nu), Nk(_Nk), Nm(_Nm), + Nstop(_Nstop), Nconv_test_interval(_Nconv_test_interval), + Nu(_Nu), Nk(_Nk), Nm(_Nm), Nblock_m(_Nm/_Nu), Nblock_k(_Nk/_Nu), //eresid(_eresid), MaxIter(10), eresid(_eresid), MaxIter(_MaxIter), @@ -117,30 +122,45 @@ public: } normalize(w); } - + void calc(std::vector& eval, std::vector& evec, - const std::vector& src, int& Nconv) + const std::vector& src, int& Nconv, LanczosType Impl) { - std::string fname = std::string(cname+"::calc()"); + switch (Impl) { + case LanczosType::irbl: + calc_irbl(eval,evec,src,Nconv); + break; + + case LanczosType::rbl: + calc_rbl(eval,evec,src,Nconv); + break; + } + } + + void calc_irbl(std::vector& eval, + std::vector& evec, + const std::vector& src, int& Nconv) + { + std::string fname = std::string(cname+"::calc_irbl()"); 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; + Glog << std::string(74,'*') << std::endl; + Glog << fname + " starting iteration 0 / "<< MaxIter<< std::endl; + Glog << std::string(74,'*') << std::endl; + Glog <<" -- seek Nk = "<< Nk <<" vectors"<< std::endl; + Glog <<" -- accept Nstop = "<< Nstop <<" vectors"<< std::endl; + Glog <<" -- total Nm = "<< Nm <<" vectors"<< std::endl; + Glog <<" -- size of eval = "<< eval.size() << std::endl; + Glog <<" -- size of evec = "<< evec.size() << std::endl; if ( diagonalisation == IRLdiagonaliseWithEigen ) { - clog << "Diagonalisation is Eigen "<< std::endl; + Glog << "Diagonalisation is Eigen "<< std::endl; } else { abort(); } - clog << std::string(74,'*') << std::endl; + Glog << std::string(74,'*') << std::endl; assert(Nm == evec.size() && Nm == eval.size()); @@ -167,10 +187,10 @@ public: // set initial vector for (int i=0; iNk ) { - clog <<" #Apply shifted QR transformations "<& eval, + std::vector& evec, + const std::vector& src, int& Nconv) + { + std::string fname = std::string(cname+"::calc_rbl()"); + GridBase *grid = evec[0]._grid; + assert(grid == src[0]._grid); + assert( Nu = src.size() ); + + int Np = (Nm-Nk); + if (Np > 0 && MaxIter > 1) Np /= MaxIter; + int Nblock_p = Np/Nu; + + Glog << std::string(74,'*') << std::endl; + Glog << fname + " starting iteration 0 / "<< MaxIter<< std::endl; + Glog << std::string(74,'*') << std::endl; + Glog <<" -- seek (min) Nk = "<< Nk <<" vectors"<< std::endl; + Glog <<" -- seek (inc) Np = "<< Np <<" vectors"<< std::endl; + Glog <<" -- seek (max) Nm = "<< Nm <<" vectors"<< std::endl; + Glog <<" -- accept Nstop = "<< Nstop <<" vectors"<< std::endl; + Glog <<" -- size of eval = "<< eval.size() << std::endl; + Glog <<" -- size of evec = "<< evec.size() << std::endl; + if ( diagonalisation == IRLdiagonaliseWithEigen ) { + Glog << "Diagonalisation is Eigen "<< std::endl; + } else { + abort(); + } + Glog << 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(Nk); + std::vector resid(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 omega; // ZMobius @@ -254,10 +285,6 @@ int main (int argc, char ** argv) ZMobiusFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,JP.omega,1.,0.,params); SchurDiagTwoOperator HermOp(Ddwf); - int Nu = JP.Nu; - int Nk = JP.Nk; - int Nm = Nk+JP.Np; - //std::vector Coeffs { 0.,-1.}; // ypj [note] this may not be supported by some compilers std::vector Coeffs({ 0.,-1.}); @@ -267,23 +294,23 @@ int main (int argc, char ** argv) // Cheb.csv(std::cout); ImplicitlyRestartedBlockLanczos IRBL(HermOp, Cheb, - JP.Nstop, - Nu,Nk,Nm, + JP.Nstop, JP.Ntest, + JP.Nu, JP.Nk, JP.Nm, JP.resid, JP.MaxIter); - std::vector eval(Nm); + std::vector eval(JP.Nm); - std::vector src(Nu,FrbGrid); - for ( int i=0; i src(JP.Nu,FrbGrid); + for ( int i=0; i evec(Nm,FrbGrid); + std::vector evec(JP.Nm,FrbGrid); for(int i=0;i<1;++i){ - std::cout << GridLogMessage << i <<" / "<< Nm <<" grid pointer "<< evec[i]._grid << std::endl; + std::cout << GridLogMessage << i <<" / "<< JP.Nm <<" grid pointer "<< evec[i]._grid << std::endl; }; int Nconv; - IRBL.calc(eval,evec,src,Nconv); + IRBL.calc(eval,evec,src,Nconv,JP.Impl); Grid_finalize(); From 386b4fcb04a123c4d9b8637e71305af34816fdcd Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Mon, 5 Mar 2018 18:37:51 -0500 Subject: [PATCH 15/15] update convergence check --- .../ImplicitlyRestartedBlockLanczos.h | 66 +++++++++---------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index 2de35244..0ea4e09d 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -440,41 +440,41 @@ public: // Convergence test Glog <<" #Convergence test: "<