Merge branch 'develop' into feature/bgq-asm

.gitignore

@@ -9,6 +9,7 @@
 ################
 *~
 *#
+*.sublime-*
 
 # Precompiled Headers #
 #######################
@@ -103,4 +104,16 @@ lib/fftw/*
 # libtool macros #
 ##################
 m4/lt*
-m4/libtool.m4
+m4/libtool.m4
+
+# Buck files #
+##############
+.buck*
+buck-out
+BUCK
+make-bin-BUCK.sh
+
+# generated sources #
+#####################
+lib/qcd/spin/gamma-gen/*.h
+lib/qcd/spin/gamma-gen/*.cc

benchmarks/Benchmark_dwf.cc

@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 typedef WilsonFermion5D<DomainWallVec5dImplR> WilsonFermion5DR;
@@ -351,7 +351,7 @@ int main (int argc, char ** argv)
     ref = zero;
     for(int mu=0;mu<Nd;mu++){
 
-      //    ref =  src - Gamma(Gamma::GammaX)* src ; // 1+gamma_x
+      //    ref =  src - Gamma(Gamma::Algebra::GammaX)* src ; // 1+gamma_x
       tmp = U[mu]*Cshift(src,mu+1,1);
       for(int i=0;i<ref._odata.size();i++){
 	ref._odata[i]+= tmp._odata[i] + Gamma(Gmu[mu])*tmp._odata[i]; ;

benchmarks/Benchmark_dwf_sweep.cc

@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 void benchDw(std::vector<int> & L, int Ls, int threads, int report =0 );

benchmarks/Benchmark_wilson.cc

@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 bool overlapComms = false;
@@ -106,7 +106,7 @@ int main (int argc, char ** argv)
   { // Naive wilson implementation
     ref = zero;
     for(int mu=0;mu<Nd;mu++){
-      //    ref =  src + Gamma(Gamma::GammaX)* src ; // 1-gamma_x
+      //    ref =  src + Gamma(Gamma::Algebra::GammaX)* src ; // 1-gamma_x
       tmp = U[mu]*Cshift(src,mu,1);
       for(int i=0;i<ref._odata.size();i++){
 	ref._odata[i]+= tmp._odata[i] - Gamma(Gmu[mu])*tmp._odata[i]; ;
@@ -159,7 +159,7 @@ int main (int argc, char ** argv)
     ref = zero;
     for(int mu=0;mu<Nd;mu++){
 
-      //    ref =  src - Gamma(Gamma::GammaX)* src ; // 1+gamma_x
+      //    ref =  src - Gamma(Gamma::Algebra::GammaX)* src ; // 1+gamma_x
       tmp = U[mu]*Cshift(src,mu,1);
       for(int i=0;i<ref._odata.size();i++){
 	ref._odata[i]+= tmp._odata[i] + Gamma(Gmu[mu])*tmp._odata[i]; ;

benchmarks/Benchmark_wilson_sweep.cc

@@ -30,11 +30,11 @@ struct scal {
   d internal;
 };
 
-Gamma::GammaMatrix Gmu [] = {
+Gamma::Algebra Gmu [] = {
-  Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-  Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-  Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-  Gamma::GammaT
+    Gamma::Algebra::GammaT
 };
 
 bool overlapComms = false;

configure.ac

@@ -99,6 +99,13 @@ case ${ac_MKL} in
         AC_DEFINE([USE_MKL], [1], [Define to 1 if you use the Intel MKL]);;
 esac
 
+############### HDF5
+AC_ARG_WITH([hdf5],
+    [AS_HELP_STRING([--with-hdf5=prefix],
+    [try this for a non-standard install prefix of the HDF5 library])],
+    [AM_CXXFLAGS="-I$with_hdf5/include $AM_CXXFLAGS"]
+    [AM_LDFLAGS="-L$with_hdf5/lib $AM_LDFLAGS"])
+
 ############### first-touch
 AC_ARG_ENABLE([numa],
     [AC_HELP_STRING([--enable-numa=yes|no|prefix], [enable first touch numa opt])], 
@@ -145,6 +152,12 @@ AC_SEARCH_LIBS([fftw_execute], [fftw3],
                [AC_DEFINE([HAVE_FFTW], [1], [Define to 1 if you have the `FFTW' library])]
                [have_fftw=true])
 
+AC_SEARCH_LIBS([H5Fopen], [hdf5_cpp],
+               [AC_DEFINE([HAVE_HDF5], [1], [Define to 1 if you have the `HDF5' library])]
+               [have_hdf5=true]
+               [LIBS="${LIBS} -lhdf5"], [], [-lhdf5])
+AM_CONDITIONAL(BUILD_HDF5, [ test "${have_hdf5}X" == "trueX" ])
+
 CXXFLAGS=$CXXFLAGS_CPY
 LDFLAGS=$LDFLAGS_CPY
 
@@ -306,7 +319,7 @@ AM_CONDITIONAL(BUILD_COMMS_MPI3L, [ test "${comms_type}X" == "mpi3lX" ] )
 AM_CONDITIONAL(BUILD_COMMS_NONE,  [ test "${comms_type}X" == "noneX" ])
 
 ############### RNG selection
-AC_ARG_ENABLE([rng],[AC_HELP_STRING([--enable-rng=ranlux48|mt19937],\
+AC_ARG_ENABLE([rng],[AC_HELP_STRING([--enable-rng=ranlux48|mt19937|sitmo],\
 	            [Select Random Number Generator to be used])],\
 	            [ac_RNG=${enable_rng}],[ac_RNG=ranlux48])
 
@@ -317,6 +330,9 @@ case ${ac_RNG} in
      mt19937)
       AC_DEFINE([RNG_MT19937],[1],[RNG_MT19937] )
      ;;
+     sitmo)
+      AC_DEFINE([RNG_SITMO],[1],[RNG_SITMO] )
+     ;;
      *)
       AC_MSG_ERROR([${ac_RNG} unsupported --enable-rng option]); 
      ;;
@@ -410,6 +426,7 @@ RNG choice                  : ${ac_RNG}
 GMP                         : `if test "x$have_gmp" = xtrue; then echo yes; else echo no; fi`
 LAPACK                      : ${ac_LAPACK}
 FFTW                        : `if test "x$have_fftw" = xtrue; then echo yes; else echo no; fi`
+HDF5                        : `if test "x$have_hdf5" = xtrue; then echo yes; else echo no; fi`
 build DOXYGEN documentation : `if test "$DX_FLAG_doc" = '1'; then echo yes; else echo no; fi`
 ----- BUILD FLAGS -------------------------------------
 CXXFLAGS:

extras/Hadrons/Application.cc

@@ -42,7 +42,6 @@ using namespace Hadrons;
  ******************************************************************************/
 // constructors ////////////////////////////////////////////////////////////////
 Application::Application(void)
-: env_(Environment::getInstance())
 {
     LOG(Message) << "Modules available:" << std::endl;
     auto list = ModuleFactory::getInstance().getBuilderList();
@@ -74,11 +73,17 @@ Application::Application(const std::string parameterFileName)
     parameterFileName_ = parameterFileName;
 }
 
+// environment shortcut ////////////////////////////////////////////////////////
+Environment & Application::env(void) const
+{
+    return Environment::getInstance();
+}
+
 // access //////////////////////////////////////////////////////////////////////
 void Application::setPar(const Application::GlobalPar &par)
 {
     par_ = par;
-    env_.setSeed(strToVec<int>(par_.seed));
+    env().setSeed(strToVec<int>(par_.seed));
 }
 
 const Application::GlobalPar & Application::getPar(void)
@@ -89,7 +94,7 @@ const Application::GlobalPar & Application::getPar(void)
 // execute /////////////////////////////////////////////////////////////////////
 void Application::run(void)
 {
-    if (!parameterFileName_.empty() and (env_.getNModule() == 0))
+    if (!parameterFileName_.empty() and (env().getNModule() == 0))
     {
         parseParameterFile(parameterFileName_);
     }
@@ -124,7 +129,7 @@ void Application::parseParameterFile(const std::string parameterFileName)
     do
     {
         read(reader, "id", id);
-        env_.createModule(id.name, id.type, reader);
+        env().createModule(id.name, id.type, reader);
     } while (reader.nextElement("module"));
     pop(reader);
     pop(reader);
@@ -134,7 +139,7 @@ void Application::saveParameterFile(const std::string parameterFileName)
 {
     XmlWriter          writer(parameterFileName);
     ObjectId           id;
-    const unsigned int nMod = env_.getNModule();
+    const unsigned int nMod = env().getNModule();
     
     LOG(Message) << "Saving application to '" << parameterFileName << "'..." << std::endl;
     write(writer, "parameters", getPar());
@@ -142,10 +147,10 @@ void Application::saveParameterFile(const std::string parameterFileName)
     for (unsigned int i = 0; i < nMod; ++i)
     {
         push(writer, "module");
-        id.name = env_.getModuleName(i);
+        id.name = env().getModuleName(i);
-        id.type = env_.getModule(i)->getRegisteredName();
+        id.type = env().getModule(i)->getRegisteredName();
         write(writer, "id", id);
-        env_.getModule(i)->saveParameters(writer, "options");
+        env().getModule(i)->saveParameters(writer, "options");
         pop(writer);
     }
     pop(writer);
@@ -164,10 +169,10 @@ auto memPeak = [this](const std::vector<unsigned int> &program)\
     \
     msg = HadronsLogMessage.isActive();\
     HadronsLogMessage.Active(false);\
-    env_.dryRun(true);\
+    env().dryRun(true);\
-    memPeak = env_.executeProgram(program);\
+    memPeak = env().executeProgram(program);\
-    env_.dryRun(false);\
+    env().dryRun(false);\
-    env_.freeAll();\
+    env().freeAll();\
     HadronsLogMessage.Active(true);\
     \
     return memPeak;\
@@ -179,7 +184,7 @@ void Application::schedule(void)
     
     // build module dependency graph
     LOG(Message) << "Building module graph..." << std::endl;
-    auto graph = env_.makeModuleGraph();
+    auto graph = env().makeModuleGraph();
     auto con = graph.getConnectedComponents();
     
     // constrained topological sort using a genetic algorithm
@@ -256,7 +261,7 @@ void Application::saveSchedule(const std::string filename)
                  << std::endl;
     for (auto address: program_)
     {
-        program.push_back(env_.getModuleName(address));
+        program.push_back(env().getModuleName(address));
     }
     write(writer, "schedule", program);
 }
@@ -274,7 +279,7 @@ void Application::loadSchedule(const std::string filename)
     program_.clear();
     for (auto &name: program)
     {
-        program_.push_back(env_.getModuleAddress(name));
+        program_.push_back(env().getModuleAddress(name));
     }
     scheduled_ = true;
     memPeak_   = memPeak(program_);
@@ -291,7 +296,7 @@ void Application::printSchedule(void)
     for (unsigned int i = 0; i < program_.size(); ++i)
     {
         LOG(Message) << std::setw(4) << i + 1 << ": "
-                     << env_.getModuleName(program_[i]) << std::endl;
+                     << env().getModuleName(program_[i]) << std::endl;
     }
 }
 
@@ -304,9 +309,9 @@ void Application::configLoop(void)
     {
         LOG(Message) << BIG_SEP << " Starting measurement for trajectory " << t
                      << " " << BIG_SEP << std::endl;
-        env_.setTrajectory(t);
+        env().setTrajectory(t);
-        env_.executeProgram(program_);
+        env().executeProgram(program_);
-        env_.freeAll();
     }
     LOG(Message) << BIG_SEP << " End of measurement " << BIG_SEP << std::endl;
+    env().freeAll();
 }

extras/Hadrons/Application.hpp

@@ -98,11 +98,13 @@ public:
     void printSchedule(void);
     // loop on configurations
     void configLoop(void);
+private:
+    // environment shortcut
+    Environment & env(void) const;
 private:
     long unsigned int         locVol_;
     std::string               parameterFileName_{""};
     GlobalPar                 par_;
-    Environment               &env_;
     std::vector<unsigned int> program_;
     Environment::Size         memPeak_;
     bool                      scheduled_{false};
@@ -115,14 +117,14 @@ private:
 template <typename M>
 void Application::createModule(const std::string name)
 {
-    env_.createModule<M>(name);
+    env().createModule<M>(name);
 }
 
 template <typename M>
 void Application::createModule(const std::string name,
                                const typename M::Par &par)
 {
-    env_.createModule<M>(name, par);
+    env().createModule<M>(name, par);
 }
 
 END_HADRONS_NAMESPACE

extras/Hadrons/Environment.cc

@@ -41,8 +41,9 @@ using namespace Hadrons;
 // constructor /////////////////////////////////////////////////////////////////
 Environment::Environment(void)
 {
+    nd_ = GridDefaultLatt().size();
     grid4d_.reset(SpaceTimeGrid::makeFourDimGrid(
-        GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
+        GridDefaultLatt(), GridDefaultSimd(nd_, vComplex::Nsimd()),
         GridDefaultMpi()));
     gridRb4d_.reset(SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_.get()));
     auto loc = getGrid()->LocalDimensions();
@@ -126,6 +127,11 @@ GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
     }
 }
 
+unsigned int Environment::getNd(void) const
+{
+    return nd_;
+}
+
 // random number generator /////////////////////////////////////////////////////
 void Environment::setSeed(const std::vector<int> &seed)
 {
@@ -410,12 +416,19 @@ Environment::Size Environment::executeProgram(const std::vector<std::string> &p)
 // general memory management ///////////////////////////////////////////////////
 void Environment::addObject(const std::string name, const int moduleAddress)
 {
-    ObjInfo info;
+    if (!hasObject(name))
-    
+    {
-    info.name   = name;
+        ObjInfo info;
-    info.module = moduleAddress;
+        
-    object_.push_back(std::move(info));
+        info.name   = name;
-    objectAddress_[name] = static_cast<unsigned int>(object_.size() - 1);
+        info.module = moduleAddress;
+        object_.push_back(std::move(info));
+        objectAddress_[name] = static_cast<unsigned int>(object_.size() - 1);
+    }
+    else
+    {
+        HADRON_ERROR("object '" + name + "' already exists");
+    }
 }
 
 void Environment::registerObject(const unsigned int address,
@@ -444,6 +457,10 @@ void Environment::registerObject(const unsigned int address,
 void Environment::registerObject(const std::string name,
                                  const unsigned int size, const unsigned int Ls)
 {
+    if (!hasObject(name))
+    {
+        addObject(name);
+    }
     registerObject(getObjectAddress(name), size, Ls);
 }
 
@@ -573,6 +590,30 @@ bool Environment::hasRegisteredObject(const std::string name) const
     }
 }
 
+bool Environment::hasCreatedObject(const unsigned int address) const
+{
+    if (hasObject(address))
+    {
+        return (object_[address].data != nullptr);
+    }
+    else
+    {
+        return false;
+    }
+}
+
+bool Environment::hasCreatedObject(const std::string name) const
+{
+    if (hasObject(name))
+    {
+        return hasCreatedObject(getObjectAddress(name));
+    }
+    else
+    {
+        return false;
+    }
+}
+
 bool Environment::isObject5d(const unsigned int address) const
 {
     return (getObjectLs(address) > 1);

extras/Hadrons/Environment.hpp

@@ -106,6 +106,7 @@ public:
     void                    createGrid(const unsigned int Ls);
     GridCartesian *         getGrid(const unsigned int Ls = 1) const;
     GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const;
+    unsigned int            getNd(void) const;
     // random number generator
     void                    setSeed(const std::vector<int> &seed);
     GridParallelRNG *       get4dRng(void) const;
@@ -137,7 +138,7 @@ public:
     Size                    executeProgram(const std::vector<std::string> &p);
     // general memory management
     void                    addObject(const std::string name,
-                                      const int moduleAddress);
+                                      const int moduleAddress = -1);
     void                    registerObject(const unsigned int address,
                                            const unsigned int size,
                                            const unsigned int Ls = 1);
@@ -176,6 +177,8 @@ public:
     bool                    hasObject(const std::string name) const;
     bool                    hasRegisteredObject(const unsigned int address) const;
     bool                    hasRegisteredObject(const std::string name) const;
+    bool                    hasCreatedObject(const unsigned int address) const;
+    bool                    hasCreatedObject(const std::string name) const;
     bool                    isObject5d(const unsigned int address) const;
     bool                    isObject5d(const std::string name) const;
     Environment::Size       getTotalSize(void) const;
@@ -198,6 +201,7 @@ private:
     std::map<unsigned int, GridPt>         grid5d_;
     GridRbPt                               gridRb4d_;
     std::map<unsigned int, GridRbPt>       gridRb5d_;
+    unsigned int                           nd_;
     // random number generator
     RngPt                                  rng4d_;
     // module and related maps

extras/Hadrons/GeneticScheduler.hpp

@@ -134,7 +134,7 @@ void GeneticScheduler<T>::nextGeneration(void)
     LOG(Debug) << "Starting population:\n" << *this << std::endl;
     
     // random mutations
-    PARALLEL_FOR_LOOP
+    //PARALLEL_FOR_LOOP
     for (unsigned int i = 0; i < par_.popSize; ++i)
     {
         doMutation();
@@ -142,7 +142,7 @@ void GeneticScheduler<T>::nextGeneration(void)
     LOG(Debug) << "After mutations:\n" << *this << std::endl;
     
     // mating
-    PARALLEL_FOR_LOOP
+    //PARALLEL_FOR_LOOP
     for (unsigned int i = 0; i < par_.popSize/2; ++i)
     {
         doCrossover();
@@ -166,7 +166,7 @@ void GeneticScheduler<T>::initPopulation(void)
     {
         auto p = graph_.topoSort(gen_);
         
-        population_.emplace(func_(p), p);
+        population_.insert(std::make_pair(func_(p), p));
     }
 }
 
@@ -180,8 +180,8 @@ void GeneticScheduler<T>::doCrossover(void)
     crossover(c1, c2, p1, p2);
     PARALLEL_CRITICAL
     {
-        population_.emplace(func_(c1), c1);
+        population_.insert(std::make_pair(func_(c1), c1));
-        population_.emplace(func_(c2), c2);
+        population_.insert(std::make_pair(func_(c2), c2));
     }
 }
 
@@ -200,7 +200,7 @@ void GeneticScheduler<T>::doMutation(void)
         mutation(m, it->second);
         PARALLEL_CRITICAL
         {
-            population_.emplace(func_(m), m);
+            population_.insert(std::make_pair(func_(m), m));
         }
     }
 }

