Merge branch 'feature/distil' of github.com:mmphys/Grid into feature/distil

Grid/qcd/utils/A2Autils.h

@@ -178,23 +178,30 @@ void A2Autils<FImpl>::NucleonFieldMom(Eigen::Tensor<ComplexD,6> &mat,
 	for(int i=0;i<oneBlock;i++){
 
 	  auto v1 = one[i]._odata[ss];
-          assert(0 && "Sorry, Felix, the next line was stopping the build");
+	  auto gv1 = Gamma(Gamma::Algebra::GammaT)*v1;
-	  //auto pv1 = 0.5*(v1 + (double)parity*Gamma(Gamma::Algebra::GammaT)*v1);
+	  auto pv1 = 0.5*(v1 + (double)parity*gv1);
-          auto pv1 = v1;
 
        	  for(int j=0;j<twoBlock;j++){
 
 	    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;  
 
 	    for(int k=0;k<threeBlock;k++){
 
 	      auto v3 = three[k]._odata[ss];
-              // C = i gamma_2 gamma_4 => C gamma_5 = - i gamma_1 gamma_3
 	      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
+                  -           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);    */
 	        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)    
@@ -208,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()()());
+		}
               }
 	    }
 	  }
@@ -222,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++){
@@ -232,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++){
 
@@ -243,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/MDistil/DistilVectors.hpp

@@ -85,7 +85,7 @@ std::vector<std::string> TDistilVectors<FImpl>::getInput(void)
 template <typename FImpl>
 std::vector<std::string> TDistilVectors<FImpl>::getOutput(void)
 {
-    std::vector<std::string> out = {getName() + "_rho", getName() + "_phi"};
+    std::vector<std::string> out = {getName() + "_rho", getName() + "_phi", getName() + "_rho_all_tsrc"};
     
     return out;
 }
@@ -106,6 +106,8 @@ void TDistilVectors<FImpl>::setup(void)
 		                    nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField));
    envCreate(std::vector<FermionField>, getName() + "_phi", 1, 
                  	            nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField)); 
+   envCreate(std::vector<FermionField>, getName() + "_rho_all_tsrc", 1, 
+		                    nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField));
 
 
   GridCartesian * grid4d = env().getGrid();
@@ -138,6 +140,7 @@ void TDistilVectors<FImpl>::execute(void)
     auto        &epack   = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
     auto        &rho       = envGet(std::vector<FermionField>, getName() + "_rho");
     auto        &phi       = envGet(std::vector<FermionField>, getName() + "_phi");
+    auto        &rho_all       = envGet(std::vector<FermionField>, getName() + "_rho_all_tsrc");
 
 
   envGetTmp(LatticeSpinColourVector, tmp2);
@@ -189,16 +192,38 @@ void TDistilVectors<FImpl>::execute(void)
               }
             }
           }
-          SpinColourVector scv0;
-          std::vector<int> siteFirst(grid4d->Nd(),0);
-          peekSite(scv0, rho[vecindex], siteFirst);
-          auto & cplx0 = scv0()(0)(0);
-          std::cout <<  "site0[rho(vecindex = " << vecindex << ")] = " << cplx0 << std::endl;
         }
       }
     }
   }
 
