From ac530636cabdb18ac9360fb6677c277bf6333305 Mon Sep 17 00:00:00 2001
From: fionnoh <fionnoh@gmail.com>
Date: Thu, 27 Jun 2019 13:44:47 +0800
Subject: [PATCH 1/8] A2Aloop bugfix

---
 Hadrons/Modules/MContraction/A2ALoop.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Hadrons/Modules/MContraction/A2ALoop.hpp b/Hadrons/Modules/MContraction/A2ALoop.hpp
index 2ef99354..7e7ffac0 100644
--- a/Hadrons/Modules/MContraction/A2ALoop.hpp
+++ b/Hadrons/Modules/MContraction/A2ALoop.hpp
@@ -112,7 +112,7 @@ void TA2ALoop<FImpl>::execute(void)
     loop = zero;
     for (unsigned int i = 0; i < left.size(); ++i)
     {
-        loop += outerProduct(adj(left[i]), right[i]);
+        loop += outerProduct(left[i], right[i]);
     }
 }
 

From 421a0a8a36015e72da769a7fa4ee647afa4ee326 Mon Sep 17 00:00:00 2001
From: fionnoh <fionnoh@gmail.com>
Date: Thu, 27 Jun 2019 13:45:20 +0800
Subject: [PATCH 2/8] Changes to A2Autils, A2AMatirx and DiskVector code that
 is needed for Hadrons 4 quark contraction module

---
 Grid/qcd/utils/A2Autils.h | 118 +++++++++++++++----
 Hadrons/A2AMatrix.hpp     |  81 ++++++++-----
 Hadrons/DiskVector.hpp    | 233 +++++++++++++++++++++++---------------
 3 files changed, 290 insertions(+), 142 deletions(-)

diff --git a/Grid/qcd/utils/A2Autils.h b/Grid/qcd/utils/A2Autils.h
index 97188ffe..966d84f7 100644
--- a/Grid/qcd/utils/A2Autils.h
+++ b/Grid/qcd/utils/A2Autils.h
@@ -68,8 +68,17 @@ public:
         const std::vector<ComplexField> &emB1,
         int orthogdim, double *t_kernel = nullptr, double *t_gsum = nullptr);
 
-  static void ContractWWVV(std::vector<PropagatorField> &WWVV,
-			   const Eigen::Tensor<ComplexD,3> &WW_sd,
+  template <typename TensorType>
+  typename std::enable_if<std::is_same<Eigen::Tensor<ComplexD,3>, 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
+  static ContractWWVV(std::vector<PropagatorField> &WWVV,
+			   const TensorType &WW_sd,
 			   const FermionField *vs,
 			   const FermionField *vd);
 
@@ -99,6 +108,11 @@ public:
 			const FermionField *vd,
 			int orthogdim);
 #endif
+private:
+  inline static void OuterProductWWVV(std::vector<PropagatorField> &WWVV,
+                               const vobj &lhs,
+                               const vobj &rhs,
+                               const int Ns, const int ss, const int t);
 };
 
 template <class FImpl>
@@ -961,12 +975,15 @@ void A2Autils<FImpl>::AslashField(TensorType &mat,
 // Take WW_sd v^dag_d (x) v_s
 // 
 
-template<class FImpl>
-void A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
-				   const Eigen::Tensor<ComplexD,3> &WW_sd,
+template <class FImpl>
+template <typename TensorType>
+typename std::enable_if<std::is_same<Eigen::Tensor<ComplexD,3>, TensorType>::value, void>::type
+A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
+				   const TensorType &WW_sd,
 				   const FermionField *vs,
 				   const FermionField *vd)
 {
+  std::cout << "Start contraction" << std::endl;
   GridBase *grid = vs[0]._grid;
 
   int nd    = grid->_ndimension;
@@ -989,33 +1006,90 @@ void A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
   vobj tmp2 = zero;
   vobj tmp3 = zero;
 
-	for(int d=d_o;d<MIN(d_o+d_unroll,N_d);d++){
+  for(int d=d_o;d<MIN(d_o+d_unroll,N_d);d++){
 	  Scalar_v coeff = WW_sd(t,s,d);
-	  tmp3 = conjugate(vd[d]._odata[ss]);
-	  mac(&tmp2, &coeff, &tmp3);
-  }
+    tmp3 = conjugate(vd[d]._odata[ss]);
+    mac(&tmp2 ,& coeff, &tmp3 );
+	}
 
 	//////////////////////////
 	// Fast outer product of tmp1 with a sum of terms suppressed by d_unroll
 	//////////////////////////
-	for(int s1=0;s1<Ns;s1++){
-	for(int s2=0;s2<Ns;s2++){
-	  WWVV[t]._odata[ss]()(s1,s2)(0,0) += tmp1()(s1)(0)*tmp2()(s2)(0);
-	  WWVV[t]._odata[ss]()(s1,s2)(0,1) += tmp1()(s1)(0)*tmp2()(s2)(1);
-	  WWVV[t]._odata[ss]()(s1,s2)(0,2) += tmp1()(s1)(0)*tmp2()(s2)(2);
-	  WWVV[t]._odata[ss]()(s1,s2)(1,0) += tmp1()(s1)(1)*tmp2()(s2)(0);
-	  WWVV[t]._odata[ss]()(s1,s2)(1,1) += tmp1()(s1)(1)*tmp2()(s2)(1);
-	  WWVV[t]._odata[ss]()(s1,s2)(1,2) += tmp1()(s1)(1)*tmp2()(s2)(2);
-	  WWVV[t]._odata[ss]()(s1,s2)(2,0) += tmp1()(s1)(2)*tmp2()(s2)(0);
-	  WWVV[t]._odata[ss]()(s1,s2)(2,1) += tmp1()(s1)(2)*tmp2()(s2)(1);
-	  WWVV[t]._odata[ss]()(s1,s2)(2,2) += tmp1()(s1)(2)*tmp2()(s2)(2);
-	}}
-
+  OuterProductWWVV(WWVV, tmp1, tmp2, Ns, ss, t);
       }}
     }
   }
 }
 
