efficient baryon test program

Grid/qcd/utils/A2Autils.h

@@ -185,21 +185,29 @@ void A2Autils<FImpl>::NucleonFieldMom(Eigen::Tensor<ComplexD,6> &mat,
 
 	    auto v2 = conjugate(two[j]._odata[ss]);
             // C = i gamma_2 gamma_4 => C gamma_5 = - i gamma_1 gamma_3
-	    auto v2g = v2*Gamma(Gamma::Algebra::SigmaXZ);
+	    //auto v2g = v2*Gamma(Gamma::Algebra::SigmaXZ);
+            auto v2g=v2;  
 
 	    for(int k=0;k<threeBlock;k++){
 
 	      auto v3 = three[k]._odata[ss];
+	      auto gv3 = Gamma(Gamma::Algebra::SigmaXZ)*v3;
 	      SpinVector_v vv;
 
               for(int s1=0;s1<Ns;s1++){
 	      for(int s2=0;s2<Ns;s2++){            
-	        vv()(s1)() =  pv1()(s1)(0) * v2g()(s2)(1) * v3()(s2)(2)   //Cross product
+	     /*   vv()(s1)() =  pv1()(s1)(0) * v2g()(s2)(1) * v3()(s2)(2)   //Cross product
                   -           pv1()(s1)(0) * v2g()(s2)(2) * v3()(s2)(1)    
                   +           pv1()(s1)(1) * v2g()(s2)(2) * v3()(s2)(0)    
                   -           pv1()(s1)(1) * v2g()(s2)(0) * v3()(s2)(2)    
                   +           pv1()(s1)(2) * v2g()(s2)(0) * v3()(s2)(1)    
-                  -           pv1()(s1)(2) * v2g()(s2)(1) * v3()(s2)(0);    
+                  -           pv1()(s1)(2) * v2g()(s2)(1) * v3()(s2)(0);    */
+	        vv()(s1)() =  pv1()(s1)(0) * v2()(s2)(1) * gv3()(s2)(2)   //Cross product
+                  -           pv1()(s1)(0) * v2()(s2)(2) * gv3()(s2)(1)    
+                  +           pv1()(s1)(1) * v2()(s2)(2) * gv3()(s2)(0)    
+                  -           pv1()(s1)(1) * v2()(s2)(0) * gv3()(s2)(2)    
+                  +           pv1()(s1)(2) * v2()(s2)(0) * gv3()(s2)(1)    
+                  -           pv1()(s1)(2) * v2()(s2)(1) * gv3()(s2)(0);    
               }}
 	    
 	      // After getting the sitewise product do the mom phase loop
@@ -207,7 +215,9 @@ void A2Autils<FImpl>::NucleonFieldMom(Eigen::Tensor<ComplexD,6> &mat,
 	      for ( int m=0;m<Nmom;m++){
 	        int idx = m+base;
 	        auto phase = mom[m]._odata[ss];
-	        mac(&lvSum[idx],&vv,&phase()()());
+                for(int is=0;is<Ns;is++){
+	          mac(&lvSum[idx]()(is)(),&vv()(is)(),&phase()()());
+		}
               }
 	    }
 	  }
@@ -221,8 +231,7 @@ void A2Autils<FImpl>::NucleonFieldMom(Eigen::Tensor<ComplexD,6> &mat,
   parallel_for(int rt=0;rt<rd;rt++){
 
     std::vector<int> icoor(nd);
-    iScalar<vector_type> temp; 
+    std::vector<SpinVector_s>  extracted(Nsimd);               
-    std::vector<iScalar<SpinVector_s> > extracted(Nsimd);               
 
     for(int i=0;i<oneBlock;i++){
     for(int j=0;j<twoBlock;j++){
@@ -231,8 +240,7 @@ void A2Autils<FImpl>::NucleonFieldMom(Eigen::Tensor<ComplexD,6> &mat,
 
       int ij_rdx = m+Nmom*i + Nmom*oneBlock * j + Nmom*oneBlock * twoBlock * k + Nmom*oneBlock * twoBlock *threeBlock * rt;
 
-      temp._internal = lvSum[ij_rdx];
+      extract(lvSum[ij_rdx],extracted);
-      extract(temp,extracted);
 
       for(int idx=0;idx<Nsimd;idx++){
 
@@ -242,7 +250,7 @@ void A2Autils<FImpl>::NucleonFieldMom(Eigen::Tensor<ComplexD,6> &mat,
 
 	int ij_ldx = m+Nmom*i + Nmom*oneBlock * j + Nmom*oneBlock * twoBlock * k + Nmom*oneBlock * twoBlock *threeBlock * ldx;
 
-	lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx]._internal;
+	lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx];
 
       }
     }}}}

Hadrons/Modules.hpp

@@ -15,6 +15,7 @@
 #include <Hadrons/Modules/MDistil/Baryon2pt.hpp>
 #include <Hadrons/Modules/MDistil/BContraction.hpp>
 #include <Hadrons/Modules/MDistil/PerambLight.hpp>
+#include <Hadrons/Modules/MDistil/BC2.hpp>
 #include <Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
 #include <Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
 #include <Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>

Hadrons/modules.inc

@@ -10,6 +10,7 @@ modules_cc =\
   Modules/MSolver/RBPrecCG.cc \
   Modules/MDistil/g5_multiply.cc \
   Modules/MDistil/BContraction.cc \
+  Modules/MDistil/BC2.cc \
   Modules/MDistil/LapEvec.cc \
   Modules/MDistil/Baryon2pt.cc \
   Modules/MDistil/PerambLight.cc \
@@ -94,6 +95,7 @@ modules_hpp =\
   Modules/MDistil/Baryon2pt.hpp \
   Modules/MDistil/BContraction.hpp \
   Modules/MDistil/PerambLight.hpp \
+  Modules/MDistil/BC2.hpp \
   Modules/MContraction/WeakHamiltonian.hpp \
   Modules/MContraction/WeakNeutral4ptDisc.hpp \
   Modules/MContraction/WeakHamiltonianEye.hpp \

tests/hadrons/Test_hadrons_distil.cc

@@ -286,6 +286,38 @@ void test_MesonFieldRhoAll(Application &application)
   application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRhoAll",A2AMesonFieldPar);
 }
 /////////////////////////////////////////////////////////////