extras/Hadrons/Modules/MContraction/Meson.hpp

@@ -45,9 +45,11 @@ class MesonPar: Serializable
 {
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(MesonPar,
-                                    std::string, q1,
+                                    std::string,    q1,
-                                    std::string, q2,
+                                    std::string,    q2,
-                                    std::string, output);
+                                    std::string,    output,
+                                    Gamma::Algebra, gammaSource,
+                                    Gamma::Algebra, gammaSink);
 };
 
 template <typename FImpl1, typename FImpl2>
@@ -59,8 +61,7 @@ public:
     class Result: Serializable
     {
     public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result, std::vector<Complex>, corr);
-                                        std::vector<std::vector<std::vector<Complex>>>, corr);
     };
 public:
     // constructor
@@ -114,29 +115,17 @@ void TMeson<FImpl1, FImpl2>::execute(void)
     PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
     PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
     LatticeComplex        c(env().getGrid());
-    SpinMatrix            g[Ns*Ns], g5;
+    Gamma                 gSrc(par().gammaSource), gSnk(par().gammaSink);
+    Gamma                 g5(Gamma::Algebra::Gamma5);
     std::vector<TComplex> buf;
     Result                result;
     
-    g5 = makeGammaProd(Ns*Ns - 1);
+    c = trace(gSnk*q1*adj(gSrc)*g5*adj(q2)*g5);
-    result.corr.resize(Ns*Ns);
+    sliceSum(c, buf, Tp);
-    for (unsigned int i = 0; i < Ns*Ns; ++i)
+    result.corr.resize(buf.size());
+    for (unsigned int t = 0; t < buf.size(); ++t)
     {
-        g[i] = makeGammaProd(i);
+        result.corr[t] = TensorRemove(buf[t]);
-    }
-    for (unsigned int iSink = 0; iSink < Ns*Ns; ++iSink)
-    {
-        result.corr[iSink].resize(Ns*Ns);
-        for (unsigned int iSrc = 0; iSrc < Ns*Ns; ++iSrc)
-        {
-            c = trace(g[iSink]*q1*g[iSrc]*g5*adj(q2)*g5);
-            sliceSum(c, buf, Tp);
-            result.corr[iSink][iSrc].resize(buf.size());
-            for (unsigned int t = 0; t < buf.size(); ++t)
-            {
-                result.corr[iSink][iSrc][t] = TensorRemove(buf[t]);
-            }
-        }
     }
     write(writer, "meson", result);
 }

extras/Hadrons/Modules/MSource/SeqGamma.hpp

@@ -60,11 +60,11 @@ class SeqGammaPar: Serializable
 {
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(SeqGammaPar,
-                                    std::string,  q,
+                                    std::string,    q,
-                                    unsigned int, tA,
+                                    unsigned int,   tA,
-                                    unsigned int, tB,
+                                    unsigned int,   tB,
-                                    unsigned int, gamma,
+                                    Gamma::Algebra, gamma,
-                                    std::string,  mom);
+                                    std::string,    mom);
 };
 
 template <typename FImpl>
@@ -140,21 +140,20 @@ void TSeqGamma<FImpl>::execute(void)
     PropagatorField &q   = *env().template getObject<PropagatorField>(par().q);
     Lattice<iScalar<vInteger>> t(env().getGrid());
     LatticeComplex             ph(env().getGrid()), coor(env().getGrid());
-    SpinMatrix                 g;
+    Gamma                      g(par().gamma);
     std::vector<Real>          p;
     Complex                    i(0.0,1.0);
     
-    g  = makeGammaProd(par().gamma);
     p  = strToVec<Real>(par().mom);
     ph = zero;
-    for(unsigned int mu = 0; mu < Nd; mu++)
+    for(unsigned int mu = 0; mu < env().getNd(); mu++)
     {
         LatticeCoordinate(coor, mu);
         ph = ph + p[mu]*coor;
     }
     ph = exp(i*ph);
     LatticeCoordinate(t, Tp);
-    src = where((t >= par().tA) and (t <= par().tB), g*ph*q, 0.*q);
+    src = where((t >= par().tA) and (t <= par().tB), ph*(g*q), 0.*q);
 }
 
 END_MODULE_NAMESPACE

lib/Grid.h

@@ -40,7 +40,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 #include <Grid/GridCore.h>
 #include <Grid/GridQCDcore.h>
-#include <Grid/qcd/action/Actions.h>
+#include <Grid/qcd/action/Action.h>
 #include <Grid/qcd/smearing/Smearing.h>
 #include <Grid/qcd/hmc/HMC_aggregate.h>
 

lib/GridCore.h

@@ -61,12 +61,12 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 ///////////////////
 #include "Config.h"
 
-#include <Grid/serialisation/Serialisation.h>
 #include <Grid/perfmon/Timer.h>
 #include <Grid/perfmon/PerfCount.h>
 #include <Grid/log/Log.h>
 #include <Grid/allocator/AlignedAllocator.h>
 #include <Grid/simd/Simd.h>
+#include <Grid/serialisation/Serialisation.h>
 #include <Grid/threads/Threads.h>
 #include <Grid/util/Util.h>
 #include <Grid/communicator/Communicator.h> 

lib/GridQCDcore.h

@@ -27,7 +27,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-
 #ifndef GRID_QCD_CORE_H
 #define GRID_QCD_CORE_H
 

lib/Makefile.am

@@ -1,4 +1,5 @@
 extra_sources=
+extra_headers=
 if BUILD_COMMS_MPI
   extra_sources+=communicator/Communicator_mpi.cc
   extra_sources+=communicator/Communicator_base.cc
@@ -24,6 +25,12 @@ if BUILD_COMMS_NONE
   extra_sources+=communicator/Communicator_base.cc
 endif
 
+if BUILD_HDF5
+  extra_sources+=serialisation/Hdf5IO.cc 
+  extra_headers+=serialisation/Hdf5IO.h
+  extra_headers+=serialisation/Hdf5Type.h
+endif
+
 #
 # Libraries
 #
@@ -32,6 +39,9 @@ include Eigen.inc
 
 lib_LIBRARIES = libGrid.a
 
-libGrid_a_SOURCES              = $(CCFILES) $(extra_sources)
+CCFILES += $(extra_sources)
+HFILES  += $(extra_headers)
+
+libGrid_a_SOURCES              = $(CCFILES)
 libGrid_adir                   = $(pkgincludedir)
 nobase_dist_pkginclude_HEADERS = $(HFILES) $(eigen_files) Config.h

lib/algorithms/approx/Remez.h

@@ -16,7 +16,7 @@
 #define INCLUDED_ALG_REMEZ_H
 
 #include <stddef.h>
-#include <Grid/Config.h>
+#include <Config.h>
 
 #ifdef HAVE_LIBGMP
 #include "bigfloat.h"

lib/algorithms/iterative/ConjugateGradient.h

@@ -45,6 +45,8 @@ class ConjugateGradient : public OperatorFunction<Field> {
                            // Defaults true.
   RealD Tolerance;
   Integer MaxIterations;
+  Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
+  
   ConjugateGradient(RealD tol, Integer maxit, bool err_on_no_conv = true)
       : Tolerance(tol),
         MaxIterations(maxit),
@@ -155,13 +157,14 @@ class ConjugateGradient : public OperatorFunction<Field> {
         std::cout << std::endl;
 
         if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
-
+	IterationsToComplete = k;	
         return;
       }
     }
     std::cout << GridLogMessage << "ConjugateGradient did NOT converge"
               << std::endl;
     if (ErrorOnNoConverge) assert(0);
+    IterationsToComplete = k;
   }
 };
 }

lib/algorithms/iterative/ConjugateGradientMixedPrec.h