+template <class FImpl>
+template <typename TensorType>
+typename std::enable_if<!std::is_same<Eigen::Tensor<ComplexD,3>, TensorType>::value, void>::type
+A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
+                                   const TensorType &WW_sd,
+                                   const FermionField *vs,
+                                   const FermionField *vd)
+{
+  GridBase *grid = vs[0]._grid;
+
+  int nd    = grid->_ndimension;
+  int Nsimd = grid->Nsimd();
+  int N_t = WW_sd.dimensions()[0];
+  int N_s = WW_sd.dimensions()[1];
+  int N_d = WW_sd.dimensions()[2];
+
+  int d_unroll = 32;// Empirical optimisation
+
+  Eigen::Matrix<Complex, -1, -1, Eigen::RowMajor> buf;
+
+  for(int t=0;t<N_t;t++){
+    WWVV[t] = zero;
+  }
+
+  for (int t = 0; t < N_t; t++){
+    std::cout << GridLogMessage << "Contraction t = " << t << std::endl;
+    buf = WW_sd[t];
+    parallel_for(int ss=0;ss<grid->oSites();ss++){
+      for(int d_o=0;d_o<N_d;d_o+=d_unroll){
+        for(int s=0;s<N_s;s++){
+    auto tmp1 = vs[s]._odata[ss];
+    vobj tmp2 = zero;
+    vobj tmp3 = zero;
+
+    for(int d=d_o;d<MIN(d_o+d_unroll,N_d);d++){
+      Scalar_v coeff = buf(s,d);
+      tmp3 = conjugate(vd[d]._odata[ss]);
+      mac(&tmp2 ,& coeff, &tmp3 );
+    }
+
+    //////////////////////////
+    // Fast outer product of tmp1 with a sum of terms suppressed by d_unroll
+    //////////////////////////
+    OuterProductWWVV(WWVV, tmp1, tmp2, Ns, ss, t);
+      }}
+    }
+  }
+}
+
+template <class FImpl>
+inline void A2Autils<FImpl>::OuterProductWWVV(std::vector<PropagatorField> &WWVV,
+                                             const vobj &lhs,
+                                             const vobj &rhs,
+                                             const int Ns, const int ss, const int t)
+{
+  for (int s1 = 0; s1 < Ns; s1++){
+    for (int s2 = 0; s2 < Ns; s2++){
+      WWVV[t]._odata[ss]()(s1, s2)(0, 0) += lhs()(s1)(0) * rhs()(s2)(0);
+      WWVV[t]._odata[ss]()(s1, s2)(0, 1) += lhs()(s1)(0) * rhs()(s2)(1);
+      WWVV[t]._odata[ss]()(s1, s2)(0, 2) += lhs()(s1)(0) * rhs()(s2)(2);
+      WWVV[t]._odata[ss]()(s1, s2)(1, 0) += lhs()(s1)(1) * rhs()(s2)(0);
+      WWVV[t]._odata[ss]()(s1, s2)(1, 1) += lhs()(s1)(1) * rhs()(s2)(1);
+      WWVV[t]._odata[ss]()(s1, s2)(1, 2) += lhs()(s1)(1) * rhs()(s2)(2);
+      WWVV[t]._odata[ss]()(s1, s2)(2, 0) += lhs()(s1)(2) * rhs()(s2)(0);
+      WWVV[t]._odata[ss]()(s1, s2)(2, 1) += lhs()(s1)(2) * rhs()(s2)(1);
+      WWVV[t]._odata[ss]()(s1, s2)(2, 2) += lhs()(s1)(2) * rhs()(s2)(2);
+    }
+  }
+}
 
 template<class FImpl>
 void A2Autils<FImpl>::ContractFourQuarkColourDiagonal(const PropagatorField &WWVV0,
diff --git a/Hadrons/A2AMatrix.hpp b/Hadrons/A2AMatrix.hpp
index 29b30a7a..50ba4848 100644
--- a/Hadrons/A2AMatrix.hpp
+++ b/Hadrons/A2AMatrix.hpp
@@ -108,7 +108,7 @@ public:
     void saveBlock(const A2AMatrixSet<T> &m, const unsigned int ext, const unsigned int str,
                    const unsigned int i, const unsigned int j);
     template <template <class> class Vec, typename VecT>
-    void load(Vec<VecT> &v, double *tRead = nullptr);
+    void load(Vec<VecT> &v, double *tRead = nullptr, GridBase *grid = nullptr);
 private:
     std::string  filename_{""}, dataname_{""};
     unsigned int nt_{0}, ni_{0}, nj_{0};
@@ -495,44 +495,53 @@ void A2AMatrixIo<T>::saveBlock(const A2AMatrixSet<T> &m,
 
 template <typename T>
 template <template <class> class Vec, typename VecT>
-void A2AMatrixIo<T>::load(Vec<VecT> &v, double *tRead)
+void A2AMatrixIo<T>::load(Vec<VecT> &v, double *tRead, GridBase *grid)
 {
 #ifdef HAVE_HDF5
-    Hdf5Reader           reader(filename_);
     std::vector<hsize_t> hdim;
     H5NS::DataSet        dataset;
     H5NS::DataSpace      dataspace;
     H5NS::CompType       datatype;
-    
-    push(reader, dataname_);
-    auto &group = reader.getGroup();
-    dataset     = group.openDataSet(HADRONS_A2AM_NAME);
-    datatype    = dataset.getCompType();
-    dataspace   = dataset.getSpace();
-    hdim.resize(dataspace.getSimpleExtentNdims());
-    dataspace.getSimpleExtentDims(hdim.data());
-    if ((nt_*ni_*nj_ != 0) and
-        ((hdim[0] != nt_) or (hdim[1] != ni_) or (hdim[2] != nj_)))
+
+    if (!(grid) || grid->IsBoss())
     {
-        HADRONS_ERROR(Size, "all-to-all matrix size mismatch (got "
-            + std::to_string(hdim[0]) + "x" + std::to_string(hdim[1]) + "x"
-            + std::to_string(hdim[2]) + ", expected "
-            + std::to_string(nt_) + "x" + std::to_string(ni_) + "x"
-            + std::to_string(nj_));
-    }
-    else if (ni_*nj_ == 0)
-    {
-        if (hdim[0] != nt_)
+        Hdf5Reader reader(filename_);
+        push(reader, dataname_);
+        auto &group = reader.getGroup();
+        dataset = group.openDataSet(HADRONS_A2AM_NAME);
+        datatype = dataset.getCompType();
+        dataspace = dataset.getSpace();
+        hdim.resize(dataspace.getSimpleExtentNdims());
+        dataspace.getSimpleExtentDims(hdim.data());
+        if ((nt_ * ni_ * nj_ != 0) and
+            ((hdim[0] != nt_) or (hdim[1] != ni_) or (hdim[2] != nj_)))
         {
-            HADRONS_ERROR(Size, "all-to-all time size mismatch (got "
-                + std::to_string(hdim[0]) + ", expected "
-                + std::to_string(nt_) + ")");
+            HADRONS_ERROR(Size, "all-to-all matrix size mismatch (got "
+                + std::to_string(hdim[0]) + "x" + std::to_string(hdim[1]) + "x"
+                + std::to_string(hdim[2]) + ", expected "
+                + std::to_string(nt_) + "x" + std::to_string(ni_) + "x"
+                + std::to_string(nj_));
         }
-        ni_ = hdim[1];
-        nj_ = hdim[2];
+        else if (ni_*nj_ == 0)
+        {
+            if (hdim[0] != nt_)
+            {
+                HADRONS_ERROR(Size, "all-to-all time size mismatch (got "
+                    + std::to_string(hdim[0]) + ", expected "
+                    + std::to_string(nt_) + ")");
+            }
+            ni_ = hdim[1];
+            nj_ = hdim[2];
+        }
+    }
+    if (grid)
+    {
+        grid->Broadcast(grid->BossRank(), &ni_, sizeof(unsigned int));
+        grid->Broadcast(grid->BossRank(), &nj_, sizeof(unsigned int));
     }
 
     A2AMatrix<T>         buf(ni_, nj_);
+    int broadcastSize =  sizeof(T) * buf.size();
     std::vector<hsize_t> count    = {1, static_cast<hsize_t>(ni_),
                                      static_cast<hsize_t>(nj_)},
                          stride   = {1, 1, 1},
@@ -554,10 +563,20 @@ void A2AMatrixIo<T>::load(Vec<VecT> &v, double *tRead)
             std::cout << " " << t;
             std::cout.flush();
         }
-        dataspace.selectHyperslab(H5S_SELECT_SET, count.data(), offset.data(),
-                                  stride.data(), block.data());
-        if (tRead) *tRead -= usecond();    
-        dataset.read(buf.data(), datatype, memspace, dataspace);
+        if (!(grid) || grid->IsBoss())
+        {
+            dataspace.selectHyperslab(H5S_SELECT_SET, count.data(), offset.data(),
+                                      stride.data(), block.data());
+        }
+        if (tRead) *tRead -= usecond();
+        if (!(grid) || grid->IsBoss())
+        {
+            dataset.read(buf.data(), datatype, memspace, dataspace);
+        }
+        if (grid)
+        {
+            grid->Broadcast(grid->BossRank(), buf.data(), broadcastSize);
+        }
         if (tRead) *tRead += usecond();
         v[t] = buf.template cast<VecT>();
     }
diff --git a/Hadrons/DiskVector.hpp b/Hadrons/DiskVector.hpp
index 3a8d0696..dd5956aa 100644
--- a/Hadrons/DiskVector.hpp
+++ b/Hadrons/DiskVector.hpp
@@ -87,13 +87,20 @@ public:
     };
 public:
     DiskVectorBase(const std::string dirname, const unsigned int size = 0,
-                   const unsigned int cacheSize = 1, const bool clean = true);
+                   const unsigned int cacheSize = 1, const bool clean = true,
+                   GridBase *grid = nullptr);
     DiskVectorBase(DiskVectorBase<T> &&v) = default;
     virtual ~DiskVectorBase(void);
     const T & operator[](const unsigned int i) const;
     RwAccessHelper operator[](const unsigned int i);
     double hitRatio(void) const;
     void resetStat(void);
+    void setSize(unsigned int size_);
+    unsigned int getSize() const;
+    unsigned int dvSize;
+    void setGrid(GridBase *grid_);
+    GridBase *getGrid() const;
+    GridBase *dvGrid;
 private:
     virtual void load(T &obj, const std::string filename) const = 0;
     virtual void save(const std::string filename, const T &obj) const = 0;
@@ -107,6 +114,7 @@ private:
     unsigned int                                          size_, cacheSize_;
     double                                                access_{0.}, hit_{0.};
     bool                                                  clean_;
+    GridBase                                              *grid_;
     // using pointers to allow modifications when class is const
     // semantic: const means data unmodified, but cache modification allowed
     std::unique_ptr<std::vector<T>>                       cachePtr_;
@@ -158,66 +166,92 @@ public:
     {
         return (*this)[i](j, k);
     }
