diff --git a/lib/Makefile.am b/lib/Makefile.am index 20c1f867..9034ce9c 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -4,18 +4,21 @@ AM_CXXFLAGS = -I$(top_srcdir)/ extra_sources= if BUILD_COMMS_MPI extra_sources+=communicator/Grid_communicator_mpi.cc - extra_sources+=stencil/Grid_stencil_common.cc endif + if BUILD_COMMS_NONE extra_sources+=communicator/Grid_communicator_fake.cc - extra_sources+=stencil/Grid_stencil_common.cc endif # # Libraries # lib_LIBRARIES = libGrid.a -libGrid_a_SOURCES = Grid_init.cc $(extra_sources) +libGrid_a_SOURCES =\ + Grid_init.cc\ + stencil/Grid_stencil_common.cc\ + qcd/Grid_qcd_dirac.cc\ + $(extra_sources) # # Include files @@ -32,7 +35,6 @@ include_HEADERS =\ Grid_math.h\ Grid_simd.h\ Grid_stencil.h\ - Grid_summation.h\ Grid_where.h nobase_include_HEADERS=\ @@ -68,7 +70,8 @@ nobase_include_HEADERS=\ math/Grid_math_trace.h\ math/Grid_math_traits.h\ math/Grid_math_transpose.h\ - qcd/Grid_QCD.h\ + qcd/Grid_qcd.h\ + qcd/Grid_qcd_dirac.h\ simd/Grid_vComplexD.h\ simd/Grid_vComplexF.h\ simd/Grid_vInteger.h\ diff --git a/lib/qcd/Grid_qcd_dirac.h b/lib/qcd/Grid_qcd_dirac.h index ef2a4cd8..945178e9 100644 --- a/lib/qcd/Grid_qcd_dirac.h +++ b/lib/qcd/Grid_qcd_dirac.h @@ -44,12 +44,14 @@ namespace QCD { // MinusSigmaYT, // MinusSigmaZT }; + + static GammaMatrix GammaMatrices[]; + static const char *GammaMatrixNames[]; Gamma (GammaMatrix g) { _g=g; } GammaMatrix _g; - }; @@ -94,13 +96,13 @@ namespace QCD { } }; - template inline void multGammaX(iVector &ret, iVector &rhs){ + template inline void multGammaX(iVector &ret, const iVector &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 inline void multMinusGammaX(iVector &ret, iVector &rhs){ + template inline void multMinusGammaX(iVector &ret, const iVector &rhs){ ret(0) = timesMinusI(rhs(3)); ret(1) = timesMinusI(rhs(2)); ret(2) = timesI(rhs(1)); @@ -146,13 +148,13 @@ namespace QCD { ret(3,i) = rhs(0,i); } }; - template inline void multGammaY(iVector &ret, iVector &rhs){ + template inline void multGammaY(iVector &ret, const iVector &rhs){ ret(0) = -rhs(3); ret(1) = rhs(2); ret(2) = rhs(1); ret(3) = -rhs(0); }; - template inline void multMinusGammaY(iVector &ret, iVector &rhs){ + template inline void multMinusGammaY(iVector &ret, const iVector &rhs){ ret(0) = rhs(3); ret(1) = -rhs(2); ret(2) = -rhs(1); @@ -196,13 +198,13 @@ namespace QCD { ret(3,i) = timesMinusI(rhs(1,i)); } }; - template inline void multGammaZ(iVector &ret, iVector &rhs){ + template inline void multGammaZ(iVector &ret, const iVector &rhs){ ret(0) = timesI(rhs(2)); ret(1) =timesMinusI(rhs(3)); ret(2) =timesMinusI(rhs(0)); ret(3) = timesI(rhs(1)); }; - template inline void multMinusGammaZ(iVector &ret, iVector &rhs){ + template inline void multMinusGammaZ(iVector &ret, const iVector &rhs){ ret(0) = timesMinusI(rhs(2)); ret(1) = timesI(rhs(3)); ret(2) = timesI(rhs(0)); @@ -246,13 +248,13 @@ namespace QCD { ret(3,i) =-rhs(1,i); } }; - template inline void multGammaT(iVector &ret, iVector &rhs){ + template inline void multGammaT(iVector &ret, const iVector &rhs){ ret(0) = rhs(2); ret(1) = rhs(3); ret(2) = rhs(0); ret(3) = rhs(1); }; - template inline void multMinusGammaT(iVector &ret, iVector &rhs){ + template inline void multMinusGammaT(iVector &ret, const iVector &rhs){ ret(0) =-rhs(2); ret(1) =-rhs(3); ret(2) =-rhs(0); @@ -298,13 +300,13 @@ namespace QCD { } }; - template inline void multGamma5(iVector &ret, iVector &rhs){ + template inline void multGamma5(iVector &ret, const iVector &rhs){ ret(0) = rhs(0); ret(1) = rhs(1); ret(2) =-rhs(2); ret(3) =-rhs(3); }; - template inline void multMinusGamma5(iVector &ret, iVector &rhs){ + template inline void multMinusGamma5(iVector &ret, const iVector &rhs){ ret(0) =-rhs(0); ret(1) =-rhs(1); ret(2) = rhs(2); @@ -313,6 +315,10 @@ namespace QCD { +#ifdef GRID_WARN_SUBOPTIMAL +#warning "Optimisation alert switch over to multGammaX early " +#endif + /////////////////////////////////////////////////////////////////////////////////////////////////// // Operator * : first case this is not a spin index, so recurse /////////////////////////////////////////////////////////////////////////////////////////////////// @@ -323,10 +329,8 @@ namespace QCD { // 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 // -#ifdef GRID_WARN_SUBOPTIMAL -#warning "Optimisation alert switch over to multGammaX early " -#endif + //left multiply template inline auto operator * ( const Gamma &G,const iScalar &arg) -> typename std::enable_if,SpinIndex>::notvalue,iScalar >::type @@ -349,6 +353,32 @@ namespace QCD { ret._internal=G*arg._internal; return ret; } + + + //right multiply + template inline auto operator * (const iScalar &arg, const Gamma &G) -> + typename std::enable_if,SpinIndex>::notvalue,iScalar >::type + + { + iScalar ret; + ret._internal=arg._internal*G; + return ret; + } + template inline auto operator * (const iVector &arg, const Gamma &G) -> + typename std::enable_if,SpinIndex>::notvalue,iVector >::type + { + iVector ret; + ret._internal=arg._internal*G; + return ret; + } + template inline auto operator * (const iMatrix &arg, const Gamma &G) -> + typename std::enable_if,SpinIndex>::notvalue,iMatrix >::type + { + iMatrix ret; + ret._internal=arg._internal*G; + return ret; + } + //////////////////////////////////////////////////////// // When we hit the spin index this matches and we stop //////////////////////////////////////////////////////// @@ -399,6 +429,104 @@ namespace QCD { } return ret; } + // Could have used type trait for Matrix/vector and then an enable if to share code + template inline auto operator * ( const Gamma &G,const iVector &arg) -> + typename std::enable_if,SpinIndex>::value,iVector >::type + { + iVector ret; + switch (G._g) { + case Gamma::Identity: + ret = arg; + break; + case Gamma::MinusIdentity: + ret = -arg; + break; + case Gamma::GammaX: + multGammaX(ret,arg); + break; + case Gamma::MinusGammaX: + multMinusGammaX(ret,arg); + break; + case Gamma::GammaY: + multGammaY(ret,arg); + break; + case Gamma::MinusGammaY: + multMinusGammaY(ret,arg); + break; + case Gamma::GammaZ: + multGammaZ(ret,arg); + break; + case Gamma::MinusGammaZ: + multMinusGammaZ(ret,arg); + break; + case Gamma::GammaT: + multGammaT(ret,arg); + break; + case Gamma::MinusGammaT: + multMinusGammaT(ret,arg); + break; + case Gamma::Gamma5: + multGamma5(ret,arg); + break; + case Gamma::MinusGamma5: + multMinusGamma5(ret,arg); + break; + default: + assert(0); + break; + } + return ret; + } + + template inline auto operator * (const iMatrix &arg, const Gamma &G) -> + typename std::enable_if,SpinIndex>::value,iMatrix >::type + { + iMatrix ret; + switch (G._g) { + case Gamma::Identity: + ret = arg; + break; + case Gamma::MinusIdentity: + ret = -arg; + break; + case Gamma::GammaX: + rmultGammaX(ret,arg); + break; + case Gamma::MinusGammaX: + rmultMinusGammaX(ret,arg); + break; + case Gamma::GammaY: + rmultGammaY(ret,arg); + break; + case Gamma::MinusGammaY: + rmultMinusGammaY(ret,arg); + break; + case Gamma::GammaZ: + rmultGammaZ(ret,arg); + break; + case Gamma::MinusGammaZ: + rmultMinusGammaZ(ret,arg); + break; + case Gamma::GammaT: + rmultGammaT(ret,arg); + break; + case Gamma::MinusGammaT: + rmultMinusGammaT(ret,arg); + break; + case Gamma::Gamma5: + rmultGamma5(ret,arg); + break; + case Gamma::MinusGamma5: + rmultMinusGamma5(ret,arg); + break; + default: + assert(0); + break; + } + return ret; + } + + /* Output from test ./Grid_gamma Identity((1,0),(0,0),(0,0),(0,0)) diff --git a/tests/Grid_gamma.cc b/tests/Grid_gamma.cc index 9c213352..d3a33039 100644 --- a/tests/Grid_gamma.cc +++ b/tests/Grid_gamma.cc @@ -23,6 +23,9 @@ int main (int argc, char ** argv) SpinMatrix rr=zero; SpinMatrix result; + SpinVector lv=zero; + SpinVector rv=zero; + for(int a=0;a