@@ -35,6 +35,7 @@ namespace Grid {
   class MixedPrecisionConjugateGradient : public LinearFunction<FieldD> {
   public:                                                
     RealD   Tolerance;
+    RealD   InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed
     Integer MaxInnerIterations;
     Integer MaxOuterIterations;
     GridBase* SinglePrecGrid; //Grid for single-precision fields
@@ -42,12 +43,16 @@ namespace Grid {
     LinearOperatorBase<FieldF> &Linop_f;
     LinearOperatorBase<FieldD> &Linop_d;
 
+    Integer TotalInnerIterations; //Number of inner CG iterations
+    Integer TotalOuterIterations; //Number of restarts
+    Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step
+
     //Option to speed up *inner single precision* solves using a LinearFunction that produces a guess
     LinearFunction<FieldF> *guesser;
     
     MixedPrecisionConjugateGradient(RealD tol, Integer maxinnerit, Integer maxouterit, GridBase* _sp_grid, LinearOperatorBase<FieldF> &_Linop_f, LinearOperatorBase<FieldD> &_Linop_d) :
       Linop_f(_Linop_f), Linop_d(_Linop_d),
-      Tolerance(tol), MaxInnerIterations(maxinnerit), MaxOuterIterations(maxouterit), SinglePrecGrid(_sp_grid),
+      Tolerance(tol), InnerTolerance(tol), MaxInnerIterations(maxinnerit), MaxOuterIterations(maxouterit), SinglePrecGrid(_sp_grid),
       OuterLoopNormMult(100.), guesser(NULL){ };
 
     void useGuesser(LinearFunction<FieldF> &g){
@@ -55,6 +60,8 @@ namespace Grid {
     }
   
     void operator() (const FieldD &src_d_in, FieldD &sol_d){
+      TotalInnerIterations = 0;
+	
       GridStopWatch TotalTimer;
       TotalTimer.Start();
     
@@ -74,7 +81,7 @@ namespace Grid {
       FieldD src_d(DoublePrecGrid);
       src_d = src_d_in; //source for next inner iteration, computed from residual during operation
     
-      RealD inner_tol = Tolerance;
+      RealD inner_tol = InnerTolerance;
     
       FieldF src_f(SinglePrecGrid);
       src_f.checkerboard = cb;
@@ -89,7 +96,9 @@ namespace Grid {
 
       GridStopWatch PrecChangeTimer;
     
-      for(Integer outer_iter = 0; outer_iter < MaxOuterIterations; outer_iter++){
+      Integer &outer_iter = TotalOuterIterations; //so it will be equal to the final iteration count
+      
+      for(outer_iter = 0; outer_iter < MaxOuterIterations; outer_iter++){
 	//Compute double precision rsd and also new RHS vector.
 	Linop_d.HermOp(sol_d, tmp_d);
 	RealD norm = axpy_norm(src_d, -1., tmp_d, src_d_in); //src_d is residual vector
@@ -117,6 +126,7 @@ namespace Grid {
 	InnerCGtimer.Start();
 	CG_f(Linop_f, src_f, sol_f);
 	InnerCGtimer.Stop();
+	TotalInnerIterations += CG_f.IterationsToComplete;
       
 	//Convert sol back to double and add to double prec solution
 	PrecChangeTimer.Start();
@@ -131,9 +141,11 @@ namespace Grid {
     
       ConjugateGradient<FieldD> CG_d(Tolerance, MaxInnerIterations);
       CG_d(Linop_d, src_d_in, sol_d);
+      TotalFinalStepIterations = CG_d.IterationsToComplete;
 
       TotalTimer.Stop();
-      std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Total " << TotalTimer.Elapsed() << " Precision change " << PrecChangeTimer.Elapsed() << " Inner CG total " << InnerCGtimer.Elapsed() << std::endl;
+      std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Inner CG iterations " << TotalInnerIterations << " Restarts " << TotalOuterIterations << " Final CG iterations " << TotalFinalStepIterations << std::endl;
+      std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Total time " << TotalTimer.Elapsed() << " Precision change " << PrecChangeTimer.Elapsed() << " Inner CG total " << InnerCGtimer.Elapsed() << std::endl;
     }
   };
 

lib/algorithms/iterative/Matrix.h

@@ -36,7 +36,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <iomanip>
 #include <complex>
 #include <typeinfo>
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 
 /** Sign function **/

lib/lattice/Lattice_rng.h

@@ -30,6 +30,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #define GRID_LATTICE_RNG_H
 
 #include <random>
+#include <Grid/sitmo_rng/sitmo_prng_engine.hpp>
 
 namespace Grid {
 
@@ -127,10 +128,14 @@ namespace Grid {
     typedef uint64_t      RngStateType;
     typedef std::ranlux48 RngEngine;
     static const int RngStateCount = 15;
-#else
+#elif RNG_MT19937 
     typedef std::mt19937 RngEngine;
     typedef uint32_t     RngStateType;
     static const int     RngStateCount = std::mt19937::state_size;
+#elif RNG_SITMO
+    typedef sitmo::prng_engine 	RngEngine;
+    typedef uint64_t    	RngStateType;
+    static const int    	RngStateCount = 4;
 #endif
     std::vector<RngEngine>                             _generators;
     std::vector<std::uniform_real_distribution<RealD>> _uniform;

lib/lattice/Lattice_transfer.h

@@ -429,7 +429,7 @@ void ExtractSlice(Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slice, in
     std::vector<int> lcoor(nl);
     std::vector<int> hcoor(nh);
     lg->LocalIndexToLocalCoor(idx,lcoor);
-    ddl=0;
+    int ddl=0;
     hcoor[orthog] = slice;
     for(int d=0;d<nh;d++){
       if ( d!=orthog ) { 

lib/perfmon/Stat.cc

@@ -2,10 +2,8 @@
 #include <Grid/perfmon/PerfCount.h>
 #include <Grid/perfmon/Stat.h>
 
-
 namespace Grid { 
 
-
 bool PmuStat::pmu_initialized=false;
 
 

lib/pugixml/pugixml.cc

@@ -14,7 +14,7 @@
 #ifndef SOURCE_PUGIXML_CPP
 #define SOURCE_PUGIXML_CPP
 
-#include <pugixml/pugixml.h>
+#include <Grid/pugixml/pugixml.h>
 
 #include <stdlib.h>
 #include <stdio.h>

lib/qcd/action/Action.h

@@ -30,23 +30,18 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#ifndef GRID_QCD_ACTIONS_H
+#ifndef GRID_QCD_ACTION_H
-#define GRID_QCD_ACTIONS_H
+#define GRID_QCD_ACTION_H
-
-// * Linear operators             (Hermitian and non-hermitian)  .. my LinearOperator
-// * System solvers               (Hermitian and non-hermitian)  .. my OperatorFunction
-// * MultiShift System solvers    (Hermitian and non-hermitian)  .. my OperatorFunction
 
 ////////////////////////////////////////////
 // Abstract base interface
 ////////////////////////////////////////////
 #include <Grid/qcd/action/ActionCore.h>
-
 ////////////////////////////////////////////////////////////////////////
 // Fermion actions; prevent coupling fermion.cc files to other headers
 ////////////////////////////////////////////////////////////////////////
+#include <Grid/qcd/action/fermion/FermionCore.h>
 #include <Grid/qcd/action/fermion/Fermion.h>
-
 ////////////////////////////////////////
 // Pseudo fermion combinations for HMC
 ////////////////////////////////////////

lib/qcd/action/fermion/Fermion.h

@@ -25,8 +25,8 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#ifndef  GRID_QCD_FERMION_ACTIONS_H
+#ifndef  GRID_QCD_FERMION_H
-#define  GRID_QCD_FERMION_ACTIONS_H
+#define  GRID_QCD_FERMION_H
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Explicit explicit template instantiation is still required in the .cc files

lib/qcd/action/fermion/PartialFractionFermion5D.cc

@@ -28,6 +28,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     /*  END LEGAL */
 #include <Grid/qcd/action/fermion/FermionCore.h>
 #include <Grid/qcd/action/fermion/PartialFractionFermion5D.h>
+
 namespace Grid {
   namespace QCD {
 

lib/qcd/action/fermion/WilsonFermion.cc

@@ -134,7 +134,6 @@ void WilsonFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out)
   out.checkerboard = in.checkerboard;
   MooeeInv(in,out);
 }
-
 template<class Impl>
 void WilsonFermion<Impl>::MomentumSpacePropagator(FermionField &out, const FermionField &in,RealD _m) 
 {  
@@ -145,11 +144,11 @@ void WilsonFermion<Impl>::MomentumSpacePropagator(FermionField &out, const Fermi
   // what type LatticeComplex 
   conformable(_grid,out._grid);
   
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
   
   std::vector<int> latt_size   = _grid->_fdimensions;

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -563,11 +563,11 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt(FermionField &out,const Fe
   typedef iSinglet<ScalComplex> Tcomplex;
   typedef Lattice<iSinglet<vector_type> > LatComplex;
   
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
   std::vector<int> latt_size   = _grid->_fdimensions;
@@ -634,11 +634,11 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt(FermionField &out,const Fe
 template<class Impl>
 void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass) 
 {
-    Gamma::GammaMatrix Gmu [] = {
+    Gamma::Algebra Gmu [] = {
-      Gamma::GammaX,
+      Gamma::Algebra::GammaX,
-      Gamma::GammaY,
+      Gamma::Algebra::GammaY,
-      Gamma::GammaZ,
+      Gamma::Algebra::GammaZ,
-      Gamma::GammaT
+      Gamma::Algebra::GammaT
     };
 
     GridBase *_grid = _FourDimGrid;

lib/qcd/action/fermion/WilsonKernels.cc

@@ -29,6 +29,7 @@ directory
 *************************************************************************************/
 /*  END LEGAL */
 #include <Grid/qcd/action/fermion/FermionCore.h>
+
 namespace Grid {
 namespace QCD {
 

lib/qcd/action/fermion/WilsonKernelsAsm.cc

@@ -32,7 +32,6 @@ Author: Guido Cossu <guido.cossu@ed.ac.uk>
 
 #include <Grid/qcd/action/fermion/FermionCore.h>
 
-
 namespace Grid {
 namespace QCD {
 

lib/qcd/action/fermion/g5HermitianLinop.h

@@ -80,7 +80,7 @@ class Gamma5HermitianLinearOperator : public LinearOperatorBase<Field> {
   Matrix &_Mat;
   Gamma g5;
 public:
-    Gamma5HermitianLinearOperator(Matrix &Mat): _Mat(Mat), g5(Gamma::Gamma5) {};
+    Gamma5HermitianLinearOperator(Matrix &Mat): _Mat(Mat), g5(Gamma::Algebra::Gamma5) {};
   void Op     (const Field &in, Field &out){
     HermOp(in,out);
   }

lib/qcd/action/gauge/Gauge.h

@@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-#ifndef GRID_QCD_GAUGE_AGGREGATE_H
+#ifndef GRID_QCD_GAUGE_H
-#define GRID_QCD_GAUGE_AGGREGATE_H
+#define GRID_QCD_GAUGE_H
 
 #include <Grid/qcd/action/gauge/GaugeImpl.h>
 #include <Grid/qcd/utils/WilsonLoops.h>

lib/qcd/spin/Dirac.cc

@@ -1,96 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/qcd/spin/Dirac.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    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 */
-#include <Grid/GridCore.h>
-#include <Grid/GridQCDcore.h>
-
-namespace Grid {
-
-  namespace QCD {
-
-    Gamma::GammaMatrix  Gamma::GammaMatrices [] = {
-      Gamma::Identity,
-      Gamma::GammaX,
-      Gamma::GammaY,
-      Gamma::GammaZ,
-      Gamma::GammaT,
-      Gamma::Gamma5,
-      Gamma::MinusIdentity,
-      Gamma::MinusGammaX,
-      Gamma::MinusGammaY,
-      Gamma::MinusGammaZ,
-      Gamma::MinusGammaT,
-      Gamma::MinusGamma5
-    };
-    const char *Gamma::GammaMatrixNames[] = { 
-      "Identity ",
-      "GammaX   ",
-      "GammaY   ",
-      "GammaZ   ",
-      "GammaT   ",
-      "Gamma5   ",
-      "-Identity",
-      "-GammaX  ",
-      "-GammaY  ",
-      "-GammaZ  ",
-      "-GammaT  ",
-      "-Gamma5  ",
-      "         "
-    };
-    
-    SpinMatrix makeGammaProd(const unsigned int i)
-    {
-      SpinMatrix g;
-      
-      g = 1.;
-      if (i & 0x1)
-      {
-        g = g*Gamma(Gamma::GammaMatrix::GammaX);
-      }
-      if (i & 0x2)
-      {
-        g = g*Gamma(Gamma::GammaMatrix::GammaY);
-      }
-      if (i & 0x4)
-      {
-        g = g*Gamma(Gamma::GammaMatrix::GammaZ);
-      }
-      if (i & 0x8)
-      {
-        g = g*Gamma(Gamma::GammaMatrix::GammaT);
-      }
-      
-      return g;
-    }
-
-    //    void sprojMul( vHalfSpinColourVector &out,vColourMatrix &u, vSpinColourVector &in){
-    //      vHalfSpinColourVector hspin;
-    //      spProjXp(hspin,in);
-    //      mult(&out,&u,&hspin);
-    //    }
-  }
-}

lib/qcd/spin/Dirac.h

@@ -1,628 +1,232 @@
-    /*************************************************************************************
+/*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+Grid physics library, www.github.com/paboyle/Grid 
 
-    Source file: ./lib/qcd/spin/Dirac.h
+Source file: lib/qcd/spin/Dirac.h
 
-    Copyright (C) 2015
+Copyright (C) 2015
+Copyright (C) 2016
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
 
-    This program is free software; you can redistribute it and/or modify
+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
+it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
+the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
+(at your option) any later version.
 
-    This program is distributed in the hope that it will be useful,
+This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
+but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+GNU General Public License for more details.
 
-    You should have received a copy of the GNU General Public License along
+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.,
+with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
+See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
+*************************************************************************************/
-    /*  END LEGAL */
+/*  END LEGAL */
 #ifndef GRID_QCD_DIRAC_H
 #define GRID_QCD_DIRAC_H
-namespace Grid{
+
+// Gamma matrices using the code generated by the Mathematica notebook 
+// gamma-gen/gamma-gen.nb in Gamma.cc & Gamma.h
+////////////////////////////////////////////////////////////////////////////////
+#include <Grid/qcd/spin/Gamma.h>
+
+namespace Grid {
+
+// Dirac algebra adjoint operator (not in QCD:: to overload other adj)
+inline QCD::Gamma adj(const QCD::Gamma &g)
+{
+   return QCD::Gamma (QCD::Gamma::adj[g.g]);
+}
 
 namespace QCD {
 
+// Dirac algebra mutliplication operator
+inline Gamma operator*(const Gamma &g1, const Gamma &g2)
+{
+   return Gamma (Gamma::mul[g1.g][g2.g]);
+}
 
-  class Gamma {
+// general left multiply
+template<class vtype> 
+inline auto operator*(const Gamma &G, const iScalar<vtype> &arg)
+->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
+{
+  iScalar<vtype> ret;
+  ret._internal=G*arg._internal;
+  return ret;
+}
 
-  public:
+template<class vtype,int N>
+inline auto operator*(const Gamma &G, const iVector<vtype, N> &arg)
+->typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N>>::type 
+{
+  iVector<vtype,N> ret;
+  for(int i=0;i<N;i++){
+    ret._internal[i]=G*arg._internal[i];
+  }
+  return ret;
+}
 
-    const int Ns=4;
+template<class vtype, int N>
-    
+inline auto operator*(const Gamma &G, const iMatrix<vtype, N> &arg) 
-    enum GammaMatrix {
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type 
-      Identity,                        
+{
-      GammaX, 
+  iMatrix<vtype,N> ret;
-      GammaY, 
+  for(int i=0;i<N;i++){
-      GammaZ, 
+  for(int j=0;j<N;j++){
-      GammaT,  
+    ret._internal[i][j]=G*arg._internal[i][j];
-      Gamma5,
+  }}
-      MinusIdentity,                        
+  return ret;
-      MinusGammaX, 
+}
-      MinusGammaY, 
-      MinusGammaZ, 
-      MinusGammaT,  
-      MinusGamma5
-      //      GammaXGamma5, // Rest are composite (willing to take hit for two calls sequentially)
-      //      GammaYGamma5, // as they are less commonly used.
-      //      GammaZGamma5, 
-      //      GammaTGamma5,  
-      //      SigmaXY,
-      //      SigmaXZ,
-      //      SigmaYZ,
-      //      SigmaXT,
-      //      SigmaYT,
-      //      SigmaZT,
-      //      MinusGammaXGamma5, easiest to form by composition
-      //      MinusGammaYGamma5, as performance is not critical for these
-      //      MinusGammaZGamma5, 
-      //      MinusGammaTGamma5,  
-      //      MinusSigmaXY,
-      //      MinusSigmaXZ,
-      //      MinusSigmaYZ,
-      //      MinusSigmaXT,
-      //      MinusSigmaYT,
-      //      MinusSigmaZT
-    };
-    
-    static GammaMatrix GammaMatrices[];
-    static const char *GammaMatrixNames[];
 
-    Gamma (GammaMatrix g) { _g=g; }
+// general right multiply
+template<class vtype>
+inline auto operator*(const iScalar<vtype> &arg, const Gamma &G)
+->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
+{
+  iScalar<vtype> ret;
+  ret._internal=arg._internal*G;
+  return ret;
+}
 
-    GammaMatrix _g;
+template<class vtype, int N>
+inline auto operator * (const iMatrix<vtype, N> &arg, const Gamma &G) 
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type 
+{
+  iMatrix<vtype,N> ret;
+  for(int i=0;i<N;i++){
+  for(int j=0;j<N;j++){
+    ret._internal[i][j]=arg._internal[i][j]*G;
+  }}
+  return ret;
+}
 
-  };
+// Gamma-left matrices gL_mu = g_mu*(1 - g5)
+////////////////////////////////////////////////////////////////////////////////
+class GammaL
+{
+public:
+  typedef Gamma::Algebra Algebra;
+  Gamma gamma;
+public:
+  GammaL(const Algebra initg): gamma(initg) {}
+  GammaL(const Gamma   initg): gamma(initg) {}
+};
+
+// vector multiply
+template<class vtype> 
+inline auto operator*(const GammaL &gl, const iVector<vtype, Ns> &arg)
+->typename std::enable_if<matchGridTensorIndex<iVector<vtype, Ns>, SpinorIndex>::value, iVector<vtype, Ns>>::type
+{
+  iVector<vtype, Ns> buf;
   
-    // Make gamma products (Chroma convention)
+  buf(0) = 0.;
-    SpinMatrix makeGammaProd(const unsigned int i);
+  buf(1) = 0.;
-    
+  buf(2) = 2.*arg(2);
-    /* Gx
+  buf(3) = 2.*arg(3);
-     *  0 0  0  i    
+  
-     *  0 0  i  0    
+  return gl.gamma*buf;
-     *  0 -i 0  0
+};
-     * -i 0  0  0
-     */
-  // right multiplication makes sense for matrix args, not for vector since there is 
-  // no concept of row versus columnar indices
-    template<class vtype> inline void rmultMinusGammaX(iMatrix<vtype,Ns> &ret,const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) = timesI(rhs(i,3));
-	ret(i,1) = timesI(rhs(i,2));
-	ret(i,2) = timesMinusI(rhs(i,1));
-	ret(i,3) = timesMinusI(rhs(i,0));
-      }
-    };
-    template<class vtype> inline void rmultGammaX(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) = timesMinusI(rhs(i,3));
-	ret(i,1) = timesMinusI(rhs(i,2));
-	ret(i,2) = timesI(rhs(i,1));
-	ret(i,3) = timesI(rhs(i,0));
-      }
-    };
-    template<class vtype> inline void multGammaX(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = timesI(rhs(3,i));
-	ret(1,i) = timesI(rhs(2,i));
-	ret(2,i) = timesMinusI(rhs(1,i));
-	ret(3,i) = timesMinusI(rhs(0,i));
-      }
-    };
-    template<class vtype> inline void multMinusGammaX(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = timesMinusI(rhs(3,i));
-	ret(1,i) = timesMinusI(rhs(2,i));
-	ret(2,i) = timesI(rhs(1,i));
-	ret(3,i) = timesI(rhs(0,i));
-      }
-    };
 
-    template<class vtype> inline void multGammaX(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
+// matrix left multiply
-      ret._internal[0] = timesI(rhs._internal[3]);
+template<class vtype> 
-      ret._internal[1] = timesI(rhs._internal[2]);
+inline auto operator*(const GammaL &gl, const iMatrix<vtype, Ns> &arg)
-      ret._internal[2] = timesMinusI(rhs._internal[1]);
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, iMatrix<vtype, Ns>>::type
-      ret._internal[3] = timesMinusI(rhs._internal[0]);
+{
-    };
+  iMatrix<vtype, Ns> buf;
-    template<class vtype> inline void multMinusGammaX(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
+  
-	ret(0) = timesMinusI(rhs(3));
+  for(unsigned int i = 0; i < Ns; ++i)
-	ret(1) = timesMinusI(rhs(2));
+  {
-	ret(2) = timesI(rhs(1));
+    buf(0, i) = 0.;
-	ret(3) = timesI(rhs(0));
+    buf(1, i) = 0.;
-    };
+    buf(2, i) = 2.*arg(2, i);
+    buf(3, i) = 2.*arg(3, i);
+  }
+  
+  return gl.gamma*buf;
+};
 
+// matrix right multiply
+template<class vtype> 
+inline auto operator*(const iMatrix<vtype, Ns> &arg, const GammaL &gl)
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, iMatrix<vtype, Ns>>::type
+{
+  iMatrix<vtype, Ns> buf;
+  
+  buf = arg*gl.gamma;
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    buf(i, 0) = 0.;
+    buf(i, 1) = 0.;
+    buf(i, 2) = 2.*buf(i, 2);
+    buf(i, 3) = 2.*buf(i, 3);
+  }
+  
+  return buf;
+};
 
-    /*Gy
+//general left multiply
-     *  0 0  0  -1  [0] -+ [3]
+template<class vtype> 
-     *  0 0  1  0   [1] +- [2]
+inline auto operator*(const GammaL &gl, const iScalar<vtype> &arg)
-     *  0 1  0  0
+->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
-     * -1 0  0  0
+{
-     */
+  iScalar<vtype> ret;
-    template<class vtype> inline void rmultGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
+  ret._internal=gl*arg._internal;
-      for(int i=0;i<Ns;i++){
+  return ret;
-	ret(i,0) = -rhs(i,3);
+}
-	ret(i,1) =  rhs(i,2);
-	ret(i,2) =  rhs(i,1);
-	ret(i,3) = -rhs(i,0);
-      }
-    };
-    template<class vtype> inline void rmultMinusGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) =  rhs(i,3);
-	ret(i,1) = -rhs(i,2);
-	ret(i,2) = -rhs(i,1);
-	ret(i,3) =  rhs(i,0);
-      }
-    };
-    template<class vtype> inline void multGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = -rhs(3,i);
-	ret(1,i) =  rhs(2,i);
-	ret(2,i) =  rhs(1,i);
-	ret(3,i) = -rhs(0,i);
-      }
-    };
-    template<class vtype> inline void multMinusGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) =  rhs(3,i);
-	ret(1,i) = -rhs(2,i);
-	ret(2,i) = -rhs(1,i);
-	ret(3,i) =  rhs(0,i);
-      }
-    };
-    template<class vtype> inline void multGammaY(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) = -rhs(3);
-      ret(1) =  rhs(2);
-      ret(2) =  rhs(1);
-      ret(3) = -rhs(0);
-    };
-    template<class vtype> inline void multMinusGammaY(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) =  rhs(3);
-      ret(1) = -rhs(2);
-      ret(2) = -rhs(1);
-      ret(3) =  rhs(0);
-    };
-    /*Gz
-     *  0 0  i  0   [0]+-i[2]
-     *  0 0  0 -i   [1]-+i[3]
-     * -i 0  0  0
-     *  0 i  0  0
-     */
-    template<class vtype> inline void rmultGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) = timesMinusI(rhs(i,2));
-	ret(i,1) =      timesI(rhs(i,3));
-	ret(i,2) =      timesI(rhs(i,0));
-	ret(i,3) = timesMinusI(rhs(i,1));
-      }
-    };
-    template<class vtype> inline void rmultMinusGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) =      timesI(rhs(i,2));
-	ret(i,1) = timesMinusI(rhs(i,3));
-	ret(i,2) = timesMinusI(rhs(i,0));
-	ret(i,3) =      timesI(rhs(i,1));
-      }
-    };
-    template<class vtype> inline void multGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = timesI(rhs(2,i));
-	ret(1,i) =timesMinusI(rhs(3,i));
-	ret(2,i) =timesMinusI(rhs(0,i));
-	ret(3,i) = timesI(rhs(1,i));
-      }
-    };
-    template<class vtype> inline void multMinusGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = timesMinusI(rhs(2,i));
-	ret(1,i) = timesI(rhs(3,i));
-	ret(2,i) = timesI(rhs(0,i));
-	ret(3,i) = timesMinusI(rhs(1,i));
-      }
-    };
-    template<class vtype> inline void multGammaZ(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) = timesI(rhs(2));
-      ret(1) =timesMinusI(rhs(3));
-      ret(2) =timesMinusI(rhs(0));
-      ret(3) = timesI(rhs(1));
-    };
-    template<class vtype> inline void multMinusGammaZ(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) = timesMinusI(rhs(2));
-      ret(1) = timesI(rhs(3));
-      ret(2) = timesI(rhs(0));
-      ret(3) = timesMinusI(rhs(1));
-    };
-    /*Gt
-     *  0 0  1  0 [0]+-[2]
-     *  0 0  0  1 [1]+-[3]
-     *  1 0  0  0
-     *  0 1  0  0
-     */
-    template<class vtype> inline void rmultGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) = rhs(i,2);
-	ret(i,1) = rhs(i,3);
-	ret(i,2) = rhs(i,0);
-	ret(i,3) = rhs(i,1);
-      }
-    };
-    template<class vtype> inline void rmultMinusGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) =- rhs(i,2);
-	ret(i,1) =- rhs(i,3);
-	ret(i,2) =- rhs(i,0);
-	ret(i,3) =- rhs(i,1);
-      }
-    };
-    template<class vtype> inline void multGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = rhs(2,i);
-	ret(1,i) = rhs(3,i);
-	ret(2,i) = rhs(0,i);
-	ret(3,i) = rhs(1,i);
-      }
-    };
-    template<class vtype> inline void multMinusGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) =-rhs(2,i);
-	ret(1,i) =-rhs(3,i);
-	ret(2,i) =-rhs(0,i);
-	ret(3,i) =-rhs(1,i);
-      }
-    };
-    template<class vtype> inline void multGammaT(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) = rhs(2);
-      ret(1) = rhs(3);
-      ret(2) = rhs(0);
-      ret(3) = rhs(1);
-    };
-    template<class vtype> inline void multMinusGammaT(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) =-rhs(2);
-      ret(1) =-rhs(3);
-      ret(2) =-rhs(0);
-      ret(3) =-rhs(1);
-    };
-    /*G5
-     *  1 0  0  0 [0]+-[2]
-     *  0 1  0  0 [1]+-[3]
-     *  0 0 -1  0
-     *  0 0  0 -1
-     */
-    template<class vtype> inline void rmultGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) = rhs(i,0);
-	ret(i,1) = rhs(i,1);
-	ret(i,2) =-rhs(i,2);
-	ret(i,3) =-rhs(i,3);
-      }
-    };
-    template<class vtype> inline void rmultMinusGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) =-rhs(i,0);
-	ret(i,1) =-rhs(i,1);
-	ret(i,2) = rhs(i,2);
-	ret(i,3) = rhs(i,3);
-      }
-    };
 