+    std::vector<int> dimensions() const
+    {
+        std::vector<int> dims(3);
+        dims[0] = (*this).getSize();
+        dims[1] = (*this)[0].rows();
+        dims[2] = (*this)[0].cols();
+        return dims;
+    }
 private:
     virtual void load(EigenDiskVectorMat<T> &obj, const std::string filename) const
     {
-        std::ifstream f(filename, std::ios::binary);
-        uint32_t      crc, check;
-        Eigen::Index  nRow, nCol;
-        size_t        matSize;
-        double        tRead, tHash;
-
-        f.read(reinterpret_cast<char *>(&crc), sizeof(crc));
-        f.read(reinterpret_cast<char *>(&nRow), sizeof(nRow));
-        f.read(reinterpret_cast<char *>(&nCol), sizeof(nCol));
-        obj.resize(nRow, nCol);
-        matSize = nRow*nCol*sizeof(T);
-        tRead  = -usecond();
-        f.read(reinterpret_cast<char *>(obj.data()), matSize);
-        tRead += usecond();
-        tHash  = -usecond();
-#ifdef USE_IPP
-        check  = GridChecksum::crc32c(obj.data(), matSize);
-#else
-        check  = GridChecksum::crc32(obj.data(), matSize);
-#endif
-        tHash += usecond();
-        DV_DEBUG_MSG(this, "Eigen read " << tRead/1.0e6 << " sec " << matSize/tRead*1.0e6/1024/1024 << " MB/s");
-        DV_DEBUG_MSG(this, "Eigen crc32 " << std::hex << check << std::dec 
-                     << " " << tHash/1.0e6 << " sec " << matSize/tHash*1.0e6/1024/1024 << " MB/s");
-        if (crc != check)
+        GridBase *loadGrid;
+        loadGrid = (*this).getGrid();
+        if (!(loadGrid) || loadGrid->IsBoss())
         {
-            HADRONS_ERROR(Io, "checksum failed")
+            std::ifstream f(filename, std::ios::binary);
+            uint32_t      crc, check;
+            Eigen::Index  nRow, nCol;
+            size_t        matSize;
+            double        tRead, tHash;
+
+            f.read(reinterpret_cast<char *>(&crc), sizeof(crc));
+            f.read(reinterpret_cast<char *>(&nRow), sizeof(nRow));
+            f.read(reinterpret_cast<char *>(&nCol), sizeof(nCol));
+            obj.resize(nRow, nCol);
+            matSize = nRow*nCol*sizeof(T);
+            tRead  = -usecond();
+            f.read(reinterpret_cast<char *>(obj.data()), matSize);
+            tRead += usecond();
+            tHash  = -usecond();
+    #ifdef USE_IPP
+            check  = GridChecksum::crc32c(obj.data(), matSize);
+    #else
+            check  = GridChecksum::crc32(obj.data(), matSize);
+    #endif
+            tHash += usecond();
+            DV_DEBUG_MSG(this, "Eigen read " << tRead/1.0e6 << " sec " << matSize/tRead*1.0e6/1024/1024 << " MB/s");
+            DV_DEBUG_MSG(this, "Eigen crc32 " << std::hex << check << std::dec 
+                        << " " << tHash/1.0e6 << " sec " << matSize/tHash*1.0e6/1024/1024 << " MB/s");
+            if (crc != check)
+            {
+                HADRONS_ERROR(Io, "checksum failed")
+            }
+        }
+        int broadcastSize;
+        broadcastSize = sizeof(T)*obj.size();
+        if (loadGrid)
+        {
+            loadGrid->Broadcast(loadGrid->BossRank(), obj.data(), broadcastSize);
+            loadGrid->Barrier();
         }
     }
 
     virtual void save(const std::string filename, const EigenDiskVectorMat<T> &obj) const
     {
-        std::ofstream f(filename, std::ios::binary);
-        uint32_t      crc;
-        Eigen::Index  nRow, nCol;
-        size_t        matSize;
-        double        tWrite, tHash;
-        
-        nRow    = obj.rows();
-        nCol    = obj.cols();
-        matSize = nRow*nCol*sizeof(T);
-        tHash   = -usecond();
-#ifdef USE_IPP
-        crc     = GridChecksum::crc32c(obj.data(), matSize);
-#else
-        crc     = GridChecksum::crc32(obj.data(), matSize);
-#endif
-        tHash  += usecond();
-        f.write(reinterpret_cast<char *>(&crc), sizeof(crc));
-        f.write(reinterpret_cast<char *>(&nRow), sizeof(nRow));
-        f.write(reinterpret_cast<char *>(&nCol), sizeof(nCol));
-        tWrite = -usecond();
-        f.write(reinterpret_cast<const char *>(obj.data()), matSize);
-        tWrite += usecond();
-        DV_DEBUG_MSG(this, "Eigen write " << tWrite/1.0e6 << " sec " << matSize/tWrite*1.0e6/1024/1024 << " MB/s");
-        DV_DEBUG_MSG(this, "Eigen crc32 " << std::hex << crc << std::dec
-                     << " " << tHash/1.0e6 << " sec " << matSize/tHash*1.0e6/1024/1024 << " MB/s");
+        GridBase *saveGrid;
+        saveGrid = (*this).getGrid();
+        if (!(saveGrid) || saveGrid->IsBoss())
+        {
+            std::ofstream f(filename, std::ios::binary);
+            uint32_t      crc;
+            Eigen::Index  nRow, nCol;
+            size_t        matSize;
+            double        tWrite, tHash;
+            
+            nRow    = obj.rows();
+            nCol    = obj.cols();
+            matSize = nRow*nCol*sizeof(T);
+            tHash   = -usecond();
+    #ifdef USE_IPP
+            crc     = GridChecksum::crc32c(obj.data(), matSize);
+    #else
+            crc     = GridChecksum::crc32(obj.data(), matSize);
+    #endif
+            tHash  += usecond();
+            f.write(reinterpret_cast<char *>(&crc), sizeof(crc));
+            f.write(reinterpret_cast<char *>(&nRow), sizeof(nRow));
+            f.write(reinterpret_cast<char *>(&nCol), sizeof(nCol));
+            tWrite = -usecond();
+            f.write(reinterpret_cast<const char *>(obj.data()), matSize);
+            tWrite += usecond();
+            DV_DEBUG_MSG(this, "Eigen write " << tWrite/1.0e6 << " sec " << matSize/tWrite*1.0e6/1024/1024 << " MB/s");
+            DV_DEBUG_MSG(this, "Eigen crc32 " << std::hex << crc << std::dec
+                        << " " << tHash/1.0e6 << " sec " << matSize/tHash*1.0e6/1024/1024 << " MB/s");
+        }
+        if (saveGrid)   saveGrid->Barrier();
     }
 };
 
@@ -228,8 +262,9 @@ template <typename T>
 DiskVectorBase<T>::DiskVectorBase(const std::string dirname, 
                                   const unsigned int size,
                                   const unsigned int cacheSize,
-                                  const bool clean)
-: dirname_(dirname), size_(size), cacheSize_(cacheSize), clean_(clean)
+                                  const bool clean,
+                                  GridBase *grid)
+: dirname_(dirname), size_(size), cacheSize_(cacheSize), clean_(clean), grid_(grid)
 , cachePtr_(new std::vector<T>(size))
 , modifiedPtr_(new std::vector<bool>(size, false))
 , indexPtr_(new std::map<unsigned int, unsigned int>())
@@ -238,15 +273,21 @@ DiskVectorBase<T>::DiskVectorBase(const std::string dirname,
 {
     struct stat s;
 
-    if(stat(dirname.c_str(), &s) == 0)
+    if (!(grid_) || grid_->IsBoss())
     {
-        HADRONS_ERROR(Io, "directory '" + dirname + "' already exists")
+        if(stat(dirname.c_str(), &s) == 0)
+        {
+            HADRONS_ERROR(Io, "directory '" + dirname + "' already exists")
+        }
+        mkdir(dirname);
     }
-    mkdir(dirname);
+    if (grid_)  grid_->Barrier();
     for (unsigned int i = 0; i < cacheSize_; ++i)
     {
         freePtr_->push(i);
     }
+    setSize(size_);
+    setGrid(grid_);
 }
 
 template <typename T>
@@ -258,6 +299,30 @@ DiskVectorBase<T>::~DiskVectorBase(void)
     }
 }
 
+template <typename T>
+void DiskVectorBase<T>::setSize(unsigned int size_)
+{
+    dvSize = size_;
+}
+
+template <typename T>
+unsigned int DiskVectorBase<T>::getSize() const
+{
+    return dvSize;
+}
+
+template <typename T>
+void DiskVectorBase<T>::setGrid(GridBase *grid_)
+{
+    dvGrid = grid_;
+}
+
+template <typename T>
+GridBase *DiskVectorBase<T>::getGrid() const
+{
+    return dvGrid;
+}
+
 template <typename T>
 const T & DiskVectorBase<T>::operator[](const unsigned int i) const
 {
@@ -299,7 +364,7 @@ const T & DiskVectorBase<T>::operator[](const unsigned int i) const
     }
     DV_DEBUG_MSG(this, "in cache: " << msg);
 #endif
-
+    if (grid_)  grid_->Barrier();
     return cache[index.at(i)];
 }
 
@@ -358,6 +423,7 @@ void DiskVectorBase<T>::evict(void) const
         index.erase(i);
         loads.pop_front();
     }
+    if (grid_)  grid_->Barrier();
 }
 
 template <typename T>
@@ -395,27 +461,14 @@ void DiskVectorBase<T>::cacheInsert(const unsigned int i, const T &obj) const
     auto &freeInd  = *freePtr_;
     auto &loads    = *loadsPtr_;
 
-    // cache miss, evict and store
-    if (index.find(i) == index.end())
-    {
-        evict();
-        index[i] = freeInd.top();
-        freeInd.pop();
-        cache[index.at(i)] = obj;
-        loads.push_back(i);
-        modified[index.at(i)] = false;
-    }
-    // cache hit, modify current value
-    else
-    {
-        auto pos = std::find(loads.begin(), loads.end(), i);
-        
-        cache[index.at(i)]    = obj;
-        modified[index.at(i)] = true;
-        loads.erase(pos);
-        loads.push_back(i);
-    }
+    evict();
+    index[i] = freeInd.top();
+    freeInd.pop();
+    cache[index.at(i)] = obj;
+    loads.push_back(i);
+    modified[index.at(i)] = false;
 
+    if (grid_)  grid_->Barrier();
 #ifdef DV_DEBUG
     std::string msg;
 
@@ -434,21 +487,23 @@ void DiskVectorBase<T>::cacheInsert(const unsigned int i, const T &obj) const
 template <typename T>
 void DiskVectorBase<T>::clean(void)
 {
-    auto unlink = [](const char *fpath, const struct stat *sb, 
-                     int typeflag, struct FTW *ftwbuf)
+    if (!(grid_) || grid_->IsBoss())
     {
-        int rv = remove(fpath);
+        auto unlink = [](const char *fpath, const struct stat *sb,
+                         int typeflag, struct FTW *ftwbuf) {
+            int rv = remove(fpath);
 
-        if (rv)
-        {
-            HADRONS_ERROR(Io, "cannot remove '" + std::string(fpath) + "': "
-                          + std::string(std::strerror(errno)));
-        }
+            if (rv)
+            {
+                HADRONS_ERROR(Io, "cannot remove '" + std::string(fpath) + "': " + std::string(std::strerror(errno)));
+            }
 
-        return rv;
-    };
+            return rv;
+        };
 
-    nftw(dirname_.c_str(), unlink, 64, FTW_DEPTH | FTW_PHYS);
+        nftw(dirname_.c_str(), unlink, 64, FTW_DEPTH | FTW_PHYS);
+    }
+    if (grid_)  grid_->Barrier();
 }
 
 END_HADRONS_NAMESPACE

From ce29b18dc9971dac36e20ddfb4ef35ea748500a0 Mon Sep 17 00:00:00 2001
From: fionnoh <fionnoh@gmail.com>
Date: Thu, 27 Jun 2019 13:46:06 +0800
Subject: [PATCH 3/8] New modules for loading in MFs as diskvectors and
 producing propagaotrs from 4 quark contractions

---
 Hadrons/Modules.hpp                           | 120 ++++-----
 .../MContraction/A2AFourQuarkContraction.cc   |   7 +
 .../MContraction/A2AFourQuarkContraction.hpp  | 105 ++++++++
 .../Modules/MIO/LoadA2AMatrixDiskVector.cc    |   7 +
 .../Modules/MIO/LoadA2AMatrixDiskVector.hpp   | 115 +++++++++
 Hadrons/modules.inc                           | 232 +++++++++---------
 6 files changed, 413 insertions(+), 173 deletions(-)
 create mode 100644 Hadrons/Modules/MContraction/A2AFourQuarkContraction.cc
 create mode 100644 Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp
 create mode 100644 Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.cc
 create mode 100644 Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.hpp

diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp
index 898cb532..4d5e9db4 100644
--- a/Hadrons/Modules.hpp
+++ b/Hadrons/Modules.hpp
@@ -1,69 +1,71 @@
-#include <Hadrons/Modules/MContraction/WeakEye3pt.hpp>
-#include <Hadrons/Modules/MContraction/Baryon.hpp>
-#include <Hadrons/Modules/MContraction/A2AAslashField.hpp>
-#include <Hadrons/Modules/MContraction/A2ALoop.hpp>
-#include <Hadrons/Modules/MContraction/A2AMesonField.hpp>
-#include <Hadrons/Modules/MContraction/Meson.hpp>
-#include <Hadrons/Modules/MContraction/DiscLoop.hpp>
-#include <Hadrons/Modules/MContraction/Gamma3pt.hpp>
-#include <Hadrons/Modules/MContraction/WeakMesonDecayKl2.hpp>
-#include <Hadrons/Modules/MContraction/WeakNonEye3pt.hpp>
-#include <Hadrons/Modules/MFermion/FreeProp.hpp>
-#include <Hadrons/Modules/MFermion/GaugeProp.hpp>
-#include <Hadrons/Modules/MFermion/EMLepton.hpp>
+#include <Hadrons/Modules/MSource/Momentum.hpp>
+#include <Hadrons/Modules/MSource/SeqConserved.hpp>
 #include <Hadrons/Modules/MSource/SeqGamma.hpp>
 #include <Hadrons/Modules/MSource/SeqAslash.hpp>
+#include <Hadrons/Modules/MSource/Z2.hpp>
 #include <Hadrons/Modules/MSource/Point.hpp>
 #include <Hadrons/Modules/MSource/Wall.hpp>
-#include <Hadrons/Modules/MSource/Z2.hpp>
-#include <Hadrons/Modules/MSource/SeqConserved.hpp>
-#include <Hadrons/Modules/MSource/Momentum.hpp>
 #include <Hadrons/Modules/MSink/Smear.hpp>
 #include <Hadrons/Modules/MSink/Point.hpp>
-#include <Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp>
-#include <Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
-#include <Hadrons/Modules/MSolver/A2AAslashVectors.hpp>
-#include <Hadrons/Modules/MSolver/Guesser.hpp>
-#include <Hadrons/Modules/MSolver/RBPrecCG.hpp>
-#include <Hadrons/Modules/MSolver/A2AVectors.hpp>
-#include <Hadrons/Modules/MGauge/UnitEm.hpp>
-#include <Hadrons/Modules/MGauge/StoutSmearing.hpp>
-#include <Hadrons/Modules/MGauge/Unit.hpp>
-#include <Hadrons/Modules/MGauge/Random.hpp>
-#include <Hadrons/Modules/MGauge/GaugeFix.hpp>
-#include <Hadrons/Modules/MGauge/FundtoHirep.hpp>
-#include <Hadrons/Modules/MGauge/StochEm.hpp>
-#include <Hadrons/Modules/MGauge/Electrify.hpp>
-#include <Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp>
-#include <Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp>
-#include <Hadrons/Modules/MUtilities/PrecisionCast.hpp>
-#include <Hadrons/Modules/MUtilities/RandomVectors.hpp>
-#include <Hadrons/Modules/MScalar/FreeProp.hpp>
-#include <Hadrons/Modules/MScalar/Scalar.hpp>
-#include <Hadrons/Modules/MScalar/ChargedProp.hpp>
-#include <Hadrons/Modules/MNPR/Bilinear.hpp>
-#include <Hadrons/Modules/MNPR/Amputate.hpp>
-#include <Hadrons/Modules/MNPR/FourQuark.hpp>
-#include <Hadrons/Modules/MAction/DWF.hpp>
-#include <Hadrons/Modules/MAction/MobiusDWF.hpp>
-#include <Hadrons/Modules/MAction/Wilson.hpp>
-#include <Hadrons/Modules/MAction/WilsonClover.hpp>
-#include <Hadrons/Modules/MAction/ZMobiusDWF.hpp>
-#include <Hadrons/Modules/MAction/ScaledDWF.hpp>
-#include <Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
+#include <Hadrons/Modules/MIO/LoadCosmHol.hpp>
+#include <Hadrons/Modules/MIO/LoadBinary.hpp>
+#include <Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.hpp>
+#include <Hadrons/Modules/MIO/LoadEigenPack.hpp>
+#include <Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
+#include <Hadrons/Modules/MIO/LoadA2AVectors.hpp>
+#include <Hadrons/Modules/MIO/LoadNersc.hpp>
+#include <Hadrons/Modules/MScalarSUN/Utils.hpp>
+#include <Hadrons/Modules/MScalarSUN/Grad.hpp>
+#include <Hadrons/Modules/MScalarSUN/TrPhi.hpp>
 #include <Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
+#include <Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
+#include <Hadrons/Modules/MScalarSUN/TransProj.hpp>
+#include <Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
+#include <Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
 #include <Hadrons/Modules/MScalarSUN/Div.hpp>
 #include <Hadrons/Modules/MScalarSUN/TrMag.hpp>
 #include <Hadrons/Modules/MScalarSUN/EMT.hpp>
