Merge remote-tracking branch 'paboyle/feature/hadrons' into feature/hadrons

# Conflicts: # extras/Hadrons/Modules/MContraction/Meson.hpp # tests/hadrons/Test_hadrons_meson_3pt.cc Updated Meson.hpp to utilise zero-flop gamma matrices.
2025-10-31 03:54:33 +00:00 · 2017-02-01 09:27:00 +00:00
parent ef31c012bf 7da7d263c4
commit f8fbe4d7a3
43 changed files with 4551 additions and 1073 deletions
--- a/.gitignore
+++ b/.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
--- a/benchmarks/Benchmark_dwf.cc
+++ b/benchmarks/Benchmark_dwf.cc
@@ -37,11 +37,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    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]; ;
--- a/benchmarks/Benchmark_dwf_sweep.cc
+++ b/benchmarks/Benchmark_dwf_sweep.cc
@@ -37,11 +37,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 void benchDw(std::vector<int> & L, int Ls, int threads, int report =0 );
--- a/benchmarks/Benchmark_wilson.cc
+++ b/benchmarks/Benchmark_wilson.cc
@@ -37,11 +37,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    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]; ;
--- a/benchmarks/Benchmark_wilson_sweep.cc
+++ b/benchmarks/Benchmark_wilson_sweep.cc
@@ -30,11 +30,11 @@ struct scal {
  d internal;
 };

-Gamma::GammaMatrix Gmu [] = {
-  Gamma::GammaX,
-  Gamma::GammaY,
-  Gamma::GammaZ,
-  Gamma::GammaT
+Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
 };

 bool overlapComms = false;
--- a/extras/Hadrons/Modules/MContraction/Meson.hpp
+++ b/extras/Hadrons/Modules/MContraction/Meson.hpp
@@ -37,17 +37,17 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/ModuleFactory.hpp>

 namespace Grid {
-    // Overload >> to extract gamma pair from "[g1 g2]" string.
+    // Overload >> to extract gamma pair from "<g1 g2>" string.
    template <typename T1, typename T2>
    inline std::istringstream &operator>>(std::istringstream &sstr,
                                          std::pair<T1, T2> &buf)
    {
-        T1 buf1;
-        T2 buf2;
+        unsigned int buf1;
+        unsigned int buf2;
        char c;
        sstr >> c >> buf1 >> buf2 >> c;
        sstr.peek();
-        buf = std::make_pair(buf1, buf2);
+        buf = std::make_pair((T1)buf1, (T2)buf2);
        return sstr;
    }
 }
@@ -76,7 +76,7 @@ BEGIN_HADRONS_NAMESPACE
 ******************************************************************************/
 BEGIN_MODULE_NAMESPACE(MContraction)

-typedef std::pair<unsigned int, unsigned int> GammaPair;
+typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;

 class MesonPar: Serializable
 {
@@ -99,8 +99,8 @@ public:
    {
    public:
        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        unsigned int, gamma_snk,
-                                        unsigned int, gamma_src,
+                                        Gamma::Algebra, gamma_snk,
+                                        Gamma::Algebra, gamma_src,
                                        std::vector<Complex>, corr);
    };
 public:
@@ -151,13 +151,14 @@ void TMeson<FImpl1, FImpl2>::parseGammaString(std::vector<GammaPair> &gammaList)
    if (par().gammas.compare("all") == 0)
    {
        // Do all contractions.
-        unsigned int n_gam = Ns*Ns;
+        unsigned int n_gam = Ns * Ns;
        gammaList.resize(n_gam*n_gam);
-        for (unsigned int i = 0; i < n_gam; ++i)
+        for (unsigned int i = 1; i < Gamma::nGamma; i += 2)
        {
-            for (unsigned int j = 0; j < n_gam; ++j)
+            for (unsigned int j = 1; j < Gamma::nGamma; j += 2)
            {
-                gammaList.push_back(std::make_pair(i, j));
+                gammaList.push_back(std::make_pair((Gamma::Algebra)i, 
+                                                   (Gamma::Algebra)j));
            }
        }
    }
