Merge branch 'feature/hadrons' into feature/qed-fvol

.gitignore

@@ -123,6 +123,7 @@ make-bin-BUCK.sh
 #####################
 lib/qcd/spin/gamma-gen/*.h
 lib/qcd/spin/gamma-gen/*.cc
+lib/version.h
 
 # vs code editor files #
 ########################

Makefile.am

@@ -5,6 +5,10 @@ include $(top_srcdir)/doxygen.inc
 
 bin_SCRIPTS=grid-config
 
+BUILT_SOURCES = version.h
+
+version.h:
+	echo "`git log -n 1 --format=format:"#define GITHASH \\"%H:%d\\"%n" HEAD`" > $(srcdir)/lib/version.h
 
 .PHONY: bench check tests doxygen-run doxygen-doc $(DX_PS_GOAL) $(DX_PDF_GOAL)
 

configure.ac

@@ -340,7 +340,7 @@ case ${ac_PRECISION} in
 esac
 
 ######################  Shared memory allocation technique under MPI3
-AC_ARG_ENABLE([shm],[AC_HELP_STRING([--enable-shm=shmopen|hugetlbfs],
+AC_ARG_ENABLE([shm],[AC_HELP_STRING([--enable-shm=shmopen|hugetlbfs|shmnone],
               [Select SHM allocation technique])],[ac_SHM=${enable_shm}],[ac_SHM=shmopen])
 
 case ${ac_SHM} in
@@ -349,6 +349,10 @@ case ${ac_SHM} in
      AC_DEFINE([GRID_MPI3_SHMOPEN],[1],[GRID_MPI3_SHMOPEN] )
      ;;
 
+     shmnone)
+     AC_DEFINE([GRID_MPI3_SHM_NONE],[1],[GRID_MPI3_SHM_NONE] )
+     ;;
+
      hugetlbfs)
      AC_DEFINE([GRID_MPI3_SHMMMAP],[1],[GRID_MPI3_SHMMMAP] )
      ;;

extras/Hadrons/EigenPack.hpp

@@ -0,0 +1,218 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/EigenPack.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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_EigenPack_hpp_
+#define Hadrons_EigenPack_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/algorithms/iterative/Deflation.h>
+#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 F>
+class EigenPack
+{
+public:
+    typedef F Field;
+public:
+    std::vector<RealD> eval;
+    std::vector<F>     evec;
+public:
+    EigenPack(void)          = default;
+    virtual ~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);
+    }
+
+    virtual void read(const std::string fileStem, const int traj = -1)
+    {
+        std::string evecFilename, evalFilename;
+
+        makeFilenames(evecFilename, evalFilename, fileStem, traj);
+        XmlReader xmlReader(evalFilename);
+        basicRead(evec, evecFilename, evec.size());
+        LOG(Message) << "Reading " << eval.size() << " eigenvalues from '" 
+                     << evalFilename << "'" << std::endl;
+        Grid::read(xmlReader, "evals", eval);
+    }
+
+    virtual void write(const std::string fileStem, const int traj = -1)
+    {
+        std::string evecFilename, evalFilename;
+
+        makeFilenames(evecFilename, evalFilename, fileStem, traj);
+        XmlWriter xmlWriter(evalFilename);
+        basicWrite(evecFilename, evec, evec.size());
+        LOG(Message) << "Writing " << eval.size() << " eigenvalues to '" 
+                     << evalFilename << "'" << std::endl;
+        Grid::write(xmlWriter, "evals", eval);
+    }
+protected:
+    void makeFilenames(std::string &evecFilename, std::string &evalFilename,
+                       const std::string stem, const int traj = -1)
+    {
+        std::string t = (traj < 0) ? "" : ("." + std::to_string(traj));
+
+        evecFilename = stem + "_evec" + t + ".bin";
+        evalFilename = stem + "_eval" + t + ".xml";
+    }
+
+    template <typename T>
+    static void basicRead(std::vector<T> &evec, const std::string filename, 
+                          const unsigned int size)
+    {
+        emptyUserRecord record;
+        ScidacReader    binReader;
+
+        binReader.open(filename);
+        for(int k = 0; k < size; ++k) 
+        {
+            binReader.readScidacFieldRecord(evec[k], record);
+        }
+        binReader.close();
+    }
+
+    template <typename T>
+    static void basicWrite(const std::string filename, std::vector<T> &evec, 
+                           const unsigned int size)
+    {
+        emptyUserRecord record;
+        ScidacWriter    binWriter;
+
+        binWriter.open(filename);
+        for(int k = 0; k < size; ++k) 
+        {
+            binWriter.writeScidacFieldRecord(evec[k], record);
+        }
+        binWriter.close();
+    }
+};
+
+template <typename FineF, typename CoarseF>
+class CoarseEigenPack: public EigenPack<FineF>
+{
+public:
+    typedef CoarseF CoarseField;
+public:
+    std::vector<RealD>   evalCoarse;
+    std::vector<CoarseF> evecCoarse;
+public:
+    CoarseEigenPack(void)          = default;
+    virtual ~CoarseEigenPack(void) = default;
+
+    CoarseEigenPack(const size_t sizeFine, const size_t sizeCoarse, 
+                    GridBase *gridFine, GridBase *gridCoarse)
+    {
+        resize(sizeFine, sizeCoarse, gridFine, gridCoarse);
+    }
+
+    void resize(const size_t sizeFine, const size_t sizeCoarse, 
+                GridBase *gridFine, GridBase *gridCoarse)
+    {
+        EigenPack<FineF>::resize(sizeFine, gridFine);
+        evalCoarse.resize(sizeCoarse);
+        evecCoarse.resize(sizeCoarse, gridCoarse);
+    }
+
+    virtual void read(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);
+        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)
+    {
+        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);
+        LOG(Message) << "Writing " << evalCoarse.size() << " coarse eigenvalues to '" 
+                     << evalCoarseFilename << "'" << std::endl;
+        Grid::write(xmlCoarseWriter, "evals", evalCoarse);
+    }
+};
+
+template <typename FImpl>
+using FermionEigenPack = EigenPack<typename FImpl::FermionField>;
+
+template <typename FImpl, int nBasis>
+using CoarseFermionEigenPack = CoarseEigenPack<
+    typename FImpl::FermionField,
+    typename LocalCoherenceLanczos<typename FImpl::SiteSpinor, 
+                                   typename FImpl::SiteComplex, 
+                                   nBasis>::CoarseField>;
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_EigenPack_hpp_

extras/Hadrons/Environment.cc

@@ -61,7 +61,7 @@ Environment::Environment(void)
 // grids ///////////////////////////////////////////////////////////////////////
 void Environment::createGrid(const unsigned int Ls)
 {
-    if (grid5d_.find(Ls) == grid5d_.end())
+    if ((Ls > 1) and (grid5d_.find(Ls) == grid5d_.end()))
     {
         auto g = getGrid();
         
@@ -70,6 +70,53 @@ void Environment::createGrid(const unsigned int Ls)
     }
 }
 
+void Environment::createCoarseGrid(const std::vector<int> &blockSize, 
+                                   const unsigned int Ls)
+{
+    int              nd      = getNd();
+    std::vector<int> fineDim = getDim(), coarseDim;
+    unsigned int     cLs;
+    auto             key4d = blockSize, key5d = blockSize;
+
+    createGrid(Ls);
+    coarseDim.resize(nd);
+    for (int d = 0; d < coarseDim.size(); d++)
+    {
+        coarseDim[d] = fineDim[d]/blockSize[d];
+        if (coarseDim[d]*blockSize[d] != fineDim[d])
+        {
+            HADRON_ERROR(Size, "Fine dimension " + std::to_string(d) 
+                         + " (" + std::to_string(fineDim[d]) 
+                         + ") not divisible by coarse dimension ("
+                         + std::to_string(coarseDim[d]) + ")"); 
+        }
+    }
+    if (blockSize.size() > nd)
+    {
+        cLs = Ls/blockSize[nd];
+        if (cLs*blockSize[nd] != Ls)
+        {
+            HADRON_ERROR(Size, "Fine Ls (" + std::to_string(Ls) 
+                         + ") not divisible by coarse Ls ("
+                         + std::to_string(cLs) + ")");
+        }
+        key4d.resize(nd);
+        key5d.push_back(Ls);
+    }
+    gridCoarse4d_[key4d].reset(
+        SpaceTimeGrid::makeFourDimGrid(coarseDim, 
+            GridDefaultSimd(nd, vComplex::Nsimd()), GridDefaultMpi()));
+    gridCoarseRb4d_[key4d].reset(
+        SpaceTimeGrid::makeFourDimRedBlackGrid(gridCoarse4d_[key4d].get()));
+    if (Ls > 1)
+    {
+        gridCoarse5d_[key5d].reset(
+            SpaceTimeGrid::makeFiveDimGrid(cLs, gridCoarse4d_[key4d].get()));
+        gridCoarseRb5d_[key5d].reset(
+            SpaceTimeGrid::makeFiveDimRedBlackGrid(cLs, gridCoarse4d_[key4d].get()));
+    }
+}
+
 GridCartesian * Environment::getGrid(const unsigned int Ls) const
 {
     try
@@ -104,7 +151,55 @@ GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
     }
     catch(std::out_of_range &)
     {
-        HADRON_ERROR(Definition, "no red-black 5D grid with Ls= " + std::to_string(Ls));
+        HADRON_ERROR(Definition, "no red-black grid with Ls= " + std::to_string(Ls));
+    }
+}
+
+GridCartesian * Environment::getCoarseGrid(
+    const std::vector<int> &blockSize, const unsigned int Ls) const
+{
+    auto key = blockSize;
+
+    try
+    {
+        if (Ls == 1)
+        {
+            key.resize(getNd());
+            return gridCoarse4d_.at(key).get();
+        }
+        else
+        {
+            key.push_back(Ls);
+            return gridCoarse5d_.at(key).get();
+        }
+    }
+    catch(std::out_of_range &)
+    {
+        HADRON_ERROR(Definition, "no coarse grid with Ls= " + std::to_string(Ls));
+    }
+}
+
+GridRedBlackCartesian * Environment::getRbCoarseGrid(
+    const std::vector<int> &blockSize, const unsigned int Ls) const
+{
+    auto key = blockSize;
+
+    try
+    {
+        if (Ls == 1)
+        {
+            key.resize(getNd());
+            return gridCoarseRb4d_.at(key).get();
+        }
+        else
+        {
+            key.push_back(Ls);
+            return gridCoarseRb5d_.at(key).get();
+        }
+    }
+    catch(std::out_of_range &)
+    {
+        HADRON_ERROR(Definition, "no coarse red-black grid with Ls= " + std::to_string(Ls));
     }
 }
 

extras/Hadrons/Environment.hpp

@@ -78,7 +78,7 @@ private:
         Size                    size{0};
         Storage                 storage{Storage::object};
         unsigned int            Ls{0};
-        const std::type_info    *type{nullptr};
+        const std::type_info    *type{nullptr}, *derivedType{nullptr};
         std::string             name;
         int                     module{-1};
         std::unique_ptr<Object> data{nullptr};
@@ -86,8 +86,14 @@ private:
 public:
     // grids
     void                    createGrid(const unsigned int Ls);
+    void                    createCoarseGrid(const std::vector<int> &blockSize,
+                                             const unsigned int Ls = 1);
     GridCartesian *         getGrid(const unsigned int Ls = 1) const;
     GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const;
+    GridCartesian *         getCoarseGrid(const std::vector<int> &blockSize,
+                                          const unsigned int Ls = 1) const;
+    GridRedBlackCartesian * getRbCoarseGrid(const std::vector<int> &blockSize,
+                                            const unsigned int Ls = 1) const;
     std::vector<int>        getDim(void) const;
     int                     getDim(const unsigned int mu) const;
     unsigned long int       getLocalVolume(void) const;
@@ -110,6 +116,10 @@ public:
                                          Ts && ... args);
     void                    setObjectModule(const unsigned int objAddress,
                                             const int modAddress);
+    template <typename B, typename T>
+    T *                     getDerivedObject(const unsigned int address) const;
+    template <typename B, typename T>
+    T *                     getDerivedObject(const std::string name) const;
     template <typename T>
     T *                     getObject(const unsigned int address) const;
     template <typename T>
@@ -155,6 +165,10 @@ private:
     std::map<unsigned int, GridPt>         grid5d_;
     GridRbPt                               gridRb4d_;
     std::map<unsigned int, GridRbPt>       gridRb5d_;
+    std::map<std::vector<int>, GridPt>     gridCoarse4d_;
+    std::map<std::vector<int>, GridRbPt>   gridCoarseRb4d_;
+    std::map<std::vector<int>, GridPt>     gridCoarse5d_;
+    std::map<std::vector<int>, GridRbPt>   gridCoarseRb5d_;
     unsigned int                           nd_;
     // random number generator
     RngPt                                  rng4d_;
@@ -176,7 +190,7 @@ Holder<T>::Holder(T *pt)
 template <typename T>
 T & Holder<T>::get(void) const
 {
-    return &objPt_.get();
+    return *objPt_.get();
 }
 
 template <typename T>
@@ -216,22 +230,24 @@ void Environment::createDerivedObject(const std::string name,
         {
             MemoryProfiler::stats = &memStats;
         }
-        size_t initMem           = MemoryProfiler::stats->currentlyAllocated;
+        size_t initMem               = MemoryProfiler::stats->currentlyAllocated;
-        object_[address].storage = storage;
+        object_[address].storage     = storage;
-        object_[address].Ls      = Ls;
+        object_[address].Ls          = Ls;
         object_[address].data.reset(new Holder<B>(new T(std::forward<Ts>(args)...)));
-        object_[address].size    = MemoryProfiler::stats->maxAllocated - initMem;
+        object_[address].size        = MemoryProfiler::stats->maxAllocated - initMem;
-        object_[address].type    = &typeid(T);
+        object_[address].type        = &typeid(B);
+        object_[address].derivedType = &typeid(T);
         if (MemoryProfiler::stats == &memStats)
         {
             MemoryProfiler::stats = nullptr;
         }
     }
     // object already exists, no error if it is a cache, error otherwise
-    else if ((object_[address].storage != Storage::cache) or 
+    else if ((object_[address].storage     != Storage::cache) or 
-             (object_[address].storage != storage)        or
+             (object_[address].storage     != storage)        or
-             (object_[address].name    != name)           or
+             (object_[address].name        != name)           or
-             (object_[address].type    != &typeid(T)))
+             (object_[address].type        != &typeid(B))     or
+             (object_[address].derivedType != &typeid(T)))
     {
         HADRON_ERROR(Definition, "object '" + name + "' already allocated");
     }
@@ -246,21 +262,37 @@ void Environment::createObject(const std::string name,
     createDerivedObject<T, T>(name, storage, Ls, std::forward<Ts>(args)...);
 }
 
-template <typename T>
+template <typename B, typename T>
-T * Environment::getObject(const unsigned int address) const
+T * Environment::getDerivedObject(const unsigned int address) const
 {
     if (hasObject(address))
     {
         if (hasCreatedObject(address))
         {
-            if (auto h = dynamic_cast<Holder<T> *>(object_[address].data.get()))
+            if (auto h = dynamic_cast<Holder<B> *>(object_[address].data.get()))
             {
-                return h->getPt();
+                if (&typeid(T) == &typeid(B))
+                {
+                    return dynamic_cast<T *>(h->getPt());
+                }
+                else
+                {
+                    if (auto hder = dynamic_cast<T *>(h->getPt()))
+                    {
+                        return hder;
+                    }
+                    else
+                    {
+                        HADRON_ERROR(Definition, "object with address " + std::to_string(address) +
+                            " cannot be casted to '" + typeName(&typeid(T)) +
+                            "' (has type '" + typeName(&typeid(h->get())) + "')");
+                    }
+                }
             }
             else
             {
                 HADRON_ERROR(Definition, "object with address " + std::to_string(address) +
-                            " does not have type '" + typeName(&typeid(T)) +
+                            " does not have type '" + typeName(&typeid(B)) +
                             "' (has type '" + getObjectType(address) + "')");
             }
         }
@@ -276,6 +308,18 @@ T * Environment::getObject(const unsigned int address) const
     }
 }
 
+template <typename B, typename T>
+T * Environment::getDerivedObject(const std::string name) const
+{
+    return getDerivedObject<B, T>(getObjectAddress(name));
+}
+
+template <typename T>
+T * Environment::getObject(const unsigned int address) const
+{
+    return getDerivedObject<T, T>(address);
+}
+
 template <typename T>
 T * Environment::getObject(const std::string name) const
 {

extras/Hadrons/Global.hpp

@@ -43,12 +43,15 @@ See the full license in the file "LICENSE" in the top level distribution directo
 namespace Grid {\
 using namespace QCD;\
 namespace Hadrons {\
-using Grid::operator<<;
+using Grid::operator<<;\
+using Grid::operator>>;
 #define END_HADRONS_NAMESPACE }}
 
 #define BEGIN_MODULE_NAMESPACE(name)\
 namespace name {\
-using Grid::operator<<;
+using Grid::operator<<;\
+using Grid::operator>>;
+
 #define END_MODULE_NAMESPACE }
 
 /* the 'using Grid::operator<<;' statement prevents a very nasty compilation
@@ -187,7 +190,7 @@ name + "." + std::to_string(vm().getTrajectory()) + "." + resultFileExt
 // default Schur convention
 
 #ifndef HADRONS_DEFAULT_SCHUR 
-#define HADRONS_DEFAULT_SCHUR DiagMooee
+#define HADRONS_DEFAULT_SCHUR DiagTwo
 #endif
 #define _HADRONS_SCHUR_OP_(conv) Schur##conv##Operator
 #define HADRONS_SCHUR_OP(conv) _HADRONS_SCHUR_OP_(conv)

extras/Hadrons/LanczosUtils.hpp

@@ -1,115 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: extras/Hadrons/LanczosUtils.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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_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 T>
-struct EigenPack
-{
-    typedef T VectorType;
-    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 LocalCoherenceLanczos<typename FImpl::SiteSpinor, 
-                                   typename FImpl::SiteComplex, 
-                                   nBasis>::CoarseField>;
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_LanczosUtils_hpp_

extras/Hadrons/Makefile.am

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

extras/Hadrons/Module.hpp

@@ -91,6 +91,9 @@ static ns##mod##ModuleRegistrar ns##mod##ModuleRegistrarInstance;
 #define envGet(type, name)\
 *env().template getObject<type>(name)
 
+#define envGetDerived(base, type, name)\
+*env().template getDerivedObject<base, type>(name)
+
 #define envGetTmp(type, var)\
 type &var = *env().template getObject<type>(getName() + "_tmp_" + #var)
 

extras/Hadrons/Modules.hpp

@@ -64,6 +64,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
 #include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
 #include <Grid/Hadrons/Modules/MAction/ZMobiusDWF.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/EMT.hpp>
@@ -72,5 +73,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TransProj.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
+#include <Grid/Hadrons/Modules/MIO/LoadEigenPack.hpp>
 #include <Grid/Hadrons/Modules/MIO/LoadNersc.hpp>
+#include <Grid/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
 #include <Grid/Hadrons/Modules/MIO/LoadBinary.hpp>

extras/Hadrons/Modules/MGauge/FundtoHirep.cc

@@ -42,7 +42,8 @@ TFundtoHirep<Rep>::TFundtoHirep(const std::string name)
 template <class Rep>
 std::vector<std::string> TFundtoHirep<Rep>::getInput(void)
 {
-    std::vector<std::string> in;
+    std::vector<std::string> in = {par().gaugeconf};
+
     return in;
 }
 
@@ -50,6 +51,7 @@ template <class Rep>
 std::vector<std::string> TFundtoHirep<Rep>::getOutput(void)
 {
     std::vector<std::string> out = {getName()};
+
     return out;
 }
 
@@ -57,19 +59,19 @@ std::vector<std::string> TFundtoHirep<Rep>::getOutput(void)
 template <typename Rep>
 void TFundtoHirep<Rep>::setup(void)
 {
-    envCreateLat(typename Rep::LatticeField, getName());
+    envCreateLat(Rep::LatticeField, getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <class Rep>
 void TFundtoHirep<Rep>::execute(void)
 {
-    auto &U      = *env().template getObject<LatticeGaugeField>(par().gaugeconf);
     LOG(Message) << "Transforming Representation" << std::endl;
 
+    auto &U    = envGet(LatticeGaugeField, par().gaugeconf);
+    auto &URep = envGet(Rep::LatticeField, getName());
+
     Rep TargetRepresentation(U._grid);
     TargetRepresentation.update_representation(U);
-
-    auto &URep = envGet(typename Rep::LatticeField, getName());
     URep = TargetRepresentation.U;
 }

extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp

@@ -0,0 +1,126 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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_MIO_LoadCoarseEigenPack_hpp_
+#define Hadrons_MIO_LoadCoarseEigenPack_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/EigenPack.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *              Load local coherence eigen vectors/values package             *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MIO)
+
+class LoadCoarseEigenPackPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(LoadCoarseEigenPackPar,
+                                    std::string, filestem,
+                                    unsigned int, sizeFine,
+                                    unsigned int, sizeCoarse,
+                                    unsigned int, Ls,
+                                    std::vector<int>, blockSize);
+};
+
+template <typename Pack>
+class TLoadCoarseEigenPack: public Module<LoadCoarseEigenPackPar>
+{
+public:
+    typedef CoarseEigenPack<typename Pack::Field, typename Pack::CoarseField> BasePack;
+public:
+    // constructor
+    TLoadCoarseEigenPack(const std::string name);
+    // destructor
+    virtual ~TLoadCoarseEigenPack(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(LoadCoarseFermionEigenPack, 
+    ARG(TLoadCoarseEigenPack<CoarseFermionEigenPack<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>>), MIO);
+
+/******************************************************************************
+ *                 TLoadCoarseEigenPack implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename Pack>
+TLoadCoarseEigenPack<Pack>::TLoadCoarseEigenPack(const std::string name)
+: Module<LoadCoarseEigenPackPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename Pack>
+std::vector<std::string> TLoadCoarseEigenPack<Pack>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename Pack>
+std::vector<std::string> TLoadCoarseEigenPack<Pack>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename Pack>
+void TLoadCoarseEigenPack<Pack>::setup(void)
+{
+    env().createGrid(par().Ls);
+    env().createCoarseGrid(par().blockSize, par().Ls);
+    envCreateDerived(BasePack, Pack, getName(), par().Ls, par().sizeFine,
+                     par().sizeCoarse, env().getRbGrid(par().Ls), 
+                     env().getCoarseGrid(par().blockSize, par().Ls));
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename Pack>
+void TLoadCoarseEigenPack<Pack>::execute(void)
+{
+    auto &epack = envGetDerived(BasePack, Pack, getName());
+
+    epack.read(par().filestem, vm().getTrajectory());
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MIO_LoadCoarseEigenPack_hpp_

extras/Hadrons/Modules/MIO/LoadEigenPack.hpp

@@ -0,0 +1,121 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MIO/LoadEigenPack.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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_MIO_LoadEigenPack_hpp_
+#define Hadrons_MIO_LoadEigenPack_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/EigenPack.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                   Load eigen vectors/values package                        *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MIO)
+
+class LoadEigenPackPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(LoadEigenPackPar,
+                                    std::string, filestem,
+                                    unsigned int, size,
+                                    unsigned int, Ls);
+};
+
+template <typename Pack>
+class TLoadEigenPack: public Module<LoadEigenPackPar>
+{
+public:
+    typedef EigenPack<typename Pack::Field> BasePack;
+public:
+    // constructor
+    TLoadEigenPack(const std::string name);
+    // destructor
+    virtual ~TLoadEigenPack(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(LoadFermionEigenPack, TLoadEigenPack<FermionEigenPack<FIMPL>>, MIO);
+
+/******************************************************************************
+ *                    TLoadEigenPack implementation                           *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename Pack>
+TLoadEigenPack<Pack>::TLoadEigenPack(const std::string name)
+: Module<LoadEigenPackPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename Pack>
+std::vector<std::string> TLoadEigenPack<Pack>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename Pack>
+std::vector<std::string> TLoadEigenPack<Pack>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename Pack>
+void TLoadEigenPack<Pack>::setup(void)
+{
+    env().createGrid(par().Ls);
+    envCreateDerived(BasePack, Pack, getName(), par().Ls, par().size, 
+                     env().getRbGrid(par().Ls));
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename Pack>
+void TLoadEigenPack<Pack>::execute(void)
+{
+    auto &epack = envGetDerived(BasePack, Pack, getName());
+
+    epack.read(par().filestem, vm().getTrajectory());
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MIO_LoadEigenPack_hpp_

extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp

@@ -0,0 +1,169 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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_MScalarSUN_ShiftProbe_hpp_
+#define Hadrons_MScalarSUN_ShiftProbe_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *         Ward identity phi^n probe with fields at different positions       *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MScalarSUN)
+
+typedef std::pair<unsigned int, unsigned int> ShiftPair;
+
+class ShiftProbePar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ShiftProbePar,
+                                    std::string, field,
+                                    std::string, shifts,
+                                    std::string, output);
+};
+
+template <typename SImpl>
+class TShiftProbe: public Module<ShiftProbePar>
+{
+public:
+    
+    typedef typename SImpl::Field                          Field;
+    typedef typename SImpl::ComplexField                   ComplexField;
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        std::string, op,
+                                        Complex    , value);
+    };
+public:
+    // constructor
+    TShiftProbe(const std::string name);
+    // destructor
+    virtual ~TShiftProbe(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(ShiftProbeSU2, TShiftProbe<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_NS(ShiftProbeSU3, TShiftProbe<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_NS(ShiftProbeSU4, TShiftProbe<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_NS(ShiftProbeSU5, TShiftProbe<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_NS(ShiftProbeSU6, TShiftProbe<ScalarNxNAdjImplR<6>>, MScalarSUN);
+
+/******************************************************************************
+ *                        TShiftProbe implementation                          *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename SImpl>
+TShiftProbe<SImpl>::TShiftProbe(const std::string name)
+: Module<ShiftProbePar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename SImpl>
+std::vector<std::string> TShiftProbe<SImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().field};
+    
+    return in;
+}
+
+template <typename SImpl>
+std::vector<std::string> TShiftProbe<SImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TShiftProbe<SImpl>::setup(void)
+{
+    envTmpLat(Field, "acc");
+    envCreateLat(ComplexField, getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TShiftProbe<SImpl>::execute(void)
+{
+    LOG(Message) << "Creating shift probe for shifts " << par().shifts
+                 << std::endl;
+
+    std::vector<ShiftPair> shift;
+    unsigned int           sign;
+    auto                   &phi   = envGet(Field, par().field);
+    auto                   &probe = envGet(ComplexField, getName());
+
+    shift = strToVec<ShiftPair>(par().shifts);
+    if (shift.size() % 2 != 0)
+    {
+        HADRON_ERROR(Size, "the number of shifts is odd");
+    }
+    sign = (shift.size() % 4 == 0) ? 1 : -1;
+    for (auto &s: shift)
+    {
+        if (s.first >= env().getNd())
+        {
+            HADRON_ERROR(Size, "dimension to large for shift <" 
+                               + std::to_string(s.first) + " " 
+                               + std::to_string(s.second) + ">" );
+        }
+    }
+    envGetTmp(Field, acc);
+    acc   = 1.;
+    for (unsigned int i = 0; i < shift.size(); ++i)
+    {
+        if (shift[i].second == 0)
+        {
+            acc *= phi;
+        }
+        else
+        {
+            acc *= Cshift(phi, shift[i].first, shift[i].second);
+        }
+    }
+    probe = sign*trace(acc);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MScalarSUN_ShiftProbe_hpp_

extras/Hadrons/Modules/MSink/Point.hpp

@@ -128,7 +128,7 @@ void TPoint<FImpl>::execute(void)
         envGetTmp(LatticeComplex, coor);
         p  = strToVec<Real>(par().mom);
         ph = zero;
-        for(unsigned int mu = 0; mu < env().getNd(); mu++)
+        for(unsigned int mu = 0; mu < p.size(); mu++)
         {
             LatticeCoordinate(coor, mu);
             ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor;

extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp

@@ -31,7 +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/Hadrons/LanczosUtils.hpp>
+#include <Grid/Hadrons/EigenPack.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -45,8 +45,7 @@ class LocalCoherenceLanczosPar: Serializable
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(LocalCoherenceLanczosPar,
                                     std::string,   action,
-                                    int,           doFine,
+                                    bool,          doCoarse,
-                                    int,           doCoarse,
                                     LanczosParams, fineParams,
                                     LanczosParams, coarseParams,
                                     ChebyParams,   smoother,
@@ -63,8 +62,8 @@ public:
     typedef LocalCoherenceLanczos<typename FImpl::SiteSpinor, 
                                   typename FImpl::SiteComplex, 
                                   nBasis>                LCL;
-    typedef FineEigenPack<FImpl>                         FinePack;
+    typedef FermionEigenPack<FImpl>                      BasePack;
-    typedef CoarseEigenPack<FImpl, nBasis>               CoarsePack; 
+    typedef CoarseFermionEigenPack<FImpl, nBasis>        CoarsePack;
     typedef HADRONS_DEFAULT_SCHUR_OP<FMat, FermionField> SchurFMat;
 public:
     // constructor
@@ -78,16 +77,6 @@ public:
     virtual void setup(void);
     // execution
     virtual void execute(void);
-private:
-    void makeCoarseGrid(void);
-private:
-    std::vector<int>                       coarseDim_;
-    int                                    Ls_, cLs_{1};
-    std::unique_ptr<GridCartesian>         coarseGrid4_{nullptr};
-    std::unique_ptr<GridCartesian>         coarseGrid_{nullptr};
-    std::unique_ptr<GridRedBlackCartesian> coarseGrid4Rb_{nullptr};
-    std::unique_ptr<GridRedBlackCartesian> coarseGridRb_{nullptr};
-    std::string                            fineName_, coarseName_;
 };
 
 MODULE_REGISTER_NS(LocalCoherenceLanczos, 
@@ -104,10 +93,7 @@ MODULE_REGISTER_NS(ZLocalCoherenceLanczos,
 template <typename FImpl, int nBasis>
 TLocalCoherenceLanczos<FImpl, nBasis>::TLocalCoherenceLanczos(const std::string name)
 : Module<LocalCoherenceLanczosPar>(name)
-{
+{}
-    fineName_   = getName() + "_fine";
-    coarseName_ = getName() + "_coarse";
-}
 
 // dependencies/products ///////////////////////////////////////////////////////
 template <typename FImpl, int nBasis>
@@ -121,61 +107,12 @@ 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 = {fineName_, coarseName_};
+    std::vector<std::string> out = {getName()};
     
     return out;
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
-template <typename FImpl, int nBasis>
-void TLocalCoherenceLanczos<FImpl, nBasis>::makeCoarseGrid(void)
-{
-    int              nd        = env().getNd();
-    std::vector<int> blockSize = strToVec<int>(par().blockSize);
-    auto             fineDim   = env().getDim();
-
-    Ls_ = env().getObjectLs(par().action);
-    env().createGrid(Ls_);
-    coarseDim_.resize(nd);
-    for (int d = 0; d < coarseDim_.size(); d++)
-    {
-        coarseDim_[d] = fineDim[d]/blockSize[d];
-        if (coarseDim_[d]*blockSize[d] != fineDim[d])
-        {
-            HADRON_ERROR(Size, "Fine dimension " + std::to_string(d) 
-                         + " (" + std::to_string(fineDim[d]) 
-                         + ") not divisible by coarse dimension ("
-                         + std::to_string(coarseDim_[d]) + ")"); 
-        }
-    }
-    if (blockSize.size() > nd)
-    {
-        cLs_ = Ls_/blockSize[nd];
-        if (cLs_*blockSize[nd] != Ls_)
-        {
-            HADRON_ERROR(Size, "Fine Ls (" + std::to_string(Ls_) 
-                         + ") not divisible by coarse Ls ("
-                         + std::to_string(cLs_) + ")");
-        }
-    }
-    if (Ls_ > 1)
-    {
-        coarseGrid4_.reset(SpaceTimeGrid::makeFourDimGrid(
-            coarseDim_, GridDefaultSimd(nd, vComplex::Nsimd()),
-            GridDefaultMpi()));
-        coarseGrid4Rb_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(coarseGrid4_.get()));
-        coarseGrid_.reset(SpaceTimeGrid::makeFiveDimGrid(cLs_, coarseGrid4_.get()));
-        coarseGridRb_.reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(cLs_, coarseGrid4_.get()));
-    }
-    else
-    {
-        coarseGrid_.reset(SpaceTimeGrid::makeFourDimGrid(
-            coarseDim_, GridDefaultSimd(nd, vComplex::Nsimd()),
-            GridDefaultMpi()));
-        coarseGridRb_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(coarseGrid_.get()));
-    }
-}
-
 template <typename FImpl, int nBasis>
 void TLocalCoherenceLanczos<FImpl, nBasis>::setup(void)
 {
@@ -183,19 +120,25 @@ void TLocalCoherenceLanczos<FImpl, nBasis>::setup(void)
                  << " action '" << par().action << "' (" << nBasis
                  << " eigenvectors)..." << std::endl;
     
-    if (!coarseGrid_)
+    unsigned int Ls        = env().getObjectLs(par().action);
-    {
+    auto         blockSize = strToVec<int>(par().blockSize);
-        makeCoarseGrid();
+
-    }
+    env().createCoarseGrid(blockSize, Ls);
-    LOG(Message) << "Coarse grid: " << coarseGrid_->GlobalDimensions() << std::endl;
+
-    envCreate(FinePack, fineName_, Ls_, par().fineParams.Nm, env().getRbGrid(Ls_));
+    auto cg   = env().getCoarseGrid(blockSize, Ls);
-    envCreate(CoarsePack, coarseName_, Ls_, par().coarseParams.Nm, coarseGridRb_.get());
+    auto cgrb = env().getRbCoarseGrid(blockSize, Ls);
-    auto &fine   = envGet(FinePack, fineName_);
+    int  cNm  = (par().doCoarse) ? par().coarseParams.Nm : 0;
-    auto &coarse = envGet(CoarsePack, coarseName_);
+
-    envTmp(SchurFMat, "mat", Ls_, envGet(FMat, par().action));
+    LOG(Message) << "Coarse grid: " << cg->GlobalDimensions() << std::endl;
+    envCreateDerived(BasePack, CoarsePack, getName(), Ls,
+                     par().fineParams.Nm, cNm, env().getRbGrid(Ls), cgrb);
+
+    auto &epack = envGetDerived(BasePack, CoarsePack, getName());
+
+    envTmp(SchurFMat, "mat", Ls, envGet(FMat, par().action));
     envGetTmp(SchurFMat, mat);
-    envTmp(LCL, "solver", Ls_, env().getRbGrid(Ls_), coarseGridRb_.get(), mat, 
+    envTmp(LCL, "solver", Ls, env().getRbGrid(Ls), cgrb, mat, 
-           Odd, fine.evec, coarse.evec, fine.eval, coarse.eval);
+           Odd, epack.evec, epack.evecCoarse, epack.eval, epack.evalCoarse);
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -204,41 +147,33 @@ void TLocalCoherenceLanczos<FImpl, nBasis>::execute(void)
 {
     auto &finePar   = par().fineParams;
     auto &coarsePar = par().coarseParams;
-    auto &fine      = envGet(FinePack, fineName_);
+    auto &epack     = envGetDerived(BasePack, CoarsePack, getName());
-    auto &coarse    = envGet(CoarsePack, coarseName_);
 
     envGetTmp(LCL, solver);
-    if (par().doFine)
+    LOG(Message) << "Performing fine grid IRL -- Nstop= " 
-    {
+                 << finePar.Nstop << ", Nk= " << finePar.Nk << ", Nm= " 
-        LOG(Message) << "Performing fine grid IRL -- Nstop= " 
+                 << finePar.Nm << std::endl;
-                     << finePar.Nstop << ", Nk= " << finePar.Nk << ", Nm= " 
+    solver.calcFine(finePar.Cheby, finePar.Nstop, finePar.Nk, finePar.Nm,
-                     << finePar.Nm << std::endl;
+                    finePar.resid,finePar.MaxIt, finePar.betastp, 
-        solver.calcFine(finePar.Cheby, finePar.Nstop, finePar.Nk, finePar.Nm,
+                    finePar.MinRes);
-                        finePar.resid,finePar.MaxIt, finePar.betastp, 
+    solver.testFine(finePar.resid*100.0);
-                        finePar.MinRes);
-        solver.testFine(finePar.resid*100.0);
-        LOG(Message) << "Orthogonalising" << std::endl;
-        solver.Orthogonalise();
-        if (!par().output.empty())
-        {
-            fine.write(par().output + "_fine");
-        }
-    }
     if (par().doCoarse)
     {
+        LOG(Message) << "Orthogonalising" << std::endl;
+        solver.Orthogonalise();
         LOG(Message) << "Performing coarse grid IRL -- Nstop= " 
-                     << coarsePar.Nstop << ", Nk= " << coarsePar.Nk << ", Nm= " 
+                    << coarsePar.Nstop << ", Nk= " << coarsePar.Nk << ", Nm= " 
-                     << coarsePar.Nm << std::endl;
+                    << coarsePar.Nm << std::endl;
         solver.calcCoarse(coarsePar.Cheby, par().smoother, par().coarseRelaxTol,
-			              coarsePar.Nstop, coarsePar.Nk, coarsePar.Nm, 
+                          coarsePar.Nstop, coarsePar.Nk, coarsePar.Nm, 
                           coarsePar.resid, coarsePar.MaxIt, coarsePar.betastp, 
                           coarsePar.MinRes);
         solver.testCoarse(coarsePar.resid*100.0, par().smoother, 
-                          par().coarseRelaxTol);
+                        par().coarseRelaxTol);
-        if (!par().output.empty())
+    }
-        {
+    if (!par().output.empty())
-            coarse.write(par().output + "_coarse");
+    {
-        }
+        epack.write(par().output, vm().getTrajectory());
     }
 }
 

extras/Hadrons/Modules/MSolver/RBPrecCG.hpp

@@ -32,6 +32,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/Hadrons/EigenPack.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -44,16 +45,24 @@ class RBPrecCGPar: Serializable
 {
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(RBPrecCGPar ,
-                                    std::string    , action,
+                                    std::string , action,
-                                    unsigned int   , maxIteration,
+                                    unsigned int, maxIteration,
-                                    double         , residual);
+                                    double      , residual,
+                                    std::string , eigenPack);
 };
 
-template <typename FImpl>
+template <typename FImpl, int nBasis>
 class TRBPrecCG: public Module<RBPrecCGPar>
 {
 public:
     FGS_TYPE_ALIASES(FImpl,);
+    typedef FermionEigenPack<FImpl>                       EPack;
+    typedef CoarseFermionEigenPack<FImpl, nBasis>         CoarseEPack;
+    typedef std::shared_ptr<Guesser<FermionField>>        GuesserPt;
+    typedef DeflatedGuesser<typename FImpl::FermionField> FineGuesser;
+    typedef LocalCoherenceDeflatedGuesser<
+        typename FImpl::FermionField,
+        typename CoarseEPack::CoarseField> CoarseGuesser;
 public:
     // constructor
     TRBPrecCG(const std::string name);
@@ -70,37 +79,39 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(RBPrecCG,  TRBPrecCG<FIMPL>, MSolver);
+MODULE_REGISTER_NS(RBPrecCG,  
-MODULE_REGISTER_NS(ZRBPrecCG, TRBPrecCG<ZFIMPL>, MSolver);
+    ARG(TRBPrecCG<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
+MODULE_REGISTER_NS(ZRBPrecCG, 
+    ARG(TRBPrecCG<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
 
 /******************************************************************************
  *                      TRBPrecCG template implementation                     *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-template <typename FImpl>
+template <typename FImpl, int nBasis>
-TRBPrecCG<FImpl>::TRBPrecCG(const std::string name)
+TRBPrecCG<FImpl, nBasis>::TRBPrecCG(const std::string name)
 : Module(name)
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
-template <typename FImpl>
+template <typename FImpl, int nBasis>
-std::vector<std::string> TRBPrecCG<FImpl>::getInput(void)
+std::vector<std::string> TRBPrecCG<FImpl, nBasis>::getInput(void)
 {
     std::vector<std::string> in = {};
     
     return in;
 }
 
-template <typename FImpl>
+template <typename FImpl, int nBasis>
-std::vector<std::string> TRBPrecCG<FImpl>::getReference(void)
+std::vector<std::string> TRBPrecCG<FImpl, nBasis>::getReference(void)
 {
     std::vector<std::string> ref = {par().action};
     
     return ref;
 }
 
-template <typename FImpl>
+template <typename FImpl, int nBasis>
-std::vector<std::string> TRBPrecCG<FImpl>::getOutput(void)
+std::vector<std::string> TRBPrecCG<FImpl, nBasis>::getOutput(void)
 {
     std::vector<std::string> out = {getName()};
     
@@ -108,30 +119,60 @@ std::vector<std::string> TRBPrecCG<FImpl>::getOutput(void)
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
-template <typename FImpl>
+template <typename FImpl, int nBasis>
-void TRBPrecCG<FImpl>::setup(void)
+void TRBPrecCG<FImpl, nBasis>::setup(void)
 {
+    if (par().maxIteration == 0)
+    {
+        HADRON_ERROR(Argument, "zero maximum iteration");
+    }
+
     LOG(Message) << "setting up Schur red-black preconditioned CG for"
                  << " action '" << par().action << "' with residual "
-                 << par().residual << std::endl;
+                 << par().residual << ", maximum iteration " 
+                 << par().maxIteration << std::endl;
 
-    auto Ls     = env().getObjectLs(par().action);
+    auto        Ls          = env().getObjectLs(par().action);
-    auto &mat   = envGet(FMat, par().action);
+    auto        &mat        = envGet(FMat, par().action);
-    auto solver = [&mat, this](FermionField &sol, const FermionField &source)
+    std::string guesserName = getName() + "_guesser";
+    GuesserPt   guesser{nullptr};
+
+    if (par().eigenPack.empty())
+    {
+        guesser.reset(new ZeroGuesser<FermionField>());
+    }
+    else
+    {
+        try
+        {
+            auto &epack = envGetDerived(EPack, CoarseEPack, par().eigenPack);
+            
+            guesser.reset(new CoarseGuesser(epack.evec, epack.evecCoarse,
+                                            epack.evalCoarse));
+        }
+        catch (Exceptions::Definition &e)
+        {
+            auto &epack = envGet(EPack, par().eigenPack);
+
+            guesser.reset(new FineGuesser(epack.evec, epack.eval));
+        }
+    }
+    auto solver = [&mat, guesser, this](FermionField &sol, 
+                                        const FermionField &source)
     {
         ConjugateGradient<FermionField>           cg(par().residual, 
                                                      par().maxIteration);
         HADRONS_DEFAULT_SCHUR_SOLVE<FermionField> schurSolver(cg);
         
-        schurSolver(mat, source, sol);
+        schurSolver(mat, source, sol, *guesser);
     };
     envCreate(SolverFn, getName(), Ls, solver);
 }
 
 
 // execution ///////////////////////////////////////////////////////////////////
-template <typename FImpl>
+template <typename FImpl, int nBasis>
-void TRBPrecCG<FImpl>::execute(void)
+void TRBPrecCG<FImpl, nBasis>::execute(void)
 {}
 
 END_MODULE_NAMESPACE

extras/Hadrons/VirtualMachine.cc

@@ -381,7 +381,6 @@ void VirtualMachine::makeMemoryProfile(void)
     env().protectObjects(false);
     GridLogMessage.Active(false);
     HadronsLogMessage.Active(false);
-    HadronsLogError.Active(false);
     for (auto it = program.rbegin(); it != program.rend(); ++it) 
     {
         auto a = *it;
@@ -397,7 +396,6 @@ void VirtualMachine::makeMemoryProfile(void)
     env().protectObjects(protect);
     GridLogMessage.Active(gmsg);
     HadronsLogMessage.Active(hmsg);
-    HadronsLogError.Active(err);
     LOG(Debug) << "Memory profile:" << std::endl;
     LOG(Debug) << "----------------" << std::endl;
     for (unsigned int a = 0; a < profile_.module.size(); ++a)
@@ -632,6 +630,17 @@ void VirtualMachine::executeProgram(const Program &p) const
     // build garbage collection schedule
     LOG(Debug) << "Building garbage collection schedule..." << std::endl;
     freeProg = makeGarbageSchedule(p);
+    for (unsigned int i = 0; i < freeProg.size(); ++i)
+    {
+        std::string msg = "";
+
+        for (auto &a: freeProg[i])
+        {
+            msg += env().getObjectName(a) + " ";
+        }
+        msg += "]";
+        LOG(Debug) << std::setw(4) << i + 1 << ": [" << msg << std::endl;
+    }
 
     // program execution
     LOG(Debug) << "Executing program..." << std::endl;

extras/Hadrons/modules.inc

@@ -50,6 +50,7 @@ modules_hpp =\
   Modules/MAction/Wilson.hpp \
   Modules/MAction/WilsonClover.hpp \
   Modules/MAction/ZMobiusDWF.hpp \
+  Modules/MScalarSUN/ShiftProbe.hpp \
   Modules/MScalarSUN/Div.hpp \
   Modules/MScalarSUN/TrMag.hpp \
   Modules/MScalarSUN/EMT.hpp \
@@ -58,6 +59,8 @@ modules_hpp =\
   Modules/MScalarSUN/Utils.hpp \
   Modules/MScalarSUN/TransProj.hpp \
   Modules/MScalarSUN/TrKinetic.hpp \
+  Modules/MIO/LoadEigenPack.hpp \
   Modules/MIO/LoadNersc.hpp \
+  Modules/MIO/LoadCoarseEigenPack.hpp \
   Modules/MIO/LoadBinary.hpp
 

lib/algorithms/iterative/Deflation.h

@@ -30,22 +30,31 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 namespace Grid { 
 
-struct ZeroGuesser {
+template<class Field>
+class Guesser {
 public:
-  template<class Field> 
+  Guesser(void) = default;
-  void operator()(const Field &src,Field &guess) { guess = Zero(); };
+  virtual ~Guesser(void) = default;
+  virtual void operator()(const Field &src, Field &guess) = 0;
 };
-struct SourceGuesser {
+
+template<class Field>
+class ZeroGuesser: public Guesser<Field> {
 public:
-  template<class Field> 
+  virtual void operator()(const Field &src, Field &guess) { guess = zero; };
-  void operator()(const Field &src,Field &guess) { guess = src; };
+};
+
+template<class Field>
+class SourceGuesser: public Guesser<Field> {
+public:
+  virtual void operator()(const Field &src, Field &guess) { guess = src; };
 };
 
 ////////////////////////////////
 // Fine grid deflation
 ////////////////////////////////
 template<class Field>
-struct DeflatedGuesser {
+class DeflatedGuesser: public Guesser<Field> {
 private:
   const std::vector<Field> &evec;
   const std::vector<RealD> &eval;
@@ -54,7 +63,7 @@ public:
 
   DeflatedGuesser(const std::vector<Field> & _evec,const std::vector<RealD> & _eval) : evec(_evec), eval(_eval) {};
 
-  void operator()(const Field &src,Field &guess) { 
+  virtual void operator()(const Field &src,Field &guess) { 
     guess = zero;
     assert(evec.size()==eval.size());
     auto N = evec.size();
@@ -66,7 +75,7 @@ public:
 };
 
 template<class FineField, class CoarseField>
-class LocalCoherenceDeflatedGuesser {
+class LocalCoherenceDeflatedGuesser: public Guesser<FineField> {
 private:
   const std::vector<FineField>   &subspace;
   const std::vector<CoarseField> &evec_coarse;

lib/algorithms/iterative/SchurRedBlack.h

@@ -108,7 +108,7 @@ namespace Grid {
 
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
-      ZeroGuesser guess;
+      ZeroGuesser<Field> guess;
       (*this)(_Matrix,in,out,guess);
     }
     template<class Matrix, class Guesser>
@@ -195,7 +195,7 @@ namespace Grid {
   };
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
-      ZeroGuesser guess;
+      ZeroGuesser<Field> guess;
       (*this)(_Matrix,in,out,guess);
     }
     template<class Matrix, class Guesser>
@@ -280,7 +280,7 @@ namespace Grid {
 
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
-      ZeroGuesser guess;
+      ZeroGuesser<Field> guess;
       (*this)(_Matrix,in,out,guess);
     }
     template<class Matrix,class Guesser>
@@ -365,7 +365,7 @@ namespace Grid {
 
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
-      ZeroGuesser guess;
+      ZeroGuesser<Field> guess;
       (*this)(_Matrix,in,out,guess);
     }
     template<class Matrix, class Guesser>

lib/communicator/SharedMemoryMPI.cc

@@ -226,6 +226,48 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 };
 #endif // MMAP
 
+#ifdef GRID_MPI3_SHM_NONE
+void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
+{
+  std::cout << "SharedMemoryAllocate "<< bytes<< " MMAP anonymous implementation "<<std::endl;
+  assert(_ShmSetup==1);
+  assert(_ShmAlloc==0);
+  //////////////////////////////////////////////////////////////////////////////////////////////////////////
+  // allocate the shared windows for our group
+  //////////////////////////////////////////////////////////////////////////////////////////////////////////
+  MPI_Barrier(WorldShmComm);
+  WorldShmCommBufs.resize(WorldShmSize);
+  
+  ////////////////////////////////////////////////////////////////////////////////////////////
+  // Hugetlbf and others map filesystems as mappable huge pages
+  ////////////////////////////////////////////////////////////////////////////////////////////
+  char shm_name [NAME_MAX];
+  assert(WorldShmSize == 1);
+  for(int r=0;r<WorldShmSize;r++){
+    
+    int fd=-1;
+    int mmap_flag = MAP_SHARED |MAP_ANONYMOUS ;
+#ifdef MAP_POPULATE    
+    mmap_flag|=MAP_POPULATE;
+#endif
+#ifdef MAP_HUGETLB
+    if ( flags ) mmap_flag |= MAP_HUGETLB;
+#endif
+    void *ptr = (void *) mmap(NULL, bytes, PROT_READ | PROT_WRITE, mmap_flag,fd, 0); 
+    if ( ptr == (void *)MAP_FAILED ) {    
+      printf("mmap %s failed\n",shm_name);
+      perror("failed mmap");      assert(0);    
+    }
+    assert(((uint64_t)ptr&0x3F)==0);
+    close(fd);
+    WorldShmCommBufs[r] =ptr;
+    std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< bytes<< "bytes)"<<std::endl;
+  }
+  _ShmAlloc=1;
+  _ShmAllocBytes  = bytes;
+};
+#endif // MMAP
+
 #ifdef GRID_MPI3_SHMOPEN
 ////////////////////////////////////////////////////////////////////////////////////////////
 // POSIX SHMOPEN ; as far as I know Linux does not allow EXPLICIT HugePages with this case
@@ -246,7 +288,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 	
       size_t size = bytes;
       
-      sprintf(shm_name,"/Grid_mpi3_shm_%d_%d",WorldNode,r);
+      sprintf(shm_name,"/myGrid_mpi3_shm_%d_%d",WorldNode,r);
       
       shm_unlink(shm_name);
       int fd=shm_open(shm_name,O_RDWR|O_CREAT,0666);
@@ -401,7 +443,10 @@ void *SharedMemory::ShmBufferTranslate(int rank,void * local_p)
 }
 SharedMemory::~SharedMemory()
 {
-  MPI_Comm_free(&ShmComm);
+  int MPI_is_finalised;  MPI_Finalized(&MPI_is_finalised);
+  if ( !MPI_is_finalised ) { 
+    MPI_Comm_free(&ShmComm);
+  }
 };
 
 }

lib/lattice/Lattice_comparison_utils.h

@@ -198,7 +198,7 @@ namespace Grid {
       typedef typename vsimd::scalar_type scalar;\
       return Comparison(functor<scalar,scalar>(),lhs,rhs);\
     }\
-  template<class vsimd>\
+  template<class vsimd,IfSimd<vsimd> = 0>\
     inline vInteger operator op(const iScalar<vsimd> &lhs,const iScalar<vsimd> &rhs)\
     {									\
       return lhs._internal op rhs._internal;				\

lib/parallelIO/BinaryIO.h

@@ -91,7 +91,7 @@ class BinaryIO {
     typedef typename vobj::scalar_object sobj;
 
     GridBase *grid = lat._grid;
-    int lsites = grid->lSites();
+    uint64_t lsites = grid->lSites();
 
     std::vector<sobj> scalardata(lsites); 
     unvectorizeToLexOrdArray(scalardata,lat);    
@@ -160,7 +160,9 @@ class BinaryIO {
 
 	/* 
 	 * Scidac csum  is rather more heavyweight
+	 * FIXME -- 128^3 x 256 x 16 will overflow.
 	 */