-#include <Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
-#include <Hadrons/Modules/MScalarSUN/TrPhi.hpp>
-#include <Hadrons/Modules/MScalarSUN/Utils.hpp>
-#include <Hadrons/Modules/MScalarSUN/TransProj.hpp>
-#include <Hadrons/Modules/MScalarSUN/Grad.hpp>
-#include <Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
-#include <Hadrons/Modules/MIO/LoadEigenPack.hpp>
-#include <Hadrons/Modules/MIO/LoadNersc.hpp>
-#include <Hadrons/Modules/MIO/LoadA2AVectors.hpp>
-#include <Hadrons/Modules/MIO/LoadCosmHol.hpp>
-#include <Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
-#include <Hadrons/Modules/MIO/LoadBinary.hpp>
+#include <Hadrons/Modules/MAction/ZMobiusDWF.hpp>
+#include <Hadrons/Modules/MAction/ScaledDWF.hpp>
+#include <Hadrons/Modules/MAction/Wilson.hpp>
+#include <Hadrons/Modules/MAction/WilsonClover.hpp>
+#include <Hadrons/Modules/MAction/MobiusDWF.hpp>
+#include <Hadrons/Modules/MAction/DWF.hpp>
+#include <Hadrons/Modules/MContraction/DiscLoop.hpp>
+#include <Hadrons/Modules/MContraction/WeakEye3pt.hpp>
+#include <Hadrons/Modules/MContraction/Meson.hpp>
+#include <Hadrons/Modules/MContraction/A2AMesonField.hpp>
+#include <Hadrons/Modules/MContraction/Gamma3pt.hpp>
+#include <Hadrons/Modules/MContraction/WeakNonEye3pt.hpp>
+#include <Hadrons/Modules/MContraction/Baryon.hpp>
+#include <Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp>
+#include <Hadrons/Modules/MContraction/A2ALoop.hpp>
+#include <Hadrons/Modules/MContraction/WeakMesonDecayKl2.hpp>
+#include <Hadrons/Modules/MContraction/A2AAslashField.hpp>
+#include <Hadrons/Modules/MScalar/Scalar.hpp>
+#include <Hadrons/Modules/MScalar/FreeProp.hpp>
+#include <Hadrons/Modules/MScalar/ChargedProp.hpp>
+#include <Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp>
+#include <Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp>
+#include <Hadrons/Modules/MUtilities/RandomVectors.hpp>
+#include <Hadrons/Modules/MUtilities/PrecisionCast.hpp>
+#include <Hadrons/Modules/MFermion/FreeProp.hpp>
+#include <Hadrons/Modules/MFermion/EMLepton.hpp>
+#include <Hadrons/Modules/MFermion/GaugeProp.hpp>
+#include <Hadrons/Modules/MSolver/A2AAslashVectors.hpp>
+#include <Hadrons/Modules/MSolver/Guesser.hpp>
+#include <Hadrons/Modules/MSolver/A2AVectors.hpp>
+#include <Hadrons/Modules/MSolver/RBPrecCG.hpp>
+#include <Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
+#include <Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp>
+#include <Hadrons/Modules/MNPR/FourQuark.hpp>
+#include <Hadrons/Modules/MNPR/Bilinear.hpp>
+#include <Hadrons/Modules/MNPR/Amputate.hpp>
+#include <Hadrons/Modules/MGauge/StoutSmearing.hpp>
+#include <Hadrons/Modules/MGauge/GaugeFix.hpp>
+#include <Hadrons/Modules/MGauge/StochEm.hpp>
+#include <Hadrons/Modules/MGauge/FundtoHirep.hpp>
+#include <Hadrons/Modules/MGauge/Unit.hpp>
+#include <Hadrons/Modules/MGauge/Random.hpp>
+#include <Hadrons/Modules/MGauge/UnitEm.hpp>
+#include <Hadrons/Modules/MGauge/Electrify.hpp>
diff --git a/Hadrons/Modules/MContraction/A2AFourQuarkContraction.cc b/Hadrons/Modules/MContraction/A2AFourQuarkContraction.cc
new file mode 100644
index 00000000..9afe12ba
--- /dev/null
+++ b/Hadrons/Modules/MContraction/A2AFourQuarkContraction.cc
@@ -0,0 +1,7 @@
+#include <Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+template class Grid::Hadrons::MContraction::TA2AFourQuarkContraction<FIMPL>;
diff --git a/Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp b/Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp
new file mode 100644
index 00000000..f977caab
--- /dev/null
+++ b/Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp
@@ -0,0 +1,105 @@
+#ifndef Hadrons_MContraction_A2AFourQuarkContraction_hpp_
+#define Hadrons_MContraction_A2AFourQuarkContraction_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Hadrons/DiskVector.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         A2AFourQuarkContraction                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class A2AFourQuarkContractionPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(A2AFourQuarkContractionPar,
+                                    std::string, v1,
+                                    std::string, v2,
+                                    std::string, mf12);
+};
+
+template <typename FImpl>
+class TA2AFourQuarkContraction: public Module<A2AFourQuarkContractionPar>
+{
+  public:
+    FERM_TYPE_ALIASES(FImpl, );
+    // constructor
+    TA2AFourQuarkContraction(const std::string name);
+    // destructor
+    virtual ~TA2AFourQuarkContraction(void) {};
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(A2AFourQuarkContraction, TA2AFourQuarkContraction<FIMPL>, MContraction);
+
+/******************************************************************************
+ *                 TA2AFourQuarkContraction implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TA2AFourQuarkContraction<FImpl>::TA2AFourQuarkContraction(const std::string name)
+: Module<A2AFourQuarkContractionPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TA2AFourQuarkContraction<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().v1, par().v2, par().mf12};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TA2AFourQuarkContraction<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TA2AFourQuarkContraction<FImpl>::setup(void)
+{
+    int nt = env().getDim(Tp);
+
+    envCreate(std::vector<PropagatorField>, getName(), 1, nt, 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());
+
+    for (int t = 0; t < nt; t++)
+    {
+        wwvv[t] = zero;
+    }
+
+    LOG(Message) << "Computing 4 quark contraction for: " << getName() << std::endl;
+    A2Autils<FImpl>::ContractWWVV(wwvv, mf12, &v1[0], &v2[0]);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MContraction_A2AFourQuarkContraction_hpp_
diff --git a/Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.cc b/Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.cc
new file mode 100644
index 00000000..696e159c
--- /dev/null
+++ b/Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.cc
@@ -0,0 +1,7 @@
+#include <Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MIO;
+
+template class Grid::Hadrons::MIO::TLoadA2AMatrixDiskVector<FIMPL>;
diff --git a/Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.hpp b/Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.hpp
new file mode 100644
index 00000000..d7799a37
--- /dev/null
+++ b/Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.hpp
@@ -0,0 +1,115 @@
+#ifndef Hadrons_MIO_LoadA2AMatrixDiskVector_hpp_
+#define Hadrons_MIO_LoadA2AMatrixDiskVector_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Hadrons/A2AMatrix.hpp>
+#include <Hadrons/DiskVector.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         LoadA2AMatrixDiskVector                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MIO)
+
+class LoadA2AMatrixDiskVectorPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(LoadA2AMatrixDiskVectorPar,
+                                    std::string,  file,
+                                    std::string,  dataset,
+                                    std::string,  diskVectorDir,
+                                    int,  cacheSize);
+};
+
+template <typename FImpl>
+class TLoadA2AMatrixDiskVector: public Module<LoadA2AMatrixDiskVectorPar>
+{
+  public:
+    FERM_TYPE_ALIASES(FImpl, );
+    // constructor
+    TLoadA2AMatrixDiskVector(const std::string name);
+    // destructor
+    virtual ~TLoadA2AMatrixDiskVector(void) {};
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(LoadA2AMatrixDiskVector, TLoadA2AMatrixDiskVector<FIMPL>, MIO);
+
+/******************************************************************************
+ *                 TLoadA2AMatrixDiskVector implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TLoadA2AMatrixDiskVector<FImpl>::TLoadA2AMatrixDiskVector(const std::string name)
+: Module<LoadA2AMatrixDiskVectorPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TLoadA2AMatrixDiskVector<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TLoadA2AMatrixDiskVector<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TLoadA2AMatrixDiskVector<FImpl>::setup(void)
+{
+    int Ls = 1;
+
+    std::string dvDir = par().diskVectorDir;
+    std::string dataset = par().dataset;
+    std::string dvFile = dvDir + "/" + getName() + "." + std::to_string(vm().getTrajectory());
+
+    int nt = env().getDim(Tp);
+    int cacheSize = par().cacheSize;
+    bool clean = true;
+    GridBase *grid = envGetGrid(FermionField);
+
+    envCreate(EigenDiskVector<ComplexD>, getName(), Ls, dvFile, nt, cacheSize, clean, grid);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TLoadA2AMatrixDiskVector<FImpl>::execute(void)
+{
+    int nt = env().getDim(Tp);
+    std::string file  = par().file;
+    std::string dataset  = par().dataset;
+    GridBase *grid = envGetGrid(FermionField);
+
+    auto &mesonFieldDV = envGet(EigenDiskVector<ComplexD>, getName());
+
+    int traj = vm().getTrajectory();
+    tokenReplace(file, "traj", traj);
+    LOG(Message) << "-- Loading '" << file << "'-- " << std::endl;
+    double t;
+    A2AMatrixIo<HADRONS_A2AM_IO_TYPE> mfIO(file, dataset, nt);
+    mfIO.load(mesonFieldDV, &t, grid);
+    LOG(Message) << "Read " << mfIO.getSize() << " bytes in " << t << " usec, " << mfIO.getSize() / t * 1.0e6 / 1024 / 1024 << " MB/s" << std::endl;
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MIO_LoadA2AMatrixDiskVector_hpp_
diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc
index 3b9e6563..d607c4bc 100644
--- a/Hadrons/modules.inc
+++ b/Hadrons/modules.inc
@@ -1,139 +1,143 @@
 modules_cc =\
-  Modules/MContraction/Baryon.cc \
-  Modules/MContraction/Meson.cc \
-  Modules/MContraction/WeakMesonDecayKl2.cc \
-  Modules/MContraction/WeakEye3pt.cc \
-  Modules/MContraction/A2ALoop.cc \
-  Modules/MContraction/WeakNonEye3pt.cc \
-  Modules/MContraction/A2AAslashField.cc \
-  Modules/MContraction/A2AMesonField.cc \
-  Modules/MContraction/DiscLoop.cc \
-  Modules/MContraction/Gamma3pt.cc \
-  Modules/MFermion/FreeProp.cc \
-  Modules/MFermion/GaugeProp.cc \
-  Modules/MFermion/EMLepton.cc \
-  Modules/MSource/Momentum.cc \
-  Modules/MSource/Point.cc \
-  Modules/MSource/Wall.cc \
   Modules/MSource/SeqConserved.cc \
   Modules/MSource/SeqGamma.cc \
-  Modules/MSource/SeqAslash.cc \
+  Modules/MSource/Momentum.cc \
+  Modules/MSource/Point.cc \
   Modules/MSource/Z2.cc \
+  Modules/MSource/SeqAslash.cc \
+  Modules/MSource/Wall.cc \
   Modules/MSink/Point.cc \
   Modules/MSink/Smear.cc \
-  Modules/MSolver/A2AVectors.cc \
-  Modules/MSolver/A2AAslashVectors.cc \
-  Modules/MSolver/RBPrecCG.cc \
-  Modules/MSolver/MixedPrecisionRBPrecCG.cc \
-  Modules/MSolver/LocalCoherenceLanczos.cc \
-  Modules/MGauge/StoutSmearing.cc \
-  Modules/MGauge/Unit.cc \
-  Modules/MGauge/UnitEm.cc \
-  Modules/MGauge/StochEm.cc \
-  Modules/MGauge/Random.cc \
-  Modules/MGauge/Electrify.cc \
-  Modules/MGauge/FundtoHirep.cc \
-  Modules/MGauge/GaugeFix.cc \
-  Modules/MNoise/FullVolumeSpinColorDiagonal.cc \
+  Modules/MIO/LoadA2AMatrixDiskVector.cc \
+  Modules/MIO/LoadNersc.cc \
+  Modules/MIO/LoadEigenPack.cc \
+  Modules/MIO/LoadCosmHol.cc \
+  Modules/MIO/LoadA2AVectors.cc \
+  Modules/MIO/LoadCoarseEigenPack.cc \
+  Modules/MIO/LoadBinary.cc \
+  Modules/MScalarSUN/TwoPointNPR.cc \
+  Modules/MScalarSUN/TrPhi.cc \
+  Modules/MScalarSUN/StochFreeField.cc \
+  Modules/MScalarSUN/TrMag.cc \
+  Modules/MScalarSUN/Div.cc \
+  Modules/MScalarSUN/Grad.cc \
+  Modules/MScalarSUN/TwoPoint.cc \
+  Modules/MScalarSUN/EMT.cc \
+  Modules/MScalarSUN/TransProj.cc \
+  Modules/MScalarSUN/TrKinetic.cc \
+  Modules/MAction/DWF.cc \
+  Modules/MAction/MobiusDWF.cc \
+  Modules/MAction/ZMobiusDWF.cc \
+  Modules/MAction/Wilson.cc \
+  Modules/MAction/ScaledDWF.cc \
+  Modules/MAction/WilsonClover.cc \
+  Modules/MContraction/WeakEye3pt.cc \
+  Modules/MContraction/DiscLoop.cc \
+  Modules/MContraction/Baryon.cc \
+  Modules/MContraction/A2AMesonField.cc \
+  Modules/MContraction/Gamma3pt.cc \
+  Modules/MContraction/A2ALoop.cc \
+  Modules/MContraction/A2AAslashField.cc \
+  Modules/MContraction/Meson.cc \
+  Modules/MContraction/A2AFourQuarkContraction.cc \
+  Modules/MContraction/WeakMesonDecayKl2.cc \
+  Modules/MContraction/WeakNonEye3pt.cc \
+  Modules/MScalar/ChargedProp.cc \
+  Modules/MScalar/FreeProp.cc \
   Modules/MNoise/TimeDilutedSpinColorDiagonal.cc \
+  Modules/MNoise/FullVolumeSpinColorDiagonal.cc \
   Modules/MUtilities/RandomVectors.cc \
   Modules/MUtilities/PrecisionCast.cc \
-  Modules/MScalar/FreeProp.cc \
-  Modules/MScalar/ChargedProp.cc \
+  Modules/MFermion/EMLepton.cc \
+  Modules/MFermion/FreeProp.cc \
+  Modules/MFermion/GaugeProp.cc \
+  Modules/MSolver/RBPrecCG.cc \
+  Modules/MSolver/LocalCoherenceLanczos.cc \
+  Modules/MSolver/MixedPrecisionRBPrecCG.cc \
+  Modules/MSolver/A2AVectors.cc \
+  Modules/MSolver/A2AAslashVectors.cc \
   Modules/MNPR/Amputate.cc \
   Modules/MNPR/Bilinear.cc \
   Modules/MNPR/FourQuark.cc \
-  Modules/MAction/Wilson.cc \
-  Modules/MAction/MobiusDWF.cc \
-  Modules/MAction/ZMobiusDWF.cc \
-  Modules/MAction/WilsonClover.cc \
-  Modules/MAction/DWF.cc \
-  Modules/MAction/ScaledDWF.cc \
-  Modules/MScalarSUN/TrPhi.cc \
-  Modules/MScalarSUN/Grad.cc \
-  Modules/MScalarSUN/TrMag.cc \
-  Modules/MScalarSUN/TrKinetic.cc \
-  Modules/MScalarSUN/EMT.cc \
-  Modules/MScalarSUN/TransProj.cc \
-  Modules/MScalarSUN/StochFreeField.cc \
-  Modules/MScalarSUN/TwoPoint.cc \
-  Modules/MScalarSUN/TwoPointNPR.cc \
-  Modules/MScalarSUN/Div.cc \
-  Modules/MIO/LoadEigenPack.cc \
-  Modules/MIO/LoadBinary.cc \
-  Modules/MIO/LoadNersc.cc \
-  Modules/MIO/LoadCoarseEigenPack.cc \
-  Modules/MIO/LoadCosmHol.cc \
-  Modules/MIO/LoadA2AVectors.cc
+  Modules/MGauge/Unit.cc \
+  Modules/MGauge/Electrify.cc \
+  Modules/MGauge/Random.cc \
+  Modules/MGauge/GaugeFix.cc \
+  Modules/MGauge/UnitEm.cc \
+  Modules/MGauge/StochEm.cc \
+  Modules/MGauge/FundtoHirep.cc \
+  Modules/MGauge/StoutSmearing.cc
 
 modules_hpp =\
-  Modules/MContraction/WeakEye3pt.hpp \
-  Modules/MContraction/Baryon.hpp \
-  Modules/MContraction/A2AAslashField.hpp \
-  Modules/MContraction/A2ALoop.hpp \
-  Modules/MContraction/A2AMesonField.hpp \
-  Modules/MContraction/Meson.hpp \
-  Modules/MContraction/DiscLoop.hpp \
-  Modules/MContraction/Gamma3pt.hpp \
-  Modules/MContraction/WeakMesonDecayKl2.hpp \
-  Modules/MContraction/WeakNonEye3pt.hpp \
-  Modules/MFermion/FreeProp.hpp \
-  Modules/MFermion/GaugeProp.hpp \
-  Modules/MFermion/EMLepton.hpp \
+  Modules/MSource/Momentum.hpp \
+  Modules/MSource/SeqConserved.hpp \
   Modules/MSource/SeqGamma.hpp \
   Modules/MSource/SeqAslash.hpp \
+  Modules/MSource/Z2.hpp \
   Modules/MSource/Point.hpp \
   Modules/MSource/Wall.hpp \
-  Modules/MSource/Z2.hpp \
-  Modules/MSource/SeqConserved.hpp \
-  Modules/MSource/Momentum.hpp \
   Modules/MSink/Smear.hpp \
   Modules/MSink/Point.hpp \
-  Modules/MSolver/MixedPrecisionRBPrecCG.hpp \
-  Modules/MSolver/LocalCoherenceLanczos.hpp \
-  Modules/MSolver/A2AAslashVectors.hpp \
-  Modules/MSolver/Guesser.hpp \
-  Modules/MSolver/RBPrecCG.hpp \
-  Modules/MSolver/A2AVectors.hpp \
-  Modules/MGauge/UnitEm.hpp \
-  Modules/MGauge/StoutSmearing.hpp \
-  Modules/MGauge/Unit.hpp \
-  Modules/MGauge/Random.hpp \
-  Modules/MGauge/GaugeFix.hpp \
-  Modules/MGauge/FundtoHirep.hpp \
-  Modules/MGauge/StochEm.hpp \
-  Modules/MGauge/Electrify.hpp \
-  Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp \
-  Modules/MNoise/FullVolumeSpinColorDiagonal.hpp \
-  Modules/MUtilities/PrecisionCast.hpp \
-  Modules/MUtilities/RandomVectors.hpp \
-  Modules/MScalar/FreeProp.hpp \
-  Modules/MScalar/Scalar.hpp \
-  Modules/MScalar/ChargedProp.hpp \
-  Modules/MNPR/Bilinear.hpp \
-  Modules/MNPR/Amputate.hpp \
-  Modules/MNPR/FourQuark.hpp \
-  Modules/MAction/DWF.hpp \
-  Modules/MAction/MobiusDWF.hpp \
-  Modules/MAction/Wilson.hpp \
-  Modules/MAction/WilsonClover.hpp \
-  Modules/MAction/ZMobiusDWF.hpp \
-  Modules/MAction/ScaledDWF.hpp \
-  Modules/MScalarSUN/StochFreeField.hpp \
+  Modules/MIO/LoadCosmHol.hpp \
+  Modules/MIO/LoadBinary.hpp \
+  Modules/MIO/LoadA2AMatrixDiskVector.hpp \
+  Modules/MIO/LoadEigenPack.hpp \
+  Modules/MIO/LoadCoarseEigenPack.hpp \
+  Modules/MIO/LoadA2AVectors.hpp \
+  Modules/MIO/LoadNersc.hpp \
+  Modules/MScalarSUN/Utils.hpp \
+  Modules/MScalarSUN/Grad.hpp \
+  Modules/MScalarSUN/TrPhi.hpp \
   Modules/MScalarSUN/TwoPointNPR.hpp \
+  Modules/MScalarSUN/TwoPoint.hpp \
+  Modules/MScalarSUN/TransProj.hpp \
+  Modules/MScalarSUN/TrKinetic.hpp \
+  Modules/MScalarSUN/StochFreeField.hpp \
   Modules/MScalarSUN/Div.hpp \
   Modules/MScalarSUN/TrMag.hpp \
   Modules/MScalarSUN/EMT.hpp \
-  Modules/MScalarSUN/TwoPoint.hpp \
-  Modules/MScalarSUN/TrPhi.hpp \
-  Modules/MScalarSUN/Utils.hpp \
-  Modules/MScalarSUN/TransProj.hpp \
-  Modules/MScalarSUN/Grad.hpp \
-  Modules/MScalarSUN/TrKinetic.hpp \
-  Modules/MIO/LoadEigenPack.hpp \
-  Modules/MIO/LoadNersc.hpp \
-  Modules/MIO/LoadA2AVectors.hpp \
-  Modules/MIO/LoadCosmHol.hpp \
-  Modules/MIO/LoadCoarseEigenPack.hpp \
-  Modules/MIO/LoadBinary.hpp
+  Modules/MAction/ZMobiusDWF.hpp \
+  Modules/MAction/ScaledDWF.hpp \
+  Modules/MAction/Wilson.hpp \
+  Modules/MAction/WilsonClover.hpp \
+  Modules/MAction/MobiusDWF.hpp \
+  Modules/MAction/DWF.hpp \
+  Modules/MContraction/DiscLoop.hpp \
+  Modules/MContraction/WeakEye3pt.hpp \
+  Modules/MContraction/Meson.hpp \
+  Modules/MContraction/A2AMesonField.hpp \
+  Modules/MContraction/Gamma3pt.hpp \
+  Modules/MContraction/WeakNonEye3pt.hpp \
+  Modules/MContraction/Baryon.hpp \
+  Modules/MContraction/A2AFourQuarkContraction.hpp \
+  Modules/MContraction/A2ALoop.hpp \
+  Modules/MContraction/WeakMesonDecayKl2.hpp \
+  Modules/MContraction/A2AAslashField.hpp \
+  Modules/MScalar/Scalar.hpp \
+  Modules/MScalar/FreeProp.hpp \
+  Modules/MScalar/ChargedProp.hpp \
+  Modules/MNoise/FullVolumeSpinColorDiagonal.hpp \
+  Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp \
+  Modules/MUtilities/RandomVectors.hpp \
+  Modules/MUtilities/PrecisionCast.hpp \
+  Modules/MFermion/FreeProp.hpp \
+  Modules/MFermion/EMLepton.hpp \
+  Modules/MFermion/GaugeProp.hpp \
+  Modules/MSolver/A2AAslashVectors.hpp \
+  Modules/MSolver/Guesser.hpp \
+  Modules/MSolver/A2AVectors.hpp \
+  Modules/MSolver/RBPrecCG.hpp \
+  Modules/MSolver/LocalCoherenceLanczos.hpp \
+  Modules/MSolver/MixedPrecisionRBPrecCG.hpp \
+  Modules/MNPR/FourQuark.hpp \
+  Modules/MNPR/Bilinear.hpp \
+  Modules/MNPR/Amputate.hpp \
+  Modules/MGauge/StoutSmearing.hpp \
+  Modules/MGauge/GaugeFix.hpp \
+  Modules/MGauge/StochEm.hpp \
+  Modules/MGauge/FundtoHirep.hpp \
+  Modules/MGauge/Unit.hpp \
+  Modules/MGauge/Random.hpp \
+  Modules/MGauge/UnitEm.hpp \
+  Modules/MGauge/Electrify.hpp
 

From 67690df3bdaa6c32f082798c99272c102ffa77f1 Mon Sep 17 00:00:00 2001
From: fionnoh <fionnoh@gmail.com>
Date: Fri, 28 Jun 2019 15:18:28 +0800
Subject: [PATCH 4/8] Changes nedded to have a current insertion on every
 second time slice - avoids unnecessary contractions

---
 Grid/qcd/utils/A2Autils.h                     | 22 ++++---
 .../MContraction/A2AFourQuarkContraction.hpp  | 61 ++++++++++++++-----
 2 files changed, 62 insertions(+), 21 deletions(-)

diff --git a/Grid/qcd/utils/A2Autils.h b/Grid/qcd/utils/A2Autils.h
index 966d84f7..e163b9c6 100644
--- a/Grid/qcd/utils/A2Autils.h
+++ b/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,11 +1029,13 @@ 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,
-                                   const FermionField *vd)
+                              const TensorType &WW_sd,
+                              const FermionField *vs,
+                              const FermionField *vd)
 {
   GridBase *grid = vs[0]._grid;
 
diff --git a/Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp b/Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp
index f977caab..7a8112e1 100644
--- a/Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp
+++ b/Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp
@@ -17,9 +17,11 @@ class A2AFourQuarkContractionPar: Serializable
 {
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(A2AFourQuarkContractionPar,
-                                    std::string, v1,
-                                    std::string, v2,
-                                    std::string, mf12);
+                                    std::string,  v1,
+                                    std::string,  v2,
+                                    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);
-
-    envCreate(std::vector<PropagatorField>, getName(), 1, nt, envGetGrid(PropagatorField));
+    if (par().allContr)
+    {
+        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;
-    A2Autils<FImpl>::ContractWWVV(wwvv, mf12, &v1[0], &v2[0]);
+        auto &WWVV = envGet(std::vector<PropagatorField>, getName());
+        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

From 25150eb2e056303fc7668d69c287e3731c570c1a Mon Sep 17 00:00:00 2001
From: Fionn O hOgain <fionn.o.hogain@ed.ac.uk>
Date: Thu, 15 Aug 2019 12:09:30 +0100
Subject: [PATCH 5/8] 3pt contraction now takes a list of gammas

---
 Hadrons/Modules/MContraction/Gamma3pt.hpp | 76 ++++++++++++++++-------
 1 file changed, 54 insertions(+), 22 deletions(-)

diff --git a/Hadrons/Modules/MContraction/Gamma3pt.hpp b/Hadrons/Modules/MContraction/Gamma3pt.hpp
index d36e184b..006d7413 100644
--- a/Hadrons/Modules/MContraction/Gamma3pt.hpp
+++ b/Hadrons/Modules/MContraction/Gamma3pt.hpp
@@ -57,7 +57,8 @@ BEGIN_HADRONS_NAMESPACE
  *   - q1: sink smeared propagator, source at i
  *   - q2: propagator, source at i
  *   - q3: propagator, source at f
- *   - gamma: gamma matrix to insert
+ *   - gammas: gamma matrices to insert
+ *             (space-separated strings e.g. "GammaT GammaX GammaY") 
  *   - tSnk: sink position for propagator q1.
  *
  */
@@ -71,12 +72,12 @@ class Gamma3ptPar: Serializable
 {
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(Gamma3ptPar,
-                                    std::string,    q1,
-                                    std::string,    q2,
-                                    std::string,    q3,
-                                    Gamma::Algebra, gamma,
-                                    unsigned int,   tSnk,
-                                    std::string,    output);
+                                    std::string,  q1,
+                                    std::string,  q2,
+                                    std::string,  q3,
+                                    std::string,  gammas,
+                                    unsigned int, tSnk,
+                                    std::string,  output);
 };
 
 template <typename FImpl1, typename FImpl2, typename FImpl3>