-    template<class vtype> inline void multGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
+template<class vtype,int N>
-      for(int i=0;i<Ns;i++){
+inline auto operator*(const GammaL &gl, const iVector<vtype, N> &arg)
-	ret(0,i) = rhs(0,i);
+->typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N>>::type 
-	ret(1,i) = rhs(1,i);
+{
-	ret(2,i) =-rhs(2,i);
+  iVector<vtype,N> ret;
-	ret(3,i) =-rhs(3,i);
+  for(int i=0;i<N;i++){
-      }
+    ret._internal[i]=gl*arg._internal[i];
-    };
+  }
-    template<class vtype> inline void multMinusGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
+  return ret;
-      for(int i=0;i<Ns;i++){
+}
-	ret(0,i) =-rhs(0,i);
-	ret(1,i) =-rhs(1,i);
-	ret(2,i) = rhs(2,i);
-	ret(3,i) = rhs(3,i);
-      }
-    };
 
-    template<class vtype> inline void multGamma5(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
+template<class vtype, int N>
-      ret(0) = rhs(0);
+inline auto operator*(const GammaL &gl, const iMatrix<vtype, N> &arg) 
-      ret(1) = rhs(1);
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type 
-      ret(2) =-rhs(2);
+{
-      ret(3) =-rhs(3);
+  iMatrix<vtype,N> ret;
-    };
+  for(int i=0;i<N;i++){
-    template<class vtype> inline void multMinusGamma5(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
+  for(int j=0;j<N;j++){
-      ret(0) =-rhs(0);
+    ret._internal[i][j]=gl*arg._internal[i][j];
-      ret(1) =-rhs(1);
+  }}
-      ret(2) = rhs(2);
+  return ret;
-      ret(3) = rhs(3);
+}
-    };
 
+//general right multiply
+template<class vtype>
+inline auto operator*(const iScalar<vtype> &arg, const GammaL &gl)
+->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
+{
+  iScalar<vtype> ret;
+  ret._internal=arg._internal*gl;
+  return ret;
+}
 
+template<class vtype, int N>
+inline auto operator * (const iMatrix<vtype, N> &arg, const GammaL &gl) 
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type 
+{
+  iMatrix<vtype,N> ret;
+  for(int i=0;i<N;i++){
+  for(int j=0;j<N;j++){
+    ret._internal[i][j]=arg._internal[i][j]*gl;
+  }}
+  return ret;
+}
 
-#ifdef GRID_WARN_SUBOPTIMAL
+}}
-#warning "Optimisation alert switch over to multGammaX early "
-#endif     
 
-    ///////////////////////////////////////////////////////////////////////////////////////////////////
-    // Operator * : first case this is not a spin index, so recurse
-    ///////////////////////////////////////////////////////////////////////////////////////////////////
-    
-    // FIXME
-    //
-    // Optimisation; switch over to a "multGammaX(ret._internal,arg._internal)" style early and
-    // note that doing so from the lattice operator will avoid copy back and case switch overhead, as
-    // was done for the tensor math operator to remove operator * notation early
-    //
-
-    //left multiply
-    template<class vtype> inline auto operator * ( const Gamma &G,const iScalar<vtype> &arg) ->
-      typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype> >::type 
-
-	{
-	  iScalar<vtype> ret;
-	  ret._internal=G*arg._internal;
-	  return ret;
-	}
-    template<class vtype,int N> inline auto operator * ( const Gamma &G,const iVector<vtype,N> &arg) ->
-      typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N> >::type 
-	{
-	  iVector<vtype,N> ret;
-	  for(int i=0;i<N;i++){
-	    ret._internal[i]=G*arg._internal[i];
-	  }
-	  return ret;
-	}
-    template<class vtype,int N> inline auto operator * ( const Gamma &G,const iMatrix<vtype,N> &arg) ->
-      typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N> >::type 
-	{
-	  iMatrix<vtype,N> ret;
-	  for(int i=0;i<N;i++){
-	  for(int j=0;j<N;j++){
-	    ret._internal[i][j]=G*arg._internal[i][j];
-	  }}
-	  return ret;
-	}
-
-
-    //right multiply
-    template<class vtype> inline auto operator * (const iScalar<vtype> &arg, const Gamma &G) ->
-      typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype> >::type 
-
-	{
-	  iScalar<vtype> ret;
-	  ret._internal=arg._internal*G;
-	  return ret;
-	}
-    template<class vtype,int N> inline auto operator * (const iVector<vtype,N> &arg, const Gamma &G) ->
-      typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N> >::type 
-	{
-	  iVector<vtype,N> ret;
-	  for(int i=0;i<N;i++){
-	    ret._internal=arg._internal[i]*G;
-	  }
-	  return ret;
-	}
-    template<class vtype,int N> inline auto operator * (const iMatrix<vtype,N> &arg, const Gamma &G) ->
-      typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N> >::type 
-	{
-	  iMatrix<vtype,N> ret;
-	  for(int i=0;i<N;i++){
-	  for(int j=0;j<N;j++){
-	    ret._internal[i][j]=arg._internal[i][j]*G;
-	  }}
-	  return ret;
-	}
-
-    ////////////////////////////////////////////////////////
-    // When we hit the spin index this matches and we stop
-    ////////////////////////////////////////////////////////
-    template<class vtype> inline auto operator * ( const Gamma &G,const iMatrix<vtype,Ns> &arg) ->
-      typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,Ns>,SpinorIndex>::value,iMatrix<vtype,Ns> >::type 
-      {
-	iMatrix<vtype,Ns> ret;
-	switch (G._g) {
-	case Gamma::Identity:
-	  ret = arg;
-	  break;
-	case Gamma::MinusIdentity:
-	  ret = -arg;
-	  break;
-	case Gamma::GammaX:
-	  multGammaX(ret,arg);
-	  break;
-	case Gamma::MinusGammaX:
-	  multMinusGammaX(ret,arg);
-	  break;
-	case Gamma::GammaY:
-	  multGammaY(ret,arg);
-	  break;
-	case Gamma::MinusGammaY:
-	  multMinusGammaY(ret,arg);
-	  break;
-	case Gamma::GammaZ:
-	  multGammaZ(ret,arg);
-	  break;
-	case Gamma::MinusGammaZ:
-	  multMinusGammaZ(ret,arg);
-	  break;
-	case Gamma::GammaT:
-	  multGammaT(ret,arg);
-	  break;
-	case Gamma::MinusGammaT:
-	  multMinusGammaT(ret,arg);
-	  break;
-	case Gamma::Gamma5:
-	  multGamma5(ret,arg);
-	  break;
-	case Gamma::MinusGamma5:
-	  multMinusGamma5(ret,arg);
-	  break;
-	default:
-	  assert(0);
-	  break;
-	}
-	return ret;
-      }
-    // Could have used type trait for Matrix/vector and then an enable if to share code
-    template<class vtype> inline auto operator * ( const Gamma &G,const iVector<vtype,Ns> &arg) ->
-      typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type 
-      {
-	iVector<vtype,Ns> ret;
-	switch (G._g) {
-	case Gamma::Identity:
-	  ret = arg;
-	  break;
-	case Gamma::MinusIdentity:
-	  ret = -arg;
-	  break;
-	case Gamma::GammaX:
-	  multGammaX(ret,arg);
-	  break;
-	case Gamma::MinusGammaX:
-	  multMinusGammaX(ret,arg);
-	  break;
-	case Gamma::GammaY:
-	  multGammaY(ret,arg);
-	  break;
-	case Gamma::MinusGammaY:
-	  multMinusGammaY(ret,arg);
-	  break;
-	case Gamma::GammaZ:
-	  multGammaZ(ret,arg);
-	  break;
-	case Gamma::MinusGammaZ:
-	  multMinusGammaZ(ret,arg);
-	  break;
-	case Gamma::GammaT:
-	  multGammaT(ret,arg);
-	  break;
-	case Gamma::MinusGammaT:
-	  multMinusGammaT(ret,arg);
-	  break;
-	case Gamma::Gamma5:
-	  multGamma5(ret,arg);
-	  break;
-	case Gamma::MinusGamma5:
-	  multMinusGamma5(ret,arg);
-	  break;
-	default:
-	  assert(0);
-	  break;
-	}
-	return ret;
-      }
-
-    template<class vtype> inline auto operator * (const iMatrix<vtype,Ns> &arg, const Gamma &G) ->
-      typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,Ns>,SpinorIndex>::value,iMatrix<vtype,Ns> >::type 
-      {
-	iMatrix<vtype,Ns> ret;
-	switch (G._g) {
-	case Gamma::Identity:
-	  ret = arg;
-	  break;
-	case Gamma::MinusIdentity:
-	  ret = -arg;
-	  break;
-	case Gamma::GammaX:
-	  rmultGammaX(ret,arg);
-	  break;
-	case Gamma::MinusGammaX:
-	  rmultMinusGammaX(ret,arg);
-	  break;
-	case Gamma::GammaY:
-	  rmultGammaY(ret,arg);
-	  break;
-	case Gamma::MinusGammaY:
-	  rmultMinusGammaY(ret,arg);
-	  break;
-	case Gamma::GammaZ:
-	  rmultGammaZ(ret,arg);
-	  break;
-	case Gamma::MinusGammaZ:
-	  rmultMinusGammaZ(ret,arg);
-	  break;
-	case Gamma::GammaT:
-	  rmultGammaT(ret,arg);
-	  break;
-	case Gamma::MinusGammaT:
-	  rmultMinusGammaT(ret,arg);
-	  break;
-	case Gamma::Gamma5:
-	  rmultGamma5(ret,arg);
-	  break;
-	case Gamma::MinusGamma5:
-	  rmultMinusGamma5(ret,arg);
-	  break;
-	default:
-	  assert(0);
-	  break;
-	}
-	return ret;
-      }
-
-
-    /* Output from test
-./Grid_gamma 
-Identity((1,0),(0,0),(0,0),(0,0))
-        ((0,0),(1,0),(0,0),(0,0))
-        ((0,0),(0,0),(1,0),(0,0))
-        ((0,0),(0,0),(0,0),(1,0))   OK
-
-GammaX  ((0,0),(0,0),(0,0),(0,1))
-        ((0,0),(0,0),(0,1),(0,0))
-        ((0,0),(0,-1),(0,0),(0,0))
-        ((0,-1),(0,0),(0,0),(0,0))  OK
-    * Gx
-    *  0 0  0  i    
-    *  0 0  i  0    
-    *  0 -i 0  0
-    * -i 0  0  0
-
-GammaY  ((-0,-0),(-0,-0),(-0,-0),(-1,-0))
-        ((0,0),(0,0),(1,0),(0,0))
-        ((0,0),(1,0),(0,0),(0,0))          OK
-        ((-1,-0),(-0,-0),(-0,-0),(-0,-0))
-     *Gy
-     *  0 0  0  -1  [0] -+ [3]
-     *  0 0  1  0   [1] +- [2]
-     *  0 1  0  0
-     * -1 0  0  0
-
-GammaZ  ((0,0),(0,0),(0,1),(0,0))
-        ((0,0),(0,0),(0,0),(0,-1))
-        ((0,-1),(0,0),(0,0),(0,0))
-        ((0,0),(0,1),(0,0),(0,0))   OK
-     *  0 0  i  0   [0]+-i[2]
-     *  0 0  0 -i   [1]-+i[3]
-     * -i 0  0  0
-     *  0 i  0  0
-
-GammaT  ((0,0),(0,0),(1,0),(0,0))
-        ((0,0),(0,0),(0,0),(1,0))  OK
-        ((1,0),(0,0),(0,0),(0,0))
-        ((0,0),(1,0),(0,0),(0,0))
-     *  0 0  1  0 [0]+-[2]
-     *  0 0  0  1 [1]+-[3]
-     *  1 0  0  0
-     *  0 1  0  0
-
-Gamma5  ((1,0),(0,0),(0,0),(0,0))
-        ((0,0),(1,0),(0,0),(0,0))
-        ((-0,-0),(-0,-0),(-1,-0),(-0,-0))
-        ((-0,-0),(-0,-0),(-0,-0),(-1,-0))
-     *  1 0  0  0 [0]+-[2]
-     *  0 1  0  0 [1]+-[3]  OK
-     *  0 0 -1  0
-     *  0 0  0 -1
-     */
-
-}   //namespace QCD
-} // Grid
 #endif

lib/qcd/utils/LinalgUtils.h

@@ -47,7 +47,7 @@ void axpibg5x(Lattice<vobj> &z,const Lattice<vobj> &x,Coeff a,Coeff b)
 
   GridBase *grid=x._grid;
 