+	
 	int global_site;
 
 	Lexicographic::CoorFromIndex(coor,local_site,local_vol);
@@ -261,7 +263,7 @@ class BinaryIO {
 			      GridBase *grid,
 			      std::vector<fobj> &iodata,
 			      std::string file,
-			      Integer offset,
+			      uint64_t offset,
 			      const std::string &format, int control,
 			      uint32_t &nersc_csum,
 			      uint32_t &scidac_csuma,
@@ -523,7 +525,7 @@ class BinaryIO {
   static inline void readLatticeObject(Lattice<vobj> &Umu,
 				       std::string file,
 				       munger munge,
-				       Integer offset,
+				       uint64_t offset,
 				       const std::string &format,
 				       uint32_t &nersc_csum,
 				       uint32_t &scidac_csuma,
@@ -533,7 +535,7 @@ class BinaryIO {
     typedef typename vobj::Realified::scalar_type word;    word w=0;
 
     GridBase *grid = Umu._grid;
-    int lsites = grid->lSites();
+    uint64_t lsites = grid->lSites();
 
     std::vector<sobj> scalardata(lsites); 
     std::vector<fobj>     iodata(lsites); // Munge, checksum, byte order in here
@@ -544,7 +546,7 @@ class BinaryIO {
     GridStopWatch timer; 
     timer.Start();
 
-    parallel_for(int x=0;x<lsites;x++) munge(iodata[x], scalardata[x]);
+    parallel_for(uint64_t x=0;x<lsites;x++) munge(iodata[x], scalardata[x]);
 
     vectorizeFromLexOrdArray(scalardata,Umu);    
     grid->Barrier();
@@ -560,7 +562,7 @@ class BinaryIO {
     static inline void writeLatticeObject(Lattice<vobj> &Umu,
 					  std::string file,
 					  munger munge,
-					  Integer offset,
+					  uint64_t offset,
 					  const std::string &format,
 					  uint32_t &nersc_csum,
 					  uint32_t &scidac_csuma,
@@ -569,7 +571,7 @@ class BinaryIO {
     typedef typename vobj::scalar_object sobj;
     typedef typename vobj::Realified::scalar_type word;    word w=0;
     GridBase *grid = Umu._grid;
-    int lsites = grid->lSites();
+    uint64_t lsites = grid->lSites();
 
     std::vector<sobj> scalardata(lsites); 
     std::vector<fobj>     iodata(lsites); // Munge, checksum, byte order in here
@@ -580,7 +582,7 @@ class BinaryIO {
     GridStopWatch timer; timer.Start();
     unvectorizeToLexOrdArray(scalardata,Umu);    
 
-    parallel_for(int x=0;x<lsites;x++) munge(scalardata[x],iodata[x]);
+    parallel_for(uint64_t x=0;x<lsites;x++) munge(scalardata[x],iodata[x]);
 
     grid->Barrier();
     timer.Stop();
@@ -597,7 +599,7 @@ class BinaryIO {
   static inline void readRNG(GridSerialRNG &serial,
 			     GridParallelRNG &parallel,
 			     std::string file,
-			     Integer offset,
+			     uint64_t offset,
 			     uint32_t &nersc_csum,
 			     uint32_t &scidac_csuma,
 			     uint32_t &scidac_csumb)
@@ -610,8 +612,8 @@ class BinaryIO {
     std::string format = "IEEE32BIG";
 
     GridBase *grid = parallel._grid;
-    int gsites = grid->gSites();
+    uint64_t gsites = grid->gSites();
-    int lsites = grid->lSites();
+    uint64_t lsites = grid->lSites();
 
     uint32_t nersc_csum_tmp   = 0;
     uint32_t scidac_csuma_tmp = 0;
@@ -626,7 +628,7 @@ class BinaryIO {
 	     nersc_csum,scidac_csuma,scidac_csumb);
 
     timer.Start();
-    parallel_for(int lidx=0;lidx<lsites;lidx++){
+    parallel_for(uint64_t lidx=0;lidx<lsites;lidx++){
       std::vector<RngStateType> tmp(RngStateCount);
       std::copy(iodata[lidx].begin(),iodata[lidx].end(),tmp.begin());
       parallel.SetState(tmp,lidx);
@@ -659,7 +661,7 @@ class BinaryIO {
   static inline void writeRNG(GridSerialRNG &serial,
 			      GridParallelRNG &parallel,
 			      std::string file,
-			      Integer offset,
+			      uint64_t offset,
 			      uint32_t &nersc_csum,
 			      uint32_t &scidac_csuma,
 			      uint32_t &scidac_csumb)
@@ -670,8 +672,8 @@ class BinaryIO {
     typedef std::array<RngStateType,RngStateCount> RNGstate;
 
     GridBase *grid = parallel._grid;
-    int gsites = grid->gSites();
+    uint64_t gsites = grid->gSites();
-    int lsites = grid->lSites();
+    uint64_t lsites = grid->lSites();
 
     uint32_t nersc_csum_tmp;
     uint32_t scidac_csuma_tmp;
@@ -684,7 +686,7 @@ class BinaryIO {
 
     timer.Start();
     std::vector<RNGstate> iodata(lsites);
-    parallel_for(int lidx=0;lidx<lsites;lidx++){
+    parallel_for(uint64_t lidx=0;lidx<lsites;lidx++){
       std::vector<RngStateType> tmp(RngStateCount);
       parallel.GetState(tmp,lidx);
       std::copy(tmp.begin(),tmp.end(),iodata[lidx].begin());

lib/parallelIO/IldgIO.h

@@ -337,6 +337,20 @@ class GridLimeWriter : public BinaryIO {
   template<class vobj>
   void writeLimeLatticeBinaryObject(Lattice<vobj> &field,std::string record_name)
   {
+    ////////////////////////////////////////////////////////////////////
+    // NB: FILE and iostream are jointly writing disjoint sequences in the
+    // the same file through different file handles (integer units).
+    // 
+    // These are both buffered, so why I think this code is right is as follows.
+    //
+    // i)  write record header to FILE *File, telegraphing the size; flush
+    // ii) ftello reads the offset from FILE *File . 
+    // iii) iostream / MPI Open independently seek this offset. Write sequence direct to disk.
+    //      Closes iostream and flushes.
+    // iv) fseek on FILE * to end of this disjoint section.
+    //  v) Continue writing scidac record.
+    ////////////////////////////////////////////////////////////////////
+
     ////////////////////////////////////////////
     // Create record header
     ////////////////////////////////////////////
@@ -350,25 +364,24 @@ class GridLimeWriter : public BinaryIO {
     //    std::cout << "W Gsites "           <<field._grid->_gsites<<std::endl;
     //    std::cout << "W Payload expected " <<PayloadSize<<std::endl;
 
-    ////////////////////////////////////////////////////////////////////
+    fflush(File);
-    // NB: FILE and iostream are jointly writing disjoint sequences in the
+
-    // the same file through different file handles (integer units).
+    ///////////////////////////////////////////
-    // 
+    // Write by other means into the binary record
-    // These are both buffered, so why I think this code is right is as follows.
+    ///////////////////////////////////////////
-    //
+    uint64_t offset1 = ftello(File);    //    std::cout << " Writing to offset "<<offset1 << std::endl;
-    // i)  write record header to FILE *File, telegraphing the size. 
-    // ii) ftello reads the offset from FILE *File .
-    // iii) iostream / MPI Open independently seek this offset. Write sequence direct to disk.
-    //      Closes iostream and flushes.
-    // iv) fseek on FILE * to end of this disjoint section.
-    //  v) Continue writing scidac record.
-    ////////////////////////////////////////////////////////////////////
-    uint64_t offset = ftello(File);
-    //    std::cout << " Writing to offset "<<offset << std::endl;
     std::string format = getFormatString<vobj>();
     BinarySimpleMunger<sobj,sobj> munge;
-    BinaryIO::writeLatticeObject<vobj,sobj>(field, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb);
+    BinaryIO::writeLatticeObject<vobj,sobj>(field, filename, munge, offset1, format,nersc_csum,scidac_csuma,scidac_csumb);
-    //    fseek(File,0,SEEK_END);    offset = ftello(File);std::cout << " offset now "<<offset << std::endl;
+
+    ///////////////////////////////////////////
+    // Wind forward and close the record
+    ///////////////////////////////////////////
+    fseek(File,0,SEEK_END);             
+    uint64_t offset2 = ftello(File);     //    std::cout << " now at offset "<<offset2 << std::endl;
+
+    assert((offset2-offset1) == PayloadSize);
+
     err=limeWriterCloseRecord(LimeW);  assert(err>=0);
 
     ////////////////////////////////////////
@@ -568,7 +581,6 @@ class IldgWriter : public ScidacWriter {
     writeLimeIldgLFN(header.ildg_lfn);                                                 // rec
     writeLimeLatticeBinaryObject(Umu,std::string(ILDG_BINARY_DATA));      // Closes message with checksum
     //    limeDestroyWriter(LimeW);
-    fclose(File);
   }
 };
 

lib/parallelIO/NerscIO.h

@@ -57,7 +57,7 @@ namespace Grid {
       // for the header-reader
       static inline int readHeader(std::string file,GridBase *grid,  FieldMetaData &field)
       {
-      int offset=0;
+      uint64_t offset=0;
       std::map<std::string,std::string> header;
       std::string line;
 
@@ -139,7 +139,7 @@ namespace Grid {
       typedef Lattice<iLorentzColourMatrix<vsimd> > GaugeField;
 
       GridBase *grid = Umu._grid;
-      int offset = readHeader(file,Umu._grid,header);
+      uint64_t offset = readHeader(file,Umu._grid,header);
 
       FieldMetaData clone(header);
 
@@ -236,7 +236,7 @@ namespace Grid {
 	GaugeStatistics(Umu,header);
 	MachineCharacteristics(header);
 
-	int offset;
+	uint64_t offset;
   
 	truncate(file);
 
@@ -278,7 +278,7 @@ namespace Grid {
 	header.plaquette=0.0;
 	MachineCharacteristics(header);
 
-	int offset;
+	uint64_t offset;
   
 #ifdef RNG_RANLUX
 	header.floating_point = std::string("UINT64");
@@ -313,7 +313,7 @@ namespace Grid {
 
 	GridBase *grid = parallel._grid;
 
-	int offset = readHeader(file,grid,header);
+	uint64_t offset = readHeader(file,grid,header);
 
 	FieldMetaData clone(header);
 

lib/qcd/action/gauge/WilsonGaugeAction.h

@@ -71,18 +71,14 @@ class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
 
     RealD factor = 0.5 * beta / RealD(Nc);
 
-    //GaugeLinkField Umu(U._grid);
+    GaugeLinkField Umu(U._grid);
     GaugeLinkField dSdU_mu(U._grid);
     for (int mu = 0; mu < Nd; mu++) {
-      //Umu = PeekIndex<LorentzIndex>(U, mu);
+      Umu = PeekIndex<LorentzIndex>(U, mu);
 
       // Staple in direction mu
-      //WilsonLoops<Gimpl>::Staple(dSdU_mu, U, mu);
+      WilsonLoops<Gimpl>::Staple(dSdU_mu, U, mu);
-      //dSdU_mu = Ta(Umu * dSdU_mu) * factor;
+      dSdU_mu = Ta(Umu * dSdU_mu) * factor;
-
-  
-      WilsonLoops<Gimpl>::StapleMult(dSdU_mu, U, mu);
-      dSdU_mu = Ta(dSdU_mu) * factor;
 
       PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
     }

lib/qcd/utils/WilsonLoops.h

@@ -212,6 +212,7 @@ public:
 
 
 // For the force term
+/*
 static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
     GridBase *grid = Umu._grid;
     std::vector<GaugeMat> U(Nd, grid);
@@ -225,7 +226,7 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
 
     for (int nu = 0; nu < Nd; nu++) {
       if (nu != mu) {
-        // this is ~10% faster than the Staple
+        // this is ~10% faster than the Staple  -- PAB: so what it gives the WRONG answers for other BC's!
         tmp1 = Cshift(U[nu], mu, 1);
         tmp2 = Cshift(U[mu], nu, 1);
         staple += tmp1* adj(U[nu]*tmp2);
@@ -235,7 +236,7 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
     }
     staple = U[mu]*staple;
 }
-
+*/
   //////////////////////////////////////////////////
   // the sum over all staples on each site
   //////////////////////////////////////////////////

lib/serialisation/BaseIO.h

@@ -31,161 +31,10 @@ Author: Guido Cossu <guido.cossu@ed.ac.uk>
 #define GRID_SERIALISATION_ABSTRACT_READER_H
 
 #include <type_traits>
+#include <Grid/tensors/Tensors.h>
+#include <Grid/serialisation/VectorUtils.h>
 
 namespace Grid {
-  // Vector IO utilities ///////////////////////////////////////////////////////
-  // helper function to read space-separated values
-  template <typename T>
-  std::vector<T> strToVec(const std::string s)
-  {
-    std::istringstream sstr(s);
-    T                  buf;
-    std::vector<T>     v;
-    
-    while(!sstr.eof())
-    {
-      sstr >> buf;
-      v.push_back(buf);
-    }
-    
-    return v;
-  }
-  
-  // output to streams for vectors
-  template < class T >
-  inline std::ostream & operator<<(std::ostream &os, const std::vector<T> &v)
-  {
-    os << "[";
-    for (auto &x: v)
-    {
-      os << x << " ";
-    }
-    if (v.size() > 0)
-    {
-      os << "\b";
-    }
-    os << "]";
-    
-    return os;
-  }
-  
-  // Vector element trait //////////////////////////////////////////////////////  
-  template <typename T>
-  struct element
-  {
-    typedef T type;
-    static constexpr bool is_number = false;
-  };
-  
-  template <typename T>
-  struct element<std::vector<T>>
-  {
-    typedef typename element<T>::type type;
-    static constexpr bool is_number = std::is_arithmetic<T>::value
-                                      or is_complex<T>::value
-                                      or element<T>::is_number;
-  };
-  
-  // Vector flattening utility class ////////////////////////////////////////////
-  // Class to flatten a multidimensional std::vector
-  template <typename V>
-  class Flatten
-  {
-  public:
-    typedef typename element<V>::type Element;
-  public:
-    explicit                     Flatten(const V &vector);
-    const V &                    getVector(void);
-    const std::vector<Element> & getFlatVector(void);
-    const std::vector<size_t>  & getDim(void);
-  private:
-    void accumulate(const Element &e);
-    template <typename W>
-    void accumulate(const W &v);
-    void accumulateDim(const Element &e);
-    template <typename W>
-    void accumulateDim(const W &v);
-  private:
-    const V              &vector_;
-    std::vector<Element> flatVector_;
-    std::vector<size_t>  dim_;
-  };
-  
-  // Class to reconstruct a multidimensional std::vector
-  template <typename V>
-  class Reconstruct
-  {
-  public:
-    typedef typename element<V>::type Element;
-  public:
-    Reconstruct(const std::vector<Element> &flatVector,
-                const std::vector<size_t> &dim);
-    const V &                    getVector(void);
-    const std::vector<Element> & getFlatVector(void);
-    const std::vector<size_t>  & getDim(void);
-  private:
-    void fill(std::vector<Element> &v);
-    template <typename W>
-    void fill(W &v);
-    void resize(std::vector<Element> &v, const unsigned int dim);
-    template <typename W>
-    void resize(W &v, const unsigned int dim);
-  private:
-    V                          vector_;
-    const std::vector<Element> &flatVector_;
-    std::vector<size_t>        dim_;
-    size_t                     ind_{0};
-    unsigned int               dimInd_{0};
-  };
-  
-  // Pair IO utilities /////////////////////////////////////////////////////////
-  // helper function to parse input in the format "<obj1 obj2>"
-  template <typename T1, typename T2>
-  inline std::istream & operator>>(std::istream &is, std::pair<T1, T2> &buf)
-  {
-    T1 buf1;
-    T2 buf2;
-    char c;
-
-    // Search for "pair" delimiters.
-    do
-    {
-      is.get(c);
-    } while (c != '<' && !is.eof());
-    if (c == '<')
-    {
-      int start = is.tellg();
-      do
-      {
-        is.get(c);
-      } while (c != '>' && !is.eof());
-      if (c == '>')
-      {
-        int end = is.tellg();
-        int psize = end - start - 1;
-
-        // Only read data between pair limiters.
-        is.seekg(start);
-        std::string tmpstr(psize, ' ');
-        is.read(&tmpstr[0], psize);
-        std::istringstream temp(tmpstr);
-        temp >> buf1 >> buf2;
-        buf = std::make_pair(buf1, buf2);
-        is.seekg(end);
-      }
-    }
-    is.peek();
-    return is;
-  }
-  
-  // output to streams for pairs
-  template <class T1, class T2>
-  inline std::ostream & operator<<(std::ostream &os, const std::pair<T1, T2> &p)
-  {
-    os << "<" << p.first << " " << p.second << ">";
-    return os;
-  }
-
   // Abstract writer/reader classes ////////////////////////////////////////////
   // static polymorphism implemented using CRTP idiom
   class Serializable;
@@ -205,6 +54,12 @@ namespace Grid {
     template <typename U>
     typename std::enable_if<!std::is_base_of<Serializable, U>::value, void>::type
     write(const std::string& s, const U &output);
+    template <typename U>
+    void write(const std::string &s, const iScalar<U> &output);
+    template <typename U, int N>
+    void write(const std::string &s, const iVector<U, N> &output);
+    template <typename U, int N>
+    void write(const std::string &s, const iMatrix<U, N> &output);
   private:
     T *upcast;
   };
@@ -224,6 +79,12 @@ namespace Grid {
     template <typename U>
     typename std::enable_if<!std::is_base_of<Serializable, U>::value, void>::type
     read(const std::string& s, U &output);
+    template <typename U>
+    void read(const std::string &s, iScalar<U> &output);
+    template <typename U, int N>
+    void read(const std::string &s, iVector<U, N> &output);
+    template <typename U, int N>
+    void read(const std::string &s, iMatrix<U, N> &output);
   protected:
     template <typename U>
     void fromString(U &output, const std::string &s);
@@ -237,203 +98,9 @@ namespace Grid {
   };
   template<typename T> struct isWriter {
     static const bool value = false;
-  }; 
-
-
-
-  // Generic writer interface
-  // serializable base class
-  class Serializable
-  {
-  public:
-    template <typename T>
-    static inline void write(Writer<T> &WR,const std::string &s,
-                             const Serializable &obj)
-    {}
-    
-    template <typename T>
-    static inline void read(Reader<T> &RD,const std::string &s,
-                            Serializable &obj)
-    {}
-    
-    friend inline std::ostream & operator<<(std::ostream &os,
-                                            const Serializable &obj)
-    {
-      return os;
-    }
   };
-  
+
-  // Flatten class template implementation /////////////////////////////////////
+  // Writer template implementation
-  template <typename V>
-  void Flatten<V>::accumulate(const Element &e)
-  {
-    flatVector_.push_back(e);
-  }
-  
-  template <typename V>
-  template <typename W>
-  void Flatten<V>::accumulate(const W &v)
-  {
-    for (auto &e: v)
-    {
-      accumulate(e);
-    }
-  }
-  
-  template <typename V>
-  void Flatten<V>::accumulateDim(const Element &e) {};
-  
-  template <typename V>
-  template <typename W>
-  void Flatten<V>::accumulateDim(const W &v)
-  {
-    dim_.push_back(v.size());
-    accumulateDim(v[0]);
-  }
-  
-  template <typename V>
-  Flatten<V>::Flatten(const V &vector)
-  : vector_(vector)
-  {
-    accumulate(vector_);
-    accumulateDim(vector_);
-  }
-  
-  template <typename V>
-  const V & Flatten<V>::getVector(void)
-  {
-    return vector_;
-  }
-  
-  template <typename V>
-  const std::vector<typename Flatten<V>::Element> &
-  Flatten<V>::getFlatVector(void)
-  {
-    return flatVector_;
-  }
-  
-  template <typename V>
-  const std::vector<size_t> & Flatten<V>::getDim(void)
-  {
-    return dim_;
-  }
-  
-  // Reconstruct class template implementation /////////////////////////////////
-  template <typename V>
-  void Reconstruct<V>::fill(std::vector<Element> &v)
-  {
-    for (auto &e: v)
-    {
-      e = flatVector_[ind_++];
-    }
-  }
-  
-  template <typename V>
-  template <typename W>
-  void Reconstruct<V>::fill(W &v)
-  {
-    for (auto &e: v)
-    {
-      fill(e);
-    }
-  }
-  
-  template <typename V>
-  void Reconstruct<V>::resize(std::vector<Element> &v, const unsigned int dim)
-  {
-    v.resize(dim_[dim]);
-  }
-  
-  template <typename V>
-  template <typename W>
-  void Reconstruct<V>::resize(W &v, const unsigned int dim)
-  {
-    v.resize(dim_[dim]);
-    for (auto &e: v)
-    {
-      resize(e, dim + 1);
-    }
-  }
-  
-  template <typename V>
-  Reconstruct<V>::Reconstruct(const std::vector<Element> &flatVector,
-                              const std::vector<size_t> &dim)
-  : flatVector_(flatVector)
-  , dim_(dim)
-  {
-    resize(vector_, 0);
-    fill(vector_);
-  }
-  
-  template <typename V>
-  const V & Reconstruct<V>::getVector(void)
-  {
-    return vector_;
-  }
-  
-  template <typename V>
-  const std::vector<typename Reconstruct<V>::Element> &
-  Reconstruct<V>::getFlatVector(void)
-  {
-    return flatVector_;
-  }
-  
-  template <typename V>
-  const std::vector<size_t> & Reconstruct<V>::getDim(void)
-  {
-    return dim_;
-  }
-  
-  // Generic writer interface //////////////////////////////////////////////////
-  template <typename T>
-  inline void push(Writer<T> &w, const std::string &s) {
-    w.push(s);
-  }
-  
-  template <typename T>
-  inline void push(Writer<T> &w, const char *s)
-  {
-    w.push(std::string(s));
-  }
-  
-  template <typename T>
-  inline void pop(Writer<T> &w)
-  {
-    w.pop();
-  }
-  
-  template <typename T, typename U>
-  inline void write(Writer<T> &w, const std::string& s, const U &output)
-  {
-    w.write(s, output);
-  }
-  
-  // Generic reader interface
-  template <typename T>
-  inline bool push(Reader<T> &r, const std::string &s)
-  {
-    return r.push(s);
-  }
-  
-  template <typename T>
-  inline bool push(Reader<T> &r, const char *s)
-  {
-    return r.push(std::string(s));
-  }
-  
-  template <typename T>
-  inline void pop(Reader<T> &r)
-  {
-    r.pop();
-  }
-  
-  template <typename T, typename U>
-  inline void read(Reader<T> &r, const std::string &s, U &output)
-  {
-    r.read(s, output);
-  }
-  
-  // Writer template implementation ////////////////////////////////////////////
   template <typename T>
   Writer<T>::Writer(void)
   {
@@ -467,6 +134,27 @@ namespace Grid {
   {
     upcast->writeDefault(s, output);
   }
+
+  template <typename T>
+  template <typename U>
+  void Writer<T>::write(const std::string &s, const iScalar<U> &output)
+  {
+    upcast->writeDefault(s, tensorToVec(output));
+  }
+
+  template <typename T>
+  template <typename U, int N>
+  void Writer<T>::write(const std::string &s, const iVector<U, N> &output)
+  {
+    upcast->writeDefault(s, tensorToVec(output));
+  }
+
+  template <typename T>
+  template <typename U, int N>
+  void Writer<T>::write(const std::string &s, const iMatrix<U, N> &output)
+  {
+    upcast->writeDefault(s, tensorToVec(output));
+  }
   
   // Reader template implementation
   template <typename T>
@@ -502,7 +190,37 @@ namespace Grid {
   {
     upcast->readDefault(s, output);
   }
+
+  template <typename T>
+  template <typename U>
+  void Reader<T>::read(const std::string &s, iScalar<U> &output)
+  {
+    typename TensorToVec<iScalar<U>>::type v;
+
+    upcast->readDefault(s, v);
+    vecToTensor(output, v);
+  }
+
+  template <typename T>
+  template <typename U, int N>
+  void Reader<T>::read(const std::string &s, iVector<U, N> &output)
+  {
+    typename TensorToVec<iVector<U, N>>::type v;
+    
+    upcast->readDefault(s, v);
+    vecToTensor(output, v);
+  }
   
+  template <typename T>
+  template <typename U, int N>
+  void Reader<T>::read(const std::string &s, iMatrix<U, N> &output)
+  {
+    typename TensorToVec<iMatrix<U, N>>::type v;
+    
+    upcast->readDefault(s, v);
+    vecToTensor(output, v);
+  }
+
   template <typename T>
   template <typename U>
   void Reader<T>::fromString(U &output, const std::string &s)
@@ -514,13 +232,83 @@ namespace Grid {
     {
       is >> std::boolalpha >> output;
     }
-    catch(std::istringstream::failure &e)
+    catch(std::ios_base::failure &e)
     {
       std::cerr << "numerical conversion failure on '" << s << "' ";
       std::cerr << "(typeid: " << typeid(U).name() << ")" << std::endl;
       abort();
     }
   }
+
+  // serializable base class ///////////////////////////////////////////////////
+  class Serializable
+  {
+  public:
+    template <typename T>
+    static inline void write(Writer<T> &WR,const std::string &s,
+                             const Serializable &obj)
+    {}
+    
+    template <typename T>
+    static inline void read(Reader<T> &RD,const std::string &s,
+                            Serializable &obj)
+    {}
+    
+    friend inline std::ostream & operator<<(std::ostream &os,
+                                            const Serializable &obj)
+    {
+      return os;
+    }
+  };
+  
+  // Generic writer interface //////////////////////////////////////////////////
+  template <typename T>
+  inline void push(Writer<T> &w, const std::string &s) {
+    w.push(s);
+  }
+  
+  template <typename T>
+  inline void push(Writer<T> &w, const char *s)
+  {
+    w.push(std::string(s));
+  }
+  
+  template <typename T>
+  inline void pop(Writer<T> &w)
+  {
+    w.pop();
+  }
+  
+  template <typename T, typename U>
+  inline void write(Writer<T> &w, const std::string& s, const U &output)
+  {
+    w.write(s, output);
+  }
+  
+  // Generic reader interface //////////////////////////////////////////////////
+  template <typename T>
+  inline bool push(Reader<T> &r, const std::string &s)
+  {
+    return r.push(s);
+  }
+  
+  template <typename T>
+  inline bool push(Reader<T> &r, const char *s)
+  {
+    return r.push(std::string(s));
+  }
+  
+  template <typename T>
+  inline void pop(Reader<T> &r)
+  {
+    r.pop();
+  }
+  
+  template <typename T, typename U>
+  inline void read(Reader<T> &r, const std::string &s, U &output)
+  {
+    r.read(s, output);
+  }
 }
 
 #endif

lib/serialisation/Hdf5IO.h

@@ -5,6 +5,7 @@
 #include <string>
 #include <vector>
 #include <H5Cpp.h>
+#include <Grid/tensors/Tensors.h>
 #include "Hdf5Type.h"
 
 #ifndef H5_NO_NAMESPACE

lib/serialisation/VectorUtils.h

@@ -0,0 +1,384 @@
+#ifndef GRID_SERIALISATION_VECTORUTILS_H
+#define GRID_SERIALISATION_VECTORUTILS_H
+
+#include <type_traits>
+#include <Grid/tensors/Tensors.h>
+
+namespace Grid {
+  // Pair IO utilities /////////////////////////////////////////////////////////
+  // helper function to parse input in the format "<obj1 obj2>"
+  template <typename T1, typename T2>
+  inline std::istream & operator>>(std::istream &is, std::pair<T1, T2> &buf)
+  {
+    T1 buf1;
+    T2 buf2;
+    char c;
+
+    // Search for "pair" delimiters.
+    do
+    {
+      is.get(c);
+    } while (c != '(' && !is.eof());
+    if (c == '(')
+    {
+      int start = is.tellg();
+      do
+      {
+        is.get(c);
+      } while (c != ')' && !is.eof());
+      if (c == ')')
+      {
+        int end = is.tellg();
+        int psize = end - start - 1;
+
+        // Only read data between pair limiters.
+        is.seekg(start);
+        std::string tmpstr(psize, ' ');
+        is.read(&tmpstr[0], psize);
+        std::istringstream temp(tmpstr);
+        temp >> buf1 >> buf2;
+        buf = std::make_pair(buf1, buf2);
+        is.seekg(end);
+      }
+    }
+    is.peek();
+    return is;
+  }
+  
+  // output to streams for pairs
+  template <class T1, class T2>
+  inline std::ostream & operator<<(std::ostream &os, const std::pair<T1, T2> &p)
+  {
+    os << "(" << p.first << " " << p.second << ")";
+    return os;
+  }
+  
+  // Grid scalar tensors to nested std::vectors //////////////////////////////////
+  template <typename T>
+  struct TensorToVec
+  {
+    typedef T type;
+  };
+
+  template <typename T>
+  struct TensorToVec<iScalar<T>>
+  {
+    typedef typename TensorToVec<T>::type type;
+  };
+
+  template <typename T, int N>
+  struct TensorToVec<iVector<T, N>>
+  {
+    typedef typename std::vector<typename TensorToVec<T>::type> type;
+  };
+
+  template <typename T, int N>
+  struct TensorToVec<iMatrix<T, N>>
+  {
+    typedef typename std::vector<std::vector<typename TensorToVec<T>::type>> type;
+  };
+
+  template <typename T>
+  typename TensorToVec<T>::type tensorToVec(const T &t)
+  {
+    return t;
+  }
+
+  template <typename T>
+  typename TensorToVec<iScalar<T>>::type tensorToVec(const iScalar<T>& t)
+  {
+    return tensorToVec(t._internal);
+  }
+
+  template <typename T, int N>
+  typename TensorToVec<iVector<T, N>>::type tensorToVec(const iVector<T, N>& t)
+  {
+    typename TensorToVec<iVector<T, N>>::type v;
+
+    v.resize(N);
+    for (unsigned int i = 0; i < N; i++) 
+    {
+      v[i] = tensorToVec(t._internal[i]);
+    }
+
+    return v;
+  }
+
+  template <typename T, int N>
+  typename TensorToVec<iMatrix<T, N>>::type tensorToVec(const iMatrix<T, N>& t)
+  {
+    typename TensorToVec<iMatrix<T, N>>::type v;
+
+    v.resize(N);
+    for (unsigned int i = 0; i < N; i++)
+    {
+      v[i].resize(N);
+      for (unsigned int j = 0; j < N; j++) 
+      {
+        v[i][j] = tensorToVec(t._internal[i][j]);
+      }
+    }
+
+    return v;
+  }
+
+  template <typename T>
+  void vecToTensor(T &t, const typename TensorToVec<T>::type &v)
+  {
+    t = v;
+  }
+
+
+  template <typename T>
+  void vecToTensor(iScalar<T> &t, const typename TensorToVec<iScalar<T>>::type &v)
+  {
+    vecToTensor(t._internal, v);
+  }
+
+  template <typename T, int N>
+  void vecToTensor(iVector<T, N> &t, const typename TensorToVec<iVector<T, N>>::type &v)
+  {
+    for (unsigned int i = 0; i < N; i++) 
+    {
+      vecToTensor(t._internal[i], v[i]);
+    }
+  }
+
+  template <typename T, int N>
+  void vecToTensor(iMatrix<T, N> &t, const typename TensorToVec<iMatrix<T, N>>::type &v)
+  {
+    for (unsigned int i = 0; i < N; i++)
+    for (unsigned int j = 0; j < N; j++)
+    {
+      vecToTensor(t._internal[i][j], v[i][j]);
+    }
+  }
+
+  // Vector element trait //////////////////////////////////////////////////////  
+  template <typename T>
+  struct element
+  {
+    typedef T type;
+    static constexpr bool is_number = false;
+  };
+  
+  template <typename T>
+  struct element<std::vector<T>>
+  {
+    typedef typename element<T>::type type;
+    static constexpr bool is_number = std::is_arithmetic<T>::value
+                                      or is_complex<T>::value
+                                      or element<T>::is_number;
+  };
+  
+  // Vector flattening utility class ////////////////////////////////////////////
+  // Class to flatten a multidimensional std::vector
+  template <typename V>
+  class Flatten
+  {
+  public:
+    typedef typename element<V>::type Element;
+  public:
+    explicit                     Flatten(const V &vector);
+    const V &                    getVector(void);
+    const std::vector<Element> & getFlatVector(void);
+    const std::vector<size_t>  & getDim(void);
+  private:
+    void accumulate(const Element &e);
+    template <typename W>
+    void accumulate(const W &v);
+    void accumulateDim(const Element &e);
+    template <typename W>
+    void accumulateDim(const W &v);
+  private:
+    const V              &vector_;
+    std::vector<Element> flatVector_;
+    std::vector<size_t>  dim_;
+  };
+  
+  // Class to reconstruct a multidimensional std::vector
+  template <typename V>
+  class Reconstruct
+  {
+  public:
+    typedef typename element<V>::type Element;
+  public:
+    Reconstruct(const std::vector<Element> &flatVector,
+                const std::vector<size_t> &dim);
+    const V &                    getVector(void);
+    const std::vector<Element> & getFlatVector(void);
+    const std::vector<size_t>  & getDim(void);
+  private:
+    void fill(std::vector<Element> &v);
+    template <typename W>
+    void fill(W &v);
+    void resize(std::vector<Element> &v, const unsigned int dim);
+    template <typename W>
+    void resize(W &v, const unsigned int dim);
+  private:
+    V                          vector_;
+    const std::vector<Element> &flatVector_;
+    std::vector<size_t>        dim_;
+    size_t                     ind_{0};
+    unsigned int               dimInd_{0};
+  };
+
+  // Flatten class template implementation
+  template <typename V>
+  void Flatten<V>::accumulate(const Element &e)
+  {
+    flatVector_.push_back(e);
+  }
+  
+  template <typename V>
+  template <typename W>
+  void Flatten<V>::accumulate(const W &v)
+  {
+    for (auto &e: v)
+    {
+      accumulate(e);
+    }
+  }
+  
+  template <typename V>
+  void Flatten<V>::accumulateDim(const Element &e) {};
+  
+  template <typename V>
+  template <typename W>
+  void Flatten<V>::accumulateDim(const W &v)
+  {
+    dim_.push_back(v.size());
+    accumulateDim(v[0]);
+  }
+  
+  template <typename V>
+  Flatten<V>::Flatten(const V &vector)
+  : vector_(vector)
+  {
+    accumulate(vector_);
+    accumulateDim(vector_);
+  }
+  
+  template <typename V>
+  const V & Flatten<V>::getVector(void)
+  {
+    return vector_;
+  }
+  
+  template <typename V>
+  const std::vector<typename Flatten<V>::Element> &
+  Flatten<V>::getFlatVector(void)
+  {
+    return flatVector_;
+  }
+  
+  template <typename V>
+  const std::vector<size_t> & Flatten<V>::getDim(void)
+  {
+    return dim_;
+  }
+  
+  // Reconstruct class template implementation
+  template <typename V>
+  void Reconstruct<V>::fill(std::vector<Element> &v)
+  {
+    for (auto &e: v)
+    {
+      e = flatVector_[ind_++];
+    }
+  }
+  
+  template <typename V>
+  template <typename W>
+  void Reconstruct<V>::fill(W &v)
+  {
+    for (auto &e: v)
+    {
+      fill(e);
+    }
+  }
+  
+  template <typename V>
+  void Reconstruct<V>::resize(std::vector<Element> &v, const unsigned int dim)
+  {
+    v.resize(dim_[dim]);
+  }
+  
+  template <typename V>
+  template <typename W>
+  void Reconstruct<V>::resize(W &v, const unsigned int dim)
+  {
+    v.resize(dim_[dim]);
+    for (auto &e: v)
+    {
+      resize(e, dim + 1);
+    }
+  }
+  
+  template <typename V>
+  Reconstruct<V>::Reconstruct(const std::vector<Element> &flatVector,
+                              const std::vector<size_t> &dim)
+  : flatVector_(flatVector)
+  , dim_(dim)
+  {
+    resize(vector_, 0);
+    fill(vector_);
+  }
+  
+  template <typename V>
+  const V & Reconstruct<V>::getVector(void)
+  {
+    return vector_;
+  }
+  
+  template <typename V>
+  const std::vector<typename Reconstruct<V>::Element> &
+  Reconstruct<V>::getFlatVector(void)
+  {
+    return flatVector_;
+  }
+  
+  template <typename V>
+  const std::vector<size_t> & Reconstruct<V>::getDim(void)
+  {
+    return dim_;
+  }
+
+  // Vector IO utilities ///////////////////////////////////////////////////////
+  // helper function to read space-separated values
+  template <typename T>
+  std::vector<T> strToVec(const std::string s)
+  {
+    std::istringstream sstr(s);
+    T                  buf;
+    std::vector<T>     v;
+    
+    while(!sstr.eof())
+    {
+      sstr >> buf;
+      v.push_back(buf);
+    }
+    
+    return v;
+  }
+  
+  // output to streams for vectors
+  template < class T >
+  inline std::ostream & operator<<(std::ostream &os, const std::vector<T> &v)
+  {
+    os << "[";
+    for (auto &x: v)
+    {
+      os << x << " ";
+    }
+    if (v.size() > 0)
+    {
+      os << "\b";
+    }
+    os << "]";
+    
+    return os;
+  }
+}
+
+#endif

lib/tensors/Tensor_logical.h

@@ -55,5 +55,38 @@ LOGICAL_BINOP(&);
 LOGICAL_BINOP(||);
 LOGICAL_BINOP(&&);
 
+template <class T>
+strong_inline bool operator==(const iScalar<T> &t1, const iScalar<T> &t2)
+{
+  return (t1._internal == t2._internal);
+}
+
+template <class T, int N>
+strong_inline bool operator==(const iVector<T, N> &t1, const iVector<T, N> &t2)
+{
+  bool res = true;
+
+  for (unsigned int i = 0; i < N; ++i)
+  {
+    res = (res && (t1._internal[i] == t2._internal[i]));
+  }
+
+  return res;
+}
+
+template <class T, int N>
+strong_inline bool operator==(const iMatrix<T, N> &t1, const iMatrix<T, N> &t2)
+{
+  bool res = true;
+
+  for (unsigned int i = 0; i < N; ++i)
+  for (unsigned int j = 0; j < N; ++j)
+  {
+    res = (res && (t1._internal[i][j] == t2._internal[i][j]));
+  }
+  
+  return res;
+}
+
 }
 #endif

lib/util/Init.cc

@@ -49,6 +49,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/Grid.h>
 
 #include <Grid/util/CompilerCompatible.h>
+#include <version.h>
 
 
 #include <fenv.h>
@@ -288,6 +289,12 @@ void Grid_init(int *argc,char ***argv)
     std::cout << "but WITHOUT ANY WARRANTY; without even the implied warranty of"<<std::endl;
     std::cout << "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the"<<std::endl;
     std::cout << "GNU General Public License for more details."<<std::endl;
+#ifdef GITHASH
+    std::cout << "Current Grid git commit hash=" << GITHASH << std::endl;
+#else
+    std::cout << "Current Grid git commit hash is undefined. Check makefile." << std::endl;
+#endif
+#undef GITHASH
     std::cout << std::endl;
   }
 

tests/IO/Test_serialisation.cc

@@ -45,7 +45,8 @@ public:
                           bool , b,
                           std::vector<double>, array,
                           std::vector<std::vector<double> >, twodimarray,
-                          std::vector<std::vector<std::vector<Complex> > >, cmplx3darray
+                          std::vector<std::vector<std::vector<Complex> > >, cmplx3darray,
+                          SpinColourMatrix, scm
                           );
   myclass() {}
   myclass(int i)
@@ -59,6 +60,12 @@ public:
     y=2*i;
     b=true;
     name="bother said pooh";
+    scm()(0, 1)(2, 1) = 2.356;
+    scm()(3, 0)(1, 1) = 1.323;
+    scm()(2, 1)(0, 1) = 5.3336;
+    scm()(0, 2)(1, 1) = 6.336;
+    scm()(2, 1)(2, 2) = 7.344;
+    scm()(1, 1)(2, 0) = 8.3534;
   }
 };
 
@@ -93,8 +100,30 @@ void ioTest(const std::string &filename, const O &object, const std::string &nam
   if (!good) exit(EXIT_FAILURE);
 }
 
+template <typename T>
+void tensorConvTestFn(GridSerialRNG &rng, const std::string label)
+{
+  T    t, ft;
+  Real n;
+  bool good;
+
+  random(rng, t);
+  auto tv = tensorToVec(t);
+  vecToTensor(ft, tv);
+  n    = norm2(t - ft);
+  good = (n == 0);
+  std::cout << label << " norm 2 diff: " << n << " -- " 
+            << (good ? "success" : "failure") << std::endl;
+}
+
+#define tensorConvTest(rng, type) tensorConvTestFn<type>(rng, #type)
+
 int main(int argc,char **argv)
 {
+  GridSerialRNG    rng;
+
+  rng.SeedFixedIntegers(std::vector<int>({42,10,81,9}));
+  
   std::cout << "==== basic IO" << std::endl;
   XmlWriter WR("bother.xml");
 
@@ -120,7 +149,7 @@ int main(int argc,char **argv)
   std::cout << "-- serialisable class writing to 'bother.xml'..." << std::endl;
   write(WR,"obj",obj);
   WR.write("obj2", obj);
-  vec.push_back(myclass(1234));
+  vec.push_back(obj);
   vec.push_back(myclass(5678));
   vec.push_back(myclass(3838));
   pair = std::make_pair(myenum::red, myenum::blue);
@@ -131,8 +160,6 @@ int main(int argc,char **argv)
   std::cout << "-- serialisable class comparison:" << std::endl;
   std::cout << "vec[0] == obj: " << ((vec[0] == obj) ? "true" : "false") << std::endl;
   std::cout << "vec[1] == obj: " << ((vec[1] == obj) ? "true" : "false") << std::endl;
-
-  write(WR, "objpair", pair);
   std::cout << "-- pair writing to std::cout:" << std::endl;
   std::cout << pair << std::endl;
 
@@ -141,26 +168,20 @@ int main(int argc,char **argv)
   //// XML
   ioTest<XmlWriter, XmlReader>("iotest.xml", obj, "XML    (object)           ");
   ioTest<XmlWriter, XmlReader>("iotest.xml", vec, "XML    (vector of objects)");
-  ioTest<XmlWriter, XmlReader>("iotest.xml", pair, "XML    (pair of objects)");
   //// binary
   ioTest<BinaryWriter, BinaryReader>("iotest.bin", obj, "binary (object)           ");
   ioTest<BinaryWriter, BinaryReader>("iotest.bin", vec, "binary (vector of objects)");
-  ioTest<BinaryWriter, BinaryReader>("iotest.bin", pair, "binary (pair of objects)");
   //// text
   ioTest<TextWriter, TextReader>("iotest.dat", obj, "text   (object)           ");
   ioTest<TextWriter, TextReader>("iotest.dat", vec, "text   (vector of objects)");
-  ioTest<TextWriter, TextReader>("iotest.dat", pair, "text   (pair of objects)");
   //// text
   ioTest<JSONWriter, JSONReader>("iotest.json", obj,  "JSON   (object)           ");
   ioTest<JSONWriter, JSONReader>("iotest.json", vec,  "JSON   (vector of objects)");
-  ioTest<JSONWriter, JSONReader>("iotest.json", pair, "JSON   (pair of objects)");
 
   //// HDF5
-#undef HAVE_HDF5
 #ifdef HAVE_HDF5
   ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", obj, "HDF5   (object)           ");
   ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", vec, "HDF5   (vector of objects)");
-  ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", pair, "HDF5   (pair of objects)");
 #endif
 
   std::cout << "\n==== vector flattening/reconstruction" << std::endl;
@@ -197,68 +218,11 @@ int main(int argc,char **argv)
   std::cout << flatdv.getVector() << std::endl;
   std::cout << std::endl;
 
-
+  std::cout << "==== Grid tensor to vector test" << std::endl;
-  std::cout << ".:::::: Testing JSON classes "<< std::endl;
+  tensorConvTest(rng, SpinColourMatrix);
-
+  tensorConvTest(rng, SpinColourVector);
-
+  tensorConvTest(rng, ColourMatrix);
-  {
+  tensorConvTest(rng, ColourVector);
-    JSONWriter JW("bother.json");
+  tensorConvTest(rng, SpinMatrix);
-
+  tensorConvTest(rng, SpinVector);
-    // test basic type writing
-    myenum a = myenum::red;
-    push(JW,"BasicTypes");
-    write(JW,std::string("i16"),i16);
-    write(JW,"myenum",a);
-    write(JW,"u16",u16);
-    write(JW,"i32",i32);
-    write(JW,"u32",u32);
-    write(JW,"i64",i64);
-    write(JW,"u64",u64);
-    write(JW,"f",f);
-    write(JW,"d",d);
-    write(JW,"b",b);
-    pop(JW);
-
-
-    // test serializable class writing
-    myclass obj(1234); // non-trivial constructor
-    std::cout << obj << std::endl;
-    std::cout << "-- serialisable class writing to 'bother.json'..." << std::endl;
-    write(JW,"obj",obj);
-    JW.write("obj2", obj);
-
-
-    std::vector<myclass> vec;
-    vec.push_back(myclass(1234));
-    vec.push_back(myclass(5678));
-    vec.push_back(myclass(3838));
-    write(JW, "objvec", vec);
-
-  }
-
-
-  {
-    JSONReader RD("bother.json");
-    myclass jcopy1;
-    std::vector<myclass> jveccopy1;
-    read(RD,"obj",jcopy1);
-    read(RD,"objvec", jveccopy1);
-    std::cout << "Loaded (JSON) -----------------" << std::endl;
-    std::cout << jcopy1 << std::endl << jveccopy1 << std::endl;
-  }
- 
-
-/*
-  // This is still work in progress
-  {
-    // Testing the next element function
-    JSONReader RD("test.json");
-    RD.push("grid");
-    RD.push("Observable");
-    std::string name;
-    read(RD,"name", name);
-  }
-*/
-
-
 }

tests/forces/Test_gp_rect_force.cc

@@ -59,8 +59,8 @@ int main (int argc, char ** argv)
   double beta = 1.0;
   double c1   = 0.331;
 
-  //GparityPlaqPlusRectangleActionR Action(beta,c1);
+  ConjugatePlaqPlusRectangleActionR Action(beta,c1);
-  ConjugateWilsonGaugeActionR Action(beta);
+  //  ConjugateWilsonGaugeActionR Action(beta);
   //WilsonGaugeActionR Action(beta);
 
   ComplexD S    = Action.S(U);

tests/forces/Test_gpwilson_force.cc

@@ -91,7 +91,7 @@ int main (int argc, char ** argv)
   ////////////////////////////////////
   // Modify the gauge field a little 
   ////////////////////////////////////
-  RealD dt = 0.0001;
+  RealD dt = 0.01;
 
   LatticeColourMatrix mommu(UGrid); 
   LatticeColourMatrix forcemu(UGrid); 

tests/hadrons/Test_QED.cc

@@ -93,8 +93,9 @@ int main(int argc, char *argv[])
         
         // solvers
         MSolver::RBPrecCG::Par solverPar;
-        solverPar.action   = "DWF_" + flavour[i];
+        solverPar.action       = "DWF_" + flavour[i];
-        solverPar.residual = 1.0e-8;
+        solverPar.residual     = 1.0e-8;
+        solverPar.maxIteration = 10000;
         application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
                                                     solverPar);
         

tests/hadrons/Test_hadrons.hpp

@@ -176,8 +176,9 @@ inline void makeRBPrecCGSolver(Application &application, std::string &solverName
     if (!(VirtualMachine::getInstance().hasModule(solverName)))
     {
         MSolver::RBPrecCG::Par solverPar;
-        solverPar.action   = actionName;
+        solverPar.action       = actionName;
-        solverPar.residual = residual;
+        solverPar.residual     = residual;
+        solverPar.maxIteration = 10000;
         application.createModule<MSolver::RBPrecCG>(solverName,
                                                     solverPar);
     }

tests/hadrons/Test_hadrons_2AS_spectrum.cc

@@ -106,8 +106,9 @@ int main(int argc, char *argv[])
         
         // solvers
         MSolver::RBPrecCG2AS::Par solverPar;
-        solverPar.action   = "WilsonClover2AS_" + flavour[i];
+        solverPar.action       = "WilsonClover2AS_" + flavour[i];
-        solverPar.residual = 1.0e-8;
+        solverPar.residual     = 1.0e-8;
+        solverPar.maxIteration = 10000;
         application.createModule<MSolver::RBPrecCG2AS>("CG_" + flavour[i],
                                                     solverPar);
         

tests/hadrons/Test_hadrons_meson_3pt.cc

@@ -82,8 +82,9 @@ int main(int argc, char *argv[])
         
         // solvers
         MSolver::RBPrecCG::Par solverPar;
-        solverPar.action   = "DWF_" + flavour[i];
+        solverPar.action       = "DWF_" + flavour[i];
-        solverPar.residual = 1.0e-8;
+        solverPar.residual     = 1.0e-8;
+        solverPar.maxIteration = 10000;
         application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
                                                     solverPar);
     }

tests/hadrons/Test_hadrons_spectrum.cc

@@ -84,8 +84,9 @@ int main(int argc, char *argv[])
         
         // solvers
         MSolver::RBPrecCG::Par solverPar;
-        solverPar.action   = "DWF_" + flavour[i];
+        solverPar.action       = "DWF_" + flavour[i];
-        solverPar.residual = 1.0e-8;
+        solverPar.residual     = 1.0e-8;
+        solverPar.maxIteration = 10000;
         application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
                                                     solverPar);
         

tests/hadrons/Test_hadrons_wilsonFund.cc

@@ -96,8 +96,9 @@ int main(int argc, char *argv[])
         
         // solvers
         MSolver::RBPrecCG::Par solverPar;
-        solverPar.action   = "WilsonClover_" + flavour[i];
+        solverPar.action       = "WilsonClover_" + flavour[i];
-        solverPar.residual = 1.0e-8;
+        solverPar.residual     = 1.0e-8;
+        solverPar.maxIteration = 10000;
         application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
                                                     solverPar);
         

tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc

@@ -180,7 +180,6 @@ int main (int argc, char ** argv) {
   GridCartesian         * CoarseGrid4    = SpaceTimeGrid::makeFourDimGrid(coarseLatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
   GridRedBlackCartesian * CoarseGrid4rb  = SpaceTimeGrid::makeFourDimRedBlackGrid(CoarseGrid4);
   GridCartesian         * CoarseGrid5    = SpaceTimeGrid::makeFiveDimGrid(cLs,CoarseGrid4);
-  GridRedBlackCartesian * CoarseGrid5rb  = SpaceTimeGrid::makeFourDimRedBlackGrid(CoarseGrid5);
 
   // Gauge field
   LatticeGaugeField Umu(UGrid);
@@ -206,7 +205,7 @@ int main (int argc, char ** argv) {
 
   const int nbasis= 60;
   assert(nbasis==Ns1);
-  LocalCoherenceLanczosScidac<vSpinColourVector,vTComplex,nbasis> _LocalCoherenceLanczos(FrbGrid,CoarseGrid5rb,HermOp,Odd);
+  LocalCoherenceLanczosScidac<vSpinColourVector,vTComplex,nbasis> _LocalCoherenceLanczos(FrbGrid,CoarseGrid5,HermOp,Odd);
   std::cout << GridLogMessage << "Constructed LocalCoherenceLanczos" << std::endl;
 
   assert( (Params.doFine)||(Params.doFineRead));