@@ -100,6 +101,7 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+    virtual void parseGammaString(std::vector<Gamma::Algebra> &gammaList);
 protected:
     // setup
     virtual void setup(void);
@@ -142,37 +144,67 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::setup(void)
     envTmpLat(LatticeComplex, "c");
 }
 
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+void TGamma3pt<FImpl1, FImpl2, FImpl3>::parseGammaString(std::vector<Gamma::Algebra> &gammaList)
+{
+    gammaList.clear();
+    // Determine gamma matrices to insert at source/sink.
+    if (par().gammas.compare("all") == 0)
+    {
+        // Do all contractions.
+        for (unsigned int i = 1; i < Gamma::nGamma; i += 2)
+        {
+            gammaList.push_back((Gamma::Algebra)i);
+        }
+    }
+    else
+    {
+        // Parse individual contractions from input string.
+        gammaList = strToVec<Gamma::Algebra>(par().gammas);
+    } 
+}
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl1, typename FImpl2, typename FImpl3>
 void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
 {
     LOG(Message) << "Computing 3pt contractions '" << getName() << "' using"
                  << " quarks '" << par().q1 << "', '" << par().q2 << "' and '"
-                 << par().q3 << "', with " << par().gamma << " insertion." 
+                 << par().q3 << "', with " << par().gammas << " insertions." 
                  << std::endl;
 
     // Initialise variables. q2 and q3 are normal propagators, q1 may be 
     // sink smeared.
-    auto                  &q1 = envGet(SlicedPropagator1, par().q1);
-    auto                  &q2 = envGet(PropagatorField2, par().q2);
-    auto                  &q3 = envGet(PropagatorField2, par().q3);
-    Gamma                 g5(Gamma::Algebra::Gamma5);
-    Gamma                 gamma(par().gamma);
-    std::vector<TComplex> buf;
-    Result                result;
+    auto                        &q1 = envGet(SlicedPropagator1, par().q1);
+    auto                        &q2 = envGet(PropagatorField2, par().q2);
+    auto                        &q3 = envGet(PropagatorField2, par().q3);
+    Gamma                       g5(Gamma::Algebra::Gamma5);
+    std::vector<Gamma::Algebra> gammaList;
+    std::vector<TComplex>       buf;
+    std::vector<Result>         result;
+    int                         nt = env().getDim(Tp);
+
+
+    parseGammaString(gammaList);
+    result.resize(gammaList.size());
+    for (unsigned int i = 0; i < result.size(); ++i)
+    {
+        result[i].gamma = gammaList[i];
+        result[i].corr.resize(nt);
+    }
     
     // Extract relevant timeslice of sinked propagator q1, then contract &
     // sum over all spacial positions of gamma insertion.
     SitePropagator1 q1Snk = q1[par().tSnk];
     envGetTmp(LatticeComplex, c);
-    c = trace(g5*q1Snk*adj(q2)*(g5*gamma)*q3);
-    sliceSum(c, buf, Tp);
-
-    result.gamma = par().gamma;
-    result.corr.resize(buf.size());
-    for (unsigned int t = 0; t < buf.size(); ++t)
+    for (unsigned int i = 0; i < result.size(); ++i)
     {
-        result.corr[t] = TensorRemove(buf[t]);
+        Gamma gamma(gammaList[i]);
+        c = trace(g5*q1Snk*adj(q2)*(g5*gamma)*q3);
+        sliceSum(c, buf, Tp);
+        for (unsigned int t = 0; t < buf.size(); ++t)
+        {
+            result[i].corr[t] = TensorRemove(buf[t]);
+        }
     }
     saveResult(par().output, "gamma3pt", result);
 }

