From dc6f637e707571ad589ec7004328af521cda638e Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Sat, 27 Jan 2018 08:21:27 -0500 Subject: [PATCH] change GparityDomainWallFermion to ZMobius and add command line options to read boundary phase and omega --- .../ImplicitlyRestartedBlockLanczos.h | 337 ++---- .../ImplicitlyRestartedBlockLanczos.h.bak2 | 1041 +++++++++++++++++ tests/lanczos/Test_dwf_block_lanczos.cc | 118 +- 3 files changed, 1197 insertions(+), 299 deletions(-) create mode 100644 lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h.bak2 diff --git a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index ff8050b7..933f2882 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -203,15 +203,6 @@ until convergence // additional (Nblock_m - Nblock_k) steps for(int b=Nblock_k; b=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;