/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid Source file: Hadrons/Modules/MDistil/LapEvec.hpp Copyright (C) 2019 Author: Felix Erben Author: Michael Marshall 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 Hadrons_MDistil_LapEvec_hpp_ #define Hadrons_MDistil_LapEvec_hpp_ #include #include #include #include // These are members of Distillation #include BEGIN_HADRONS_NAMESPACE BEGIN_MODULE_NAMESPACE(MDistil) /****************************************************************************** Laplacian eigenvectors - parameters ******************************************************************************/ struct StoutParameters: Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(StoutParameters, int, steps, double, rho) StoutParameters() = default; template StoutParameters(Reader& Reader){read(Reader,"StoutSmearing",*this);} }; struct ChebyshevParameters: Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(ChebyshevParameters, int, PolyOrder, double, alpha, double, beta) ChebyshevParameters() = default; template ChebyshevParameters(Reader& Reader){read(Reader,"Chebyshev",*this);} }; struct LanczosParameters: Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters, int, Nvec, int, Nk, int, Np, int, MaxIt, double, resid, int, IRLLog) LanczosParameters() = default; template LanczosParameters(Reader& Reader){read(Reader,"Lanczos",*this);} }; // These are the actual parameters passed to the module during construction struct LapEvecPar: Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(LapEvecPar ,std::string, gauge ,StoutParameters, Stout ,ChebyshevParameters, Cheby ,LanczosParameters, Lanczos) }; /****************************************************************************** Laplacian eigenvectors - Module (class) definition ******************************************************************************/ template class TLapEvec: public Module { public: GAUGE_TYPE_ALIASES(GImpl,); // constructor TLapEvec(const std::string name); // destructor virtual ~TLapEvec(void); // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); // setup virtual void setup(void); // execution virtual void execute(void); protected: // These variables are created in setup() and freed in Cleanup() GridCartesian * gridLD; // Owned by me, so I must delete it GridCartesian * gridHD; // Owned by environment (so I won't delete it) std::string sGaugeName; protected: virtual void Cleanup(void); }; MODULE_REGISTER_TMP(LapEvec, TLapEvec, MDistil); /****************************************************************************** TLapEvec implementation ******************************************************************************/ // constructor ///////////////////////////////////////////////////////////////// template TLapEvec::TLapEvec(const std::string name) : gridLD{nullptr}, Module(name) { } // destructor ///////////////////////////////////////////////////////////////// template TLapEvec::~TLapEvec() { Cleanup(); } // dependencies/products /////////////////////////////////////////////////////// template std::vector TLapEvec::getInput(void) { sGaugeName = par().gauge; if( sGaugeName.size() == 0 ) { sGaugeName = getName(); sGaugeName.append( "_gauge" ); } return std::vector{ sGaugeName }; } template std::vector TLapEvec::getOutput(void) { std::vector out = {getName()}; // This is the higher dimensional eigenpack return out; } // setup /////////////////////////////////////////////////////////////////////// template void TLapEvec::setup(void) { Cleanup(); Environment & e{env()}; gridHD = e.getGrid(); gridLD = MakeLowerDimGrid( gridHD ); const int Ntlocal{gridHD->LocalDimensions()[Tdir]}; // Temporaries envTmpLat(GaugeField, "Umu_stout"); envTmpLat(GaugeField, "Umu_smear"); envTmp(LatticeGaugeField, "UmuNoTime",1,LatticeGaugeField(gridLD)); envTmp(LatticeColourVector, "src",1,LatticeColourVector(gridLD)); envTmp(std::vector, "eig",1,std::vector(Ntlocal)); // Output objects envCreate(LapEvecs, getName(), 1, par().Lanczos.Nvec, gridHD ); } // clean up any temporaries created by setup (that aren't stored in the environment) template void TLapEvec::Cleanup(void) { if( gridLD != nullptr ) { delete gridLD; gridLD = nullptr; } gridHD = nullptr; } /****************************************************************************** Calculate low-mode eigenvalues of the Laplacian ******************************************************************************/ // execution /////////////////////////////////////////////////////////////////// template void TLapEvec::execute(void) { const ChebyshevParameters &ChebPar{par().Cheby}; const LanczosParameters &LPar{par().Lanczos}; // Disable IRL logging if requested LOG(Message) << "IRLLog=" << LPar.IRLLog << std::endl; const int PreviousIRLLogState{GridLogIRL.isActive()}; GridLogIRL.Active( LPar.IRLLog == 0 ? 0 : 1 ); // Stout smearing envGetTmp(GaugeField, Umu_smear); Umu_smear = envGet(GaugeField, sGaugeName); // The smeared field starts off as the Gauge field LOG(Message) << "Initial plaquette: " << WilsonLoops::avgPlaquette(Umu_smear) << std::endl; const StoutParameters &Stout{par().Stout}; if( Stout.steps ) { envGetTmp(GaugeField, Umu_stout); Smear_Stout LS(Stout.rho, Tdir); // spatial smearing only for (int i = 0; i < Stout.steps; i++) { LS.smear(Umu_stout, Umu_smear); Umu_smear = Umu_stout; } LOG(Message) << "Smeared plaquette: " << WilsonLoops::avgPlaquette(Umu_smear) << std::endl; } //////////////////////////////////////////////////////////////////////// // Invert Peardon Nabla operator separately on each time-slice //////////////////////////////////////////////////////////////////////// auto & eig4d = envGet(LapEvecs, getName() ); envGetTmp(std::vector, eig); // Eigenpack for each timeslice envGetTmp(LatticeGaugeField, UmuNoTime); // Gauge field without time dimension envGetTmp(LatticeColourVector, src); const int Ntlocal{gridHD->LocalDimensions()[Tdir]}; const int Ntfirst{gridHD->LocalStarts()[Tdir]}; uint32_t ConvergenceErrors{0}; for(int t = 0; t < Ntlocal; t++ ) { LOG(Message) << "------------------------------------------------------------" << std::endl; LOG(Message) << " Compute eigenpack, local timeslice = " << t << " / " << Ntlocal << std::endl; LOG(Message) << "------------------------------------------------------------" << std::endl; eig[t].resize(LPar.Nk+LPar.Np,gridLD); // Construct smearing operator ExtractSliceLocal(UmuNoTime,Umu_smear,0,t,Tdir); // switch to 3d/4d objects LinOpPeardonNabla PeardonNabla(UmuNoTime); LOG(Debug) << "Chebyshev preconditioning to order " << ChebPar.PolyOrder << " with parameters (alpha,beta) = (" << ChebPar.alpha << "," << ChebPar.beta << ")" << std::endl; Chebyshev Cheb(ChebPar.alpha,ChebPar.beta,ChebPar.PolyOrder); // Construct source vector according to Test_dwf_compressed_lanczos.cc src = 11.0; //TODO: Why hard-coded 11? RealD nn = norm2(src); nn = Grid::sqrt(nn); src = src * (1.0/nn); LinOpPeardonNablaHerm PeardonNablaCheby(Cheb,PeardonNabla); ImplicitlyRestartedLanczos IRL(PeardonNablaCheby,PeardonNabla,LPar.Nvec,LPar.Nk,LPar.Nk+LPar.Np,LPar.resid,LPar.MaxIt); int Nconv = 0; IRL.calc(eig[t].eval,eig[t].evec,src,Nconv); if( Nconv < LPar.Nvec ) { // NB: Can't assert here since we are processing local slices - i.e. not all nodes would assert ConvergenceErrors = 1; LOG(Error) << "MDistil::LapEvec : Not enough eigenvectors converged. If this occurs in practice, we should modify the eigensolver to iterate once more to ensure the second convergence test does not take us below the requested number of eigenvectors" << std::endl; } if( Nconv != LPar.Nvec ) eig[t].resize( LPar.Nvec, gridLD ); RotateEigen( eig[t].evec ); // Rotate the eigenvectors into our phase convention for (int i=0;iGlobalSum(ConvergenceErrors); assert(ConvergenceErrors==0 && "The eingensolver failed to find enough eigenvectors on at least one node"); #if DEBUG // Now write out the 4d eigenvectors eig4d.record.operatorXml = "Distillation"; eig4d.record.solverXml = "CG"; std::string sEigenPackName(getName()); sEigenPackName.append("."); sEigenPackName.append(std::to_string(vm().getTrajectory())); eig4d.write(sEigenPackName,false); #endif } END_MODULE_NAMESPACE END_HADRONS_NAMESPACE #endif // Hadrons_MDistil_LapEvec_hpp_