@@ -181,7 +182,7 @@ 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                  g5(Gamma::Algebra::Gamma5);
    std::vector<GammaPair> gammaList;
    std::vector<TComplex>  buf;
    std::vector<Result>    result;
@@ -198,17 +199,14 @@ void TMeson<FImpl1, FImpl2>::execute(void)
    }
    ph = exp(-2*M_PI*i*ph);
    
-    g5 = makeGammaProd(Ns*Ns - 1);
-    for (int i = 0; i < Ns*Ns; ++i)
-    {
-        g[i] = makeGammaProd(i);
-    }
    parseGammaString(gammaList);

    result.resize(gammaList.size());
    for (unsigned int i = 0; i < result.size(); ++i)
    {
-        c = trace(g[gammaList[i].first]*q1*g[gammaList[i].second]*g5*adj(q2)*g5*ph);
+        Gamma gSnk(gammaList[i].first);
+        Gamma gSrc(gammaList[i].second);
+        c = trace((g5*gSnk)*q1*(gSrc*g5)*adj(q2))*ph;
        sliceSum(c, buf, Tp);

        result[i].gamma_snk = gammaList[i].first;
--- a/extras/Hadrons/Modules/MSource/SeqGamma.hpp
+++ b/extras/Hadrons/Modules/MSource/SeqGamma.hpp
@@ -64,7 +64,7 @@ public:
                                    std::string,    q,
                                    unsigned int,   tA,
                                    unsigned int,   tB,
-                                    unsigned int, gamma,
+                                    Gamma::Algebra, gamma,
                                    std::string,    mom);
 };

@@ -141,11 +141,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++)
@@ -155,7 +154,7 @@ void TSeqGamma<FImpl>::execute(void)
    }
    ph = exp(2*M_PI*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
--- a/lib/qcd/action/fermion/WilsonFermion.cc
+++ b/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::GammaX,
-      Gamma::GammaY,
-      Gamma::GammaZ,
-      Gamma::GammaT
+    Gamma::Algebra Gmu [] = {
+      Gamma::Algebra::GammaX,
+      Gamma::Algebra::GammaY,
+      Gamma::Algebra::GammaZ,
+      Gamma::Algebra::GammaT
    };

    std::vector<int> latt_size   = _grid->_fdimensions;
--- a/lib/qcd/action/fermion/WilsonFermion5D.cc
+++ b/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::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    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::GammaX,
-      Gamma::GammaY,
-      Gamma::GammaZ,
-      Gamma::GammaT
+    Gamma::Algebra Gmu [] = {
+      Gamma::Algebra::GammaX,
+      Gamma::Algebra::GammaY,
+      Gamma::Algebra::GammaZ,
+      Gamma::Algebra::GammaT
    };

    GridBase *_grid = _FourDimGrid;
--- a/lib/qcd/action/fermion/g5HermitianLinop.h
+++ b/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);
  }
