Hadrons: Lanczos code improvement

extras/Hadrons/Global.cc

@@ -37,6 +37,7 @@ HadronsLogger Hadrons::HadronsLogWarning(1,"Warning");
 HadronsLogger Hadrons::HadronsLogMessage(1,"Message");
 HadronsLogger Hadrons::HadronsLogIterative(1,"Iterative");
 HadronsLogger Hadrons::HadronsLogDebug(1,"Debug");
+HadronsLogger Hadrons::HadronsLogIRL(1,"IRL");
 
 void Hadrons::initLogger(void)
 {
@@ -46,11 +47,13 @@ void Hadrons::initLogger(void)
     GridLogMessage.setTopWidth(w);
     GridLogIterative.setTopWidth(w);
     GridLogDebug.setTopWidth(w);
+    GridLogIRL.setTopWidth(w);
     HadronsLogError.Active(GridLogError.isActive());
     HadronsLogWarning.Active(GridLogWarning.isActive());
     HadronsLogMessage.Active(GridLogMessage.isActive());
     HadronsLogIterative.Active(GridLogIterative.isActive());
     HadronsLogDebug.Active(GridLogDebug.isActive());
+    HadronsLogIRL.Active(GridLogIRL.isActive());
 }
 
 // type utilities //////////////////////////////////////////////////////////////

extras/Hadrons/Global.hpp

@@ -111,6 +111,7 @@ extern HadronsLogger HadronsLogWarning;
 extern HadronsLogger HadronsLogMessage;
 extern HadronsLogger HadronsLogIterative;
 extern HadronsLogger HadronsLogDebug;
+extern HadronsLogger HadronsLogIRL;
 
 void initLogger(void);
 

extras/Hadrons/LanczosUtils.hpp

@@ -0,0 +1,88 @@
+#ifndef Hadrons_LanczosUtils_hpp_
+#define Hadrons_LanczosUtils_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/algorithms/iterative/LocalCoherenceLanczos.h>
+
+BEGIN_HADRONS_NAMESPACE
+
+// Lanczos type
+#ifndef HADRONS_DEFAULT_LANCZOS_NBASIS
+#define HADRONS_DEFAULT_LANCZOS_NBASIS 60
+#endif
+
+template <typename FImpl, int nBasis>
+using LCL = LocalCoherenceLanczos<typename FImpl::SiteSpinor, 
+                                  typename FImpl::SiteComplex, nBasis>;
+
+template <typename T>
+struct EigenPack
+{
+    std::vector<RealD> eval;
+    std::vector<T>     evec;
+
+    EigenPack(void) = default;
+
+    EigenPack(const size_t size, GridBase *grid)
+    {
+        resize(size, grid);
+    }
+
+    void resize(const size_t size, GridBase *grid)
+    {
+        eval.resize(size);
+        evec.resize(size, grid);
+    }
+
+    void read(const std::string fileStem)
+    {
+        std::string     evecFilename = fileStem + "_evec.bin";
+        std::string     evalFilename = fileStem + "_eval.xml";
+        emptyUserRecord record;
+        ScidacReader    binReader;
+        XmlReader       xmlReader(evalFilename);
+
+        LOG(Message) << "Reading " << evec.size() << " eigenvectors from '" 
+                     << evecFilename << "'" << std::endl;
+        binReader.open(evecFilename);
+        for(int k = 0; k < evec.size(); ++k) 
+        {
+            binReader.readScidacFieldRecord(evec[k], record);
+        }
+        binReader.close();
+        LOG(Message) << "Reading " << eval.size() << " eigenvalues from '" 
+                     << evalFilename << "'" << std::endl;
+        Grid::read(xmlReader, "evals", eval);
+    }
+
+    void write(const std::string fileStem)
+    {
+        std::string     evecFilename = fileStem + "_evec.bin";
+        std::string     evalFilename = fileStem + "_eval.xml";
+        emptyUserRecord record;
+        ScidacWriter    binWriter;
+        XmlWriter       xmlWriter(evalFilename);
+
+        LOG(Message) << "Writing " << evec.size() << " eigenvectors to '" 
+                     << evecFilename << "'" << std::endl;
+        binWriter.open(fileStem + "_evec.bin");
+        for(int k = 0; k < evec.size(); ++k) 
+        {
+            binWriter.writeScidacFieldRecord(evec[k], record);
+        }
+        binWriter.close();
+        LOG(Message) << "Writing " << eval.size() << " eigenvalues to '" 
+                     << evalFilename << "'" << std::endl;
+        Grid::write(xmlWriter, "evals", eval);
+    }
+};
+
+template <typename FImpl>
+using FineEigenPack = EigenPack<typename FImpl::FermionField>;
+
+template <typename FImpl, int nBasis>
+using CoarseEigenPack = EigenPack<typename LCL<FImpl, nBasis>::CoarseField>;
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_LanczosUtils_hpp_