-  Gamma G5(Gamma::Gamma5);
+  Gamma G5(Gamma::Algebra::Gamma5);
   parallel_for(int ss=0;ss<grid->oSites();ss++){
     vobj tmp;
     tmp = a*x._odata[ss];
@@ -78,7 +78,8 @@ void ag5xpby_ssp(Lattice<vobj> &z,Coeff a,const Lattice<vobj> &x,Coeff b,const L
   conformable(x,z);
   GridBase *grid=x._grid;
   int Ls = grid->_rdimensions[0];
-  Gamma G5(Gamma::Gamma5);
+
+  Gamma G5(Gamma::Algebra::Gamma5);
   parallel_for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     vobj tmp;
     tmp = G5*x._odata[ss+s]*a;
@@ -95,7 +96,7 @@ void axpbg5y_ssp(Lattice<vobj> &z,Coeff a,const Lattice<vobj> &x,Coeff b,const L
   conformable(x,z);
   GridBase *grid=x._grid;
   int Ls = grid->_rdimensions[0];
-  Gamma G5(Gamma::Gamma5);
+  Gamma G5(Gamma::Algebra::Gamma5);
   parallel_for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     vobj tmp;
     tmp = G5*y._odata[ss+sp]*b;
@@ -112,7 +113,8 @@ void ag5xpbg5y_ssp(Lattice<vobj> &z,Coeff a,const Lattice<vobj> &x,Coeff b,const
   conformable(x,z);
   GridBase *grid=x._grid;
   int Ls = grid->_rdimensions[0];
-  Gamma G5(Gamma::Gamma5);
+
+  Gamma G5(Gamma::Algebra::Gamma5);
   parallel_for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     vobj tmp1;
     vobj tmp2;
@@ -161,8 +163,8 @@ void G5R5(Lattice<vobj> &z,const Lattice<vobj> &x)
   z.checkerboard = x.checkerboard;
   conformable(x,z);
   int Ls = grid->_rdimensions[0];
-  Gamma G5(Gamma::Gamma5);
+  Gamma G5(Gamma::Algebra::Gamma5);
-  parallel_for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+  parallel_for(int ss=0;ss<grid->oSites();ss+=Ls) {
     vobj tmp;
     for(int s=0;s<Ls;s++){
       int sp = Ls-1-s;

lib/serialisation/BaseIO.h

@@ -32,6 +32,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <type_traits>
 
 namespace Grid {
+  // Vector IO utilities ///////////////////////////////////////////////////////
   // helper function to read space-separated values
   template <typename T>
   std::vector<T> strToVec(const std::string s)
@@ -67,6 +68,77 @@ namespace Grid {
     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 flatening 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};
+  };
+  
+  // Abstract writer/reader classes ////////////////////////////////////////////
   // static polymorphism implemented using CRTP idiom
   class Serializable;
   
@@ -83,12 +155,7 @@ namespace Grid {
     typename std::enable_if<std::is_base_of<Serializable, U>::value, void>::type
     write(const std::string& s, const U &output);
     template <typename U>
-    typename std::enable_if<std::is_enum<U>::value, void>::type
+    typename std::enable_if<!std::is_base_of<Serializable, U>::value, void>::type
-    write(const std::string& s, const U &output);
-    template <typename U>
-    typename std::enable_if<
-      !(std::is_base_of<Serializable, U>::value or std::is_enum<U>::value),
-      void>::type
     write(const std::string& s, const U &output);
   private:
     T *upcast;
@@ -107,12 +174,7 @@ namespace Grid {
     typename std::enable_if<std::is_base_of<Serializable, U>::value, void>::type
     read(const std::string& s, U &output);
     template <typename U>
-    typename std::enable_if<std::is_enum<U>::value, void>::type
+    typename std::enable_if<!std::is_base_of<Serializable, U>::value, void>::type
-    read(const std::string& s, U &output);
-    template <typename U>
-    typename std::enable_if<
-      !(std::is_base_of<Serializable, U>::value or std::is_enum<U>::value),
-      void>::type
     read(const std::string& s, U &output);
   protected:
     template <typename U>
@@ -142,7 +204,128 @@ namespace Grid {
     }
   };
   
-  // Generic writer interface
+  // 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_;
+  }
+  
+  // Generic writer interface //////////////////////////////////////////////////
   template <typename T>
   inline void push(Writer<T> &w, const std::string &s)
   {
@@ -221,23 +404,13 @@ namespace Grid {
   
   template <typename T>
   template <typename U>
-  typename std::enable_if<std::is_enum<U>::value, void>::type
+  typename std::enable_if<!std::is_base_of<Serializable, U>::value, void>::type
-  Writer<T>::write(const std::string &s, const U &output)
-  {
-    EnumIO<U>::write(*this, s, output);
-  }
-  
-  template <typename T>
-  template <typename U>
-  typename std::enable_if<
-    !(std::is_base_of<Serializable, U>::value or std::is_enum<U>::value),
-    void>::type
   Writer<T>::write(const std::string &s, const U &output)
   {
     upcast->writeDefault(s, output);
   }
   
-  // Reader template implementation ////////////////////////////////////////////
+  // Reader template implementation
   template <typename T>
   Reader<T>::Reader(void)
   {
@@ -266,17 +439,7 @@ namespace Grid {
   
   template <typename T>
   template <typename U>
-  typename std::enable_if<std::is_enum<U>::value, void>::type
+  typename std::enable_if<!std::is_base_of<Serializable, U>::value, void>::type
-  Reader<T>::read(const std::string &s, U &output)
-  {
-    EnumIO<U>::read(*this, s, output);
-  }
-  
-  template <typename T>
-  template <typename U>
-  typename std::enable_if<
-    !(std::is_base_of<Serializable, U>::value or std::is_enum<U>::value),
-    void>::type
   Reader<T>::read(const std::string &s, U &output)
   {
     upcast->readDefault(s, output);
@@ -300,7 +463,6 @@ namespace Grid {
       abort();
     }
   }
-
 }
 
 #endif

lib/serialisation/BinaryIO.cc

@@ -26,7 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-
 #include <Grid/GridCore.h>
 
 using namespace Grid;

lib/serialisation/Hdf5IO.cc

@@ -0,0 +1,103 @@
+#include <Grid/Grid.h>
+
+using namespace Grid;
+#ifndef H5_NO_NAMESPACE
+using namespace H5NS;
+#endif
+
+// Writer implementation ///////////////////////////////////////////////////////
+Hdf5Writer::Hdf5Writer(const std::string &fileName)
+: fileName_(fileName)
+, file_(fileName.c_str(), H5F_ACC_TRUNC)
+{
+  group_ = file_.openGroup("/");
+  writeSingleAttribute(dataSetThres_, HDF5_GRID_GUARD "dataset_threshold",
+                       Hdf5Type<unsigned int>::type());
+}
+
+void Hdf5Writer::push(const std::string &s)
+{
+  group_ = group_.createGroup(s);
+  path_.push_back(s);
+}
+
+void Hdf5Writer::pop(void)
+{
+  path_.pop_back();
+  if (path_.empty())
+  {
+    group_ = file_.openGroup("/");
+  }
+  else
+  {
+    auto binOp = [](const std::string &a, const std::string &b)->std::string
+    {
+      return a + "/" + b;
+    };
+    
+    group_ = group_.openGroup(std::accumulate(path_.begin(), path_.end(),
+                                              std::string(""), binOp));
+  }
+}
+
+template <>
+void Hdf5Writer::writeDefault(const std::string &s, const std::string &x)
+{
+  StrType     strType(PredType::C_S1, x.size());
+  
+  writeSingleAttribute(*(x.data()), s, strType);
+}
+
+void Hdf5Writer::writeDefault(const std::string &s, const char *x)
+{
+  std::string sx(x);
+  
+  writeDefault(s, sx);
+}
+
+// Reader implementation ///////////////////////////////////////////////////////
+Hdf5Reader::Hdf5Reader(const std::string &fileName)
+: fileName_(fileName)
+, file_(fileName.c_str(), H5F_ACC_RDONLY)
+{
+  group_ = file_.openGroup("/");
+  readSingleAttribute(dataSetThres_, HDF5_GRID_GUARD "dataset_threshold",
+                      Hdf5Type<unsigned int>::type());
+}
+
+void Hdf5Reader::push(const std::string &s)
+{
+  group_ = group_.openGroup(s);
+  path_.push_back(s);
+}
+
+void Hdf5Reader::pop(void)
+{
+  path_.pop_back();
+  if (path_.empty())
+  {
+    group_ = file_.openGroup("/");
+  }
+  else
+  {
+    auto binOp = [](const std::string &a, const std::string &b)->std::string
+    {
+      return a + "/" + b;
+    };
+    
+    group_ = group_.openGroup(std::accumulate(path_.begin(), path_.end(),
+                                              std::string(""), binOp));
+  }
+}
+
+template <>
+void Hdf5Reader::readDefault(const std::string &s, std::string &x)
+{
+  Attribute attribute;
+  
+  attribute       = group_.openAttribute(s);
+  StrType strType = attribute.getStrType();
+  
+  x.resize(strType.getSize());
+  attribute.read(strType, &(x[0]));
+}

lib/serialisation/Hdf5IO.h

@@ -0,0 +1,242 @@
+#ifndef GRID_SERIALISATION_HDF5_H
+#define GRID_SERIALISATION_HDF5_H
+
+#include <stack>
+#include <string>
+#include <vector>
+#include <H5Cpp.h>
+#include "Hdf5Type.h"
+
+#ifndef H5_NO_NAMESPACE
+#define H5NS H5
+#endif
+
+// default thresold above which datasets are used instead of attributes
+#ifndef HDF5_DEF_DATASET_THRES
+#define HDF5_DEF_DATASET_THRES 6u
+#endif
+
+// name guard for Grid metadata
+#define HDF5_GRID_GUARD "_Grid_"
+
+namespace Grid
+{
+  class Hdf5Writer: public Writer<Hdf5Writer>
+  {
+  public:
+    Hdf5Writer(const std::string &fileName);
+    virtual ~Hdf5Writer(void) = default;
+    void push(const std::string &s);
+    void pop(void);
+    void writeDefault(const std::string &s, const char *x);
+    template <typename U>
+    void writeDefault(const std::string &s, const U &x);
+    template <typename U>
+    typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+    writeDefault(const std::string &s, const std::vector<U> &x);
+    template <typename U>
+    typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+    writeDefault(const std::string &s, const std::vector<U> &x);
+  private:
+    template <typename U>
+    void writeSingleAttribute(const U &x, const std::string &name,
+                              const H5NS::DataType &type);
+  private:
+    std::string              fileName_;
+    std::vector<std::string> path_;
+    H5NS::H5File             file_;
+    H5NS::Group              group_;
+    unsigned int             dataSetThres_{HDF5_DEF_DATASET_THRES};
+  };
+  
+  class Hdf5Reader: public Reader<Hdf5Reader>
+  {
+  public:
+    Hdf5Reader(const std::string &fileName);
+    virtual ~Hdf5Reader(void) = default;
+    void push(const std::string &s);
+    void pop(void);
+    template <typename U>
+    void readDefault(const std::string &s, U &output);
+    template <typename U>
+    typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+    readDefault(const std::string &s, std::vector<U> &x);
+    template <typename U>
+    typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+    readDefault(const std::string &s, std::vector<U> &x);
+  private:
+    template <typename U>
+    void readSingleAttribute(U &x, const std::string &name,
+                             const H5NS::DataType &type);
+  private:
+    std::string              fileName_;
+    std::vector<std::string> path_;
+    H5NS::H5File             file_;
+    H5NS::Group              group_;
+    unsigned int             dataSetThres_;
+  };
+  
+  // Writer template implementation ////////////////////////////////////////////
+  template <typename U>
+  void Hdf5Writer::writeSingleAttribute(const U &x, const std::string &name,
+                                        const H5NS::DataType &type)
+  {
+    H5NS::Attribute attribute;
+    hsize_t         attrDim = 1;
+    H5NS::DataSpace attrSpace(1, &attrDim);
+    
+    attribute = group_.createAttribute(name, type, attrSpace);
+    attribute.write(type, &x);
+  }
+  
+  template <typename U>
+  void Hdf5Writer::writeDefault(const std::string &s, const U &x)
+  {
+    writeSingleAttribute(x, s, Hdf5Type<U>::type());
+  }
+  
+  template <>
+  void Hdf5Writer::writeDefault(const std::string &s, const std::string &x);
+  
+  template <typename U>
+  typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+  Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
+  {
+    // alias to element type
+    typedef typename element<std::vector<U>>::type Element;
+    
+    // flatten the vector and getting dimensions
+    Flatten<std::vector<U>> flat(x);
+    std::vector<hsize_t> dim;
+    const auto           &flatx = flat.getFlatVector();
+    
+    for (auto &d: flat.getDim())
+    {
+      dim.push_back(d);
+    }
+    
+    // write to file
+    H5NS::DataSpace dataSpace(dim.size(), dim.data());
+    
+    if (flatx.size() > dataSetThres_)
+    {
+      H5NS::DataSet dataSet;
+      
+      dataSet = group_.createDataSet(s, Hdf5Type<Element>::type(), dataSpace);
+      dataSet.write(flatx.data(), Hdf5Type<Element>::type());
+    }
+    else
+    {
+      H5NS::Attribute attribute;
+      
+      attribute = group_.createAttribute(s, Hdf5Type<Element>::type(), dataSpace);
+      attribute.write(Hdf5Type<Element>::type(), flatx.data());
+    }
+  }
+  
+  template <typename U>
+  typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+  Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
+  {
+    push(s);
+    writeSingleAttribute(x.size(), HDF5_GRID_GUARD "vector_size",
+                         Hdf5Type<uint64_t>::type());
+    for (hsize_t i = 0; i < x.size(); ++i)
+    {
+      write(s + "_" + std::to_string(i), x[i]);
+    }
+    pop();
+  }
+  
+  // Reader template implementation ////////////////////////////////////////////
+  template <typename U>
+  void Hdf5Reader::readSingleAttribute(U &x, const std::string &name,
+                                       const H5NS::DataType &type)
+  {
+    H5NS::Attribute attribute;
+    
+    attribute = group_.openAttribute(name);
+    attribute.read(type, &x);
+  }
+  
+  template <typename U>
+  void Hdf5Reader::readDefault(const std::string &s, U &output)
+  {
+    readSingleAttribute(output, s, Hdf5Type<U>::type());
+  }
+  
+  template <>
+  void Hdf5Reader::readDefault(const std::string &s, std::string &x);
+  
+  template <typename U>
+  typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+  Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
+  {
+    // alias to element type
+    typedef typename element<std::vector<U>>::type Element;
+    
+    // read the dimensions
+    H5NS::DataSpace       dataSpace;
+    std::vector<hsize_t>  hdim;
+    std::vector<size_t>   dim;
+    hsize_t               size = 1;
+    
+    if (group_.attrExists(s))
+    {
+      dataSpace = group_.openAttribute(s).getSpace();
+    }
+    else
+    {
+      dataSpace = group_.openDataSet(s).getSpace();
+    }
+    hdim.resize(dataSpace.getSimpleExtentNdims());
+    dataSpace.getSimpleExtentDims(hdim.data());
+    for (auto &d: hdim)
+    {
+      dim.push_back(d);
+      size *= d;
+    }
+    
+    // read the flat vector
+    std::vector<Element> buf(size);
+
+    if (size > dataSetThres_)
+    {
+      H5NS::DataSet dataSet;
+      
+      dataSet = group_.openDataSet(s);
+      dataSet.read(buf.data(), Hdf5Type<Element>::type());
+    }
+    else
+    {
+      H5NS::Attribute attribute;
+      
+      attribute = group_.openAttribute(s);
+      attribute.read(Hdf5Type<Element>::type(), buf.data());
+    }
+    
+    // reconstruct the multidimensional vector
+    Reconstruct<std::vector<U>> r(buf, dim);
+    
+    x = r.getVector();
+  }
+  
+  template <typename U>
+  typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
+  Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
+  {
+    uint64_t size;
+    
+    push(s);
+    readSingleAttribute(size, HDF5_GRID_GUARD "vector_size",
+                        Hdf5Type<uint64_t>::type());
+    x.resize(size);
+    for (hsize_t i = 0; i < x.size(); ++i)
+    {
+      read(s + "_" + std::to_string(i), x[i]);
+    }
+    pop();
+  }
+}
+
+#endif

lib/serialisation/Hdf5Type.h

@@ -0,0 +1,77 @@
+#ifndef GRID_SERIALISATION_HDF5_TYPE_H
+#define GRID_SERIALISATION_HDF5_TYPE_H
+
+#include <H5Cpp.h>
+#include <complex>
+#include <memory>
+
+#ifndef H5_NO_NAMESPACE
+#define H5NS H5
+#endif
+
+#define HDF5_NATIVE_TYPE(predType, cType)\
+template <>\
+class Hdf5Type<cType>\
+{\
+public:\
+  static inline const H5NS::DataType & type(void)\
+  {\
+    return H5NS::PredType::predType;\
+  }\
+  static constexpr bool isNative = true;\
+};
+
+#define DEFINE_HDF5_NATIVE_TYPES \
+HDF5_NATIVE_TYPE(NATIVE_B8,      bool);\
+HDF5_NATIVE_TYPE(NATIVE_CHAR,    char);\
+HDF5_NATIVE_TYPE(NATIVE_SCHAR,   signed char);\
+HDF5_NATIVE_TYPE(NATIVE_UCHAR,   unsigned char);\
+HDF5_NATIVE_TYPE(NATIVE_SHORT,   short);\
+HDF5_NATIVE_TYPE(NATIVE_USHORT,  unsigned short);\
+HDF5_NATIVE_TYPE(NATIVE_INT,     int);\
+HDF5_NATIVE_TYPE(NATIVE_UINT,    unsigned int);\
+HDF5_NATIVE_TYPE(NATIVE_LONG,    long);\
+HDF5_NATIVE_TYPE(NATIVE_ULONG,   unsigned long);\
+HDF5_NATIVE_TYPE(NATIVE_LLONG,   long long);\
+HDF5_NATIVE_TYPE(NATIVE_ULLONG,  unsigned long long);\
+HDF5_NATIVE_TYPE(NATIVE_FLOAT,   float);\
+HDF5_NATIVE_TYPE(NATIVE_DOUBLE,  double);\
+HDF5_NATIVE_TYPE(NATIVE_LDOUBLE, long double);
+
+namespace Grid
+{
+  template <typename T> class Hdf5Type
+  {
+  public:
+    static constexpr bool isNative = false;
+  };
+  
+  DEFINE_HDF5_NATIVE_TYPES;
+  
+  template <typename R>
+  class Hdf5Type<std::complex<R>>
+  {
+  public:
+    static inline const H5NS::DataType & type(void)
+    {
+      if (typePtr_ == nullptr)
+      {
+        typePtr_.reset(new H5NS::CompType(sizeof(std::complex<R>)));
+        typePtr_->insertMember("re", 0,         Hdf5Type<R>::type());
+        typePtr_->insertMember("im", sizeof(R), Hdf5Type<R>::type());
+      }
+
+      return *typePtr_;
+    }
+    static constexpr bool isNative = false;
+  private:
+    static std::unique_ptr<H5NS::CompType> typePtr_;
+  };
+  
+  template <typename R>
+  std::unique_ptr<H5NS::CompType> Hdf5Type<std::complex<R>>::typePtr_ = nullptr;
+}
+
+#undef HDF5_NATIVE_TYPE
+
+#endif /* GRID_SERIALISATION_HDF5_TYPE_H */

lib/serialisation/MacroMagic.h

@@ -109,40 +109,36 @@ THE SOFTWARE.
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
 #define GRID_MACRO_MEMBER(A,B)        A B;
+#define GRID_MACRO_COMP_MEMBER(A,B) result = (result and (lhs. B == rhs. B));
 #define GRID_MACRO_OS_WRITE_MEMBER(A,B) os<< #A <<" "#B <<" = "<< obj. B <<" ; " <<std::endl;
 #define GRID_MACRO_READ_MEMBER(A,B) Grid::read(RD,#B,obj. B);
 #define GRID_MACRO_WRITE_MEMBER(A,B) Grid::write(WR,#B,obj. B);
 
-#define GRID_SERIALIZABLE_CLASS_MEMBERS(cname,...)		\
+#define GRID_SERIALIZABLE_CLASS_MEMBERS(cname,...)\
-  \
+GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_MEMBER,__VA_ARGS__))\
-  \
+template <typename T>\
-  GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_MEMBER,__VA_ARGS__))		\
+static inline void write(Writer<T> &WR,const std::string &s, const cname &obj){ \
-  \
+  push(WR,s);\
-  \
+  GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_WRITE_MEMBER,__VA_ARGS__))	\
-  template <typename T>\
+  pop(WR);\
-  static inline void write(Writer<T> &WR,const std::string &s, const cname &obj){ \
+}\
-    push(WR,s);\
+template <typename T>\
-    GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_WRITE_MEMBER,__VA_ARGS__))	\
+static inline void read(Reader<T> &RD,const std::string &s, cname &obj){	\
-    pop(WR);\
+  push(RD,s);\
-  } \
+  GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_READ_MEMBER,__VA_ARGS__))	\
-  \
+  pop(RD);\
-  \
+}\
-  template <typename T>\
+friend inline std::ostream & operator << (std::ostream &os, const cname &obj ) { \
-  static inline void read(Reader<T> &RD,const std::string &s, cname &obj){	\
+  os<<"class "<<#cname<<" {"<<std::endl;\
-    push(RD,s);\
+  GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_OS_WRITE_MEMBER,__VA_ARGS__))	\
-    GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_READ_MEMBER,__VA_ARGS__))	\
+    os<<"}";								\
-    pop(RD);\
+  return os;\
-  } \
+}\
-  \
+friend inline bool operator==(const cname &lhs, const cname &rhs) {\
-  \
+  bool result = true;\
-  friend inline std::ostream & operator << (std::ostream &os, const cname &obj ) { \
+  GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_COMP_MEMBER,__VA_ARGS__))\
-    os<<"class "<<#cname<<" {"<<std::endl;\
+  return result;\
-    GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_OS_WRITE_MEMBER,__VA_ARGS__))	\
+}
-      os<<"}";								\
-    return os;\
-  };  
-
-
 
 #define GRID_ENUM_TYPE(obj) std::remove_reference<decltype(obj)>::type
 #define GRID_MACRO_ENUMVAL(A,B) A = B,