From 10de4bfc23372bbab2e4a9295428f7bc551b357e Mon Sep 17 00:00:00 2001
From: Fionn O hOgain <fionn.o.hogain@ed.ac.uk>
Date: Tue, 8 Oct 2019 13:24:56 +0100
Subject: [PATCH 6/8] Changes to A2Autils

---
 Grid/qcd/utils/A2Autils.h | 32 ++++++++++++++++----------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/Grid/qcd/utils/A2Autils.h b/Grid/qcd/utils/A2Autils.h
index f2f26aa5..b603eede 100644
--- a/Grid/qcd/utils/A2Autils.h
+++ b/Grid/qcd/utils/A2Autils.h
@@ -112,7 +112,7 @@ public:
 			int orthogdim);
 #endif
 private:
-  inline static void OuterProductWWVV(std::vector<PropagatorField> &WWVV,
+  inline static void OuterProductWWVV(PropagatorField &WWVV,
                                const vobj &lhs,
                                const vobj &rhs,
                                const int Ns, const int ss);
@@ -1028,8 +1028,7 @@ A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
   //////////////////////////
   // Fast outer product of tmp1 with a sum of terms suppressed by d_unroll
   //////////////////////////
-  auto WWVV_v = WWVV[t].View();
-  OuterProductWWVV(WWVV_v, tmp1, tmp2, Ns, ss);
+  OuterProductWWVV(WWVV[t], tmp1, tmp2, Ns, ss);
 
       }}
     }
@@ -1046,7 +1045,7 @@ A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
                               const FermionField *vs,
                               const FermionField *vd)
 {
-  GridBase *grid = vs[0]._grid;
+  GridBase *grid = vs[0].Grid();
 
   int nd    = grid->_ndimension;
   int Nsimd = grid->Nsimd();
@@ -1059,7 +1058,7 @@ A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
   Eigen::Matrix<Complex, -1, -1, Eigen::RowMajor> buf;
 
   for(int t=0;t<N_t;t++){
-    WWVV[t] = zero;
+    WWVV[t] = Zero();
   }
 
   for (int t = 0; t < N_t; t++){
@@ -1082,7 +1081,7 @@ A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
     //////////////////////////
     // Fast outer product of tmp1 with a sum of terms suppressed by d_unroll
     //////////////////////////
-    OuterProductWWVV(WWVV, tmp1, tmp2, Ns, ss);
+    OuterProductWWVV(WWVV[t], tmp1, tmp2, Ns, ss);
       }}
     });
   }
