From 68e6a58f123a6da1f540b156f61e8d5b2f8c61d7 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 3 Apr 2018 17:42:21 +0100 Subject: [PATCH] Hadrons: several Lanczos fixes and improvements --- extras/Hadrons/EigenPack.hpp | 53 +++++++++++++------ .../Modules/MIO/LoadCoarseEigenPack.hpp | 11 +++- .../Modules/MSolver/LocalCoherenceLanczos.hpp | 15 +++--- .../iterative/LocalCoherenceLanczos.h | 8 +-- 4 files changed, 62 insertions(+), 25 deletions(-) diff --git a/extras/Hadrons/EigenPack.hpp b/extras/Hadrons/EigenPack.hpp index d27e35b9..08811107 100644 --- a/extras/Hadrons/EigenPack.hpp +++ b/extras/Hadrons/EigenPack.hpp @@ -152,55 +152,78 @@ public: evecCoarse.resize(sizeCoarse, gridCoarse); } - virtual void read(const std::string fileStem, const int traj = -1) + void readFine(const std::string fileStem, const int traj = -1) { std::string evecFineFilename, evalFineFilename; std::string evecCoarseFilename, evalCoarseFilename; this->makeFilenames(evecFineFilename, evalFineFilename, fileStem + "_fine", traj); - this->makeFilenames(evecCoarseFilename, evalCoarseFilename, - fileStem + "_coarse", traj); XmlReader xmlFineReader(evalFineFilename); - XmlReader xmlCoarseReader(evalCoarseFilename); LOG(Message) << "Reading " << this->evec.size() << " fine eigenvectors from '" << evecFineFilename << "'" << std::endl; this->basicRead(this->evec, evecFineFilename, this->evec.size()); - LOG(Message) << "Reading " << evecCoarse.size() << " coarse eigenvectors from '" - << evecCoarseFilename << "'" << std::endl; - this->basicRead(evecCoarse, evecCoarseFilename, evecCoarse.size()); LOG(Message) << "Reading " << this->eval.size() << " fine eigenvalues from '" << evalFineFilename << "'" << std::endl; Grid::read(xmlFineReader, "evals", this->eval); + } + + void readCoarse(const std::string fileStem, const int traj = -1) + { + std::string evecCoarseFilename, evalCoarseFilename; + + this->makeFilenames(evecCoarseFilename, evalCoarseFilename, + fileStem + "_coarse", traj); + XmlReader xmlCoarseReader(evalCoarseFilename); + LOG(Message) << "Reading " << evecCoarse.size() << " coarse eigenvectors from '" + << evecCoarseFilename << "'" << std::endl; + this->basicRead(evecCoarse, evecCoarseFilename, evecCoarse.size()); LOG(Message) << "Reading " << evalCoarse.size() << " coarse eigenvalues from '" << evalCoarseFilename << "'" << std::endl; Grid::read(xmlCoarseReader, "evals", evalCoarse); } - virtual void write(const std::string fileStem, const int traj = -1) + virtual void read(const std::string fileStem, const int traj = -1) + { + readFine(fileStem, traj); + readCoarse(fileStem, traj); + } + + void writeFine(const std::string fileStem, const int traj = -1) { std::string evecFineFilename, evalFineFilename; - std::string evecCoarseFilename, evalCoarseFilename; this->makeFilenames(evecFineFilename, evalFineFilename, fileStem + "_fine", traj); - this->makeFilenames(evecCoarseFilename, evalCoarseFilename, - fileStem + "_coarse", traj); XmlWriter xmlFineWriter(evalFineFilename); - XmlWriter xmlCoarseWriter(evalCoarseFilename); LOG(Message) << "Writing " << this->evec.size() << " fine eigenvectors to '" << evecFineFilename << "'" << std::endl; this->basicWrite(evecFineFilename, this->evec, this->evec.size()); - LOG(Message) << "Writing " << evecCoarse.size() << " coarse eigenvectors to '" - << evecCoarseFilename << "'" << std::endl; - this->basicWrite(evecCoarseFilename, evecCoarse, evecCoarse.size()); LOG(Message) << "Writing " << this->eval.size() << " fine eigenvalues to '" << evalFineFilename << "'" << std::endl; Grid::write(xmlFineWriter, "evals", this->eval); + } + + void writeCoarse(const std::string fileStem, const int traj = -1) + { + std::string evecCoarseFilename, evalCoarseFilename; + + this->makeFilenames(evecCoarseFilename, evalCoarseFilename, + fileStem + "_coarse", traj); + XmlWriter xmlCoarseWriter(evalCoarseFilename); + LOG(Message) << "Writing " << evecCoarse.size() << " coarse eigenvectors to '" + << evecCoarseFilename << "'" << std::endl; + this->basicWrite(evecCoarseFilename, evecCoarse, evecCoarse.size()); LOG(Message) << "Writing " << evalCoarse.size() << " coarse eigenvalues to '" << evalCoarseFilename << "'" << std::endl; Grid::write(xmlCoarseWriter, "evals", evalCoarse); } + + virtual void write(const std::string fileStem, const int traj = -1) + { + writeFine(fileStem, traj); + writeCoarse(fileStem, traj); + } }; template diff --git a/extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp b/extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp index 77c3a7ee..b88844f4 100644 --- a/extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp +++ b/extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp @@ -56,6 +56,9 @@ class TLoadCoarseEigenPack: public Module { public: typedef CoarseEigenPack BasePack; + template + using iImplScalar = iScalar>>; + typedef iImplScalar SiteComplex; public: // constructor TLoadCoarseEigenPack(const std::string name); @@ -114,9 +117,15 @@ void TLoadCoarseEigenPack::setup(void) template void TLoadCoarseEigenPack::execute(void) { - auto &epack = envGetDerived(BasePack, Pack, getName()); + auto cg = env().getCoarseGrid(par().blockSize, par().Ls); + auto &epack = envGetDerived(BasePack, Pack, getName()); + Lattice dummy(cg); epack.read(par().filestem, vm().getTrajectory()); + LOG(Message) << "Block Gramm-Schmidt pass 1"<< std::endl; + blockOrthogonalise(dummy, epack.evec); + LOG(Message) << "Block Gramm-Schmidt pass 2"<< std::endl; + blockOrthogonalise(dummy, epack.evec); } END_MODULE_NAMESPACE diff --git a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp index 80717f32..60b897c8 100644 --- a/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp +++ b/extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp @@ -125,19 +125,18 @@ void TLocalCoherenceLanczos::setup(void) env().createCoarseGrid(blockSize, Ls); - auto cg = env().getCoarseGrid(blockSize, Ls); - auto cgrb = env().getRbCoarseGrid(blockSize, Ls); - int cNm = (par().doCoarse) ? par().coarseParams.Nm : 0; + auto cg = env().getCoarseGrid(blockSize, Ls); + int cNm = (par().doCoarse) ? par().coarseParams.Nm : 0; LOG(Message) << "Coarse grid: " << cg->GlobalDimensions() << std::endl; envCreateDerived(BasePack, CoarsePack, getName(), Ls, - par().fineParams.Nm, cNm, env().getRbGrid(Ls), cgrb); + par().fineParams.Nm, cNm, env().getRbGrid(Ls), cg); auto &epack = envGetDerived(BasePack, CoarsePack, getName()); envTmp(SchurFMat, "mat", Ls, envGet(FMat, par().action)); envGetTmp(SchurFMat, mat); - envTmp(LCL, "solver", Ls, env().getRbGrid(Ls), cgrb, mat, + envTmp(LCL, "solver", Ls, env().getRbGrid(Ls), cg, mat, Odd, epack.evec, epack.evecCoarse, epack.eval, epack.evalCoarse); } @@ -157,6 +156,10 @@ void TLocalCoherenceLanczos::execute(void) finePar.resid,finePar.MaxIt, finePar.betastp, finePar.MinRes); solver.testFine(finePar.resid*100.0); + if (!par().output.empty()) + { + epack.writeFine(par().output, vm().getTrajectory()); + } if (par().doCoarse) { LOG(Message) << "Orthogonalising" << std::endl; @@ -173,7 +176,7 @@ void TLocalCoherenceLanczos::execute(void) } if (!par().output.empty()) { - epack.write(par().output, vm().getTrajectory()); + epack.writeCoarse(par().output, vm().getTrajectory()); } } diff --git a/lib/algorithms/iterative/LocalCoherenceLanczos.h b/lib/algorithms/iterative/LocalCoherenceLanczos.h index 0c7f4f65..66ded6eb 100644 --- a/lib/algorithms/iterative/LocalCoherenceLanczos.h +++ b/lib/algorithms/iterative/LocalCoherenceLanczos.h @@ -284,9 +284,11 @@ public: }; void Orthogonalise(void ) { - CoarseScalar InnerProd(_CoarseGrid); - blockOrthogonalise(InnerProd,subspace);std::cout << GridLogMessage <<" Gramm-Schmidt pass 1"< static RealD normalise(T& v)