New automatic implementation of gamma matrices, Meson and SeqGamma are broken

2025-11-04 05:54:32 +00:00 · 2017-01-23 19:13:43 -08:00
parent d6401e6d2c
commit a37e71f362
41 changed files with 4129 additions and 882 deletions
--- 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
@@ -118,11 +118,11 @@ void TMeson<FImpl1, FImpl2>::execute(void)
    std::vector<TComplex> buf;
    Result                result;
    
-    g5 = makeGammaProd(Ns*Ns - 1);
+    //g5 = makeGammaProd(Ns*Ns - 1);
    result.corr.resize(Ns*Ns);
    for (unsigned int i = 0; i < Ns*Ns; ++i)
    {
-        g[i] = makeGammaProd(i);
+        //g[i] = makeGammaProd(i);
    }
    for (unsigned int iSink = 0; iSink < Ns*Ns; ++iSink)
    {
--- a/extras/Hadrons/Modules/MSource/SeqGamma.hpp
+++ b/extras/Hadrons/Modules/MSource/SeqGamma.hpp
@@ -144,7 +144,7 @@ void TSeqGamma<FImpl>::execute(void)
    std::vector<Real>          p;
    Complex                    i(0.0,1.0);
    
-    g  = makeGammaProd(par().gamma);
+    //g  = makeGammaProd(par().gamma);
    p  = strToVec<Real>(par().mom);
    ph = zero;
    for(unsigned int mu = 0; mu < env().getNd(); mu++)
--- 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,42 @@
-    /*************************************************************************************
-
-    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.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  ",
-      "-Gamma5  ",
-      "         "
-    };
+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/scripts/filelist
+++ b/scripts/filelist
@@ -4,9 +4,9 @@ home=`pwd`

 # library Make.inc
 cd $home/lib
-HFILES=`find . -type f -name '*.h' -not -path '*/Old/*' -not -path '*/Eigen/*'`
+HFILES=`find . -type f -name '*.h' -not -path '*/gamma-gen/*' -not -path '*/Old/*' -not -path '*/Eigen/*'`
 HFILES="$HFILES"
-CCFILES=`find . -type f -name '*.cc' -not  -name '*ommunicator*.cc'`
+CCFILES=`find . -type f -name '*.cc' -not -path '*/gamma-gen/*' -not -name '*ommunicator*.cc'`
 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,12 +36,12 @@ 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;
 typedef ZMobiusFermion<ZDomainWallVec5dImplR>                        ZMobiusVecFermionR;
@@ -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
@@ -81,7 +81,7 @@ int main (int argc, char ** argv)
    ident()(a,a) = ComplexF(1.0);
  }

-  const Gamma::GammaMatrix *g = Gamma::GammaMatrices;
+  const Gamma::Algebra *g = Gamma::GammaMatrices;
  const char **list           = Gamma::GammaMatrixNames;

  result =ll*Gamma(g[0])*rr;
@@ -152,70 +152,70 @@ int main (int argc, char ** argv)
  double mag;
  spProjXp(hsm,rv);
  spReconXp(recon,hsm);
-  full = rv + Gamma(Gamma::GammaX) *rv;
+  full = rv + Gamma(Gamma::Algebra::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;
+  full = rv - Gamma(Gamma::Algebra::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;
+  full = rv + Gamma(Gamma::Algebra::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;
+  full = rv - Gamma(Gamma::Algebra::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;
+  full = rv + Gamma(Gamma::Algebra::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;
+  full = rv - Gamma(Gamma::Algebra::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;
+  full = rv + Gamma(Gamma::Algebra::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;
+  full = rv - Gamma(Gamma::Algebra::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;
+  full = rv + Gamma(Gamma::Algebra::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;
+  full = rv - Gamma(Gamma::Algebra::Gamma5) *rv;
  mag = TensorRemove(norm2(full-recon));
  std::cout << "5m "<< mag<<std::endl;

--- 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)