@@ -1092,19 +1091,20 @@ template <class FImpl>
 inline void A2Autils<FImpl>::OuterProductWWVV(PropagatorField &WWVV,
                                              const vobj &lhs,
                                              const vobj &rhs,
-                                             const int Ns, const int ss, const int t)
+                                             const int Ns, const int ss)
 {
+  auto WWVV_v = WWVV.View();
   for (int s1 = 0; s1 < Ns; s1++){
     for (int s2 = 0; s2 < Ns; s2++){
-      WWVV[ss]()(s1,s2)(0, 0) += lhs()(s1)(0) * rhs()(s2)(0);
-      WWVV[ss]()(s1,s2)(0, 1) += lhs()(s1)(0) * rhs()(s2)(1);
-      WWVV[ss]()(s1,s2)(0, 2) += lhs()(s1)(0) * rhs()(s2)(2);
-      WWVV[ss]()(s1,s2)(1, 0) += lhs()(s1)(1) * rhs()(s2)(0);
-      WWVV[ss]()(s1,s2)(1, 1) += lhs()(s1)(1) * rhs()(s2)(1);
-      WWVV[ss]()(s1,s2)(1, 2) += lhs()(s1)(1) * rhs()(s2)(2);
-      WWVV[ss]()(s1,s2)(2, 0) += lhs()(s1)(2) * rhs()(s2)(0);
-      WWVV[ss]()(s1,s2)(2, 1) += lhs()(s1)(2) * rhs()(s2)(1);
-      WWVV[ss]()(s1,s2)(2, 2) += lhs()(s1)(2) * rhs()(s2)(2);
+      WWVV_v[ss]()(s1,s2)(0, 0) += lhs()(s1)(0) * rhs()(s2)(0);
+      WWVV_v[ss]()(s1,s2)(0, 1) += lhs()(s1)(0) * rhs()(s2)(1);
+      WWVV_v[ss]()(s1,s2)(0, 2) += lhs()(s1)(0) * rhs()(s2)(2);
+      WWVV_v[ss]()(s1,s2)(1, 0) += lhs()(s1)(1) * rhs()(s2)(0);
+      WWVV_v[ss]()(s1,s2)(1, 1) += lhs()(s1)(1) * rhs()(s2)(1);
+      WWVV_v[ss]()(s1,s2)(1, 2) += lhs()(s1)(1) * rhs()(s2)(2);
+      WWVV_v[ss]()(s1,s2)(2, 0) += lhs()(s1)(2) * rhs()(s2)(0);
+      WWVV_v[ss]()(s1,s2)(2, 1) += lhs()(s1)(2) * rhs()(s2)(1);
+      WWVV_v[ss]()(s1,s2)(2, 2) += lhs()(s1)(2) * rhs()(s2)(2);
     }
   }
 }

From 5de9547db52ee8530eb20aa24105ca97800b0385 Mon Sep 17 00:00:00 2001
From: Fionn O hOgain <fionn.o.hogain@ed.ac.uk>
Date: Tue, 8 Oct 2019 15:51:28 +0100
Subject: [PATCH 7/8] Removing old debug code

---
 Grid/lattice/Lattice_reduction.h | 110 -------------------------------
 1 file changed, 110 deletions(-)

diff --git a/Grid/lattice/Lattice_reduction.h b/Grid/lattice/Lattice_reduction.h
index 3a076038..3c5b03e5 100644
--- a/Grid/lattice/Lattice_reduction.h
+++ b/Grid/lattice/Lattice_reduction.h
@@ -317,116 +317,6 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
   }
 }
 
-template<class vobj>
-static void mySliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim)
-{
-  // std::cout << GridLogMessage << "Start mySliceInnerProductVector" << std::endl;
-
-  typedef typename vobj::scalar_type scalar_type;
-  std::vector<scalar_type> lsSum;
-  localSliceInnerProductVector(result, lhs, rhs, lsSum, orthogdim);
-  globalSliceInnerProductVector(result, lhs, lsSum, orthogdim);
-  // std::cout << GridLogMessage << "End mySliceInnerProductVector" << std::endl;
-}
-
-template <class vobj>
-static void localSliceInnerProductVector(std::vector<ComplexD> &result, const Lattice<vobj> &lhs, const Lattice<vobj> &rhs, std::vector<typename vobj::scalar_type> &lsSum, int orthogdim)
-{
-  // std::cout << GridLogMessage << "Start prep" << std::endl;
-  typedef typename vobj::vector_type   vector_type;
-  typedef typename vobj::scalar_type   scalar_type;
-  GridBase  *grid = lhs.Grid();
-  assert(grid!=NULL);
-  conformable(grid,rhs.Grid());
-
-  const int    Nd = grid->_ndimension;
-  const int Nsimd = grid->Nsimd();
-
-  assert(orthogdim >= 0);
-  assert(orthogdim < Nd);
-
-  int fd=grid->_fdimensions[orthogdim];
-  int ld=grid->_ldimensions[orthogdim];
-  int rd=grid->_rdimensions[orthogdim];
-  // std::cout << GridLogMessage << "Start alloc" << std::endl;
-
-  Vector<vector_type> lvSum(rd); // will locally sum vectors first
-  lsSum.resize(ld,scalar_type(0.0));                    // sum across these down to scalars
-  ExtractBuffer<iScalar<scalar_type> > extracted(Nsimd);   // splitting the SIMD  
-  // std::cout << GridLogMessage << "End alloc" << std::endl;
-
-  result.resize(fd); // And then global sum to return the same vector to every node for IO to file
-  for(int r=0;r<rd;r++){
-    lvSum[r]=Zero();
-  }
-
-  int e1=    grid->_slice_nblock[orthogdim];
-  int e2=    grid->_slice_block [orthogdim];
-  int stride=grid->_slice_stride[orthogdim];
-  // std::cout << GridLogMessage << "End prep" << std::endl;
-  // std::cout << GridLogMessage << "Start parallel inner product, _rd = " << rd << std::endl;
-  vector_type vv;
-  auto l_v=lhs.View();
-  auto r_v=rhs.View();
-  thread_for( r,rd,{
-
-    int so=r*grid->_ostride[orthogdim]; // base offset for start of plane 
-
-    for(int n=0;n<e1;n++){
-      for(int b=0;b<e2;b++){
-        int ss = so + n * stride + b;
-        vv = TensorRemove(innerProduct(l_v[ss], r_v[ss]));
-        lvSum[r] = lvSum[r] + vv;
-      }
-    }
-  });
-  // std::cout << GridLogMessage << "End parallel inner product" << std::endl;
-
-  // Sum across simd lanes in the plane, breaking out orthog dir.
-  Coordinate icoor(Nd);
-  for(int rt=0;rt<rd;rt++){
-
-    iScalar<vector_type> temp; 
-    temp._internal = lvSum[rt];
-    extract(temp,extracted);
-
-    for(int idx=0;idx<Nsimd;idx++){
-
-      grid->iCoorFromIindex(icoor,idx);
-
-      int ldx =rt+icoor[orthogdim]*rd;
-
-      lsSum[ldx]=lsSum[ldx]+extracted[idx]._internal;
-
-    }
-  }
-  // std::cout << GridLogMessage << "End sum over simd lanes" << std::endl;
-}
-template <class vobj>
-static void globalSliceInnerProductVector(std::vector<ComplexD> &result, const Lattice<vobj> &lhs, std::vector<typename vobj::scalar_type> &lsSum, int orthogdim)
-{
-  typedef typename vobj::scalar_type scalar_type;
-  GridBase *grid = lhs.Grid();
-  int fd = result.size();
-  int ld = lsSum.size();
-  // sum over nodes.
-  std::vector<scalar_type> gsum;
-  gsum.resize(fd, scalar_type(0.0));
-  // std::cout << GridLogMessage << "Start of gsum[t] creation:" << std::endl;
-  for(int t=0;t<fd;t++){
-    int pt = t/ld; // processor plane
-    int lt = t%ld;
-    if ( pt == grid->_processor_coor[orthogdim] ) {
-      gsum[t]=lsSum[lt];
-    }
-  }
-  // std::cout << GridLogMessage << "End of gsum[t] creation:" << std::endl;
-  // std::cout << GridLogMessage << "Start of GlobalSumVector:" << std::endl;
-  grid->GlobalSumVector(&gsum[0], fd);
-  // std::cout << GridLogMessage << "End of GlobalSumVector:" << std::endl;
-
-  result = gsum;
-}
 template<class vobj>
 static void sliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim) 
 {

From a55d0ba8fe91f2ba1ef6d06676790d3fa7ebe670 Mon Sep 17 00:00:00 2001
From: Fionn O hOgain <fionn.o.hogain@ed.ac.uk>
Date: Tue, 8 Oct 2019 15:52:01 +0100
Subject: [PATCH 8/8] Changing back interface for Gamma3pt

---
 Hadrons/Modules/MContraction/Gamma3pt.hpp | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/Hadrons/Modules/MContraction/Gamma3pt.hpp b/Hadrons/Modules/MContraction/Gamma3pt.hpp
index 006d7413..cd5c89e6 100644
--- a/Hadrons/Modules/MContraction/Gamma3pt.hpp
+++ b/Hadrons/Modules/MContraction/Gamma3pt.hpp
@@ -75,7 +75,7 @@ public:
                                     std::string,  q1,
                                     std::string,  q2,
                                     std::string,  q3,
-                                    std::string,  gammas,
+                                    std::string,  gamma,
                                     unsigned int, tSnk,
                                     std::string,  output);
 };
@@ -149,7 +149,7 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::parseGammaString(std::vector<Gamma::Alge
 {
     gammaList.clear();
     // Determine gamma matrices to insert at source/sink.
-    if (par().gammas.compare("all") == 0)
+    if (par().gamma.compare("all") == 0)
     {
         // Do all contractions.
         for (unsigned int i = 1; i < Gamma::nGamma; i += 2)
@@ -160,7 +160,7 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::parseGammaString(std::vector<Gamma::Alge
     else
     {
         // Parse individual contractions from input string.
-        gammaList = strToVec<Gamma::Algebra>(par().gammas);
+        gammaList = strToVec<Gamma::Algebra>(par().gamma);
     } 
 }
 // execution ///////////////////////////////////////////////////////////////////
@@ -169,7 +169,7 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
 {
     LOG(Message) << "Computing 3pt contractions '" << getName() << "' using"
                  << " quarks '" << par().q1 << "', '" << par().q2 << "' and '"
-                 << par().q3 << "', with " << par().gammas << " insertions." 
+                 << par().q3 << "', with " << par().gamma << " insertions." 
                  << std::endl;
 
     // Initialise variables. q2 and q3 are normal propagators, q1 may be