--- a/lib/qcd/spin/Dirac.cc
+++ b/lib/qcd/spin/Dirac.cc
@@ -1,95 +1,72 @@
-    /*************************************************************************************
+/*************************************************************************************

-    Grid physics library, www.github.com/paboyle/Grid 
+Grid physics library, www.github.com/paboyle/Grid 

-    Source file: ./lib/qcd/spin/Dirac.cc
+Source file: lib/qcd/spin/Dirac.cc

-    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: paboyle <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 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.
+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.
+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 */
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
 #include <Grid/Grid.h>

 namespace Grid {
+namespace QCD {

-  namespace QCD {
+#include "GammaMulTable.h"

-    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  ",
+const std::array<const char *, Gamma::nGamma> Gamma::name = {{
  "-Gamma5      ",
-      "         "
-    };
+  "Gamma5       ",
+  "-GammaT      ",
+  "GammaT       ",
+  "-GammaTGamma5",
+  "GammaTGamma5 ",
+  "-GammaX      ",
+  "GammaX       ",
+  "-GammaXGamma5",
+  "GammaXGamma5 ",
+  "-GammaY      ",
+  "GammaY       ",
+  "-GammaYGamma5",
+  "GammaYGamma5 ",
+  "-GammaZ      ",
+  "GammaZ       ",
+  "-GammaZGamma5",
+  "GammaZGamma5 ",
+  "-Identity    ",
+  "Identity     ",
+  "-SigmaXT     ",
+  "SigmaXT      ",
+  "-SigmaXY     ",
+  "SigmaXY      ",
+  "-SigmaXZ     ",
+  "SigmaXZ      ",
+  "-SigmaYT     ",
+  "SigmaYT      ",
+  "-SigmaYZ     ",
+  "SigmaYZ      ",
+  "-SigmaZT     ",
+  "SigmaZT      "}};

-    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);
-    //    }
-  }
-}
+}}
--- a/lib/qcd/spin/Dirac.h
+++ b/lib/qcd/spin/Dirac.h
--- a/lib/qcd/spin/GammaMulTable.h
+++ b/lib/qcd/spin/GammaMulTable.h
--- a/lib/qcd/spin/gamma-gen/gamma-gen.nb
+++ b/lib/qcd/spin/gamma-gen/gamma-gen.nb
--- a/lib/qcd/utils/LinalgUtils.h
+++ b/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;
--- a/lib/serialisation/MacroMagic.h
+++ b/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_;\
  }\
@@ -190,7 +190,7 @@ public:\
    return os;\
  }\
 private:\
-  EnumType value_;\
+  int value_;\
 };


--- a/scripts/filelist
+++ b/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/*'`
-CCFILES=`find . -type f -name '*.cc' -not -name '*Communicator*.cc' -not -name '*Hdf5*'`
+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 -path '*/gamma-gen/*' -not -name '*Communicator*.cc' -not -name '*Hdf5*'`
 echo HFILES=$HFILES > Make.inc
 echo >> Make.inc
 echo CCFILES=$CCFILES >> Make.inc
--- a/tests/Test_cayley_even_odd_vec.cc
+++ b/tests/Test_cayley_even_odd_vec.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+    Gamma::Algebra Gmu [] = {
+        Gamma::Algebra::GammaX,
+        Gamma::Algebra::GammaY,
+        Gamma::Algebra::GammaZ,
+        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);  
--- a/tests/Test_dwf_mixedcg_prec.cc
+++ b/tests/Test_dwf_mixedcg_prec.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/core/Test_cf_coarsen_support.cc
+++ b/tests/core/Test_cf_coarsen_support.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/core/Test_checker.cc
+++ b/tests/core/Test_checker.cc
@@ -37,11 +37,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int toint(const char* str){
--- a/tests/core/Test_contfrac_even_odd.cc
+++ b/tests/core/Test_contfrac_even_odd.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };


--- a/tests/core/Test_dwf_even_odd.cc
+++ b/tests/core/Test_dwf_even_odd.cc
@@ -37,11 +37,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };


--- a/tests/core/Test_dwf_rb5d.cc
+++ b/tests/core/Test_dwf_rb5d.cc
@@ -37,11 +37,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 typedef WilsonFermion5D<DomainWallVec5dImplR> WilsonFermion5DR;
--- a/tests/core/Test_fft.cc
+++ b/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::GammaX,
-      Gamma::GammaY,
-      Gamma::GammaZ,
-      Gamma::GammaT,
-      Gamma::Gamma5
+    Gamma::Algebra Gmu [] = {
+      Gamma::Algebra::GammaX,
+      Gamma::Algebra::GammaY,
+      Gamma::Algebra::GammaZ,
+      Gamma::Algebra::GammaT,
+      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); 
--- a/tests/core/Test_gamma.cc
+++ b/tests/core/Test_gamma.cc
@@ -1,43 +1,218 @@
-    /*************************************************************************************
+/*************************************************************************************

-    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
-    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 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.
+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.
+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
-    *************************************************************************************/
+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 checkMat(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
+}
+
+int main(int argc, char *argv[])
 {
  Grid_init(&argc,&argv);
  
@@ -46,178 +221,36 @@ 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;
-  SpinMatrix rnd  ; random(sRNG,rnd);
-
-  SpinMatrix ll; ll=zero;
-  SpinMatrix rr; rr=zero;
-  SpinMatrix result;
-
-  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 << "======== Test algebra" << std::endl;
+  createTestAlgebra();
+  std::cout << GridLogMessage << "======== Multiplication operators check" << std::endl;
+  for (int i = 0; i < Gamma::nGamma; ++i)
+  {
+    checkMat(i, sRNG);
  }
-
-  const Gamma::GammaMatrix *g = Gamma::GammaMatrices;
-  const char **list           = Gamma::GammaMatrixNames;
-
-  result =ll*Gamma(g[0])*rr;
-  result =ll*Gamma(g[0]);
-  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 << GridLogMessage << std::endl;
+  std::cout << GridLogMessage << "======== Algebra multiplication table check" << std::endl;
+  for (int i = 0; i < Gamma::nGamma; ++i)
+  for (int j = 0; j < Gamma::nGamma; ++j)
+  {
+    checkProd(i, j);
  }
-	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 << 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;
+  std::cout << GridLogMessage << std::endl;
+  std::cout << GridLogMessage << "======== Spin projectors check" << std::endl;
+  checkProject(sRNG);
+  std::cout << GridLogMessage << std::endl;
  
  Grid_finalize();
+  
+  return EXIT_SUCCESS;
 }
