Hadrons: several Lanczos fixes and improvements

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 <typename FImpl>

extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp

@@ -56,6 +56,9 @@ class TLoadCoarseEigenPack: public Module<LoadCoarseEigenPackPar>
 {
 public:
     typedef CoarseEigenPack<typename Pack::Field, typename Pack::CoarseField> BasePack;
+    template <typename vtype> 
+    using iImplScalar = iScalar<iScalar<iScalar<vtype>>>;
+    typedef iImplScalar<typename Pack::Field::vector_type> SiteComplex;
 public:
     // constructor
     TLoadCoarseEigenPack(const std::string name);
@@ -114,9 +117,15 @@ void TLoadCoarseEigenPack<Pack>::setup(void)
 template <typename Pack>
 void TLoadCoarseEigenPack<Pack>::execute(void)
 {
-    auto &epack = envGetDerived(BasePack, Pack, getName());
+    auto                 cg     = env().getCoarseGrid(par().blockSize, par().Ls);
+    auto                 &epack = envGetDerived(BasePack, Pack, getName());
+    Lattice<SiteComplex> 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

extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp

@@ -125,19 +125,18 @@ void TLocalCoherenceLanczos<FImpl, nBasis>::setup(void)
 
     env().createCoarseGrid(blockSize, Ls);
 
-    auto cg   = env().getCoarseGrid(blockSize, Ls);
+    auto cg  = env().getCoarseGrid(blockSize, Ls);
-    auto cgrb = env().getRbCoarseGrid(blockSize, Ls);
+    int  cNm = (par().doCoarse) ? par().coarseParams.Nm : 0;
-    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<FImpl, nBasis>::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<FImpl, nBasis>::execute(void)
     }
     if (!par().output.empty())
     {
-        epack.write(par().output, vm().getTrajectory());
+        epack.writeCoarse(par().output, vm().getTrajectory());
     }
 }
 

lib/algorithms/iterative/LocalCoherenceLanczos.h

@@ -284,9 +284,11 @@ public:
   };
 
   void Orthogonalise(void ) {
-    CoarseScalar InnerProd(_CoarseGrid); 
+    CoarseScalar InnerProd(_CoarseGrid);
-    blockOrthogonalise(InnerProd,subspace);std::cout << GridLogMessage <<" Gramm-Schmidt pass 1"<<std::endl;
+    std::cout << GridLogMessage <<" Gramm-Schmidt pass 1"<<std::endl;
-    blockOrthogonalise(InnerProd,subspace);std::cout << GridLogMessage <<" Gramm-Schmidt pass 2"<<std::endl;
+    blockOrthogonalise(InnerProd,subspace);
+    std::cout << GridLogMessage <<" Gramm-Schmidt pass 2"<<std::endl;
+    blockOrthogonalise(InnerProd,subspace);
   };
 
   template<typename T>  static RealD normalise(T& v)