Changes nedded to have a current insertion on every second time slice - avoids unnecessary contractions

Grid/qcd/utils/A2Autils.h

@@ -69,14 +69,18 @@ public:
         int orthogdim, double *t_kernel = nullptr, double *t_gsum = nullptr);
 
   template <typename TensorType>
-  typename std::enable_if<std::is_same<Eigen::Tensor<ComplexD,3>, TensorType>::value, void>::type
+  typename std::enable_if<(std::is_same<Eigen::Tensor<ComplexD,3>, TensorType>::value ||
+                           std::is_same<Eigen::TensorMap<Eigen::Tensor<Complex, 3, Eigen::RowMajor>>, TensorType>::value),
+                           void>::type
   static ContractWWVV(std::vector<PropagatorField> &WWVV,
 			   const TensorType &WW_sd,
 			   const FermionField *vs,
 			   const FermionField *vd);
 
   template <typename TensorType>
-  typename std::enable_if<!std::is_same<Eigen::Tensor<ComplexD,3>, TensorType>::value, void>::type
+  typename std::enable_if<!(std::is_same<Eigen::Tensor<ComplexD,3>, TensorType>::value ||
+                            std::is_same<Eigen::TensorMap<Eigen::Tensor<Complex, 3, Eigen::RowMajor>>, TensorType>::value),
+                            void>::type
   static ContractWWVV(std::vector<PropagatorField> &WWVV,
 			   const TensorType &WW_sd,
 			   const FermionField *vs,
@@ -977,7 +981,9 @@ void A2Autils<FImpl>::AslashField(TensorType &mat,
 
 template <class FImpl>
 template <typename TensorType>
-typename std::enable_if<std::is_same<Eigen::Tensor<ComplexD,3>, TensorType>::value, void>::type
+typename std::enable_if<(std::is_same<Eigen::Tensor<ComplexD,3>, TensorType>::value ||
+                         std::is_same<Eigen::TensorMap<Eigen::Tensor<Complex, 3, Eigen::RowMajor>>, TensorType>::value),
+                         void>::type
 A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
 				   const TensorType &WW_sd,
 				   const FermionField *vs,
@@ -1023,7 +1029,9 @@ A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
 
 template <class FImpl>
 template <typename TensorType>
-typename std::enable_if<!std::is_same<Eigen::Tensor<ComplexD,3>, TensorType>::value, void>::type
+typename std::enable_if<!(std::is_same<Eigen::Tensor<ComplexD, 3>, TensorType>::value ||
+                          std::is_same<Eigen::TensorMap<Eigen::Tensor<Complex, 3, Eigen::RowMajor>>, TensorType>::value),
+                          void>::type
 A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
                               const TensorType &WW_sd,
                               const FermionField *vs,

Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp

@@ -19,7 +19,9 @@ public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(A2AFourQuarkContractionPar,
                                     std::string,  v1,
                                     std::string,  v2,
-                                    std::string, mf12);
+                                    std::string,  mf12,
+                                    bool,         allContr,
+                                    unsigned int, dt);
 };
 
 template <typename FImpl>
@@ -38,6 +40,8 @@ class TA2AFourQuarkContraction: public Module<A2AFourQuarkContractionPar>
     virtual void setup(void);
     // execution
     virtual void execute(void);
+  private:
+    unsigned int nt_;
 };
 
 MODULE_REGISTER_TMP(A2AFourQuarkContraction, TA2AFourQuarkContraction<FIMPL>, MContraction);
@@ -72,30 +76,59 @@ std::vector<std::string> TA2AFourQuarkContraction<FImpl>::getOutput(void)
 template <typename FImpl>
 void TA2AFourQuarkContraction<FImpl>::setup(void)
 {
-    int nt = env().getDim(Tp);
+    if (par().allContr)
-
+    {
-    envCreate(std::vector<PropagatorField>, getName(), 1, nt, envGetGrid(PropagatorField));
+        nt_ = env().getDim(Tp);
+        envTmp(std::vector<PropagatorField>, "tmpWWVV", 1, nt_, envGetGrid(PropagatorField));
+        envCreate(std::vector<PropagatorField>, getName(), 1, nt_, envGetGrid(PropagatorField));
+    }
+    else
+    {
+        envTmp(std::vector<PropagatorField>, "tmpWWVV", 1, 1, envGetGrid(PropagatorField));
+        envCreate(PropagatorField, getName(), 1, envGetGrid(PropagatorField));
+    }
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
 void TA2AFourQuarkContraction<FImpl>::execute(void)
 {
-    int nt = env().getDim(Tp);
-
     auto &v1   = envGet(std::vector<FermionField>, par().v1);
     auto &v2   = envGet(std::vector<FermionField>, par().v2);
     auto &mf12 = envGet(EigenDiskVector<Complex>, par().mf12);
 
-    auto &wwvv = envGet(std::vector<PropagatorField>, getName());
+    envGetTmp(std::vector<PropagatorField>, tmpWWVV);
 
-    for (int t = 0; t < nt; t++)
+    unsigned int dt = par().dt;
+    unsigned int nt = env().getDim(Tp);
+
+    if (par().allContr)
     {
-        wwvv[t] = zero;
+        LOG(Message) << "Computing 4 quark contraction for " << getName()
-    }
+                     << " for all t0 time translations "
+                     << "with nt = " << nt_ << " and dt = " << dt << std::endl;
 
-    LOG(Message) << "Computing 4 quark contraction for: " << getName() << std::endl;
+        auto &WWVV = envGet(std::vector<PropagatorField>, getName());
-    A2Autils<FImpl>::ContractWWVV(wwvv, mf12, &v1[0], &v2[0]);
+        A2Autils<FImpl>::ContractWWVV(tmpWWVV, mf12, &v1[0], &v2[0]);
+        for(unsigned int t = 0; t < nt_; t++){
+            unsigned int t0 = (t + dt) % nt_;
+            WWVV[t] = tmpWWVV[t0];
+        }
+    }
+    else
+    {
+        LOG(Message) << "Computing 4 quark contraction for: " << getName()
+                     << " for time dt = " << dt << std::endl;
+
+        auto &WWVV = envGet(PropagatorField, getName());
+        int ni = v1.size();
+        int nj = v2.size();
+        Eigen::Matrix<Complex, -1, -1, Eigen::RowMajor> mf;
+        mf = mf12[dt];
+        Eigen::TensorMap<Eigen::Tensor<Complex, 3, Eigen::RowMajor>> mfT(mf.data(), 1, ni, nj);
+        A2Autils<FImpl>::ContractWWVV(tmpWWVV, mfT, &v1[0], &v2[0]);
+        WWVV = tmpWWVV[0];
+    }
 }
 
 END_MODULE_NAMESPACE