+// BaryonFields - phiphiphi - efficient
+/////////////////////////////////////////////////////////////
+
+void test_BaryonFieldPhi2(Application &application)
+{
+  // DistilVectors parameters
+  MDistil::BC2::Par BC2Par;
+  BC2Par.one="DistilVecs_phi";
+  BC2Par.two="DistilVecs_phi";
+  BC2Par.three="DistilVecs_phi";
+  BC2Par.output="BaryonFieldPhi2";
+  BC2Par.parity=1;
+  BC2Par.mom={"0 0 0"};
+  application.createModule<MDistil::BC2>("BaryonFieldPhi2",BC2Par);
+}
+/////////////////////////////////////////////////////////////
+// BaryonFields - rhorhorho - efficient
+/////////////////////////////////////////////////////////////
+
+void test_BaryonFieldRho2(Application &application)
+{
+  // DistilVectors parameters
+  MDistil::BC2::Par BC2Par;
+  BC2Par.one="DistilVecs_rho";
+  BC2Par.two="DistilVecs_rho";
+  BC2Par.three="DistilVecs_rho";
+  BC2Par.output="BaryonFieldRho2";
+  BC2Par.parity=1;
+  BC2Par.mom={"0 0 0"};
+  application.createModule<MDistil::BC2>("BaryonFieldRho2",BC2Par);
+}
+/////////////////////////////////////////////////////////////
 // BaryonFields - phiphiphi
 /////////////////////////////////////////////////////////////
 
@@ -836,6 +868,14 @@ int main(int argc, char *argv[])
       test_em( application );
       test_Aslash( application );
       break;
+    case 11: // 3
+      test_Global( application );
+      test_LapEvec( application );
+      test_Perambulators( application );
+      test_DistilVectors( application );
+      test_BaryonFieldPhi2( application );
+      test_BaryonFieldRho2( application );
+      break;
   }
   LOG(Message) << "====== XML creation for test " << iTestNum << " complete ======" << std::endl;
 

tests/hadrons/Test_tesseract.cc

@@ -285,6 +285,38 @@ void test_MesonFieldRhoAll(Application &application)
   application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRhoAll",A2AMesonFieldPar);
 }
 /////////////////////////////////////////////////////////////
+// BaryonFields - phiphiphi - efficient
+/////////////////////////////////////////////////////////////
+
+void test_BaryonFieldPhi2(Application &application)
+{
+  // DistilVectors parameters
+  MDistil::BC2::Par BC2Par;
+  BC2Par.one="DistilVecs_phi";
+  BC2Par.two="DistilVecs_phi";
+  BC2Par.three="DistilVecs_phi";
+  BC2Par.output="BaryonFieldPhi2";
+  BC2Par.parity=1;
+  BC2Par.mom={"0 0 0"};
+  application.createModule<MDistil::BC2>("BaryonFieldPhi2",BC2Par);
+}
+/////////////////////////////////////////////////////////////
+// BaryonFields - rhorhorho - efficient
+/////////////////////////////////////////////////////////////
+
+void test_BaryonFieldRho2(Application &application)
+{
+  // DistilVectors parameters
+  MDistil::BC2::Par BC2Par;
+  BC2Par.one="DistilVecs_rho";
+  BC2Par.two="DistilVecs_rho";
+  BC2Par.three="DistilVecs_rho";
+  BC2Par.output="BaryonFieldRho2";
+  BC2Par.parity=1;
+  BC2Par.mom={"0 0 0"};
+  application.createModule<MDistil::BC2>("BaryonFieldRho2",BC2Par);
+}
+/////////////////////////////////////////////////////////////
 // BaryonFields - phiphiphi
 /////////////////////////////////////////////////////////////