extras/Hadrons/Makefile.am

@@ -21,6 +21,7 @@ nobase_libHadrons_a_HEADERS = \
 	GeneticScheduler.hpp      \
 	Global.hpp                \
 	Graph.hpp                 \
+	LanczosUtils.hpp          \
 	Module.hpp                \
 	Modules.hpp               \
 	ModuleFactory.hpp         \

extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp

@@ -31,11 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
 #include <Grid/Hadrons/ModuleFactory.hpp>
-#include <Grid/algorithms/iterative/LocalCoherenceLanczos.h>
+#include <Grid/Hadrons/LanczosUtils.hpp>
-
-#ifndef DEFAULT_LANCZOS_NBASIS
-#define DEFAULT_LANCZOS_NBASIS 60
-#endif
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -55,7 +51,8 @@ public:
                                     LanczosParams,    coarseParams,
                                     ChebyParams,      smoother,
                                     RealD,            coarseRelaxTol,
-                                    std::string,      blockSize);
+                                    std::string,      blockSize,
+                                    std::string,      output);
 };
 
 template <typename FImpl, int nBasis>
@@ -63,9 +60,9 @@ class TLocalCoherenceLanczos: public Module<LocalCoherenceLanczosPar>
 {
 public:
     FERM_TYPE_ALIASES(FImpl,);
-    typedef LocalCoherenceLanczos<vSpinColourVector, vTComplex, nBasis> LCL;
+    typedef LCL<FImpl, nBasis>                         LCL;
-    typedef typename LCL::FineField                             FineField;
+    typedef FineEigenPack<FImpl>                       FineEigenPack;
-    typedef typename LCL::CoarseField                           CoarseField;
+    typedef CoarseEigenPack<FImpl, nBasis>             CoarseEigenPack; 
     typedef SchurDiagMooeeOperator<FMat, FermionField> SchurFMat;
 public:
     // constructor
@@ -88,12 +85,12 @@ private:
     std::unique_ptr<GridCartesian>         coarseGrid_{nullptr};
     std::unique_ptr<GridRedBlackCartesian> coarseGrid4Rb_{nullptr};
     std::unique_ptr<GridRedBlackCartesian> coarseGridRb_{nullptr};
-    std::string fevecName_, cevecName_, fevalName_, cevalName_;
+    std::string                            fineName_, coarseName_;
 };
 
 MODULE_REGISTER_NS(LocalCoherenceLanczos, 
-                   ARG(TLocalCoherenceLanczos<FIMPL, DEFAULT_LANCZOS_NBASIS>), 
+    ARG(TLocalCoherenceLanczos<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), 
-                   MSolver);
+    MSolver);
 
 /******************************************************************************
  *                 TLocalCoherenceLanczos implementation                             *
@@ -103,10 +100,8 @@ template <typename FImpl, int nBasis>
 TLocalCoherenceLanczos<FImpl, nBasis>::TLocalCoherenceLanczos(const std::string name)
 : Module<LocalCoherenceLanczosPar>(name)
 {
-    fevecName_ = getName() + "_fineEvec";
+    fineName_   = getName() + "_fine";
-    cevecName_ = getName() + "_coarseEvec";
+    coarseName_ = getName() + "_coarse";
-    fevalName_ = getName() + "_fineEval";
-    cevalName_ = getName() + "_coarseEval";
 }
 
 // dependencies/products ///////////////////////////////////////////////////////
@@ -121,7 +116,7 @@ std::vector<std::string> TLocalCoherenceLanczos<FImpl, nBasis>::getInput(void)
 template <typename FImpl, int nBasis>
 std::vector<std::string> TLocalCoherenceLanczos<FImpl, nBasis>::getOutput(void)
 {
-    std::vector<std::string> out = {fevecName_, cevecName_, fevalName_, cevalName_};
+    std::vector<std::string> out = {fineName_, coarseName_};
     
     return out;
 }
@@ -187,47 +182,56 @@ void TLocalCoherenceLanczos<FImpl, nBasis>::setup(void)
     {
         makeCoarseGrid();
     }
-    envCreate(std::vector<FineField>, fevecName_, Ls_, par().fineParams.Nm,
+    envCreate(FineEigenPack, fineName_, Ls_, par().fineParams.Nm,
               env().getRbGrid(Ls_));
-    envCreate(std::vector<CoarseField>, cevecName_, Ls_, par().coarseParams.Nm, 
+    envCreate(CoarseEigenPack, coarseName_, Ls_, par().coarseParams.Nm, 
               coarseGridRb_.get());
-    envCreate(std::vector<RealD>, fevalName_, Ls_, par().fineParams.Nm);
+    auto &fine   = envGet(FineEigenPack, fineName_);
-    envCreate(std::vector<RealD>, cevalName_, Ls_, par().coarseParams.Nm);
+    auto &coarse = envGet(CoarseEigenPack, coarseName_);
     envTmp(SchurFMat, "mat", Ls_, envGet(FMat, par().action));
     envGetTmp(SchurFMat, mat);
-    envTmp(LCL, "solver", Ls_, env().getRbGrid(Ls_), coarseGridRb_.get(), mat, Odd, 
+    envTmp(LCL, "solver", Ls_, env().getRbGrid(Ls_), coarseGridRb_.get(), mat, 
-           envGet(std::vector<FineField>, fevecName_),
+           Odd, fine.evec, coarse.evec, fine.eval, coarse.eval);
-           envGet(std::vector<CoarseField>, cevecName_),
-           envGet(std::vector<RealD>, fevalName_),
-           envGet(std::vector<RealD>, cevalName_));
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl, int nBasis>
 void TLocalCoherenceLanczos<FImpl, nBasis>::execute(void)
 {
-    auto &fine   = par().fineParams;
+    auto &finePar   = par().fineParams;
-    auto &coarse = par().coarseParams;
+    auto &coarsePar = par().coarseParams;
-    
+    auto &fine      = envGet(FineEigenPack, fineName_);
+    auto &coarse    = envGet(CoarseEigenPack, coarseName_);
+
     envGetTmp(LCL, solver);
     if (par().doFine)
     {
         LOG(Message) << "Performing fine grid IRL -- Nstop= " 
-                     << fine.Nstop << ", Nk= " << fine.Nk << ", Nm= " 
+                     << finePar.Nstop << ", Nk= " << finePar.Nk << ", Nm= " 
-                     << fine.Nm << std::endl;
+                     << finePar.Nm << std::endl;
-        solver.calcFine(fine.Cheby, fine.Nstop, fine.Nk, fine.Nm, fine.resid,
+        solver.calcFine(finePar.Cheby, finePar.Nstop, finePar.Nk, finePar.Nm,
-                        fine.MaxIt, fine.betastp, fine.MinRes);
+                        finePar.resid,finePar.MaxIt, finePar.betastp, 
+                        finePar.MinRes);
         LOG(Message) << "Orthogonalising" << std::endl;
         solver.Orthogonalise();
+        if (!par().output.empty())
+        {
+            fine.write(par().output + "_fine");
+        }
     }
     if (par().doCoarse)
     {
         LOG(Message) << "Performing coarse grid IRL -- Nstop= " 
-                     << fine.Nstop << ", Nk= " << fine.Nk << ", Nm= " 
+                     << coarsePar.Nstop << ", Nk= " << coarsePar.Nk << ", Nm= " 
-                     << fine.Nm << std::endl;
+                     << coarsePar.Nm << std::endl;
-        solver.calcCoarse(coarse.Cheby, par().smoother, par().coarseRelaxTol,
+        solver.calcCoarse(coarsePar.Cheby, par().smoother, par().coarseRelaxTol,
-			              coarse.Nstop, coarse.Nk, coarse.Nm, coarse.resid, 
+			              coarsePar.Nstop, coarsePar.Nk, coarsePar.Nm, 
-                          coarse.MaxIt, coarse.betastp, coarse.MinRes);
+                          coarsePar.resid, coarsePar.MaxIt, coarsePar.betastp, 
+                          coarsePar.MinRes);
+        if (!par().output.empty())
+        {
+            coarse.write(par().output + "_coarse");
+        }
     }
 }