Gammas: code cleaning and gamma_L implementation & test

lib/qcd/spin/Dirac.h

@@ -32,32 +32,28 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #ifndef GRID_QCD_DIRAC_H
 #define GRID_QCD_DIRAC_H
 
-// code generated by the Mathematica notebook gamma-gen/gamma-gen.nb
+// Gamma matrices using the code generated by the Mathematica notebook 
-#include <array>
+// gamma-gen/gamma-gen.nb in Gamma.cc & Gamma.h
+////////////////////////////////////////////////////////////////////////////////
 #include <Grid/qcd/spin/Gamma.h>
 
 namespace Grid {
-  
+
-inline QCD::Gamma adj (const QCD::Gamma & g)
+// 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 {
 
-inline Gamma operator*(const Gamma & g1, const Gamma & g2)
+// Dirac algebra mutliplication operator
+inline Gamma operator*(const Gamma &g1, const Gamma &g2)
 {
    return Gamma (Gamma::mul[g1.g][g2.g]);
 }
 
-// FIXME
+// general left multiply
-//
-// 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 
@@ -90,7 +86,7 @@ inline auto operator*(const Gamma &G, const iMatrix<vtype, N> &arg)
   return ret;
 }
 
-//right multiply
+// general right multiply
 template<class vtype>
 inline auto operator*(const iScalar<vtype> &arg, const Gamma &G)
 ->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
@@ -112,6 +108,125 @@ inline auto operator * (const iMatrix<vtype, N> &arg, const Gamma &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;
+};
+
+// matrix right multiply
+template<class vtype> 
+inline auto operator*(const iMatrix<vtype, Ns> &arg, const GammaL &gl)
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, iMatrix<vtype, Ns>>::type
+{
+  iMatrix<vtype, Ns> buf;
+  
+  buf = arg*gl.gamma;
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    buf(i, 0) = 0.;
+    buf(i, 1) = 0.;
+    buf(i, 2) = 2.*buf(i, 2);
+    buf(i, 3) = 2.*buf(i, 3);
+  }
+  
+  return buf;
+};
+
+//general left multiply
+template<class vtype> 
+inline auto operator*(const GammaL &gl, const iScalar<vtype> &arg)
+->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
+{
+  iScalar<vtype> ret;
+  ret._internal=gl*arg._internal;
+  return ret;
+}
+
+template<class vtype,int N>
+inline auto operator*(const GammaL &gl, 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]=gl*arg._internal[i];
+  }
+  return ret;
+}
+
+template<class vtype, int N>
+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 
+{
+  iMatrix<vtype,N> ret;
+  for(int i=0;i<N;i++){
+  for(int j=0;j<N;j++){
+    ret._internal[i][j]=gl*arg._internal[i][j];
+  }}
+  return ret;
+}
+
+//general right multiply
+template<class vtype>
+inline auto operator*(const iScalar<vtype> &arg, const GammaL &gl)
+->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
+{
+  iScalar<vtype> ret;
+  ret._internal=arg._internal*gl;
+  return ret;
+}
+
+template<class vtype, int N>
+inline auto operator * (const iMatrix<vtype, N> &arg, const GammaL &gl) 
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type 
+{
+  iMatrix<vtype,N> ret;
+  for(int i=0;i<N;i++){
+  for(int j=0;j<N;j++){
+    ret._internal[i][j]=arg._internal[i][j]*gl;
+  }}
+  return ret;
+}
+
 }}
 
 #endif

tests/core/Test_gamma.cc

@@ -137,7 +137,7 @@ void test(const Expr &a, const Expr &b)
   }
 }
 
-void checkMat(const Gamma::Algebra a, GridSerialRNG &rng)
+void checkGamma(const Gamma::Algebra a, GridSerialRNG &rng)
 {
   SpinVector v;
   SpinMatrix m, &testg = testAlgebra[a];
@@ -212,6 +212,28 @@ std::cout << std::endl;
 #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);
@@ -230,7 +252,7 @@ int main(int argc, char *argv[])
   std::cout << GridLogMessage << "======== Multiplication operators check" << std::endl;
   for (int i = 0; i < Gamma::nGamma; ++i)
   {
-    checkMat(i, sRNG);
+    checkGamma(i, sRNG);
   }
   std::cout << GridLogMessage << std::endl;
   std::cout << GridLogMessage << "======== Algebra multiplication table check" << std::endl;
@@ -249,6 +271,11 @@ int main(int argc, char *argv[])
   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);
+  }
   
   Grid_finalize();