@@ -150,44 +146,61 @@ THE SOFTWARE.
 #define GRID_MACRO_ENUMTEST(A,B) else if (buf == #A) {obj = GRID_ENUM_TYPE(obj)::A;}
 #define GRID_MACRO_ENUMCASEIO(A,B) case GRID_ENUM_TYPE(obj)::A: os << #A; break;
 
-namespace Grid {
-  template <typename U>
-  class EnumIO {};
-}
-
 #define GRID_SERIALIZABLE_ENUM(name,undefname,...)\
-  enum class name {\
+class name: public Grid::Serializable\
-      GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMVAL,__VA_ARGS__))\
+{\
-      undefname = -1\
+public:\
+  enum\
+  {\
+    GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMVAL,__VA_ARGS__))\
+    undefname = -1\
   };\
+public:\
+  name(void): value_(undefname) {};\
+  name(int value): value_(value) {};\
+  template <typename T>\
+  static inline void write(Grid::Writer<T> &WR,const std::string &s, const name &obj)\
+  {\
+    switch (obj.value_)\
+    {\
+      GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMCASE,__VA_ARGS__))\
+      default: Grid::write(WR,s,#undefname); break;\
+    }\
+  }\
   \
-  template<>\
+  template <typename T>\
-  class EnumIO<name> {\
+  static inline void read(Grid::Reader<T> &RD,const std::string &s, name &obj)\
-    public:\
+  {\
-      template <typename T>\
+    std::string buf;\
-      static inline void write(Writer<T> &WR,const std::string &s, const name &obj){ \
+    Grid::read(RD, s, buf);\
-        switch (obj) {\
+    if (buf == #undefname) {obj = name::undefname;}\
-          GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMCASE,__VA_ARGS__))\
+    GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMTEST,__VA_ARGS__))\
-          default: Grid::write(WR,s,#undefname); break;\
+    else {obj = name::undefname;}\
-        }\
+  }\
-      }\
+  inline operator int(void) const\
-      \
+  {\
-      template <typename T>\
+    return value_;\
-      static inline void read(Reader<T> &RD,const std::string &s, name &obj){ \
+  }\
-        std::string buf;\
+  inline friend std::ostream & operator<<(std::ostream &os, const name &obj)\
-        Grid::read(RD, s, buf);\
+  {\
-        if (buf == #undefname) {obj = name::undefname;}\
-        GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMTEST,__VA_ARGS__))\
-        else {obj = name::undefname;}\
-      }\
-  };\
-  \
-  inline std::ostream & operator << (std::ostream &os, const name &obj ) { \
     switch (obj) {\
-        GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMCASEIO,__VA_ARGS__))\
+      GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMCASEIO,__VA_ARGS__))\
-        default: os << #undefname; break;\
+      default: os << #undefname; break;\
     }\
     return os;\
-  };
+  }\
+  inline friend std::istream & operator>>(std::istream &is, name &obj)\
+  {\
+    std::string buf;\
+    is >> buf;\
+    if (buf == #undefname) {obj = name::undefname;}\
+    GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMTEST,__VA_ARGS__))\
+    else {obj = name::undefname;}\
+    return is;\
+  }\
+private:\
+  int value_;\
+};
+
 
 #endif

lib/serialisation/Serialisation.h

@@ -36,6 +36,9 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include "BinaryIO.h"
 #include "TextIO.h"
 #include "XmlIO.h"
+#ifdef HAVE_HDF5
+#include "Hdf5IO.h"
+#endif
 //////////////////////////////////////////
 // Todo:
 //////////////////////////////////////////

lib/sitmo_rng/sitmo_prng_engine.hpp

@@ -0,0 +1,390 @@
+//  Copyright (c) 2012-2016 M.A. (Thijs) van den Berg, http://sitmo.com/
+//
+//  Use, modification and distribution are subject to the MIT Software License. 
+//  
+//  The MIT License (MIT)
+//  Permission is hereby granted, free of charge, to any person obtaining a copy
+//  of this software and associated documentation files (the "Software"), to deal
+//  in the Software without restriction, including without limitation the rights
+//  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+//  copies of the Software, and to permit persons to whom the Software is
+//  furnished to do so, subject to the following conditions:
+//  
+//  The above copyright notice and this permission notice shall be included in
+//  all copies or substantial portions of the Software.
+//  
+//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+//  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+//  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+//  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+//  THE SOFTWARE.
+
+// version history:
+// version 1,  6 Sep 2012
+// version 2, 10 Dec 2013
+//      bug fix in the discard() routine, it was discarding to many elements
+//      added the version() method
+// version 3...5, 13 Dec 2013
+//      fixed type-conversion earning
+//      fixed potential issues with constructor template matching
+// version 6, 4 March 2016
+//      made min() max() constexpr for C+11 compiler (thanks to James Joseph Balamuta)
+
+#ifndef SITMO_PRNG_ENGINE_HPP
+#define SITMO_PRNG_ENGINE_HPP
+#include <iostream>
+
+#ifdef __GNUC__
+    #include <stdint.h>                 // respecting the C99 standard.
+#endif
+#ifdef _MSC_VER
+    typedef unsigned __int64 uint64_t;  // Visual Studio 6.0(VC6) and newer..
+    typedef unsigned __int32 uint32_t;
+#endif
+
+// Double mixing function
+#define MIX2(x0,x1,rx,z0,z1,rz) \
+    x0 += x1; \
+    z0 += z1; \
+    x1 = (x1 << rx) | (x1 >> (64-rx)); \
+    z1 = (z1 << rz) | (z1 >> (64-rz)); \
+    x1 ^= x0; \
+    z1 ^= z0; 
+
+
+// Double mixing function with key adition
+#define MIXK(x0,x1,rx,z0,z1,rz,k0,k1,l0,l1) \
+    x1 += k1; \
+    z1 += l1; \
+    x0 += x1+k0; \
+    z0 += z1+l0; \
+    x1 = (x1 << rx) | (x1 >> (64-rx)); \
+    z1 = (z1 << rz) | (z1 >> (64-rz)); \
+    x1 ^= x0; \
+    z1 ^= z0; \
+
+
+namespace sitmo {
+
+    // enable_if for C__98 compilers
+    template<bool C, typename T = void>
+    struct sitmo_enable_if { typedef T type; };
+
+    template<typename T>
+    struct sitmo_enable_if<false, T> { };
+
+    // SFINAE check for the existence of a "void generate(int*,int*)"member function
+    template<typename T>
+    struct has_generate_template
+    {
+        typedef char (&Two)[2];;
+        template<typename F, void (F::*)(int *, int *)> struct helper {};
+        template<typename C> static char test(helper<C, &C::template generate<int*> >*);
+        template<typename C> static Two test(...);
+        static bool const value = sizeof(test<T>(0)) == sizeof(char);
+    };
+
+
+class prng_engine
+{
+public:
+    // "req" are requirements as stated in the C++ 11 draft n3242=11-0012
+    //
+    // req: 26.5.1.3 Uniform random number generator requirements, p.906, table 116, row 1
+    typedef uint32_t result_type;
+
+    // req: 26.5.1.3 Uniform random number generator requirements, p.906, table 116, row 3 & 4
+#if __cplusplus <= 199711L
+    static result_type (min)() { return 0; }
+    static result_type (max)() { return 0xFFFFFFFF; }
+#else
+    static constexpr result_type (min)() { return 0; }
+    static constexpr result_type (max)() { return 0xFFFFFFFF; }
+#endif    
+
+    // -------------------------------------------------
+    // Constructors
+    // -------------------------------------------------
+
+    // req: 26.5.1.4 Random number engine requirements, p.907 table 117, row 1
+    // Creates an engine with the same initial state as all other
+    // default-constructed engines of type E.
+    prng_engine() 
+    {
+        seed();
+    }
+
+    // req: 26.5.1.4 Random number engine requirements, p.907 table 117, row 2
+    // Creates an engine that compares equal to x.
+    prng_engine(const prng_engine& x)
+    {
+        for (unsigned short i=0; i<4; ++i) {
+            _s[i] = x._s[i];
+            _k[i] = x._k[i];
+            _o[i] = x._o[i];
+        }
+        _o_counter = x._o_counter;
+    }
+    
+    
+    // req: 26.5.1.4 Random number engine requirements, p.907 table 117, row 3
+    // Creates an engine with initial O(size of state) state determined by s.
+    prng_engine(uint32_t s) 
+    {
+        seed(s);
+    }
+    
+        // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 4
+    // Creates an engine with an initial state that depends on a sequence
+    // produced by one call to q.generate.
+    template<class Seq> 
+    prng_engine(Seq& q, typename sitmo_enable_if< has_generate_template<Seq>::value >::type* = 0 )
+    {
+        seed(q);
+    }
+
+    // -------------------------------------------------
+    // Seeding
+    // -------------------------------------------------
+    
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 5
+    void seed()
+    {
+        for (unsigned short i=0; i<4; ++i) {
+            _k[i] = 0;
+            _s[i] = 0;
+        }
+        _o_counter = 0;
+
+        _o[0] = 0x09218ebde6c85537;
+        _o[1] = 0x55941f5266d86105;
+        _o[2] = 0x4bd25e16282434dc;
+        _o[3] = 0xee29ec846bd2e40b;
+    }
+
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 6
+    // s needs to be of return_type, which is uint32_t
+    void seed(uint32_t s)
+    { 
+        for (unsigned short i=0; i<4; ++i) {
+            _k[i] = 0;
+            _s[i] = 0;
+        }
+        _k[0] = s; 
+        _o_counter = 0;
+        
+        encrypt_counter();
+    }
+    
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 7
+    template<class Seq> 
+    void seed(Seq& q, typename sitmo_enable_if< has_generate_template<Seq>::value >::type* = 0 )
+    {
+        typename Seq::result_type w[8];
+        q.generate(&w[0], &w[8]);
+
+        for (unsigned short i=0; i<4; ++i) {
+            _k[i] = ( static_cast<uint64_t>(w[2*i]) << 32) | w[2*i+1];
+            _s[i] = 0;
+        }
+        _o_counter = 0;
+        
+        encrypt_counter();
+    }
+    
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 8
+    // Advances e’s state ei to ei+1 = TA(ei) and returns GA(ei).
+    uint32_t operator()()
+    {
+        // can we return a value from the current block?
+        if (_o_counter < 8) {
+            unsigned short _o_index = _o_counter >> 1;
+            _o_counter++;
+            if (_o_counter&1) 
+                return _o[_o_index] & 0xFFFFFFFF; 
+            else
+                return _o[_o_index] >> 32;
+        }
+        
+        // generate a new block and return the first 32 bits
+        inc_counter();
+        encrypt_counter();
+        _o_counter = 1; // the next call
+        return _o[0] & 0xFFFFFFFF;   // this call
+    }
+    
+    // -------------------------------------------------
+    // misc
+    // -------------------------------------------------
+    
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 9
+    // Advances e’s state ei to ei+z by any means equivalent to z
+    // consecutive calls e().
+    void discard(uint64_t z)
+    {
+        // check if we stay in the current block
+        if (z < 8 - _o_counter) {
+            _o_counter += static_cast<unsigned short>(z);
+            return;
+        }
+
+        // we will have to generate a new block...
+        z -= (8 - _o_counter);  // discard the remainder of the current blok
+        _o_counter = z % 8;     // set the pointer in the correct element in the new block
+        z -= _o_counter;        // update z
+        z >>= 3;                // the number of buffers is elements/8
+        ++z;                    // and one more because we crossed the buffer line
+        inc_counter(z);
+        encrypt_counter();
+    }
+
+   // -------------------------------------------------
+   // IO
+   // -------------------------------------------------
+   template<class CharT, class Traits>
+   friend std::basic_ostream<CharT,Traits>&
+   operator<<(std::basic_ostream<CharT,Traits>& os, const prng_engine& s) {
+       for (unsigned short i=0; i<4; ++i)
+           os << s._k[i] << ' ' << s._s[i] << ' ' << s._o[i] << ' ';
+       os << s._o_counter;
+       return os;
+   }
+
+   template<class CharT, class Traits>
+   friend std::basic_istream<CharT,Traits>&
+   operator>>(std::basic_istream<CharT,Traits>& is, prng_engine& s) {
+       for (unsigned short i=0; i<4; ++i)
+           is >> s._k[i] >> s._s[i] >> s._o[i];
+       is >> s._o_counter;
+       return is;
+   } 
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 10
+    // This operator is an equivalence relation. With Sx and Sy as the infinite 
+    // sequences of values that would be generated by repeated future calls to 
+    // x() and y(), respectively, returns true if Sx = Sy; else returns false.
+    bool operator==(const prng_engine& y) 
+    {
+        if (_o_counter != y._o_counter) return false;
+        for (unsigned short i=0; i<4; ++i) {
+            if (_s[i] != y._s[i]) return false;
+            if (_k[i] != y._k[i]) return false;
+            if (_o[i] != y._o[i]) return false;
+        }
+        return true;
+    }
+    
+    // req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 11
+    bool operator!=(const prng_engine& y) 
+    {
+        return !(*this == y);
+    }
+    
+    // Extra function to set the key
+    void set_key(uint64_t k0=0, uint64_t k1=0, uint64_t k2=0, uint64_t k3=0)
+    {
+        _k[0] = k0; _k[1] = k1; _k[2] = k2; _k[3] = k3;
+        encrypt_counter();
+    }
+    
+    // set the counter
+    void set_counter(uint64_t s0=0, uint64_t s1=0, uint64_t s2=0, uint64_t s3=0, unsigned short o_counter=0)
+    {
+        _s[0] = s0; 
+        _s[1] = s1; 
+        _s[2] = s2; 
+        _s[3] = s3;
+        _o_counter = o_counter % 8;
+        encrypt_counter();
+    }
+    
+    
+    // versioning
+    uint32_t version()
+    {
+        return 5;
+    }
+    
+private:
+    void encrypt_counter()
+    {
+        uint64_t b[4];
+        uint64_t k[5];
+
+        for (unsigned short i=0; i<4; ++i) b[i] = _s[i];
+        for (unsigned short i=0; i<4; ++i) k[i] = _k[i];
+
+        k[4] = 0x1BD11BDAA9FC1A22 ^ k[0] ^ k[1] ^ k[2] ^ k[3];
+
+        MIXK(b[0], b[1], 14,   b[2], b[3], 16,   k[0], k[1], k[2], k[3]);
+        MIX2(b[0], b[3], 52,   b[2], b[1], 57);
+        MIX2(b[0], b[1], 23,   b[2], b[3], 40);
+        MIX2(b[0], b[3],  5,   b[2], b[1], 37);
+        MIXK(b[0], b[1], 25,   b[2], b[3], 33,   k[1], k[2], k[3], k[4]+1);
+        MIX2(b[0], b[3], 46,   b[2], b[1], 12);
+        MIX2(b[0], b[1], 58,   b[2], b[3], 22);
+        MIX2(b[0], b[3], 32,   b[2], b[1], 32);
+
+        MIXK(b[0], b[1], 14,   b[2], b[3], 16,   k[2], k[3], k[4], k[0]+2);
+        MIX2(b[0], b[3], 52,   b[2], b[1], 57);
+        MIX2(b[0], b[1], 23,   b[2], b[3], 40);
+        MIX2(b[0], b[3],  5,   b[2], b[1], 37);
+        MIXK(b[0], b[1], 25,   b[2], b[3], 33,   k[3], k[4], k[0], k[1]+3);
+
+        MIX2(b[0], b[3], 46,   b[2], b[1], 12);
+        MIX2(b[0], b[1], 58,   b[2], b[3], 22);
+        MIX2(b[0], b[3], 32,   b[2], b[1], 32);
+
+        MIXK(b[0], b[1], 14,   b[2], b[3], 16,   k[4], k[0], k[1], k[2]+4);
+        MIX2(b[0], b[3], 52,   b[2], b[1], 57);
+        MIX2(b[0], b[1], 23,   b[2], b[3], 40);
+        MIX2(b[0], b[3],  5,   b[2], b[1], 37);
+
+        for (unsigned int i=0; i<4; ++i) _o[i] = b[i] + k[i];
+        _o[3] += 5;
+    }
+    
+    void inc_counter()
+    {
+        ++_s[0]; 
+        if (_s[0] != 0) return;
+        
+        ++_s[1]; 
+        if (_s[1] != 0) return;
+        
+        ++_s[2]; 
+        if (_s[2] != 0) return;
+        
+        ++_s[3];
+    }
+    
+    void inc_counter(uint64_t z)
+    {
+        if (z > 0xFFFFFFFFFFFFFFFF - _s[0]) {   // check if we will overflow the first 64 bit int
+            ++_s[1];
+            if (_s[1] == 0) {
+                ++_s[2];
+                if (_s[2] == 0) {
+                    ++_s[3];
+                }
+            }
+        }
+        _s[0] += z;
+    }
+    
+private:
+    uint64_t _k[4];             // key
+    uint64_t _s[4];             // state (counter)
+    uint64_t _o[4];             // cipher output    4 * 64 bit = 256 bit output
+    unsigned short _o_counter;  // output chunk counter, the 256 random bits in _o 
+                                // are returned in eight 32 bit chunks
+};
+
+
+} // namespace sitmo
+
+#undef MIXK
+#undef MIX2
+
+#endif