--- a/tests/core/Test_wilson_even_odd.cc
+++ b/tests/core/Test_wilson_even_odd.cc
@@ -37,11 +37,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/core/Test_wilson_tm_even_odd.cc
+++ b/tests/core/Test_wilson_tm_even_odd.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/debug/Test_cayley_cg.cc
+++ b/tests/debug/Test_cayley_cg.cc
@@ -37,11 +37,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 template<class What> 
--- a/tests/debug/Test_cayley_coarsen_support.cc
+++ b/tests/debug/Test_cayley_coarsen_support.cc
@@ -38,11 +38,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/debug/Test_cayley_even_odd.cc
+++ b/tests/debug/Test_cayley_even_odd.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };


--- a/tests/solver/Test_cf_cr_unprec.cc
+++ b/tests/solver/Test_cf_cr_unprec.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/solver/Test_contfrac_cg.cc
+++ b/tests/solver/Test_contfrac_cg.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };


--- a/tests/solver/Test_dwf_cg_prec.cc
+++ b/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::GammaT};
+Gamma::Algebra Gmu[] = {Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::Algebra::GammaZ,
+                            Gamma::Algebra::GammaT};

 int main(int argc, char** argv) {
  Grid_init(&argc, &argv);
--- a/tests/solver/Test_dwf_cg_schur.cc
+++ b/tests/solver/Test_dwf_cg_schur.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/solver/Test_dwf_cg_unprec.cc
+++ b/tests/solver/Test_dwf_cg_unprec.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/solver/Test_dwf_cr_unprec.cc
+++ b/tests/solver/Test_dwf_cr_unprec.cc
@@ -37,11 +37,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/solver/Test_dwf_fpgcr.cc
+++ b/tests/solver/Test_dwf_fpgcr.cc
@@ -38,11 +38,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/solver/Test_dwf_hdcr.cc
+++ b/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;
--- a/tests/solver/Test_wilson_cg_prec.cc
+++ b/tests/solver/Test_wilson_cg_prec.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/solver/Test_wilson_cg_schur.cc
+++ b/tests/solver/Test_wilson_cg_schur.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/solver/Test_wilson_cg_unprec.cc
+++ b/tests/solver/Test_wilson_cg_unprec.cc
@@ -37,11 +37,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)
--- a/tests/solver/Test_wilson_cr_unprec.cc
+++ b/tests/solver/Test_wilson_cr_unprec.cc
@@ -36,11 +36,11 @@ struct scal {
  d internal;
 };

-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
  };

 int main (int argc, char ** argv)