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

.gitignore

@@ -112,3 +112,8 @@ m4/libtool.m4
 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;
@@ -321,7 +321,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

@@ -319,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])
 
@@ -330,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]); 
      ;;

extras/Hadrons/Modules/MContraction/Meson.hpp

@@ -47,7 +47,9 @@ public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(MesonPar,
                                     std::string,    q1,
                                     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);
-    for (unsigned int i = 0; i < Ns*Ns; ++i)
-    {
-        g[i] = makeGammaProd(i);
-    }
-    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());
+    result.corr.resize(buf.size());
     for (unsigned int t = 0; t < buf.size(); ++t)
     {
-                result.corr[iSink][iSrc][t] = TensorRemove(buf[t]);
+        result.corr[t] = TensorRemove(buf[t]);
-            }
-        }
     }
     write(writer, "meson", result);
 }

extras/Hadrons/Modules/MSource/SeqGamma.hpp

@@ -63,7 +63,7 @@ public:
                                     std::string,    q,
                                     unsigned int,   tA,
                                     unsigned int,   tB,
-                                    unsigned int, gamma,
+                                    Gamma::Algebra, gamma,
                                     std::string,    mom);
 };
 
@@ -140,11 +140,10 @@ 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 < env().getNd(); mu++)
@@ -154,7 +153,7 @@ void TSeqGamma<FImpl>::execute(void)
     }
     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/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 {
 
@@ -114,10 +115,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/qcd/action/fermion/WilsonFermion.cc

@@ -149,11 +149,11 @@ void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
 
     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;

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -503,11 +503,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;
@@ -574,11 +574,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/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/spin/Dirac.cc

@@ -1,95 +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/Grid.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> 
-  public:
+inline auto operator*(const Gamma &G, const iScalar<vtype> &arg)
-
+->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
-    const int Ns=4;
+{
-    
-    enum GammaMatrix {
-      Identity,                        
-      GammaX, 
-      GammaY, 
-      GammaZ, 
-      GammaT,  
-      Gamma5,
-      MinusIdentity,                        
-      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; }
-
-    GammaMatrix _g;
-
-  };
-  
-    // Make gamma products (Chroma convention)
-    SpinMatrix makeGammaProd(const unsigned int i);
-    
-    /* Gx
-     *  0 0  0  i    
-     *  0 0  i  0    
-     *  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){
-      ret._internal[0] = timesI(rhs._internal[3]);
-      ret._internal[1] = timesI(rhs._internal[2]);
-      ret._internal[2] = timesMinusI(rhs._internal[1]);
-      ret._internal[3] = timesMinusI(rhs._internal[0]);
-    };
-    template<class vtype> inline void multMinusGammaX(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-	ret(0) = timesMinusI(rhs(3));
-	ret(1) = timesMinusI(rhs(2));
-	ret(2) = timesI(rhs(1));
-	ret(3) = timesI(rhs(0));
-    };
-
-
-    /*Gy
-     *  0 0  0  -1  [0] -+ [3]
-     *  0 0  1  0   [1] +- [2]
-     *  0 1  0  0
-     * -1 0  0  0
-     */
-    template<class vtype> inline void rmultGammaY(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 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){
-      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 multMinusGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      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){
-      ret(0) = rhs(0);
-      ret(1) = rhs(1);
-      ret(2) =-rhs(2);
-      ret(3) =-rhs(3);
-    };
-    template<class vtype> inline void multMinusGamma5(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) =-rhs(0);
-      ret(1) =-rhs(1);
-      ret(2) = rhs(2);
-      ret(3) = rhs(3);
-    };
-
-
-
-#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 
+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 
+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;
-	}
+}
 