scripts/filelist

@@ -4,13 +4,16 @@ home=`pwd`
 
 # library Make.inc
 cd $home/lib
-HFILES=`find . -type f -name '*.h' -not -path '*/Old/*' -not -path '*/Eigen/*'`
+HFILES=`find . -type f -name '*.h' -not -name '*Hdf5*' -not -path '*/gamma-gen/*' -not -path '*/Old/*' -not -path '*/Eigen/*'`
 HFILES="$HFILES"
-CCFILES=`find . -type f -name '*.cc' -not  -name '*ommunicator*.cc'`
+CCFILES=`find . -type f -name '*.cc' -not -path '*/gamma-gen/*' -not -name '*Communicator*.cc' -not -name '*Hdf5*'`
-echo HFILES=$HFILES > Make.inc
+HPPFILES=`find . -type f -name '*.hpp'`
+echo HFILES=$HFILES $HPPFILES > Make.inc
 echo >> Make.inc
 echo CCFILES=$CCFILES >> Make.inc
 
+
+
 # tests Make.inc
 cd $home/tests
 dirs=`find . -type d -not -path '*/\.*'`

tests/IO/Test_serialisation.cc

@@ -28,130 +28,152 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     /*  END LEGAL */
 #include <Grid/Grid.h>
 
-namespace Grid {
-  
-  GRID_SERIALIZABLE_ENUM(myenum, undef, red, 1, blue, 2, green, 3);
-    
-  class myclass: Serializable {
-  public:
-    
-    GRID_SERIALIZABLE_CLASS_MEMBERS(myclass,
-                            myenum, e,
-                            std::vector<myenum>, ve,
-                            std::string, name,
-                            int, x,
-                            double, y,
-                            bool , b,
-                            std::vector<double>, array,
-                            std::vector<std::vector<double>>, twodimarray,
-                            );
-    
-    myclass() {}
-    myclass(int i)
-    : array(4,5.1), twodimarray(3,std::vector<double>(2,1.23456)), ve(2, myenum::blue)
-    {
-      e=myenum::red;
-      x=i;
-      y=2*i;
-      b=true;
-      name="bother said pooh";
-    }
-  };
-  
-}
-
 using namespace Grid;
 
-int16_t i16 = 1;
+GRID_SERIALIZABLE_ENUM(myenum, undef, red, 1, blue, 2, green, 3);
+  
+class myclass: Serializable {
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(myclass,
+                          myenum, e,
+                          std::vector<myenum>, ve,
+                          std::string, name,
+                          int, x,
+                          double, y,
+                          bool , b,
+                          std::vector<double>, array,
+                          std::vector<std::vector<double>>, twodimarray,
+                          std::vector<std::vector<std::vector<Complex>>>, cmplx3darray
+                          );
+  myclass() {}
+  myclass(int i)
+  : array(4,5.1)
+  , twodimarray(3,std::vector<double>(5, 1.23456))
+  , cmplx3darray(3,std::vector<std::vector<Complex>>(5, std::vector<Complex>(7, Complex(1.2, 3.4))))
+  , ve(2, myenum::blue)
+  {
+    e=myenum::red;
+    x=i;
+    y=2*i;
+    b=true;
+    name="bother said pooh";
+  }
+};
+
+int16_t  i16 = 1;
 uint16_t u16 = 2;
-int32_t i32 = 3;
+int32_t  i32 = 3;
 uint32_t u32 = 4;
-int64_t i64 = 5;
+int64_t  i64 = 5;
 uint64_t u64 = 6;
-float    f = M_PI;
+float    f   = M_PI;
-double   d = 2*M_PI;
+double   d   = 2*M_PI;
-bool     b = false;
+bool     b   = false;
+
+template <typename W, typename R, typename O>
+void ioTest(const std::string &filename, const O &object, const std::string &name)
+{
+  // writer needs to be destroyed so that writing physically happens
+  {
+    W writer(filename);
+    
+    write(writer, "testobject", object);
+  }
+  
+  R    reader(filename);
+  O    buf;
+  bool good;
+  
+  read(reader, "testobject", buf);
+  good = (object == buf);
+  std::cout << name << " IO test: " << (good ? "success" : "failure");
+  std::cout << std::endl;
+  if (!good) exit(EXIT_FAILURE);
+}
 
 int main(int argc,char **argv)
 {
-  {
+  std::cout << "==== basic IO" << std::endl;
-    XmlWriter WR("bother.xml");
+  XmlWriter WR("bother.xml");
-    
+  
-    // test basic type writing
+  // test basic type writing
-    push(WR,"BasicTypes");
+  std::cout << "-- basic writing to 'bother.xml'..." << std::endl;
-    write(WR,std::string("i16"),i16);
+  push(WR,"BasicTypes");
-    write(WR,"u16",u16);
+  write(WR,std::string("i16"),i16);
-    write(WR,"i32",i32);
+  write(WR,"u16",u16);
-    write(WR,"u32",u32);
+  write(WR,"i32",i32);
-    write(WR,"i64",i64);
+  write(WR,"u32",u32);
-    write(WR,"u64",u64);
+  write(WR,"i64",i64);
-    write(WR,"f",f);
+  write(WR,"u64",u64);
-    write(WR,"d",d);
+  write(WR,"f",f);
-    write(WR,"b",b);
+  write(WR,"d",d);
-    pop(WR);
+  write(WR,"b",b);
-    
+  pop(WR);
-    // test serializable class writing
+  
-    myclass obj(1234); // non-trivial constructor
+  // test serializable class writing
-    write(WR,"obj",obj);
+  myclass              obj(1234); // non-trivial constructor
-    WR.write("obj2", obj);
+  std::vector<myclass> vec;
-    std::cout << obj << std::endl;
+  
-    
+  std::cout << "-- serialisable class writing to 'bother.xml'..." << std::endl;
-    std::vector<myclass> vec;
+  write(WR,"obj",obj);
-    vec.push_back(myclass(1234));
+  WR.write("obj2", obj);
-    vec.push_back(myclass(5678));
+  vec.push_back(myclass(1234));
-    vec.push_back(myclass(3838));
+  vec.push_back(myclass(5678));
-    write(WR, "objvec", vec);
+  vec.push_back(myclass(3838));
-  };
+  write(WR, "objvec", vec);
+  std::cout << "-- serialisable class writing to std::cout:" << std::endl;
+  std::cout << obj << std::endl;
+  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;
   
   // read tests
-  myclass copy1, copy2, copy3;
+  std::cout << "\n==== IO self-consistency tests" << std::endl;
-  std::vector<myclass> veccopy1, veccopy2, veccopy3;
   //// XML
-  {
+  ioTest<XmlWriter, XmlReader>("iotest.xml", obj, "XML    (object)           ");
-    XmlReader RD("bother.xml");
+  ioTest<XmlWriter, XmlReader>("iotest.xml", vec, "XML    (vector of objects)");
-    read(RD,"obj",copy1);
-    read(RD,"objvec", veccopy1);
-    std::cout << "Loaded (XML) -----------------" << std::endl;
-    std::cout << copy1 << std::endl << veccopy1 << std::endl;
-  }
   //// binary
-  {
+  ioTest<BinaryWriter, BinaryReader>("iotest.bin", obj, "binary (object)           ");
-    BinaryWriter BWR("bother.bin");
+  ioTest<BinaryWriter, BinaryReader>("iotest.bin", vec, "binary (vector of objects)");
-    write(BWR,"discard",copy1 );
-    write(BWR,"discard",veccopy1 );
-  }
-  {
-    BinaryReader BRD("bother.bin");
-    read (BRD,"discard",copy2 );
-    read (BRD,"discard",veccopy2 );
-    std::cout << "Loaded (bin) -----------------" << std::endl;
-    std::cout << copy2 << std::endl << veccopy2 << std::endl;
-  }
   //// text
-  {
+  ioTest<TextWriter, TextReader>("iotest.dat", obj, "text   (object)           ");
-    TextWriter TWR("bother.txt");
+  ioTest<TextWriter, TextReader>("iotest.dat", vec, "text   (vector of objects)");
-    write(TWR,"discard",copy1 );
+  //// HDF5
-    write(TWR,"discard",veccopy1 );
+#ifdef HAVE_HDF5
-  }
+  ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", obj, "HDF5   (object)           ");
-  {
+  ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", vec, "HDF5   (vector of objects)");
-    TextReader TRD("bother.txt");
+#endif
-    read (TRD,"discard",copy3 );
-    read (TRD,"discard",veccopy3 );
-    std::cout << "Loaded (txt) -----------------" << std::endl;
-    std::cout << copy3 << std::endl << veccopy3 << std::endl;
-  }
   
-  std::vector<int> iv = strToVec<int>("1 2 2 4");
+  std::cout << "\n==== vector flattening/reconstruction" << std::endl;
-  std::vector<std::string> sv = strToVec<std::string>("bli bla blu");
+  typedef std::vector<std::vector<std::vector<double>>> vec3d;
   
-  for (auto &e: iv)
+  vec3d dv, buf;
+  double d = 0.;
+  
+  dv.resize(4);
+  for (auto &v1: dv)
   {
-    std::cout << e << " ";
+    v1.resize(3);
+    for (auto &v2: v1)
+    {
+      v2.resize(5);
+      for (auto &x: v2)
+      {
+        x = d++;
+      }
+    }
   }
-  std::cout << std::endl;
+  std::cout << "original 3D vector:" << std::endl;
-  for (auto &e: sv)
+  std::cout << dv << std::endl;
-  {
+  
-    std::cout << e << " ";
+  Flatten<vec3d> flatdv(dv);
-  }
+  
-  std::cout << std::endl;
+  std::cout << "\ndimensions:" << std::endl;
+  std::cout << flatdv.getDim() << std::endl;
+  std::cout << "\nflattened vector:" << std::endl;
+  std::cout << flatdv.getFlatVector() << std::endl;
+  
+  Reconstruct<vec3d> rec(flatdv.getFlatVector(), flatdv.getDim());
+  std::cout << "\nreconstructed vector:" << std::endl;
+  std::cout << flatdv.getVector() << std::endl;
 }

tests/Test_cayley_even_odd_vec.cc

@@ -36,12 +36,12 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+    Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+        Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+        Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+        Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+        Gamma::Algebra::GammaT
-  };
+    };
 
 typedef DomainWallFermion<DomainWallVec5dImplR>                      DomainWallVecFermionR;
 typedef ZMobiusFermion<ZDomainWallVec5dImplR>                        ZMobiusVecFermionR;
@@ -340,7 +340,7 @@ void  TestMoo(This & Dw, That &sDw)
   LatticeFermion ndiff(ngrid);
   LatticeFermion sdiff(sgrid);
 
-  Gamma g5( Gamma::Gamma5 );
+    Gamma g5( Gamma::Algebra::Gamma5 );
 
   std::vector<int> seeds({1,2,3,4,5,7,8});
   GridParallelRNG    RNG5(ngrid);  

tests/Test_dwf_mixedcg_prec.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/core/Test_cf_coarsen_support.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/core/Test_checker.cc

@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int toint(const char* str){

tests/core/Test_contfrac_even_odd.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 

tests/core/Test_dwf_even_odd.cc

@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 

tests/core/Test_dwf_rb5d.cc

@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 typedef WilsonFermion5D<DomainWallVec5dImplR> WilsonFermion5DR;

tests/core/Test_fft.cc

@@ -177,7 +177,7 @@ int main (int argc, char ** argv)
     const int sdir=0;
     RealD mass=0.01;
     RealD M5  =1.0;
-    Gamma G5(Gamma::Gamma5);
+    Gamma G5(Gamma::Algebra::Gamma5);
 
     GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID);
     GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID);
@@ -218,12 +218,12 @@ int main (int argc, char ** argv)
     /////////////////////////////////////////////////////////////////
     // work out the predicted from Fourier
     /////////////////////////////////////////////////////////////////