+  for (int inoise = 0; inoise < nnoise; inoise++) {
+    for (int dk = 0; dk < LI; dk++) {
+      for (int dt = 0; dt < Nt_inv; dt++) {
+        for (int ds = 0; ds < Ns; ds++) {
+          vecindex = inoise + nnoise * dk + nnoise * LI * ds + nnoise *LI * Ns*dt;
+          rho_all[vecindex] = zero;
+          tmp3d_nospin = zero;
+          for (int it = 0; it < Nt; it++){
+            t_inv = it;
+            if( t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal ) {
+              for (int ik = dk; ik < nvec; ik += LI){
+                for (int is = ds; is < Ns; is += Ns){ //at the moment, full spin dilution is enforced
+                  ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv,3);
+                  tmp3d_nospin = evec3d * noise[inoise + nnoise*(t_inv + Nt*(ik+nvec*is))];
+                  tmp3d=zero;
+                  pokeSpin(tmp3d,tmp3d_nospin,is);
+                  tmp2=zero;
+                  InsertSliceLocal(tmp3d,tmp2,0,t_inv-Ntfirst,Grid::QCD::Tdir);
+                  rho_all[vecindex] += tmp2;
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
 
   for (int inoise = 0; inoise < nnoise; inoise++) {
     for (int dk = 0; dk < LI; dk++) {
@@ -214,11 +239,6 @@ void TDistilVectors<FImpl>::execute(void)
             }
             InsertSliceLocal(sink_tslice,phi[vecindex],0,t-Ntfirst,Grid::QCD::Tdir);
           }
-          SpinColourVector scv0;
-          std::vector<int> siteFirst(grid4d->Nd(),0);
-          peekSite(scv0, phi[vecindex], siteFirst);
-          auto & cplx0 = scv0()(0)(0);
-          std::cout <<  "site0[phi(vecindex = " << vecindex << ")] = " << cplx0 << std::endl;
         }
       }
     }

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

@@ -269,6 +269,55 @@ void test_MesonFieldRho(Application &application)
   application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRho",A2AMesonFieldPar);
 }
 /////////////////////////////////////////////////////////////
+// MesonFields - rhorhoall
+/////////////////////////////////////////////////////////////
+
+void test_MesonFieldRhoAll(Application &application)
+{
+  // DistilVectors parameters
+  MContraction::A2AMesonField::Par A2AMesonFieldPar;
+  A2AMesonFieldPar.left="DistilVecs_rho_all_tsrc";
+  A2AMesonFieldPar.right="DistilVecs_rho_all_tsrc";
+  A2AMesonFieldPar.output="MesonSinksRhoAll";
+  A2AMesonFieldPar.gammas="all";
+  A2AMesonFieldPar.mom={"0 0 0"};
+  A2AMesonFieldPar.cacheBlock=2;
+  A2AMesonFieldPar.block=4;
+  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
 /////////////////////////////////////////////////////////////
 
@@ -774,6 +823,7 @@ int main(int argc, char *argv[])
       test_Perambulators( application );
       test_DistilVectors( application );
       test_MesonField( application );
+      test_MesonFieldRhoAll( application );
       break;
     case 5: // 3
       test_Global( application );
@@ -818,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

@@ -268,6 +268,55 @@ void test_MesonFieldRho(Application &application)
   application.createModule<MContraction::A2AMesonField>("DistilMesonFieldRho",A2AMesonFieldPar);
 }
 /////////////////////////////////////////////////////////////
+// MesonFields - rhorhoAll
+/////////////////////////////////////////////////////////////
+
+void test_MesonFieldRhoAll(Application &application)
+{
+  // DistilVectors parameters
+  MContraction::A2AMesonField::Par A2AMesonFieldPar;
+  A2AMesonFieldPar.left="DistilVecs_rho_all_tsrc";
+  A2AMesonFieldPar.right="DistilVecs_rho_all_tsrc";
+  A2AMesonFieldPar.output="MesonSinksRhoAll";
+  A2AMesonFieldPar.gammas="all";
+  A2AMesonFieldPar.mom={"0 0 0"};
+  A2AMesonFieldPar.cacheBlock=2;
+  A2AMesonFieldPar.block=4;
+  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
 /////////////////////////////////////////////////////////////
 
@@ -757,6 +806,7 @@ int main(int argc, char *argv[])
       test_DistilVectors( application );
       test_MesonField( application );
       test_MesonFieldRho( application );
+      test_MesonFieldRhoAll( application );
       break;
   }
   LOG(Message) << "====== XML creation for test " << iTestNum << " complete ======" << std::endl;