-
+// general right multiply
-    //right multiply
+template<class vtype>
-    template<class vtype> inline auto operator * (const iScalar<vtype> &arg, const Gamma &G) ->
+inline auto operator*(const iScalar<vtype> &arg, const Gamma &G)
-      typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype> >::type 
+->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 
+template<class vtype, int N>
-	{
+inline auto operator * (const iMatrix<vtype, N> &arg, const Gamma &G) 
-	  iVector<vtype,N> ret;
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type 
-	  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;
+}
+
+// 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;
+  
+  buf(0) = 0.;
+  buf(1) = 0.;
+  buf(2) = 2.*arg(2);
+  buf(3) = 2.*arg(3);
+  
+  return gl.gamma*buf;
+};
+
+// matrix left multiply
+template<class vtype> 
+inline auto operator*(const GammaL &gl, const iMatrix<vtype, Ns> &arg)
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, iMatrix<vtype, Ns>>::type
+{
+  iMatrix<vtype, Ns> buf;
+  
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    buf(0, i) = 0.;
+    buf(1, i) = 0.;
+    buf(2, i) = 2.*arg(2, i);
+    buf(3, i) = 2.*arg(3, i);
   }
   
-    ////////////////////////////////////////////////////////
+  return gl.gamma*buf;
-    // 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) ->
+// matrix right multiply
-      typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,Ns>,SpinorIndex>::value,iMatrix<vtype,Ns> >::type 
+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)
   {
-	iMatrix<vtype,Ns> ret;
+    buf(i, 0) = 0.;
-	switch (G._g) {
+    buf(i, 1) = 0.;
-	case Gamma::Identity:
+    buf(i, 2) = 2.*buf(i, 2);
-	  ret = arg;
+    buf(i, 3) = 2.*buf(i, 3);
-	  break;
+  }
-	case Gamma::MinusIdentity:
+  
-	  ret = -arg;
+  return buf;
-	  break;
+};
-	case Gamma::GammaX:
+
-	  multGammaX(ret,arg);
+//general left multiply
-	  break;
+template<class vtype> 
-	case Gamma::MinusGammaX:
+inline auto operator*(const GammaL &gl, const iScalar<vtype> &arg)
-	  multMinusGammaX(ret,arg);
+->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
-	  break;
+{
-	case Gamma::GammaY:
+  iScalar<vtype> ret;
-	  multGammaY(ret,arg);
+  ret._internal=gl*arg._internal;
-	  break;
+  return ret;
-	case Gamma::MinusGammaY:
+}
-	  multMinusGammaY(ret,arg);
+
-	  break;
+template<class vtype,int N>
-	case Gamma::GammaZ:
+inline auto operator*(const GammaL &gl, const iVector<vtype, N> &arg)
-	  multGammaZ(ret,arg);
+->typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N>>::type 
-	  break;
+{
-	case Gamma::MinusGammaZ:
+  iVector<vtype,N> ret;
-	  multMinusGammaZ(ret,arg);
+  for(int i=0;i<N;i++){
-	  break;
+    ret._internal[i]=gl*arg._internal[i];
-	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) ->
+template<class vtype, int N>
-      typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type 
+inline auto operator*(const GammaL &gl, const iMatrix<vtype, N> &arg) 
-      {
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type 
-	iVector<vtype,Ns> ret;
+{
-	switch (G._g) {
+  iMatrix<vtype,N> ret;
-	case Gamma::Identity:
+  for(int i=0;i<N;i++){
-	  ret = arg;
+  for(int j=0;j<N;j++){
-	  break;
+    ret._internal[i][j]=gl*arg._internal[i][j];
-	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) ->
+//general right multiply
-      typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,Ns>,SpinorIndex>::value,iMatrix<vtype,Ns> >::type 
+template<class vtype>
-      {
+inline auto operator*(const iScalar<vtype> &arg, const GammaL &gl)
-	iMatrix<vtype,Ns> ret;
+->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
-	switch (G._g) {
+{
-	case Gamma::Identity:
+  iScalar<vtype> ret;
-	  ret = arg;
+  ret._internal=arg._internal*gl;
-	  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;
-      }
+}
 
+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;
+}
 
-    /* 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_LOOP
   for(int ss=0;ss<grid->oSites();ss++){
     vobj tmp;
@@ -80,7 +80,7 @@ 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_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     vobj tmp;
@@ -98,7 +98,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_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     vobj tmp;
@@ -116,7 +116,7 @@ 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_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     vobj tmp1;
@@ -168,7 +168,7 @@ 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_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     vobj tmp;

lib/serialisation/MacroMagic.h

@@ -150,14 +150,14 @@ friend inline bool operator==(const cname &lhs, const cname &rhs) {\
 class name: public Grid::Serializable\
 {\
 public:\
-  enum EnumType\
+  enum\
   {\
     GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMVAL,__VA_ARGS__))\
     undefname = -1\
   };\
 public:\
   name(void): value_(undefname) {};\
-  name(EnumType value): value_(value) {};\
+  name(int value): value_(value) {};\
   template <typename T>\
   static inline void write(Grid::Writer<T> &WR,const std::string &s, const name &obj)\
   {\
@@ -177,7 +177,7 @@ public:\
     GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMTEST,__VA_ARGS__))\
     else {obj = name::undefname;}\
   }\
-  inline operator EnumType(void) const\
+  inline operator int(void) const\
   {\
     return value_;\
   }\
@@ -189,8 +189,17 @@ public:\
     }\
     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:\
-  EnumType value_;\
+  int value_;\
 };
 
 

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,8 +4,9 @@ home=`pwd`
 
 # library Make.inc
 cd $home/lib
-HFILES=`find . -type f -name '*.h' -not -name '*Hdf5*' -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/*'`
-CCFILES=`find . -type f -name '*.cc' -not -name '*Communicator*.cc' -not -name '*Hdf5*'`
+HFILES="$HFILES"
+CCFILES=`find . -type f -name '*.cc' -not -path '*/gamma-gen/*' -not -name '*Communicator*.cc' -not -name '*Hdf5*'`
 echo HFILES=$HFILES > Make.inc
 echo >> Make.inc
 echo CCFILES=$CCFILES >> Make.inc

tests/Test_cayley_even_odd_vec.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
     };
 
 typedef DomainWallFermion<DomainWallVec5dImplR>                      DomainWallVecFermionR;
@@ -339,7 +339,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,43 +1,240 @@
-    /*************************************************************************************
+/*************************************************************************************
 
-    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);
   
@@ -46,178 +243,41 @@ int main (int argc, char ** argv)
   std::vector<int> mpi_layout  = GridDefaultMpi();
   
   GridCartesian Grid(latt_size,simd_layout,mpi_layout);
-
-  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 << GridLogMessage << std::endl;
+  std::cout << GridLogMessage << "======== Adjoints check" << std::endl;
+  for (int i = 0; i < Gamma::nGamma; ++i)
+  {
+    checkAdj(i);
   }
-      std::cout << std::endl;
+  std::cout << GridLogMessage << std::endl;
+  std::cout << GridLogMessage << "======== Spin projectors check" << std::endl;
+  checkProject(sRNG);
+  std::cout << GridLogMessage << std::endl;
+  std::cout << GridLogMessage << "======== Gamma-left matrices check" << std::endl;
+  for (int i = 0; i < Gamma::nGamma; ++i)
+  {
+    checkGammaL(i, sRNG);
   }
   
-    std::cout << std::endl;
-
-  }
-
-  std::cout << "Testing Gamma^2 - 1 = 0"<<std::endl;
-  for(int mu=0;mu<6;mu++){
-    result =  Gamma(g[mu])* ident * Gamma(g[mu]);
-    result = result - ident;
-    RealD mag = norm2(result);
-    std::cout << list[mu]<<" " << mag<<std::endl;
-  }
-
-  std::cout << "Testing (MinusGamma + G )M = 0"<<std::endl;
-  for(int mu=0;mu<6;mu++){
-    result =          rnd * Gamma(g[mu]);
-    result = result + rnd * Gamma(g[mu+6]);
-    RealD mag = norm2(result);
-    std::cout << list[mu]<<" " << mag<<std::endl;
-  }
-
-  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/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];
@@ -130,6 +138,8 @@ int main(int argc, char *argv[])
             mesPar.output      = "mesons/Z2_" + flavour[i] + flavour[j];
             mesPar.q1          = qName[i];
             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/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)