-    Gamma::GammaMatrix Gmu [] = {
+    Gamma::Algebra Gmu [] = {
-      Gamma::GammaX,
+      Gamma::Algebra::GammaX,
-      Gamma::GammaY,
+      Gamma::Algebra::GammaY,
-      Gamma::GammaZ,
+      Gamma::Algebra::GammaZ,
-      Gamma::GammaT,
+      Gamma::Algebra::GammaT,
-      Gamma::Gamma5
+      Gamma::Algebra::Gamma5
     };
     LatticeFermionD    Kinetic(FGrid); Kinetic = zero;
     LatticeComplexD    kmu(FGrid); 
@@ -311,7 +311,7 @@ int main (int argc, char ** argv)
     std::cout << " Solving by FFT and Feynman rules" <<std::endl;
     Ddwf.FreePropagator(src,ref,mass) ;
 
-    Gamma G5(Gamma::Gamma5);
+    Gamma G5(Gamma::Algebra::Gamma5);
 
     LatticeFermionD    src5(FGrid); src5=zero;
     LatticeFermionD    tmp5(FGrid); 
@@ -391,7 +391,7 @@ int main (int argc, char ** argv)
     std::cout << " Solving by FFT and Feynman rules" <<std::endl;
     Dov.FreePropagator(src,ref,mass) ;
 
-    Gamma G5(Gamma::Gamma5);
+    Gamma G5(Gamma::Algebra::Gamma5);
 
     LatticeFermionD    src5(FGrid); src5=zero;
     LatticeFermionD    tmp5(FGrid); 

tests/core/Test_gamma.cc

@@ -1,223 +1,283 @@
-    /*************************************************************************************
+/*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+Grid physics library, www.github.com/paboyle/Grid 
 
-    Source file: ./tests/Test_gamma.cc
+Source file: ./tests/Test_gamma.cc
 
-    Copyright (C) 2015
+Copyright (C) 2015-2017
 
 Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Antonin Portelli <antonin.portelli@ed.ac.uk>
 
-    This program is free software; you can redistribute it and/or modify
+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
+it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
+the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
+(at your option) any later version.
 
-    This program is distributed in the hope that it will be useful,
+This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
+but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+GNU General Public License for more details.
 
-    You should have received a copy of the GNU General Public License along
+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.,
+with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
+See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
+*************************************************************************************/
     /*  END LEGAL */
 #include <Grid/Grid.h>
 
 using namespace std;
 using namespace Grid;
-using namespace Grid::QCD;
+using namespace QCD;
 
-//template<class vobj> class is_pod< iScalar<vobj> >
+static constexpr double                      tolerance = 1.0e-6;
-//{
+static std::array<SpinMatrix, Gamma::nGamma> testAlgebra;
-//
-//};
 
-int main (int argc, char ** argv)
+void print(const SpinMatrix &g)
+{
+  for(int i = 0; i < Ns; i++)
+  {
+    std::cout << GridLogMessage << "(";
+    for(int j=0;j<Ns;j++){
+      if ( abs(g()(i,j)()) == 0 ) {
+        std::cout<< " 0";
+      } else if ( abs(g()(i,j)() - Complex(0,1)) == 0){
+        std::cout<< " i";
+      } else if ( abs(g()(i,j)() + Complex(0,1)) == 0){
+        std::cout<< "-i";
+      } else if ( abs(g()(i,j)() - Complex(1,0)) == 0){
+        std::cout<< " 1";
+      } else if ( abs(g()(i,j)() + Complex(1,0)) == 0){
+        std::cout<< "-1";
+      }
+      std::cout<<((j == Ns-1) ? ")" : "," );
+    }
+    std::cout << std::endl;
+  }
+  std::cout << GridLogMessage << std::endl;
+}
+
+void createTestAlgebra(void)
+{
+  std::array<SpinMatrix, 4> testg;
+  SpinMatrix                testg5;
+  const Complex             I(0., 1.), mI(0., -1.);
+  
+  testg[0] = zero;
+  testg[0]()(0, 3) = I;
+  testg[0]()(1, 2) = I;
+  testg[0]()(2, 1) = mI;
+  testg[0]()(3, 0) = mI;
+  std::cout << GridLogMessage << "test GammaX= " << std::endl;
+  print(testg[0]);
+  testg[1] = zero;
+  testg[1]()(0, 3) = -1.;
+  testg[1]()(1, 2) = 1.;
+  testg[1]()(2, 1) = 1.;
+  testg[1]()(3, 0) = -1.;
+  std::cout << GridLogMessage << "test GammaY= " << std::endl;
+  print(testg[1]);
+  testg[2] = zero;
+  testg[2]()(0, 2) = I;
+  testg[2]()(1, 3) = mI;
+  testg[2]()(2, 0) = mI;
+  testg[2]()(3, 1) = I;
+  std::cout << GridLogMessage << "test GammaZ= " << std::endl;
+  print(testg[2]);
+  testg[3] = zero;
+  testg[3]()(0, 2) = 1.;
+  testg[3]()(1, 3) = 1.;
+  testg[3]()(2, 0) = 1.;
+  testg[3]()(3, 1) = 1.;
+  std::cout << GridLogMessage << "test GammaT= " << std::endl;
+  print(testg[3]);
+  testg5 = testg[0]*testg[1]*testg[2]*testg[3];
+  
+#define DEFINE_TEST_G(g, exp)\
+testAlgebra[Gamma::Algebra::g]        = exp;\
+testAlgebra[Gamma::Algebra::Minus##g] = -exp;\
+  
+  DEFINE_TEST_G(Identity    , 1.);
+  DEFINE_TEST_G(Gamma5      , testg5);
+  DEFINE_TEST_G(GammaX      , testg[0]);
+  DEFINE_TEST_G(GammaY      , testg[1]);
+  DEFINE_TEST_G(GammaZ      , testg[2]);
+  DEFINE_TEST_G(GammaT      , testg[3]);
+  DEFINE_TEST_G(GammaXGamma5, testg[0]*testg5);
+  DEFINE_TEST_G(GammaYGamma5, testg[1]*testg5);
+  DEFINE_TEST_G(GammaZGamma5, testg[2]*testg5);
+  DEFINE_TEST_G(GammaTGamma5, testg[3]*testg5);
+  DEFINE_TEST_G(SigmaXY     , .5*(testg[0]*testg[1] - testg[1]*testg[0]));
+  DEFINE_TEST_G(SigmaXZ     , .5*(testg[0]*testg[2] - testg[2]*testg[0]));
+  DEFINE_TEST_G(SigmaXT     , .5*(testg[0]*testg[3] - testg[3]*testg[0]));
+  DEFINE_TEST_G(SigmaYZ     , .5*(testg[1]*testg[2] - testg[2]*testg[1]));
+  DEFINE_TEST_G(SigmaYT     , .5*(testg[1]*testg[3] - testg[3]*testg[1]));
+  DEFINE_TEST_G(SigmaZT     , .5*(testg[2]*testg[3] - testg[3]*testg[2]));
+  
+#undef DEFINE_TEST_G
+}
+
+template <typename Expr>
+void test(const Expr &a, const Expr &b)
+{
+  if (norm2(a - b) < tolerance)
+  {
+    std::cout << "[OK] ";
+  }
+  else
+  {
+    std::cout << "[fail]" << std::endl;
+    std::cout << GridLogError << "a= " << a << std::endl;
+    std::cout << GridLogError << "is different (tolerance= " << tolerance << ") from " << std::endl;
+    std::cout << GridLogError << "b= " << b << std::endl;
+    exit(EXIT_FAILURE);
+  }
+}
+
+void checkGamma(const Gamma::Algebra a, GridSerialRNG &rng)
+{
+  SpinVector v;
+  SpinMatrix m, &testg = testAlgebra[a];
+  Gamma      g(a);
+  bool       pass = true;
+  
+  random(rng, v);
+  random(rng, m);
+  
+  std::cout << GridLogMessage << "Checking " << Gamma::name[a] << ": ";
+  std::cout << "vecmul ";
+  test(g*v, testg*v);
+  std::cout << "matlmul ";
+  test(g*m, testg*m);
+  std::cout << "matrmul ";
+  test(m*g, m*testg);
+  std::cout << std::endl;
+}
+
+void checkProd(const Gamma::Algebra a, const Gamma::Algebra b)
+{
+  SpinMatrix gm, testg = testAlgebra[a]*testAlgebra[b];
+  Gamma      g = Gamma(a)*Gamma(b);
+  bool       pass = true;
+  
+  std::cout << GridLogMessage << "Checking " << Gamma::name[a] << " * "
+            << Gamma::name[b] << ": ";
+  gm = 1.0;
+  gm = g*gm;
+  test(gm, testg);
+  std::cout << "(= " << Gamma::name[g.g] << ")" << std::endl;
+}
+
+void checkAdj(const Gamma::Algebra a)
+{
+  SpinMatrix gm, testg = adj(testAlgebra[a]);
+  Gamma      g(adj(Gamma(a)));
+  bool       pass = true;
+  
+  std::cout << GridLogMessage << "Checking adj(" << Gamma::name[a] << "): ";
+  gm = 1.0;
+  gm = g*gm;
+  test(gm, testg);
+  std::cout << "(= " << Gamma::name[g.g] << ")" << std::endl;
+}
+
+void checkProject(GridSerialRNG &rng)
+{
+  SpinVector     rv, recon, full;
+  HalfSpinVector hsp, hsm;
+  
+  random(rng, rv);
+  
+#define CHECK_PROJ(dir, gamma)\
+std::cout << GridLogMessage << "Checking " << #dir << " projector: ";\
+spProj##dir(hsm,rv);\
+spRecon##dir(recon,hsm);\
+test(recon, rv + Gamma(Gamma::Algebra::gamma)*rv);\
+std::cout << std::endl;
+  
+  CHECK_PROJ(Xp, GammaX);
+  CHECK_PROJ(Yp, GammaY);
+  CHECK_PROJ(Zp, GammaZ);
+  CHECK_PROJ(Tp, GammaT);
+  CHECK_PROJ(5p, Gamma5);
+  CHECK_PROJ(Xm, MinusGammaX);
+  CHECK_PROJ(Ym, MinusGammaY);
+  CHECK_PROJ(Zm, MinusGammaZ);
+  CHECK_PROJ(Tm, MinusGammaT);
+  CHECK_PROJ(5m, MinusGamma5);
+  
+#undef CHECK_PROJ
+}
+
+void checkGammaL(const Gamma::Algebra a, GridSerialRNG &rng)
+{
+  SpinVector v;
+  SpinMatrix m, &testg = testAlgebra[a], pl;
+  GammaL     gl(a);
+  bool       pass = true;
+  
+  random(rng, v);
+  random(rng, m);
+  
+  pl = testAlgebra[Gamma::Algebra::Identity]
+       - testAlgebra[Gamma::Algebra::Gamma5];
+  std::cout << GridLogMessage << "Checking left-projected " << Gamma::name[a] << ": ";
+  std::cout << "vecmul ";
+  test(gl*v, testg*pl*v);
+  std::cout << "matlmul ";
+  test(gl*m, testg*pl*m);
+  std::cout << "matrmul ";
+  test(m*gl, m*testg*pl);
+  std::cout << std::endl;
+}
+
+int main(int argc, char *argv[])
 {
   Grid_init(&argc,&argv);
-
+  
   std::vector<int> latt_size   = GridDefaultLatt();
   std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
-     
+  
-  GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
+  GridCartesian Grid(latt_size,simd_layout,mpi_layout);
-
+  GridSerialRNG sRNG;
-  GridParallelRNG          pRNG(&Grid);
+  
-  pRNG.SeedRandomDevice();
-
-  GridSerialRNG            sRNG;
   sRNG.SeedRandomDevice();
-
+  
-  SpinMatrix ident; ident=zero;
+  std::cout << GridLogMessage << "======== Test algebra" << std::endl;
-  SpinMatrix rnd  ; random(sRNG,rnd);
+  createTestAlgebra();
-
+  std::cout << GridLogMessage << "======== Multiplication operators check" << std::endl;
-  SpinMatrix ll; ll=zero;
+  for (int i = 0; i < Gamma::nGamma; ++i)
-  SpinMatrix rr; rr=zero;
+  {
-  SpinMatrix result;
+    checkGamma(i, sRNG);
-
-  SpinVector lv; random(sRNG,lv);
-  SpinVector rv; random(sRNG,rv);
-
-  //  std::cout<<GridLogMessage << " Is pod " << std::is_pod<SpinVector>::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod double   " << std::is_pod<double>::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod ComplexF " << std::is_pod<ComplexF>::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is triv double " << std::has_trivial_default_constructor<double>::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is triv ComplexF " << std::has_trivial_default_constructor<ComplexF>::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<double> " << std::is_pod<iScalar<double> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<ComplexF> " << std::is_pod<iScalar<ComplexF> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<vComplexF> " << std::is_pod<iScalar<vComplexF> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<vComplexD> " << std::is_pod<iScalar<vComplexD> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<vRealF> " << std::is_pod<iScalar<vRealF> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<vRealD> " << std::is_pod<iScalar<vRealD> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is triv Scalar<double> " <<std::has_trivial_default_constructor<iScalar<double> >::value << std::endl;
-  //  std::cout<<GridLogMessage << " Is triv Scalar<vComplexD> "<<std::has_trivial_default_constructor<iScalar<vComplexD> >::value  << std::endl;
-
-  for(int a=0;a<Ns;a++){
-    ident()(a,a) = ComplexF(1.0);
   }
-
+  std::cout << GridLogMessage << std::endl;
-  const Gamma::GammaMatrix *g = Gamma::GammaMatrices;
+  std::cout << GridLogMessage << "======== Algebra multiplication table check" << std::endl;
-  const char **list           = Gamma::GammaMatrixNames;
+  for (int i = 0; i < Gamma::nGamma; ++i)
-
+  for (int j = 0; j < Gamma::nGamma; ++j)
-  result =ll*Gamma(g[0])*rr;
+  {
-  result =ll*Gamma(g[0]);
+    checkProd(i, j);
-  rv = Gamma(g[0])*lv;
-
-  for(int mu=0;mu<12;mu++){
-
-    result = Gamma(g[mu])* ident;
-
-    for(int i=0;i<Ns;i++){
-
-      if(i==0) std::cout<<GridLogMessage << list[mu];
-      else     std::cout<<GridLogMessage << list[12];
-
-      std::cout<<"(";
-      for(int j=0;j<Ns;j++){
-	if ( abs(result()(i,j)())==0 ) {
-	  std::cout<< " 0";
-	} else if ( abs(result()(i,j)() - Complex(0,1))==0){
-	  std::cout<< " i";
-	} else if ( abs(result()(i,j)() + Complex(0,1))==0){
-	  std::cout<< "-i";
-	} else if ( abs(result()(i,j)() - Complex(1,0))==0){
-	  std::cout<< " 1";
-	} else if ( abs(result()(i,j)() + Complex(1,0))==0){
-	  std::cout<< "-1";
-	}
-	std::cout<<((j==Ns-1) ? ")" : "," );
-      }
-      std::cout << std::endl;
-    }
-
-    std::cout << std::endl;
-
   }
-
+  std::cout << GridLogMessage << std::endl;
-  std::cout << "Testing Gamma^2 - 1 = 0"<<std::endl;
+  std::cout << GridLogMessage << "======== Adjoints check" << std::endl;
-  for(int mu=0;mu<6;mu++){
+  for (int i = 0; i < Gamma::nGamma; ++i)
-    result =  Gamma(g[mu])* ident * Gamma(g[mu]);
+  {
-    result = result - ident;
+    checkAdj(i);
-    RealD mag = norm2(result);
-    std::cout << list[mu]<<" " << mag<<std::endl;
   }
-
+  std::cout << GridLogMessage << std::endl;
-  std::cout << "Testing (MinusGamma + G )M = 0"<<std::endl;
+  std::cout << GridLogMessage << "======== Spin projectors check" << std::endl;
-  for(int mu=0;mu<6;mu++){
+  checkProject(sRNG);
-    result =          rnd * Gamma(g[mu]);
+  std::cout << GridLogMessage << std::endl;
-    result = result + rnd * Gamma(g[mu+6]);
+  std::cout << GridLogMessage << "======== Gamma-left matrices check" << std::endl;
-    RealD mag = norm2(result);
+  for (int i = 0; i < Gamma::nGamma; ++i)
-    std::cout << list[mu]<<" " << mag<<std::endl;
+  {
+    checkGammaL(i, sRNG);
   }
-
+  
-  std::cout << "Testing M(MinusGamma + G )  = 0"<<std::endl;
-  for(int mu=0;mu<6;mu++){
-    result =           Gamma(g[mu])  *rnd;
-    result = result +  Gamma(g[mu+6])*rnd;
-    RealD mag = norm2(result);
-    std::cout << list[mu]<<" " << mag<<std::endl;
-  }
-
-  // Testing spins and reconstructs
-  SpinVector     recon; random(sRNG,rv);
-  SpinVector      full;
-  HalfSpinVector hsp,hsm;
-
-  // Xp
-  double mag;
-  spProjXp(hsm,rv);
-  spReconXp(recon,hsm);
-  full = rv + Gamma(Gamma::GammaX) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Xp "<< mag<<std::endl;
-
-  // Xm
-  spProjXm(hsm,rv);
-  spReconXm(recon,hsm);
-  full = rv - Gamma(Gamma::GammaX) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Xm "<< mag<<std::endl;
-
-  // Yp
-  spProjYp(hsm,rv);
-  spReconYp(recon,hsm);
-  full = rv + Gamma(Gamma::GammaY) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Yp "<< mag<<std::endl;
-
-  // Ym
-  spProjYm(hsm,rv);
-  spReconYm(recon,hsm);
-  full = rv - Gamma(Gamma::GammaY) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Ym "<< mag<<std::endl;
-
-  // Zp
-  spProjZp(hsm,rv);
-  spReconZp(recon,hsm);
-  full = rv + Gamma(Gamma::GammaZ) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Zp "<< mag<<std::endl;
-
-  // Zm
-  spProjZm(hsm,rv);
-  spReconZm(recon,hsm);
-  full = rv - Gamma(Gamma::GammaZ) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Zm "<< mag<<std::endl;
-
-  // Tp
-  spProjTp(hsm,rv);
-  spReconTp(recon,hsm);
-  full = rv + Gamma(Gamma::GammaT) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Tp "<< mag<<std::endl;
-
-  // Tm
-  spProjTm(hsm,rv);
-  spReconTm(recon,hsm);
-  full = rv - Gamma(Gamma::GammaT) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Tm "<< mag<<std::endl;
-
-  // 5p
-  spProj5p(hsm,rv);
-  spRecon5p(recon,hsm);
-  full = rv + Gamma(Gamma::Gamma5) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "5p "<< mag<<std::endl;
-
-  // 5m
-  spProj5m(hsm,rv);
-  spRecon5m(recon,hsm);
-  full = rv - Gamma(Gamma::Gamma5) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "5m "<< mag<<std::endl;
-
   Grid_finalize();
+  
+  return EXIT_SUCCESS;
 }

tests/core/Test_wilson_even_odd.cc

@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/core/Test_wilson_tm_even_odd.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/debug/Test_cayley_cg.cc

@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 template<class What> 

tests/debug/Test_cayley_coarsen_support.cc

@@ -38,11 +38,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/debug/Test_cayley_even_odd.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 

tests/hadrons/Makefile.am

@@ -1,3 +1,3 @@
-AM_LDFLAGS += -L../../extras/hadrons
+AM_LDFLAGS += -L../../extras/Hadrons
 
 include Make.inc

tests/hadrons/Test_hadrons_meson_3pt.cc

@@ -30,6 +30,14 @@
 using namespace Grid;
 using namespace Hadrons;
 
+static Gamma::Algebra gmu[4] =
+{
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
+};
+
 int main(int argc, char *argv[])
 {
     // initialization //////////////////////////////////////////////////////////
@@ -102,7 +110,7 @@ int main(int argc, char *argv[])
             seqName.push_back(std::vector<std::string>(Nd));
             for (unsigned int mu = 0; mu < Nd; ++mu)
             {
-                seqPar.gamma   = 0x1 << mu;
+                seqPar.gamma   = gmu[mu];
                 seqName[i][mu] = "G" + std::to_string(seqPar.gamma)
                                  + "_" + std::to_string(seqPar.tA) + "-"
                                  + qName[i];
@@ -127,9 +135,11 @@ int main(int argc, char *argv[])
         for (unsigned int i = 0; i < flavour.size(); ++i)
         for (unsigned int j = i; j < flavour.size(); ++j)
         {
-            mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
+            mesPar.output      = "mesons/Z2_" + flavour[i] + flavour[j];
-            mesPar.q1     = qName[i];
+            mesPar.q1          = qName[i];
-            mesPar.q2     = qName[j];
+            mesPar.q2          = qName[j];
+            mesPar.gammaSource = Gamma::Algebra::Gamma5;
+            mesPar.gammaSink   = Gamma::Algebra::Gamma5;
             application.createModule<MContraction::Meson>("meson_Z2_"
                                                           + std::to_string(t)
                                                           + "_"

tests/qdpxx/Test_qdpxx_loops_staples.cc

@@ -25,7 +25,7 @@ Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 double calc_grid_p      (Grid::QCD::LatticeGaugeField & lat);
 double calc_chroma_p    (Grid::QCD::LatticeGaugeField & lat);

tests/qdpxx/Test_qdpxx_munprec.cc

@@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 int    Ls=8;
 double M5=1.6;

tests/solver/Test_cf_cr_unprec.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/solver/Test_contfrac_cg.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 

tests/solver/Test_dwf_cg_prec.cc

@@ -37,8 +37,8 @@ struct scal {
   d internal;
 };
 
-Gamma::GammaMatrix Gmu[] = {Gamma::GammaX, Gamma::GammaY, Gamma::GammaZ,
+Gamma::Algebra Gmu[] = {Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::Algebra::GammaZ,
-                            Gamma::GammaT};
+                            Gamma::Algebra::GammaT};
 
 int main(int argc, char** argv) {
   Grid_init(&argc, &argv);

tests/solver/Test_dwf_cg_schur.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/solver/Test_dwf_cg_unprec.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/solver/Test_dwf_cr_unprec.cc

@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/solver/Test_dwf_fpgcr.cc

@@ -38,11 +38,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/solver/Test_dwf_hdcr.cc

@@ -504,7 +504,7 @@ int main (int argc, char ** argv)
   GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
   GridParallelRNG          CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
 
-  Gamma g5(Gamma::Gamma5);
+  Gamma g5(Gamma::Algebra::Gamma5);
 
   LatticeFermion    src(FGrid); gaussian(RNG5,src);// src=src+g5*src;
   LatticeFermion result(FGrid); result=zero;

tests/solver/Test_wilson_cg_prec.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/solver/Test_wilson_cg_schur.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/solver/Test_wilson_cg_unprec.cc

@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)

tests/solver/Test_wilson_cr_unprec.cc

@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
+  Gamma::Algebra Gmu [] = {
-    Gamma::GammaX,
+    Gamma::Algebra::GammaX,
-    Gamma::GammaY,
+    Gamma::Algebra::GammaY,
-    Gamma::GammaZ,
+    Gamma::Algebra::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)