From 3e1d268fa33869a4b18ccfe52e99ab271fe342cc Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Wed, 17 Oct 2018 20:25:32 +0100
Subject: [PATCH 01/33] Hadrons: DiskVector optimisation

---
 Hadrons/DiskVector.hpp | 65 ++++++++++++++++++++++++++++--------------
 1 file changed, 44 insertions(+), 21 deletions(-)

diff --git a/Hadrons/DiskVector.hpp b/Hadrons/DiskVector.hpp
index 94c3e597..eb29755f 100644
--- a/Hadrons/DiskVector.hpp
+++ b/Hadrons/DiskVector.hpp
@@ -83,6 +83,7 @@ public:
 public:
     DiskVectorBase(const std::string dirname, const unsigned int size = 0,
                    const unsigned int cacheSize = 1, const bool clean = true);
+    DiskVectorBase(DiskVectorBase<T> &&v) = default;
     virtual ~DiskVectorBase(void);
     const T & operator[](const unsigned int i) const;
     RwAccessHelper operator[](const unsigned int i);
@@ -97,14 +98,16 @@ private:
     void cacheInsert(const unsigned int i, const T &obj) const;
     void clean(void);
 private:
-    std::string                                dirname_;
-    unsigned int                               size_, cacheSize_;
-    double                                     access_{0.}, hit_{0.};
-    bool                                       clean_;
+    std::string                               dirname_;
+    unsigned int                              size_, cacheSize_;
+    double                                    access_{0.}, hit_{0.};
+    bool                                      clean_;
     // using pointers to allow modifications when class is const
     // semantic: const means data unmodified, but cache modification allowed
-    std::unique_ptr<std::map<unsigned int, T>> cachePtr_;
-    std::unique_ptr<std::deque<unsigned int>>  loadsPtr_;                
+    std::unique_ptr<std::vector<T>>                       cachePtr_;
+    std::unique_ptr<std::map<unsigned int, unsigned int>> indexPtr_;
+    std::unique_ptr<std::stack<unsigned int>>             freePtr_;
+    std::unique_ptr<std::deque<unsigned int>>             loadsPtr_;                
 };
 
 /******************************************************************************
@@ -201,7 +204,9 @@ DiskVectorBase<T>::DiskVectorBase(const std::string dirname,
                                   const unsigned int cacheSize,
                                   const bool clean)
 : dirname_(dirname), size_(size), cacheSize_(cacheSize), clean_(clean)
-, cachePtr_(new std::map<unsigned int, T>())
+, cachePtr_(new std::vector<T>(size))
+, indexPtr_(new std::map<unsigned int, unsigned int>())
+, freePtr_(new std::stack<unsigned int>)
 , loadsPtr_(new std::deque<unsigned int>())
 {
     struct stat s;
@@ -211,6 +216,10 @@ DiskVectorBase<T>::DiskVectorBase(const std::string dirname,
         HADRONS_ERROR(Io, "directory '" + dirname + "' already exists")
     }
     mkdir(dirname);
+    for (unsigned int i = 0; i < cacheSize_; ++i)
+    {
+        freePtr_->push(i);
+    }
 }
 
 template <typename T>
@@ -225,8 +234,10 @@ DiskVectorBase<T>::~DiskVectorBase(void)
 template <typename T>
 const T & DiskVectorBase<T>::operator[](const unsigned int i) const
 {
-    auto &cache  = *cachePtr_;
-    auto &loads  = *loadsPtr_;
+    auto &cache   = *cachePtr_;
+    auto &index   = *indexPtr_;
+    auto &freeInd = *freePtr_;
+    auto &loads   = *loadsPtr_;
 
     DV_DEBUG_MSG(this, "accessing " << i << " (RO)");
 
@@ -235,7 +246,7 @@ const T & DiskVectorBase<T>::operator[](const unsigned int i) const
         HADRONS_ERROR(Size, "index out of range");
     }
     const_cast<double &>(access_)++;
-    if (cache.find(i) == cache.end())
+    if (index.find(i) == index.end())
     {
         // cache miss
         DV_DEBUG_MSG(this, "cache miss");
@@ -262,7 +273,7 @@ const T & DiskVectorBase<T>::operator[](const unsigned int i) const
     DV_DEBUG_MSG(this, "in cache: " << msg);
 #endif
 
-    return cache.at(i);
+    return cache[index.at(i)];
 }
 
 template <typename T>
@@ -300,13 +311,15 @@ std::string DiskVectorBase<T>::filename(const unsigned int i) const
 template <typename T>
 void DiskVectorBase<T>::evict(void) const
 {
-    auto &cache = *cachePtr_;
-    auto &loads = *loadsPtr_;
+    auto &index   = *indexPtr_;
+    auto &freeInd = *freePtr_;
+    auto &loads   = *loadsPtr_;
 
-    if (cache.size() >= cacheSize_)
+    if (index.size() >= cacheSize_)
     {
         DV_DEBUG_MSG(this, "evicting " << loads.front());
-        cache.erase(loads.front());
+        freeInd.push(index.at(loads.front()));
+        index.erase(loads.front());
         loads.pop_front();
     }
 }
@@ -314,29 +327,39 @@ void DiskVectorBase<T>::evict(void) const
 template <typename T>
 void DiskVectorBase<T>::fetch(const unsigned int i) const
 {
-    auto &cache = *cachePtr_;
-    auto &loads = *loadsPtr_;
+    auto &cache   = *cachePtr_;
+    auto &index   = *indexPtr_;
+    auto &freeInd = *freePtr_;
+    auto &loads   = *loadsPtr_;
+
     struct stat s;
 
     DV_DEBUG_MSG(this, "loading " << i << " from disk");
 
     evict();
+    
     if(stat(filename(i).c_str(), &s) != 0)
     {
         HADRONS_ERROR(Io, "disk vector element " + std::to_string(i) + " uninitialised");
     }
-    load(cache[i], filename(i));
+    index[i] = freeInd.top();
+    freeInd.pop();
+    load(cache[index.at(i)], filename(i));
     loads.push_back(i);
 }
 
 template <typename T>
 void DiskVectorBase<T>::cacheInsert(const unsigned int i, const T &obj) const
 {
-    auto &cache = *cachePtr_;
-    auto &loads = *loadsPtr_;
+    auto &cache   = *cachePtr_;
+    auto &index   = *indexPtr_;
+    auto &freeInd = *freePtr_;
+    auto &loads   = *loadsPtr_;
 
     evict();
-    cache[i] = obj;
+    index[i] = freeInd.top();
+    freeInd.pop();
+    cache[index.at(i)] = obj;
     loads.push_back(i);
 
 #ifdef DV_DEBUG

From f05b25dae4566fe7f5ed8c887e39bf759c753ade Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Wed, 17 Oct 2018 20:26:26 +0100
Subject: [PATCH 02/33] Hadrons: A2AMatrix load

---
 Hadrons/A2AMatrix.hpp                         | 122 +++++++++++++++++-
 .../Modules/MContraction/A2AAslashField.hpp   |   8 +-
 .../Modules/MContraction/A2AMesonField.hpp    |   8 +-
 3 files changed, 122 insertions(+), 16 deletions(-)

diff --git a/Hadrons/A2AMatrix.hpp b/Hadrons/A2AMatrix.hpp
index 95133f30..480d327f 100644
--- a/Hadrons/A2AMatrix.hpp
+++ b/Hadrons/A2AMatrix.hpp
@@ -37,6 +37,10 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #define HADRONS_A2AM_NAME "a2aMatrix"
 #endif
 
+#ifndef HADRONS_A2AM_IO_TYPE
+#define HADRONS_A2AM_IO_TYPE ComplexF
+#endif
+
 #define HADRONS_A2AM_PARALLEL_IO
 
 BEGIN_HADRONS_NAMESPACE
@@ -51,6 +55,9 @@ BEGIN_HADRONS_NAMESPACE
 template <typename T>
 using A2AMatrixSet = Eigen::TensorMap<Eigen::Tensor<T, 5, Eigen::RowMajor>>;
 
+template <typename T>
+using A2AMatrix = Eigen::Matrix<T, -1, -1, Eigen::RowMajor>;
+
 /******************************************************************************
  *                      Abstract class for A2A kernels                        *
  ******************************************************************************/
@@ -76,10 +83,15 @@ public:
     // constructors
     A2AMatrixIo(void) = default;
     A2AMatrixIo(std::string filename, std::string dataname, 
-                const unsigned int nt, const unsigned int ni,
-                const unsigned int nj);
+                const unsigned int nt, const unsigned int ni = 0,
+                const unsigned int nj = 0);
     // destructor
     ~A2AMatrixIo(void) = default;
+    // access
+    unsigned int getNi(void) const;
+    unsigned int getNj(void) const;
+    unsigned int getNt(void) const;
+    size_t       getSize(void) const;
     // file allocation
     template <typename MetadataType>
     void initFile(const MetadataType &d, const unsigned int chunkSize);
@@ -88,9 +100,11 @@ public:
                    const unsigned int blockSizei, const unsigned int blockSizej);
     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);
 private:
-    std::string  filename_, dataname_;
-    unsigned int nt_, ni_, nj_;
+    std::string  filename_{""}, dataname_{""};
+    unsigned int nt_{0}, ni_{0}, nj_{0};
 };
 
 /******************************************************************************
@@ -148,6 +162,31 @@ A2AMatrixIo<T>::A2AMatrixIo(std::string filename, std::string dataname,
 , nt_(nt), ni_(ni), nj_(nj)
 {}
 
+// access //////////////////////////////////////////////////////////////////////
+template <typename T>
+unsigned int A2AMatrixIo<T>::getNt(void) const
+{
+    return nt_;
+}
+
+template <typename T>
+unsigned int A2AMatrixIo<T>::getNi(void) const
+{
+    return ni_;
+}
+
+template <typename T>
+unsigned int A2AMatrixIo<T>::getNj(void) const
+{
+    return nj_;
+}
+
+template <typename T>
+size_t A2AMatrixIo<T>::getSize(void) const
+{
+    return nt_*ni_*nj_*sizeof(T);
+}
+
 // file allocation /////////////////////////////////////////////////////////////
 template <typename T>
 template <typename MetadataType>
@@ -226,6 +265,81 @@ void A2AMatrixIo<T>::saveBlock(const A2AMatrixSet<T> &m,
     saveBlock(m.data() + offset, i, j, blockSizei, blockSizej);
 }
 
+template <typename T>
+template <template <class> class Vec, typename VecT>
+void A2AMatrixIo<T>::load(Vec<VecT> &v, double *tRead)
+{
+#ifdef HAVE_HDF5
+    Hdf5Reader           reader(filename_);
+    std::vector<hsize_t> hdim;
+    H5NS::DataSet        dataset;
+    H5NS::DataSpace      dataspace;
+    H5NS::CompType       datatype;
+    H5NS::DSetCreatPropList plist;
+    
+    push(reader, dataname_);
+    auto &group = reader.getGroup();
+    dataset     = group.openDataSet(HADRONS_A2AM_NAME);
+    datatype    = dataset.getCompType();
+    dataspace   = dataset.getSpace();
+    plist       = dataset.getCreatePlist();
+    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 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_)
+        {
+            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];
+    }
+
+    A2AMatrix<T>         buf(ni_, nj_);
+    std::vector<hsize_t> count    = {1, static_cast<hsize_t>(ni_),
+                                     static_cast<hsize_t>(nj_)},
+                         stride   = {1, 1, 1},
+                         block    = {1, 1, 1},
+                         memCount = {static_cast<hsize_t>(ni_),
+                                     static_cast<hsize_t>(nj_)};
+    H5NS::DataSpace      memspace(memCount.size(), memCount.data());
+
+    std::cout << "Loading timeslice";
+    std::cout.flush();
+    *tRead = 0.;
+    for (unsigned int t = 0; t < nt_; ++t)
+    {
+        std::vector<hsize_t> offset = {static_cast<hsize_t>(t), 0, 0};
+        
+        if (t % 10 == 0)
+        {
+            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 (tRead) *tRead += usecond();
+        v[t] = buf.template cast<VecT>();
+    }
+    std::cout << std::endl;
+#else
+    HADRONS_ERROR(Implementation, "all-to-all matrix I/O needs HDF5 library");
+#endif
+}
+
 /******************************************************************************
  *               A2AMatrixBlockComputation template implementation            *
  ******************************************************************************/
diff --git a/Hadrons/Modules/MContraction/A2AAslashField.hpp b/Hadrons/Modules/MContraction/A2AAslashField.hpp
index 792ff6ee..a1d64d3d 100644
--- a/Hadrons/Modules/MContraction/A2AAslashField.hpp
+++ b/Hadrons/Modules/MContraction/A2AAslashField.hpp
@@ -6,10 +6,6 @@
 #include <Hadrons/ModuleFactory.hpp>
 #include <Hadrons/A2AMatrix.hpp>
 
-#ifndef ASF_IO_TYPE
-#define ASF_IO_TYPE ComplexF
-#endif
-
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
@@ -86,7 +82,7 @@ public:
     typedef A2AMatrixBlockComputation<Complex, 
                                       FermionField, 
                                       A2AAslashFieldMetadata, 
-                                      ASF_IO_TYPE> Computation;
+                                      HADRONS_A2AM_IO_TYPE> Computation;
     typedef AslashFieldKernel<Complex, FImpl> Kernel;
 public:
     // constructor
@@ -169,7 +165,7 @@ void TA2AAslashField<FImpl, PhotonImpl>::execute(void)
         LOG(Message) << "  " << name << std::endl;
     }
     LOG(Message) << "A-slash field size: " << nt << "*" << N_i << "*" << N_j 
-                 << " (filesize " << sizeString(nt*N_i*N_j*sizeof(ASF_IO_TYPE)) 
+                 << " (filesize " << sizeString(nt*N_i*N_j*sizeof(HADRONS_A2AM_IO_TYPE)) 
                  << "/EM field)" << std::endl;
     
     // preparing "B" complexified fields
diff --git a/Hadrons/Modules/MContraction/A2AMesonField.hpp b/Hadrons/Modules/MContraction/A2AMesonField.hpp
index 28759b26..1b3e7d8f 100644
--- a/Hadrons/Modules/MContraction/A2AMesonField.hpp
+++ b/Hadrons/Modules/MContraction/A2AMesonField.hpp
@@ -35,10 +35,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Hadrons/ModuleFactory.hpp>
 #include <Hadrons/A2AMatrix.hpp>
 
-#ifndef MF_IO_TYPE
-#define MF_IO_TYPE ComplexF
-#endif
-
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
@@ -118,7 +114,7 @@ public:
     typedef A2AMatrixBlockComputation<Complex, 
                                       FermionField, 
                                       A2AMesonFieldMetadata, 
-                                      MF_IO_TYPE> Computation;
+                                      HADRONS_A2AM_IO_TYPE> Computation;
     typedef MesonFieldKernel<Complex, FImpl> Kernel;
 public:
     // constructor
@@ -248,7 +244,7 @@ void TA2AMesonField<FImpl>::execute(void)
         LOG(Message) << "  " << g << std::endl;
     }
     LOG(Message) << "Meson field size: " << nt << "*" << N_i << "*" << N_j 
-                 << " (filesize " << sizeString(nt*N_i*N_j*sizeof(MF_IO_TYPE)) 
+                 << " (filesize " << sizeString(nt*N_i*N_j*sizeof(HADRONS_A2AM_IO_TYPE)) 
                  << "/momentum/bilinear)" << std::endl;
 
     auto &ph = envGet(std::vector<ComplexField>, momphName_);

From f709329d969cc832fce71e8861d3f345e3f2a2c3 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Wed, 17 Oct 2018 20:26:48 +0100
Subject: [PATCH 03/33] Hadrons: first version of a contractor utility

---
 Hadrons/Global.hpp               |   1 +
 Hadrons/Utilities/Contractor.cc  | 215 +++++++++++++++++++++++++++++++
 Hadrons/Utilities/Contractor.hpp |  12 ++
 Hadrons/Utilities/Makefile.am    |   5 +-
 4 files changed, 232 insertions(+), 1 deletion(-)
 create mode 100644 Hadrons/Utilities/Contractor.cc
 create mode 100644 Hadrons/Utilities/Contractor.hpp

diff --git a/Hadrons/Global.hpp b/Hadrons/Global.hpp
index e4003e73..d8047f09 100644
--- a/Hadrons/Global.hpp
+++ b/Hadrons/Global.hpp
@@ -32,6 +32,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 
 #include <set>
 #include <stack>
+#include <regex>
 #include <Grid/Grid.h>
 #include <cxxabi.h>
 
diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
new file mode 100644
index 00000000..ee7b6f2b
--- /dev/null
+++ b/Hadrons/Utilities/Contractor.cc
@@ -0,0 +1,215 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Utilities/Contractor.cc
+
+Copyright (C) 2015-2018
+
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#define DV_DEBUG
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/A2AMatrix.hpp>
+#include <Hadrons/DiskVector.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+namespace Contractor
+{
+    class GlobalPar: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar,
+                                        unsigned int, nt,
+                                        std::string, diskVectorDir,
+                                        std::string, output);
+    };
+
+    class A2AMatrixPar: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMatrixPar,
+                                        std::string, file,
+                                        std::string, dataset,
+                                        unsigned int, cacheSize,
+                                        std::string, name);
+    };
+
+    class ProductPar: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(ProductPar,
+                                        std::string, terms,
+                                        std::vector<std::string>, timeRange);
+    };
+}
+
+struct ContractorPar
+{
+    Contractor::GlobalPar                  global;
+    std::vector<Contractor::A2AMatrixPar>  a2aMatrix;
+    std::vector<Contractor::ProductPar>    product;
+};
+
+std::set<unsigned int> parseTimeRange(const std::string str)
+{
+    std::regex               rex("([0-9]+)|(([0-9]+)\\.\\.([0-9]+))");
+    std::smatch              sm;
+    std::vector<std::string> rstr = strToVec<std::string>(str);
+    std::set<unsigned int>   t;
+
+    for (auto &s: rstr)
+    {
+        std::regex_match(s, sm, rex);
+        if (sm[1].matched)
+        {
+            t.insert(std::stoi(sm[1].str()));
+        }
+        else if (sm[2].matched)
+        {
+            unsigned int ta, tb;
+
+            ta = std::stoi(sm[3].str());
+            tb = std::stoi(sm[4].str());
+            for (unsigned int ti = ta; ti <= tb; ++ti)
+            {
+                t.insert(ti);
+            }
+        }
+    }
+
+    return t;
+}
+
+int main(int argc, char* argv[])
+{
+    // parse command line
+    std::string   parFilename;
+
+    if (argc != 2)
+    {
+        std::cerr << "usage: " << argv[0] << " <parameter file>";
+        std::cerr << std::endl;
+        
+        return EXIT_FAILURE;
+    }
+    parFilename = argv[1];
+
+    // parse parameter file
+    ContractorPar par;
+    
+    unsigned int                        nMat, nCont;
+    {
+        XmlReader                           reader(parFilename);
+        read(reader, "global",    par.global);
+        read(reader, "a2aMatrix", par.a2aMatrix);
+        read(reader, "product",   par.product);
+    }
+    nMat  = par.a2aMatrix.size();
+    nCont = par.product.size();
+
+    // create diskvectors
+    std::map<std::string, EigenDiskVector<ComplexD>> a2aMat;
+    unsigned int                                     cacheSize;
+
+    for (auto &p: par.a2aMatrix)
+    {
+        std::string dirName = par.global.diskVectorDir + "/" + p.name;
+
+        a2aMat.emplace(p.name, EigenDiskVector<ComplexD>(dirName, par.global.nt, p.cacheSize));
+    }
+
+    // load data
+    for (unsigned int i = 0; i < a2aMat.size(); ++i)
+    {
+        auto   &p = par.a2aMatrix[i];
+        double t, size;
+
+        std::cout << "-- Loading '" << p.file << "'..." << std::endl;
+
+        A2AMatrixIo<HADRONS_A2AM_IO_TYPE> a2aIo(p.file, p.dataset, par.global.nt);
+
+        a2aIo.load(a2aMat.at(p.name), &t);
+        std::cout << "Read " << a2aIo.getSize() << " bytes in " << t << " usec, " << a2aIo.getSize()/t*1.0e6/1024/1024 << " MB/s" << std::endl;
+    }
+
+    // contract
+    EigenDiskVector<ComplexD>::Matrix buf;
+
+    for (auto &p: par.product)
+    {
+        std::vector<std::string>            term = strToVec<std::string>(p.terms);
+        std::vector<std::set<unsigned int>> times;
+
+        if (term.size() != p.timeRange.size())
+        {
+            HADRONS_ERROR(Size, "number of terms (" + std::to_string(term.size()) 
+                          + ") different from number of time ranges (" 
+                          + std::to_string(p.timeRange.size()) + ")");
+        }
+        for (auto &s: p.timeRange)
+        {
+            times.push_back(parseTimeRange(s));
+        }
+        
+        // test: just 2-pt function for now
+        if (term.size() != 2)
+        {
+            HADRONS_ERROR(Implementation, "only 2-pt function implemented");
+        }
+
+        std::vector<ComplexD> corr(par.global.nt);
+        double                tTrace = 0., flops, bytes;
+
+#define TIME_MOD(t) (((t) + par.global.nt) % par.global.nt)
+        for (unsigned int t = 0; t < par.global.nt; ++t)
+        {
+            corr[t] = ComplexD(0., 0.);
+        }
+        for (auto &trans: times[1])
+        for (auto &t0: times[0])
+        {
+            unsigned int tt0 = TIME_MOD(t0 + trans);
+            EigenDiskVector<ComplexD>::Matrix m0 = a2aMat.at(term[0])[tt0];
+            EigenDiskVector<ComplexD>::Matrix m1 = a2aMat.at(term[1])[trans];
+
+            std::cout << "tr(" << term[0] << "[" << tt0 << "]*"
+                      << term[1] << "[" << trans << "])";
+            tTrace  = -usecond();
+            corr[t0] += (m0*m1).trace();
+            tTrace += usecond();
+            flops = 6.*m0.rows()*m0.cols() + 2.*(m0.rows()*m0.cols() - 1.);
+            bytes = 2.*m0.rows()*m0.cols()*sizeof(ComplexD);
+            std::cout << " -- Perf " << flops/tTrace/1.0e3 << " GFlop/s " 
+                      << bytes/tTrace*1.0e6/1024/1024/1024 << " GB/s "<< std::endl;
+        }
+        for (unsigned int t = 0; t < par.global.nt; ++t)
+        {
+            std::cout << t << " " << corr[t] << std::endl;
+        }
+#undef TIME_MOD
+   }
+
+    
+    return EXIT_SUCCESS;
+}
diff --git a/Hadrons/Utilities/Contractor.hpp b/Hadrons/Utilities/Contractor.hpp
new file mode 100644
index 00000000..9640c7c8
--- /dev/null
+++ b/Hadrons/Utilities/Contractor.hpp
@@ -0,0 +1,12 @@
+#ifndef  Hadrons_Contractor_hpp_
+#define Hadrons_Contractor_hpp_
+
+#include <Hadrons/Global.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Contractor_hpp_
diff --git a/Hadrons/Utilities/Makefile.am b/Hadrons/Utilities/Makefile.am
index 529def24..1b257c32 100644
--- a/Hadrons/Utilities/Makefile.am
+++ b/Hadrons/Utilities/Makefile.am
@@ -1,4 +1,4 @@
-bin_PROGRAMS = HadronsXmlRun HadronsFermionEP64To32
+bin_PROGRAMS = HadronsXmlRun HadronsFermionEP64To32 HadronsContractor
 
 HadronsXmlRun_SOURCES = HadronsXmlRun.cc
 HadronsXmlRun_LDADD   = ../libHadrons.a ../../Grid/libGrid.a
@@ -6,3 +6,6 @@ HadronsXmlRun_LDADD   = ../libHadrons.a ../../Grid/libGrid.a
 HadronsFermionEP64To32_SOURCES  = EigenPackCast.cc
 HadronsFermionEP64To32_CXXFLAGS = $(AM_CXXFLAGS) -DFIN=WilsonImplD::FermionField -DFOUT=WilsonImplF::FermionField
 HadronsFermionEP64To32_LDADD    = ../libHadrons.a ../../Grid/libGrid.a
+
+HadronsContractor_SOURCES = Contractor.cc
+HadronsContractor_LDADD   = ../libHadrons.a ../../Grid/libGrid.a

From f333f3e5755c12a9613e4213816bd582db66ffde Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Thu, 18 Oct 2018 17:48:25 +0100
Subject: [PATCH 04/33] Hadrons: DiskVector save-on-eviction and faster CRC32
 for Eigen matrices

---
 Grid/util/Sha.h        |   3 +-
 Hadrons/DiskVector.hpp | 103 ++++++++++++++++++++++++-----------------
 2 files changed, 63 insertions(+), 43 deletions(-)

diff --git a/Grid/util/Sha.h b/Grid/util/Sha.h
index b0a8cc10..15cee338 100644
--- a/Grid/util/Sha.h
+++ b/Grid/util/Sha.h
@@ -34,10 +34,11 @@ extern "C" {
 class GridChecksum
 {
 public:
-  static inline uint32_t crc32(void *data,size_t bytes)
+  static inline uint32_t crc32(const void *data,size_t bytes)
   {
     return ::crc32(0L,(unsigned char *)data,bytes);
   }
+
   template <typename T>
   static inline std::string sha256_string(const std::vector<T> &hash)
   {
diff --git a/Hadrons/DiskVector.hpp b/Hadrons/DiskVector.hpp
index eb29755f..f3023a6e 100644
--- a/Hadrons/DiskVector.hpp
+++ b/Hadrons/DiskVector.hpp
@@ -59,14 +59,18 @@ public:
         : master_(master), cmaster_(master), i_(i) {}
 
         // operator=: somebody is trying to store a vector element
-        // write to disk and cache
+        // write to cache and tag as modified
         T &operator=(const T &obj) const
         {
+            auto &cache    = *master_.cachePtr_;
+            auto &modified = *master_.modifiedPtr_;
+            auto &index    = *master_.indexPtr_;
+
             DV_DEBUG_MSG(&master_, "writing to " << i_);
             master_.cacheInsert(i_, obj);
-            master_.save(master_.filename(i_), obj);
+            modified[index.at(i_)] = true;
             
-            return master_.cachePtr_->at(i_);
+            return cache[index.at(i_)];
         }
 
         // implicit cast to const object reference and redirection
@@ -98,13 +102,14 @@ private:
     void cacheInsert(const unsigned int i, const T &obj) const;
     void clean(void);
 private:
-    std::string                               dirname_;
-    unsigned int                              size_, cacheSize_;
-    double                                    access_{0.}, hit_{0.};
-    bool                                      clean_;
+    std::string                                           dirname_;
+    unsigned int                                          size_, cacheSize_;
+    double                                                access_{0.}, hit_{0.};
+    bool                                                  clean_;
     // 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_;
+    std::unique_ptr<std::vector<bool>>                    modifiedPtr_;
     std::unique_ptr<std::map<unsigned int, unsigned int>> indexPtr_;
     std::unique_ptr<std::stack<unsigned int>>             freePtr_;
     std::unique_ptr<std::deque<unsigned int>>             loadsPtr_;                
@@ -149,24 +154,24 @@ public:
 private:
     virtual void load(EigenDiskVectorMat<T> &obj, const std::string filename) const
     {
-        std::ifstream              f(filename, std::ios::binary);
-        std::vector<unsigned char> hash(SHA256_DIGEST_LENGTH);
-        Eigen::Index               nRow, nCol;
-        size_t                     matSize;
-        double                     t;
+        std::ifstream f(filename, std::ios::binary);
+        uint32_t      crc, check;
+        Eigen::Index  nRow, nCol;
+        size_t        matSize;
+        double        t;
 
-        f.read(reinterpret_cast<char *>(hash.data()), hash.size()*sizeof(unsigned char));
-        f.read(reinterpret_cast<char *>(&nRow), sizeof(Eigen::Index));
-        f.read(reinterpret_cast<char *>(&nCol), sizeof(Eigen::Index));
+        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);
         t  = -usecond();
         f.read(reinterpret_cast<char *>(obj.data()), matSize);
         t += usecond();
         DV_DEBUG_MSG(this, "Eigen read " << matSize/t*1.0e6/1024/1024 << " MB/s");
-        auto check = GridChecksum::sha256(obj.data(), matSize);
-        DV_DEBUG_MSG(this, "Eigen sha256 " << GridChecksum::sha256_string(check));
-        if (hash != check)
+        check = GridChecksum::crc32(obj.data(), matSize);
+        DV_DEBUG_MSG(this, "Eigen crc32 " << std::hex << check << std::dec);
+        if (crc != check)
         {
             HADRONS_ERROR(Io, "checksum failed")
         }
@@ -174,24 +179,24 @@ private:
 
     virtual void save(const std::string filename, const EigenDiskVectorMat<T> &obj) const
     {
-        std::ofstream              f(filename, std::ios::binary);
-        std::vector<unsigned char> hash(SHA256_DIGEST_LENGTH);
-        Eigen::Index               nRow, nCol;
-        size_t                     matSize;
-        double                     t;
+        std::ofstream f(filename, std::ios::binary);
+        uint32_t      crc;
+        Eigen::Index  nRow, nCol;
+        size_t        matSize;
+        double        t;
         
         nRow    = obj.rows();
         nCol    = obj.cols();
         matSize = nRow*nCol*sizeof(T);
-        hash    = GridChecksum::sha256(obj.data(), matSize);
-        DV_DEBUG_MSG(this, "Eigen sha256 " << GridChecksum::sha256_string(hash));
-        f.write(reinterpret_cast<char *>(hash.data()), hash.size()*sizeof(unsigned char));
-        f.write(reinterpret_cast<char *>(&nRow), sizeof(Eigen::Index));
-        f.write(reinterpret_cast<char *>(&nCol), sizeof(Eigen::Index));
+        crc     = GridChecksum::crc32(obj.data(), matSize);
+        f.write(reinterpret_cast<char *>(&crc), sizeof(crc));
+        f.write(reinterpret_cast<char *>(&nRow), sizeof(nRow));
+        f.write(reinterpret_cast<char *>(&nCol), sizeof(nCol));
         t  = -usecond();
         f.write(reinterpret_cast<const char *>(obj.data()), matSize);
         t += usecond();
         DV_DEBUG_MSG(this, "Eigen write " << matSize/t*1.0e6/1024/1024 << " MB/s");
+        DV_DEBUG_MSG(this, "Eigen crc32 " << std::hex << crc << std::dec);
     }
 };
 
@@ -205,6 +210,7 @@ DiskVectorBase<T>::DiskVectorBase(const std::string dirname,
                                   const bool clean)
 : dirname_(dirname), size_(size), cacheSize_(cacheSize), clean_(clean)
 , cachePtr_(new std::vector<T>(size))
+, modifiedPtr_(new std::vector<bool>(size, false))
 , indexPtr_(new std::map<unsigned int, unsigned int>())
 , freePtr_(new std::stack<unsigned int>)
 , loadsPtr_(new std::deque<unsigned int>())
@@ -311,15 +317,24 @@ std::string DiskVectorBase<T>::filename(const unsigned int i) const
 template <typename T>
 void DiskVectorBase<T>::evict(void) const
 {
-    auto &index   = *indexPtr_;
-    auto &freeInd = *freePtr_;
-    auto &loads   = *loadsPtr_;
+    auto &cache    = *cachePtr_;
+    auto &modified = *modifiedPtr_;
+    auto &index    = *indexPtr_;
+    auto &freeInd  = *freePtr_;
+    auto &loads    = *loadsPtr_;
 
     if (index.size() >= cacheSize_)
     {
-        DV_DEBUG_MSG(this, "evicting " << loads.front());
-        freeInd.push(index.at(loads.front()));
-        index.erase(loads.front());
+        unsigned int i = loads.front();
+        
+        DV_DEBUG_MSG(this, "evicting " << i);
+        if (modified[index.at(i)])
+        {
+            DV_DEBUG_MSG(this, "element " << i << " modified, saving to disk");
+            save(filename(i), cache[index.at(i)]);
+        }
+        freeInd.push(index.at(i));
+        index.erase(i);
         loads.pop_front();
     }
 }
@@ -327,10 +342,11 @@ void DiskVectorBase<T>::evict(void) const
 template <typename T>
 void DiskVectorBase<T>::fetch(const unsigned int i) const
 {
-    auto &cache   = *cachePtr_;
-    auto &index   = *indexPtr_;
-    auto &freeInd = *freePtr_;
-    auto &loads   = *loadsPtr_;
+    auto &cache    = *cachePtr_;
+    auto &modified = *modifiedPtr_;
+    auto &index    = *indexPtr_;
+    auto &freeInd  = *freePtr_;
+    auto &loads    = *loadsPtr_;
 
     struct stat s;
 
@@ -346,21 +362,24 @@ void DiskVectorBase<T>::fetch(const unsigned int i) const
     freeInd.pop();
     load(cache[index.at(i)], filename(i));
     loads.push_back(i);
+    modified[index.at(i)] = false;
 }
 
 template <typename T>
 void DiskVectorBase<T>::cacheInsert(const unsigned int i, const T &obj) const
 {
-    auto &cache   = *cachePtr_;
-    auto &index   = *indexPtr_;
-    auto &freeInd = *freePtr_;
-    auto &loads   = *loadsPtr_;
+    auto &cache    = *cachePtr_;
+    auto &modified = *modifiedPtr_;
+    auto &index    = *indexPtr_;
+    auto &freeInd  = *freePtr_;
+    auto &loads    = *loadsPtr_;
 
     evict();
     index[i] = freeInd.top();
     freeInd.pop();
     cache[index.at(i)] = obj;
     loads.push_back(i);
+    modified[index.at(i)] = false;
 
 #ifdef DV_DEBUG
     std::string msg;

From a7cfa26901b83b1825a7eee91f3ce1372efcfd78 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Thu, 18 Oct 2018 17:50:16 +0100
Subject: [PATCH 05/33] Hadrons: reverse A2A matrix load for better DiskVector
 cache reuse

---
 Hadrons/A2AMatrix.hpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/Hadrons/A2AMatrix.hpp b/Hadrons/A2AMatrix.hpp
index 480d327f..9842e15b 100644
--- a/Hadrons/A2AMatrix.hpp
+++ b/Hadrons/A2AMatrix.hpp
@@ -318,8 +318,9 @@ void A2AMatrixIo<T>::load(Vec<VecT> &v, double *tRead)
     std::cout << "Loading timeslice";
     std::cout.flush();
     *tRead = 0.;
-    for (unsigned int t = 0; t < nt_; ++t)
+    for (unsigned int tp1 = nt_; tp1 > 0; --tp1)
     {
+        unsigned int         t      = tp1 - 1;
         std::vector<hsize_t> offset = {static_cast<hsize_t>(t), 0, 0};
         
         if (t % 10 == 0)

From f31d6bfec2894804ef667c6722429a6d4016db84 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Thu, 18 Oct 2018 17:50:35 +0100
Subject: [PATCH 06/33] Hadrons: contractor cleaning and better error check

---
 Hadrons/Utilities/Contractor.cc | 40 ++++++++++++++++++++-------------
 1 file changed, 24 insertions(+), 16 deletions(-)

diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index ee7b6f2b..a66ad373 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -24,8 +24,6 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-#define DV_DEBUG
-
 #include <Hadrons/Global.hpp>
 #include <Hadrons/A2AMatrix.hpp>
 #include <Hadrons/DiskVector.hpp>
@@ -71,19 +69,27 @@ struct ContractorPar
     std::vector<Contractor::ProductPar>    product;
 };
 
-std::set<unsigned int> parseTimeRange(const std::string str)
+std::set<unsigned int> parseTimeRange(const std::string str, const unsigned int nt)
 {
     std::regex               rex("([0-9]+)|(([0-9]+)\\.\\.([0-9]+))");
     std::smatch              sm;
     std::vector<std::string> rstr = strToVec<std::string>(str);
-    std::set<unsigned int>   t;
+    std::set<unsigned int>   tSet;
+    
 
     for (auto &s: rstr)
     {
         std::regex_match(s, sm, rex);
         if (sm[1].matched)
         {
-            t.insert(std::stoi(sm[1].str()));
+            unsigned int t;
+            
+            t = std::stoi(sm[1].str());
+            if (t >= nt)
+            {
+                HADRONS_ERROR(Range, "time out of range (from expression '" + str + "')");
+            }
+            tSet.insert(t);
         }
         else if (sm[2].matched)
         {
@@ -91,14 +97,18 @@ std::set<unsigned int> parseTimeRange(const std::string str)
 
             ta = std::stoi(sm[3].str());
             tb = std::stoi(sm[4].str());
+            if ((ta >= nt) or (tb >= nt))
+            {
+                HADRONS_ERROR(Range, "time out of range (from expression '" + str + "')");
+            }
             for (unsigned int ti = ta; ti <= tb; ++ti)
             {
-                t.insert(ti);
+                tSet.insert(ti);
             }
         }
     }
 
-    return t;
+    return tSet;
 }
 
 int main(int argc, char* argv[])
@@ -117,14 +127,12 @@ int main(int argc, char* argv[])
 
     // parse parameter file
     ContractorPar par;
-    
-    unsigned int                        nMat, nCont;
-    {
-        XmlReader                           reader(parFilename);
-        read(reader, "global",    par.global);
-        read(reader, "a2aMatrix", par.a2aMatrix);
-        read(reader, "product",   par.product);
-    }
+    unsigned int  nMat, nCont;
+    XmlReader     reader(parFilename);
+
+    read(reader, "global",    par.global);
+    read(reader, "a2aMatrix", par.a2aMatrix);
+    read(reader, "product",   par.product);
     nMat  = par.a2aMatrix.size();
     nCont = par.product.size();
 
@@ -169,7 +177,7 @@ int main(int argc, char* argv[])
         }
         for (auto &s: p.timeRange)
         {
-            times.push_back(parseTimeRange(s));
+            times.push_back(parseTimeRange(s, par.global.nt));
         }
         
         // test: just 2-pt function for now

From 0ffcfea7240a2645ad58c77499055b956c5c8279 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Tue, 23 Oct 2018 17:08:16 +0100
Subject: [PATCH 07/33] Hadrons: contractor benchmark

---
 Hadrons/Utilities/ContractorBenchmark.cc | 124 +++++++++++++++++++++++
 Hadrons/Utilities/Makefile.am            |   5 +-
 2 files changed, 128 insertions(+), 1 deletion(-)
 create mode 100644 Hadrons/Utilities/ContractorBenchmark.cc

diff --git a/Hadrons/Utilities/ContractorBenchmark.cc b/Hadrons/Utilities/ContractorBenchmark.cc
new file mode 100644
index 00000000..9a8c0f4e
--- /dev/null
+++ b/Hadrons/Utilities/ContractorBenchmark.cc
@@ -0,0 +1,124 @@
+#include <Hadrons/Global.hpp>
+#include <Hadrons/DiskVector.hpp>
+
+using namespace Grid;
+
+#ifndef EIGEN_ORDER 
+#define EIGEN_ORDER RowMajor
+#endif
+
+typedef Eigen::Matrix<ComplexD, -1, -1, Eigen::EIGEN_ORDER> CMat;
+typedef std::vector<std::vector<CMat>>                      CMatSet;
+
+template <typename TwoMatFn>
+inline void trBenchmark(const std::string name, const CMatSet &mat, TwoMatFn fn)
+{
+    double       t, flops, bytes, n = mat[0][0].rows()*mat[0][0].cols();
+    unsigned int nMat = mat[0].size();
+    ComplexD     buf;
+
+    t = 0.;
+    for (unsigned int i = 0; i < nMat; ++i)
+    {
+        t -= usecond();
+        fn(buf, mat[0][i], mat[1][i]);
+        t += usecond();
+    }
+    flops = nMat*(6.*n + 2.*(n - 1.));
+    bytes = nMat*(2.*n*sizeof(ComplexD));
+
+    std::cout << std::setw(30) << name << ": result= "
+              << std::setw(10) << buf
+              << std::setw(10) << flops/t/1.0e3 << " GFlop/s " 
+              << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s "
+              << std::endl;
+    ::sleep(1);
+}
+
+int main(int argc, char *argv[])
+{
+    // parse command line
+    Eigen::Index ni, nj, nMat;
+
+    if (argc != 4)
+    {
+        std::cerr << "usage: " << argv[0] << " <Ni> <Nj> <#matrices>";
+        std::cerr << std::endl;
+        
+        return EXIT_FAILURE;
+    }
+    ni   = std::stoi(argv[1]);
+    nj   = std::stoi(argv[2]);
+    nMat = std::stoi(argv[3]);
+
+    std::cout << "==== generating random matrices" << std::endl;    
+    CMatSet mat(2);
+    CMat    buf;
+
+    mat[0].resize(nMat, Eigen::MatrixXcd::Random(ni, nj));
+    mat[1].resize(nMat, Eigen::MatrixXcd::Random(nj, ni));
+
+    std::cout << "==== tr(A*B) benchmark" << std::endl;
+    trBenchmark("Naive loop rows first", mat,
+    [](ComplexD &res, const CMat &a, const CMat &b)
+    { 
+        res = 0.;
+        #pragma omp parallel for schedule(static) reduction(+:res)
+        for (unsigned int i = 0; i < a.rows(); ++i)
+        for (unsigned int j = 0; j < a.cols(); ++j)
+        {
+            res += a(i, j)*b(j, i);
+        }
+    });
+    trBenchmark("Naive loop cols first", mat,
+    [](ComplexD &res, const CMat &a, const CMat &b)
+    { 
+        res = 0.;
+        for (unsigned int j = 0; j < a.cols(); ++j)
+        for (unsigned int i = 0; i < a.rows(); ++i)
+        {
+            res += a(i, j)*b(j, i);
+        }
+    });
+    trBenchmark("Eigen tr(A*B)", mat,
+    [](ComplexD &res, const CMat &a, const CMat &b)
+    { 
+        res = (a*b).trace();
+    });
+    trBenchmark("Eigen global dot", mat,
+    [&buf](ComplexD &res, const CMat &a, const CMat &b)
+    { 
+        buf = b.transpose();
+        Eigen::Map<const Eigen::VectorXcd> av(a.data(), a.rows()*a.cols());
+        Eigen::Map<const Eigen::VectorXcd> bv(buf.data(), b.rows()*b.cols());
+
+        res = av.conjugate().dot(bv);
+    });
+    trBenchmark("Eigen row-wise dot", mat,
+    [](ComplexD &res, const CMat &a, const CMat &b)
+    { 
+        res = 0.;
+        #pragma omp parallel for schedule(static) reduction(+:res)
+        for (unsigned int r = 0; r < a.rows(); ++r)
+        {
+            res += a.row(r).conjugate().dot(b.col(r));
+        }
+    });
+    trBenchmark("Eigen col-wise dot", mat,
+    [](ComplexD &res, const CMat &a, const CMat &b)
+    { 
+        res = 0.;
+        #pragma omp parallel for schedule(static) reduction(+:res)
+        for (unsigned int r = 0; r < a.cols(); ++r)
+        {
+            res += a.col(r).conjugate().dot(b.row(r));
+        }
+    });
+    trBenchmark("Eigen Hadamard", mat,
+    [](ComplexD &res, const CMat &a, const CMat &b)
+    { 
+        res = a.cwiseProduct(b.transpose()).sum();
+    });
+
+    return EXIT_SUCCESS;
+}
diff --git a/Hadrons/Utilities/Makefile.am b/Hadrons/Utilities/Makefile.am
index 1b257c32..4f324d6d 100644
--- a/Hadrons/Utilities/Makefile.am
+++ b/Hadrons/Utilities/Makefile.am
@@ -1,4 +1,4 @@
-bin_PROGRAMS = HadronsXmlRun HadronsFermionEP64To32 HadronsContractor
+bin_PROGRAMS = HadronsXmlRun HadronsFermionEP64To32 HadronsContractor HadronsContractorBenchmark
 
 HadronsXmlRun_SOURCES = HadronsXmlRun.cc
 HadronsXmlRun_LDADD   = ../libHadrons.a ../../Grid/libGrid.a
@@ -9,3 +9,6 @@ HadronsFermionEP64To32_LDADD    = ../libHadrons.a ../../Grid/libGrid.a
 
 HadronsContractor_SOURCES = Contractor.cc
 HadronsContractor_LDADD   = ../libHadrons.a ../../Grid/libGrid.a
+
+HadronsContractorBenchmark_SOURCES = ContractorBenchmark.cc
+HadronsContractorBenchmark_LDADD   = ../libHadrons.a ../../Grid/libGrid.a

From 4a7415e83c8827270390f044734063e698261b52 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Tue, 23 Oct 2018 21:00:54 +0100
Subject: [PATCH 08/33] Hadrons: contractor benchmark update

---
 Hadrons/Utilities/ContractorBenchmark.cc | 187 ++++++++++++++++++++---
 1 file changed, 165 insertions(+), 22 deletions(-)

diff --git a/Hadrons/Utilities/ContractorBenchmark.cc b/Hadrons/Utilities/ContractorBenchmark.cc
index 9a8c0f4e..d2ff7c9e 100644
--- a/Hadrons/Utilities/ContractorBenchmark.cc
+++ b/Hadrons/Utilities/ContractorBenchmark.cc
@@ -1,17 +1,30 @@
 #include <Hadrons/Global.hpp>
 #include <Hadrons/DiskVector.hpp>
+#ifdef USE_MKL
+#include "mkl.h"
+#include "mkl_cblas.h"
+#endif
 
 using namespace Grid;
 
+#define EIGEN_ROW_MAJOR 1
+#define EIGEN_COL_MAJOR 2
 #ifndef EIGEN_ORDER 
-#define EIGEN_ORDER RowMajor
+#define EIGEN_ORDER EIGEN_ROW_MAJOR
+#endif
+#if (EIGEN_ORDER == EIGEN_ROW_MAJOR)
+typedef Eigen::Matrix<ComplexD, -1, -1, Eigen::RowMajor> CMat;
+#elif (EIGEN_ORDER == EIGEN_COL_MAJOR)
+typedef Eigen::Matrix<ComplexD, -1, -1, Eigen::ColMajor> CMat;
 #endif
 
-typedef Eigen::Matrix<ComplexD, -1, -1, Eigen::EIGEN_ORDER> CMat;
 typedef std::vector<std::vector<CMat>>                      CMatSet;
 
-template <typename TwoMatFn>
-inline void trBenchmark(const std::string name, const CMatSet &mat, TwoMatFn fn)
+#pragma omp declare reduction(ComplexPlus: ComplexD: omp_out += omp_in) 
+
+template <typename Function>
+inline void trBenchmark(const std::string name, const CMatSet &mat, 
+                        const ComplexD ref, Function fn)
 {
     double       t, flops, bytes, n = mat[0][0].rows()*mat[0][0].cols();
     unsigned int nMat = mat[0].size();
@@ -27,8 +40,37 @@ inline void trBenchmark(const std::string name, const CMatSet &mat, TwoMatFn fn)
     flops = nMat*(6.*n + 2.*(n - 1.));
     bytes = nMat*(2.*n*sizeof(ComplexD));
 
-    std::cout << std::setw(30) << name << ": result= "
-              << std::setw(10) << buf
+    std::cout << std::setw(30) << name << ": diff= "
+              << std::setw(12) << std::norm(buf-ref)
+              << std::setw(10) << t/1.0e6 << " sec "
+              << std::setw(10) << flops/t/1.0e3 << " GFlop/s " 
+              << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s "
+              << std::endl;
+    ::sleep(1);
+}
+
+template <typename Function>
+inline void mulBenchmark(const std::string name, const CMatSet &mat, 
+                         const CMat ref, Function fn)
+{
+    double       t, flops, bytes;
+    double       nr = mat[0][0].rows(), nc = mat[0][0].cols(), n = nr*nc;
+    unsigned int nMat = mat[0].size();
+    CMat         buf(mat[0][0].rows(), mat[0][0].rows());
+
+    t = 0.;
+    for (unsigned int i = 0; i < nMat; ++i)
+    {
+        t -= usecond();
+        fn(buf, mat[0][i], mat[1][i]);
+        t += usecond();
+    }
+    flops = nMat*(n*(6.*nc + 2.*(nc - 1.)));
+    bytes = nMat*((2.*n+nr*nr)*sizeof(ComplexD));
+
+    std::cout << std::setw(30) << name << ": diff= "
+              << std::setw(12) << (buf-ref).squaredNorm()
+              << std::setw(10) << t/1.0e6 << " sec "
               << std::setw(10) << flops/t/1.0e3 << " GFlop/s " 
               << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s "
               << std::endl;
@@ -51,41 +93,60 @@ int main(int argc, char *argv[])
     nj   = std::stoi(argv[2]);
     nMat = std::stoi(argv[3]);
 
-    std::cout << "==== generating random matrices" << std::endl;    
-    CMatSet mat(2);
-    CMat    buf;
+    std::cout << "\n*** ALL-TO-ALL MATRIX CONTRACTION BENCHMARK ***\n" << std::endl;
+    std::cout << nMat << " " << ni << "x" << nj << " matrices\n" << std::endl;
+#ifdef EIGEN_USE_MKL_ALL
+    std::cout << "Eigen uses the MKL" << std::endl;
+#endif
+    std::cout << "Eigen uses ";
+#if (EIGEN_ORDER == EIGEN_ROW_MAJOR)
+    std::cout << "row-major ordering" << std::endl;
+#elif (EIGEN_ORDER == EIGEN_COL_MAJOR)
+    std::cout << "column-major ordering" << std::endl;
+#endif
+    std::cout << "Eigen uses " << Eigen::nbThreads() << " thread(s)" << std::endl;
+#ifdef USE_MKL
+    std::cout << "MKL   uses " << mkl_get_max_threads() << " thread(s)" << std::endl;
+#endif
+    std::cout << std::endl;
+
+    CMatSet  mat(2);
+    CMat     buf;
+    ComplexD ref;
 
     mat[0].resize(nMat, Eigen::MatrixXcd::Random(ni, nj));
     mat[1].resize(nMat, Eigen::MatrixXcd::Random(nj, ni));
 
-    std::cout << "==== tr(A*B) benchmark" << std::endl;
-    trBenchmark("Naive loop rows first", mat,
+    std::cout << "==== tr(A*B) benchmarks" << std::endl;
+    ref = (mat[0].back()*mat[1].back()).trace();
+    trBenchmark("Naive loop rows first", mat, ref,
     [](ComplexD &res, const CMat &a, const CMat &b)
     { 
         res = 0.;
-        #pragma omp parallel for schedule(static) reduction(+:res)
+        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
         for (unsigned int i = 0; i < a.rows(); ++i)
         for (unsigned int j = 0; j < a.cols(); ++j)
         {
             res += a(i, j)*b(j, i);
         }
     });
-    trBenchmark("Naive loop cols first", mat,
+    trBenchmark("Naive loop cols first", mat, ref,
     [](ComplexD &res, const CMat &a, const CMat &b)
     { 
         res = 0.;
+        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
         for (unsigned int j = 0; j < a.cols(); ++j)
         for (unsigned int i = 0; i < a.rows(); ++i)
         {
             res += a(i, j)*b(j, i);
         }
     });
-    trBenchmark("Eigen tr(A*B)", mat,
+    trBenchmark("Eigen tr(A*B)", mat, ref,
     [](ComplexD &res, const CMat &a, const CMat &b)
     { 
         res = (a*b).trace();
     });
-    trBenchmark("Eigen global dot", mat,
+    trBenchmark("Eigen global dot", mat, ref,
     [&buf](ComplexD &res, const CMat &a, const CMat &b)
     { 
         buf = b.transpose();
@@ -94,31 +155,113 @@ int main(int argc, char *argv[])
 
         res = av.conjugate().dot(bv);
     });
-    trBenchmark("Eigen row-wise dot", mat,
+    trBenchmark("Eigen row-wise dot", mat, ref,
     [](ComplexD &res, const CMat &a, const CMat &b)
     { 
         res = 0.;
-        #pragma omp parallel for schedule(static) reduction(+:res)
+        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
         for (unsigned int r = 0; r < a.rows(); ++r)
         {
             res += a.row(r).conjugate().dot(b.col(r));
         }
     });
-    trBenchmark("Eigen col-wise dot", mat,
+    trBenchmark("Eigen col-wise dot", mat, ref,
     [](ComplexD &res, const CMat &a, const CMat &b)
     { 
         res = 0.;
-        #pragma omp parallel for schedule(static) reduction(+:res)
-        for (unsigned int r = 0; r < a.cols(); ++r)
+        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
+        for (unsigned int c = 0; c < a.cols(); ++c)
         {
-            res += a.col(r).conjugate().dot(b.row(r));
+            res += a.col(c).conjugate().dot(b.row(c));
         }
     });
-    trBenchmark("Eigen Hadamard", mat,
+    trBenchmark("Eigen Hadamard", mat, ref,
     [](ComplexD &res, const CMat &a, const CMat &b)
     { 
         res = a.cwiseProduct(b.transpose()).sum();
     });
+#ifdef USE_MKL
+    trBenchmark("MKL row-wise zdotu", mat, ref,
+    [](ComplexD &res, const CMat &a, const CMat &b)
+    { 
+        res = 0.;
+        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
+        for (unsigned int r = 0; r < a.rows(); ++r)
+        {
+            ComplexD tmp;
+
+#if (EIGEN_ORDER == EIGEN_ROW_MAJOR)
+            cblas_zdotu_sub(a.cols(), a.data() + r*a.cols(), 1, b.data() + r, b.cols(), &tmp);
+#elif (EIGEN_ORDER == EIGEN_COL_MAJOR)
+            cblas_zdotu_sub(a.cols(), a.data() + r, a.rows(), b.data() + r*b.rows(), 1, &tmp);
+#endif
+            res += tmp;
+        }
+    });
+    trBenchmark("MKL col-wise zdotu", mat, ref,
+    [](ComplexD &res, const CMat &a, const CMat &b)
+    { 
+        res = 0.;
+        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
+        for (unsigned int c = 0; c < a.cols(); ++c)
+        {
+            ComplexD tmp;
+
+#if (EIGEN_ORDER == EIGEN_ROW_MAJOR)
+            cblas_zdotu_sub(a.rows(), a.data() + c, a.cols(), b.data() + c*b.cols(), 1, &tmp);
+#elif (EIGEN_ORDER == EIGEN_COL_MAJOR)
+            cblas_zdotu_sub(a.rows(), a.data() + c*a.rows(), 1, b.data() + c, b.rows(), &tmp);
+#endif
+            res += tmp;
+        }
+    });
+#endif
+
+    std::cout << std::endl;
+    std::cout << "==== A*B benchmarks" << std::endl;
+    buf = mat[0].back()*mat[1].back();
+    mulBenchmark("Naive", mat, buf,
+    [](CMat &res, const CMat &a, const CMat &b)
+    { 
+        unsigned int ni = a.rows(), nj = a.cols();
+
+        #pragma omp parallel for collapse(2)
+        for (unsigned int i = 0; i < ni; ++i)
+        for (unsigned int k = 0; k < ni; ++k)
+        {
+            res(i, k) = a(i, 0)*b(0, k);
+        }
+        #pragma omp parallel for collapse(2)
+        for (unsigned int i = 0; i < ni; ++i)
+        for (unsigned int k = 0; k < ni; ++k)
+        for (unsigned int j = 1; j < nj; ++j)
+        {
+            res(i, k) += a(i, j)*b(j, k);
+        }
+    });
+    mulBenchmark("Eigen A*B", mat, buf,
+    [](CMat &res, const CMat &a, const CMat &b)
+    { 
+        res = a*b;
+    });
+#ifdef USE_MKL
+    mulBenchmark("MKL A*B", mat, buf,
+    [](CMat &res, const CMat &a, const CMat &b)
+    {
+        const ComplexD one(1., 0.), zero(0., 0.);
+#if (EIGEN_ORDER == EIGEN_ROW_MAJOR)
+        cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),
+                    a.cols(), &one, a.data(), a.cols(), b.data(), b.cols(), &zero,
+                    res.data(), res.cols());
+#elif (EIGEN_ORDER == EIGEN_COL_MAJOR)
+        cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),
+                    a.cols(), &one, a.data(), a.rows(), b.data(), b.rows(), &zero,
+                    res.data(), res.rows());
+#endif
+    });
+#endif
+
+    std::cout << std::endl;
 
     return EXIT_SUCCESS;
 }

From 1651111d18a53afd8b848a7d61e7404f24fc12ed Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Mon, 5 Nov 2018 21:29:13 +0000
Subject: [PATCH 09/33] Hadrons: final, portable form of the contractor
 benchmark

---
 Grid/threads/Threads.h                   |  13 +
 Hadrons/Utilities/ContractorBenchmark.cc | 552 ++++++++++++++---------
 2 files changed, 352 insertions(+), 213 deletions(-)

diff --git a/Grid/threads/Threads.h b/Grid/threads/Threads.h
index 5e448c44..6403bff9 100644
--- a/Grid/threads/Threads.h
+++ b/Grid/threads/Threads.h
@@ -38,8 +38,19 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #ifdef GRID_OMP
 #include <omp.h>
 
+// complex reductions
+#pragma omp declare reduction(ComplexPlus: Grid::Complex: omp_out += omp_in)
+#pragma omp declare reduction(GridVComplexPlus: Grid::vComplex: omp_out += omp_in)
+#pragma omp declare reduction(ComplexDPlus: Grid::ComplexD: omp_out += omp_in)
+#pragma omp declare reduction(GridVComplexDPlus: Grid::vComplexD: omp_out += omp_in)
+#pragma omp declare reduction(ComplexFPlus: Grid::ComplexF: omp_out += omp_in)
+#pragma omp declare reduction(GridVComplexFPlus: Grid::vComplexF: omp_out += omp_in)
+
 #define PARALLEL_FOR_LOOP        _Pragma("omp parallel for schedule(static)")
 #define PARALLEL_FOR_LOOP_INTERN _Pragma("omp for schedule(static)")
+#define PARALLEL_FOR_REDUCE_HELPER0(x) #x
+#define PARALLEL_FOR_REDUCE_HELPER1(op, var) PARALLEL_FOR_REDUCE_HELPER0(omp parallel for schedule(static) reduction(op:var))
+#define PARALLEL_FOR_LOOP_REDUCE(op, var) _Pragma(PARALLEL_FOR_REDUCE_HELPER1(op, var))
 #define PARALLEL_NESTED_LOOP2 _Pragma("omp parallel for collapse(2)")
 #define PARALLEL_NESTED_LOOP5 _Pragma("omp parallel for collapse(5)")
 #define PARALLEL_REGION       _Pragma("omp parallel")
@@ -47,6 +58,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #else
 #define PARALLEL_FOR_LOOP
 #define PARALLEL_FOR_LOOP_INTERN
+#define PARALLEL_FOR_LOOP_REDUCE(op, var)
 #define PARALLEL_NESTED_LOOP2
 #define PARALLEL_NESTED_LOOP5
 #define PARALLEL_REGION
@@ -56,6 +68,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #define parallel_region    PARALLEL_REGION
 #define parallel_for       PARALLEL_FOR_LOOP for
 #define parallel_for_internal PARALLEL_FOR_LOOP_INTERN for
+#define parallel_for_reduce(op, var) PARALLEL_FOR_LOOP_REDUCE(op, var) for
 #define parallel_for_nest2 PARALLEL_NESTED_LOOP2 for
 #define parallel_for_nest5 PARALLEL_NESTED_LOOP5 for
 
diff --git a/Hadrons/Utilities/ContractorBenchmark.cc b/Hadrons/Utilities/ContractorBenchmark.cc
index d2ff7c9e..2a8b033b 100644
--- a/Hadrons/Utilities/ContractorBenchmark.cc
+++ b/Hadrons/Utilities/ContractorBenchmark.cc
@@ -1,86 +1,346 @@
 #include <Hadrons/Global.hpp>
-#include <Hadrons/DiskVector.hpp>
 #ifdef USE_MKL
 #include "mkl.h"
 #include "mkl_cblas.h"
 #endif
+#ifdef AVX2
+#include <immintrin.h>
+#endif
 
 using namespace Grid;
 
-#define EIGEN_ROW_MAJOR 1
-#define EIGEN_COL_MAJOR 2
-#ifndef EIGEN_ORDER 
-#define EIGEN_ORDER EIGEN_ROW_MAJOR
-#endif
-#if (EIGEN_ORDER == EIGEN_ROW_MAJOR)
-typedef Eigen::Matrix<ComplexD, -1, -1, Eigen::RowMajor> CMat;
-#elif (EIGEN_ORDER == EIGEN_COL_MAJOR)
-typedef Eigen::Matrix<ComplexD, -1, -1, Eigen::ColMajor> CMat;
+typedef Eigen::Matrix<ComplexD, -1, -1, Eigen::RowMajor> RMCMat;
+typedef Eigen::Matrix<ComplexD, -1, -1, Eigen::ColMajor> CMCMat;
+
+#ifdef GRID_COMMS_MPI3
+#define GET_RANK(rank, nMpi) \
+MPI_Comm_size(MPI_COMM_WORLD, &(nMpi));\
+MPI_Comm_rank(MPI_COMM_WORLD, &(rank))
+#define BARRIER() MPI_Barrier(MPI_COMM_WORLD)
+#define INIT() MPI_Init(NULL, NULL)
+#define FINALIZE() MPI_Finalize()
+#else
+#define GET_RANK(rank, nMpi) (nMpi) = 1; (rank) = 0
+#define BARRIER()
+#define INIT()
+#define FINALIZE()
 #endif
 
-typedef std::vector<std::vector<CMat>>                      CMatSet;
-
-#pragma omp declare reduction(ComplexPlus: ComplexD: omp_out += omp_in) 
-
-template <typename Function>
-inline void trBenchmark(const std::string name, const CMatSet &mat, 
-                        const ComplexD ref, Function fn)
+template <typename Function, typename MatLeft, typename MatRight>
+inline void trBenchmark(const std::string name, const MatLeft &left,
+                        const MatRight &right, const ComplexD ref, Function fn)
 {
-    double       t, flops, bytes, n = mat[0][0].rows()*mat[0][0].cols();
-    unsigned int nMat = mat[0].size();
+    double       t, flops, bytes, n = left[0].rows()*left[0].cols();
+    unsigned int nMat = left.size();
+    int          nMpi, rank;
     ComplexD     buf;
 
     t = 0.;
-    for (unsigned int i = 0; i < nMat; ++i)
+    GET_RANK(rank, nMpi);
+    t = -usecond();
+    BARRIER();
+    for (unsigned int i = rank*nMat/nMpi; i < (rank+1)*nMat/nMpi; ++i)
     {
-        t -= usecond();
-        fn(buf, mat[0][i], mat[1][i]);
-        t += usecond();
+        fn(buf, left[i], right[i]);      
     }
+    BARRIER();
+    t += usecond();
     flops = nMat*(6.*n + 2.*(n - 1.));
     bytes = nMat*(2.*n*sizeof(ComplexD));
 
-    std::cout << std::setw(30) << name << ": diff= "
-              << std::setw(12) << std::norm(buf-ref)
-              << std::setw(10) << t/1.0e6 << " sec "
-              << std::setw(10) << flops/t/1.0e3 << " GFlop/s " 
-              << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s "
-              << std::endl;
+    if (rank == 0)
+    {
+        std::cout << std::setw(30) << name << ": diff= "
+                  << std::setw(12) << std::norm(buf-ref)
+                  << std::setw(10) << t/1.0e6 << " sec "
+                  << std::setw(10) << flops/t/1.0e3 << " GFlop/s " 
+                  << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s "
+                  << std::endl;
+    }
     ::sleep(1);
 }
 
-template <typename Function>
-inline void mulBenchmark(const std::string name, const CMatSet &mat, 
-                         const CMat ref, Function fn)
+template <typename Function, typename MatV, typename Mat>
+inline void mulBenchmark(const std::string name, const MatV &left,
+                         const MatV &right, const Mat &ref, Function fn)
 {
     double       t, flops, bytes;
-    double       nr = mat[0][0].rows(), nc = mat[0][0].cols(), n = nr*nc;
-    unsigned int nMat = mat[0].size();
-    CMat         buf(mat[0][0].rows(), mat[0][0].rows());
+    double       nr = left[0].rows(), nc = left[0].cols(), n = nr*nc;
+    unsigned int nMat = left.size();
+    int          nMpi, rank;
+    Mat          buf(left[0].rows(), left[0].rows());
 
     t = 0.;
-    for (unsigned int i = 0; i < nMat; ++i)
+    GET_RANK(rank, nMpi);
+    t = -usecond();
+    BARRIER();
+    for (unsigned int i = rank*nMat/nMpi; i < (rank+1)*nMat/nMpi; ++i)
     {
-        t -= usecond();
-        fn(buf, mat[0][i], mat[1][i]);
-        t += usecond();
+        fn(buf, left[i], right[i]);
     }
-    flops = nMat*(n*(6.*nc + 2.*(nc - 1.)));
-    bytes = nMat*((2.*n+nr*nr)*sizeof(ComplexD));
+    BARRIER();
+    t += usecond();
+    flops = nMat*(nr*nr*(6.*nc + 2.*(nc - 1.)));
+    bytes = nMat*(2*nc*nr*sizeof(ComplexD));
 
-    std::cout << std::setw(30) << name << ": diff= "
-              << std::setw(12) << (buf-ref).squaredNorm()
-              << std::setw(10) << t/1.0e6 << " sec "
-              << std::setw(10) << flops/t/1.0e3 << " GFlop/s " 
-              << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s "
-              << std::endl;
+    if (rank == 0)
+    {
+        std::cout << std::setw(30) << name << ": diff= "
+                << std::setw(12) << (buf-ref).squaredNorm()
+                << std::setw(10) << t/1.0e6 << " sec "
+                << std::setw(10) << flops/t/1.0e3 << " GFlop/s " 
+                << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s "
+                << std::endl;
+    }
     ::sleep(1);
 }
 
+template <typename MatLeft, typename MatRight>
+static inline void zdotuRow(ComplexD &res, const unsigned int aRow,
+                            const MatLeft &a, const MatRight &b)
+{
+    const ComplexD *aPt, *bPt;
+    unsigned int   aInc, bInc;
+
+    if (MatLeft::Options == Eigen::RowMajor)
+    {
+        aPt  = a.data() + aRow*a.cols();
+        aInc = 1;
+    }
+    else if (MatLeft::Options == Eigen::ColMajor)
+    {
+        aPt  = a.data() + aRow;
+        aInc = a.rows();
+    }
+    if (MatRight::Options == Eigen::RowMajor)
+    {
+        bPt  = b.data() + aRow;
+        bInc = b.cols();
+    }
+    else if (MatRight::Options == Eigen::ColMajor)
+    {
+        bPt  = b.data() + aRow*b.rows();
+        bInc = 1;
+    }
+    cblas_zdotu_sub(a.cols(), aPt, aInc, bPt, bInc, &res);
+}
+
+template <typename MatLeft, typename MatRight>
+static inline void zdotuCol(ComplexD &res, const unsigned int aCol,
+                            const MatLeft &a, const MatRight &b)
+{
+    const ComplexD *aPt, *bPt;
+    unsigned int   aInc, bInc;
+
+    if (MatLeft::Options == Eigen::RowMajor)
+    {
+        aPt  = a.data() + aCol;
+        aInc = a.cols();
+    }
+    else if (MatLeft::Options == Eigen::ColMajor)
+    {
+        aPt  = a.data() + aCol*a.rows();
+        aInc = 1;
+    }
+    if (MatRight::Options == Eigen::RowMajor)
+    {
+        bPt  = b.data() + aCol*b.cols();
+        bInc = 1;
+    }
+    else if (MatRight::Options == Eigen::ColMajor)
+    {
+        bPt  = b.data() + aCol;
+        bInc = b.rows();
+    }
+    cblas_zdotu_sub(a.rows(), aPt, aInc, bPt, bInc, &res);
+}
+
+template <typename MatLeft, typename MatRight>
+void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigned int nMat)
+{
+    std::vector<MatLeft>  left;
+    std::vector<MatRight> right;
+    MatRight              buf;
+    ComplexD              ref;
+    int                   rank, nMpi;
+
+    left.resize(nMat, MatLeft::Random(ni, nj));
+    right.resize(nMat, MatRight::Random(nj, ni));
+    GET_RANK(rank, nMpi);
+    if (rank == 0)
+    {
+        std::cout << "==== tr(A*B) benchmarks" << std::endl;
+        std::cout << "A matrices use ";
+        if (MatLeft::Options == Eigen::RowMajor)
+        {
+            std::cout << "row-major ordering" << std::endl;
+        }
+        else if (MatLeft::Options == Eigen::ColMajor)
+        {
+            std::cout << "col-major ordering" << std::endl;
+        }
+        std::cout << "B matrices use ";
+        if (MatRight::Options == Eigen::RowMajor)
+        {
+            std::cout << "row-major ordering" << std::endl;
+        }
+        else if (MatRight::Options == Eigen::ColMajor)
+        {
+            std::cout << "col-major ordering" << std::endl;
+        }
+        std::cout << std::endl;
+    }
+    BARRIER();
+    ref = (left.back()*right.back()).trace();
+    trBenchmark("Naive loop rows first", left, right, ref,
+    [](ComplexD &res, const MatLeft &a, const MatRight &b)
+    { 
+        res = 0.;
+        auto nr = a.rows(), nc = a.cols();
+        parallel_for_reduce(ComplexDPlus, res) (unsigned int i = 0; i < nr; ++i)
+        for (unsigned int j = 0; j < nc; ++j)
+        {
+            res += a(i, j)*b(j, i);
+        }
+    });
+    trBenchmark("Naive loop cols first", left, right, ref,
+    [](ComplexD &res, const MatLeft &a, const MatRight &b)
+    { 
+        res = 0.;
+        auto nr = a.rows(), nc = a.cols();
+        parallel_for_reduce(ComplexDPlus, res) (unsigned int j = 0; j < nc; ++j)
+        for (unsigned int i = 0; i < nr; ++i)
+        {
+            res += a(i, j)*b(j, i);
+        }
+    });
+    trBenchmark("Eigen tr(A*B)", left, right, ref,
+    [](ComplexD &res, const MatLeft &a, const MatRight &b)
+    { 
+        res = (a*b).trace();
+    });
+    trBenchmark("Eigen row-wise dot", left, right, ref,
+    [](ComplexD &res, const MatLeft &a, const MatRight &b)
+    { 
+        res = 0.;
+
+        parallel_for_reduce(ComplexDPlus, res) (unsigned int r = 0; r < a.rows(); ++r)
+        {
+            res += a.row(r).conjugate().dot(b.col(r));
+        }
+    });
+    trBenchmark("Eigen col-wise dot", left, right, ref,
+    [](ComplexD &res, const MatLeft &a, const MatRight &b)
+    { 
+        res = 0.;
+
+        parallel_for_reduce(ComplexDPlus, res) (unsigned int c = 0; c < a.cols(); ++c)
+        {
+            res += a.col(c).conjugate().dot(b.row(c));
+        }
+    });
+    trBenchmark("Eigen Hadamard", left, right, ref,
+    [](ComplexD &res, const MatLeft &a, const MatRight &b)
+    { 
+        res = a.cwiseProduct(b.transpose()).sum();
+    });
+#ifdef USE_MKL
+    trBenchmark("MKL row-wise zdotu", left, right, ref,
+    [](ComplexD &res, const MatLeft &a, const MatRight &b)
+    { 
+        res = 0.;
+
+        parallel_for_reduce(ComplexDPlus, res) (unsigned int r = 0; r < a.rows(); ++r)
+        {
+            ComplexD tmp;
+
+            zdotuRow(tmp, r, a, b);
+            res += tmp;
+        }
+    });
+    trBenchmark("MKL col-wise zdotu", left, right, ref,
+    [](ComplexD &res, const MatLeft &a, const MatRight &b)
+    { 
+        res = 0.;
+
+        parallel_for_reduce(ComplexDPlus, res) (unsigned int c = 0; c < a.cols(); ++c)
+        {
+            ComplexD tmp;
+
+            zdotuCol(tmp, c, a, b);
+            res += tmp;
+        }
+    });
+#endif
+    BARRIER();
+    if (rank == 0)
+    {
+        std::cout << std::endl;
+    }
+}
+
+template <typename Mat>
+void fullMulBenchmark(const unsigned int ni, const unsigned int nj, const unsigned int nMat)
+{
+    std::vector<Mat> left, right;
+    Mat              ref;
+    int              rank, nMpi;
+
+    left.resize(nMat, Mat::Random(ni, nj));
+    right.resize(nMat, Mat::Random(nj, ni));
+    GET_RANK(rank, nMpi);
+    if (rank == 0)
+    {
+        std::cout << "==== A*B benchmarks" << std::endl;
+        std::cout << "all matrices use ";
+        if (Mat::Options == Eigen::RowMajor)
+        {
+            std::cout << "row-major ordering" << std::endl;
+        }
+        else if (Mat::Options == Eigen::ColMajor)
+        {
+            std::cout << "col-major ordering" << std::endl;
+        }
+        std::cout << std::endl;
+    }
+    BARRIER();
+    ref = left.back()*right.back();
+    mulBenchmark("Eigen A*B", left, right, ref,
+    [](Mat &res, const Mat &a, const Mat &b)
+    { 
+        res = a*b;
+    });
+#ifdef USE_MKL
+    mulBenchmark("MKL A*B", left, right, ref,
+    [](Mat &res, const Mat &a, const Mat &b)
+    {
+        const ComplexD one(1., 0.), zero(0., 0.);
+        if (Mat::Options == Eigen::RowMajor)
+        {
+            cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),
+                        a.cols(), &one, a.data(), a.cols(), b.data(), b.cols(), &zero,
+                        res.data(), res.cols());
+        }
+        else if (Mat::Options == Eigen::ColMajor)
+        {
+            cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),
+                        a.cols(), &one, a.data(), a.rows(), b.data(), b.rows(), &zero,
+                        res.data(), res.rows());
+        }
+    });
+#endif
+    BARRIER();
+    if (rank == 0)
+    {
+        std::cout << std::endl;
+    }
+}
+
 int main(int argc, char *argv[])
 {
     // parse command line
     Eigen::Index ni, nj, nMat;
+    int          nMpi, rank;
 
     if (argc != 4)
     {
@@ -93,175 +353,41 @@ int main(int argc, char *argv[])
     nj   = std::stoi(argv[2]);
     nMat = std::stoi(argv[3]);
 
-    std::cout << "\n*** ALL-TO-ALL MATRIX CONTRACTION BENCHMARK ***\n" << std::endl;
-    std::cout << nMat << " " << ni << "x" << nj << " matrices\n" << std::endl;
-#ifdef EIGEN_USE_MKL_ALL
-    std::cout << "Eigen uses the MKL" << std::endl;
-#endif
-    std::cout << "Eigen uses ";
-#if (EIGEN_ORDER == EIGEN_ROW_MAJOR)
-    std::cout << "row-major ordering" << std::endl;
-#elif (EIGEN_ORDER == EIGEN_COL_MAJOR)
-    std::cout << "column-major ordering" << std::endl;
-#endif
-    std::cout << "Eigen uses " << Eigen::nbThreads() << " thread(s)" << std::endl;
-#ifdef USE_MKL
-    std::cout << "MKL   uses " << mkl_get_max_threads() << " thread(s)" << std::endl;
-#endif
-    std::cout << std::endl;
-
-    CMatSet  mat(2);
-    CMat     buf;
-    ComplexD ref;
-
-    mat[0].resize(nMat, Eigen::MatrixXcd::Random(ni, nj));
-    mat[1].resize(nMat, Eigen::MatrixXcd::Random(nj, ni));
-
-    std::cout << "==== tr(A*B) benchmarks" << std::endl;
-    ref = (mat[0].back()*mat[1].back()).trace();
-    trBenchmark("Naive loop rows first", mat, ref,
-    [](ComplexD &res, const CMat &a, const CMat &b)
-    { 
-        res = 0.;
-        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
-        for (unsigned int i = 0; i < a.rows(); ++i)
-        for (unsigned int j = 0; j < a.cols(); ++j)
-        {
-            res += a(i, j)*b(j, i);
-        }
-    });
-    trBenchmark("Naive loop cols first", mat, ref,
-    [](ComplexD &res, const CMat &a, const CMat &b)
-    { 
-        res = 0.;
-        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
-        for (unsigned int j = 0; j < a.cols(); ++j)
-        for (unsigned int i = 0; i < a.rows(); ++i)
-        {
-            res += a(i, j)*b(j, i);
-        }
-    });
-    trBenchmark("Eigen tr(A*B)", mat, ref,
-    [](ComplexD &res, const CMat &a, const CMat &b)
-    { 
-        res = (a*b).trace();
-    });
-    trBenchmark("Eigen global dot", mat, ref,
-    [&buf](ComplexD &res, const CMat &a, const CMat &b)
-    { 
-        buf = b.transpose();
-        Eigen::Map<const Eigen::VectorXcd> av(a.data(), a.rows()*a.cols());
-        Eigen::Map<const Eigen::VectorXcd> bv(buf.data(), b.rows()*b.cols());
-
-        res = av.conjugate().dot(bv);
-    });
-    trBenchmark("Eigen row-wise dot", mat, ref,
-    [](ComplexD &res, const CMat &a, const CMat &b)
-    { 
-        res = 0.;
-        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
-        for (unsigned int r = 0; r < a.rows(); ++r)
-        {
-            res += a.row(r).conjugate().dot(b.col(r));
-        }
-    });
-    trBenchmark("Eigen col-wise dot", mat, ref,
-    [](ComplexD &res, const CMat &a, const CMat &b)
-    { 
-        res = 0.;
-        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
-        for (unsigned int c = 0; c < a.cols(); ++c)
-        {
-            res += a.col(c).conjugate().dot(b.row(c));
-        }
-    });
-    trBenchmark("Eigen Hadamard", mat, ref,
-    [](ComplexD &res, const CMat &a, const CMat &b)
-    { 
-        res = a.cwiseProduct(b.transpose()).sum();
-    });
-#ifdef USE_MKL
-    trBenchmark("MKL row-wise zdotu", mat, ref,
-    [](ComplexD &res, const CMat &a, const CMat &b)
-    { 
-        res = 0.;
-        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
-        for (unsigned int r = 0; r < a.rows(); ++r)
-        {
-            ComplexD tmp;
-
-#if (EIGEN_ORDER == EIGEN_ROW_MAJOR)
-            cblas_zdotu_sub(a.cols(), a.data() + r*a.cols(), 1, b.data() + r, b.cols(), &tmp);
-#elif (EIGEN_ORDER == EIGEN_COL_MAJOR)
-            cblas_zdotu_sub(a.cols(), a.data() + r, a.rows(), b.data() + r*b.rows(), 1, &tmp);
-#endif
-            res += tmp;
-        }
-    });
-    trBenchmark("MKL col-wise zdotu", mat, ref,
-    [](ComplexD &res, const CMat &a, const CMat &b)
-    { 
-        res = 0.;
-        #pragma omp parallel for schedule(static) reduction(ComplexPlus:res)
-        for (unsigned int c = 0; c < a.cols(); ++c)
-        {
-            ComplexD tmp;
-
-#if (EIGEN_ORDER == EIGEN_ROW_MAJOR)
-            cblas_zdotu_sub(a.rows(), a.data() + c, a.cols(), b.data() + c*b.cols(), 1, &tmp);
-#elif (EIGEN_ORDER == EIGEN_COL_MAJOR)
-            cblas_zdotu_sub(a.rows(), a.data() + c*a.rows(), 1, b.data() + c, b.rows(), &tmp);
-#endif
-            res += tmp;
-        }
-    });
-#endif
-
-    std::cout << std::endl;
-    std::cout << "==== A*B benchmarks" << std::endl;
-    buf = mat[0].back()*mat[1].back();
-    mulBenchmark("Naive", mat, buf,
-    [](CMat &res, const CMat &a, const CMat &b)
-    { 
-        unsigned int ni = a.rows(), nj = a.cols();
-
-        #pragma omp parallel for collapse(2)
-        for (unsigned int i = 0; i < ni; ++i)
-        for (unsigned int k = 0; k < ni; ++k)
-        {
-            res(i, k) = a(i, 0)*b(0, k);
-        }
-        #pragma omp parallel for collapse(2)
-        for (unsigned int i = 0; i < ni; ++i)
-        for (unsigned int k = 0; k < ni; ++k)
-        for (unsigned int j = 1; j < nj; ++j)
-        {
-            res(i, k) += a(i, j)*b(j, k);
-        }
-    });
-    mulBenchmark("Eigen A*B", mat, buf,
-    [](CMat &res, const CMat &a, const CMat &b)
-    { 
-        res = a*b;
-    });
-#ifdef USE_MKL
-    mulBenchmark("MKL A*B", mat, buf,
-    [](CMat &res, const CMat &a, const CMat &b)
+    INIT();
+    GET_RANK(rank, nMpi);
+    if (rank == 0)
     {
-        const ComplexD one(1., 0.), zero(0., 0.);
-#if (EIGEN_ORDER == EIGEN_ROW_MAJOR)
-        cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),
-                    a.cols(), &one, a.data(), a.cols(), b.data(), b.cols(), &zero,
-                    res.data(), res.cols());
-#elif (EIGEN_ORDER == EIGEN_COL_MAJOR)
-        cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),
-                    a.cols(), &one, a.data(), a.rows(), b.data(), b.rows(), &zero,
-                    res.data(), res.rows());
-#endif
-    });
+        std::cout << "\n*** ALL-TO-ALL MATRIX CONTRACTION BENCHMARK ***\n" << std::endl;
+        std::cout << nMat << " couples of " << ni << "x" << nj << " matrices\n" << std::endl;
+
+        std::cout << nMpi << " MPI processes" << std::endl;
+#ifdef GRID_OMP
+        #pragma omp parallel
+        {
+            #pragma omp single
+            std::cout << omp_get_num_threads() << " threads\n" << std::endl; 
+        }
+#else
+        std::cout << "Single-threaded\n" << std::endl; 
 #endif
 
-    std::cout << std::endl;
+#ifdef EIGEN_USE_MKL_ALL
+        std::cout << "Eigen uses the MKL" << std::endl;
+#endif
+        std::cout << "Eigen uses " << Eigen::nbThreads() << " threads" << std::endl;
+#ifdef USE_MKL
+        std::cout << "MKL   uses " << mkl_get_max_threads() << " threads" << std::endl;
+#endif
+        std::cout << std::endl;
+    }
+
+    fullTrBenchmark<RMCMat, RMCMat>(ni, nj, nMat);
+    fullTrBenchmark<RMCMat, CMCMat>(ni, nj, nMat);
+    fullTrBenchmark<CMCMat, RMCMat>(ni, nj, nMat);
+    fullTrBenchmark<CMCMat, CMCMat>(ni, nj, nMat);
+    fullMulBenchmark<RMCMat>(ni, nj, nMat);
+    fullMulBenchmark<CMCMat>(ni, nj, nMat);
+    FINALIZE();
 
     return EXIT_SUCCESS;
 }

From 9734e3ee58bbdb228a9c671a534518bb994bb076 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Tue, 6 Nov 2018 19:47:41 +0000
Subject: [PATCH 10/33] Hadrons: (somewhat) faster build

---
 Hadrons/Makefile.am | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/Hadrons/Makefile.am b/Hadrons/Makefile.am
index 206011f4..b748085a 100644
--- a/Hadrons/Makefile.am
+++ b/Hadrons/Makefile.am
@@ -5,17 +5,17 @@ lib_LIBRARIES = libHadrons.a
 include modules.inc
 
 libHadrons_a_SOURCES = \
-    $(modules_cc)      \
     Application.cc     \
     Environment.cc     \
 	Exceptions.cc      \
     Global.cc          \
     Module.cc		   \
 	TimerArray.cc      \
-	VirtualMachine.cc
+	VirtualMachine.cc  \
+	$(modules_cc)
+	
 libHadrons_adir = $(includedir)/Hadrons
 nobase_libHadrons_a_HEADERS = \
-	$(modules_hpp)            \
 	A2AVectors.hpp            \
 	A2AMatrix.hpp             \
 	Application.hpp           \
@@ -33,4 +33,5 @@ nobase_libHadrons_a_HEADERS = \
 	ModuleFactory.hpp         \
 	Solver.hpp                \
 	TimerArray.hpp            \
-	VirtualMachine.hpp
+	VirtualMachine.hpp        \
+	$(modules_hpp)

From 88d9922e4fa9bcaec11819c08b866712ee0cf60c Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Tue, 6 Nov 2018 19:49:09 +0000
Subject: [PATCH 11/33] Hadrons: fast A2A matrix contraction kernels

---
 Grid/threads/Threads.h                   |   7 +-
 Hadrons/A2AMatrix.hpp                    | 199 +++++++++++++++++++++++
 Hadrons/Utilities/ContractorBenchmark.cc |  59 ++++---
 3 files changed, 234 insertions(+), 31 deletions(-)

diff --git a/Grid/threads/Threads.h b/Grid/threads/Threads.h
index 6403bff9..580278d6 100644
--- a/Grid/threads/Threads.h
+++ b/Grid/threads/Threads.h
@@ -39,12 +39,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <omp.h>
 
 // complex reductions
-#pragma omp declare reduction(ComplexPlus: Grid::Complex: omp_out += omp_in)
-#pragma omp declare reduction(GridVComplexPlus: Grid::vComplex: omp_out += omp_in)
-#pragma omp declare reduction(ComplexDPlus: Grid::ComplexD: omp_out += omp_in)
-#pragma omp declare reduction(GridVComplexDPlus: Grid::vComplexD: omp_out += omp_in)
-#pragma omp declare reduction(ComplexFPlus: Grid::ComplexF: omp_out += omp_in)
-#pragma omp declare reduction(GridVComplexFPlus: Grid::vComplexF: omp_out += omp_in)
+#pragma omp declare reduction(ComplexPlus:Grid::ComplexD, Grid::vComplexD, Grid::ComplexF, Grid::vComplexF: omp_out += omp_in)
 
 #define PARALLEL_FOR_LOOP        _Pragma("omp parallel for schedule(static)")
 #define PARALLEL_FOR_LOOP_INTERN _Pragma("omp for schedule(static)")
diff --git a/Hadrons/A2AMatrix.hpp b/Hadrons/A2AMatrix.hpp
index 9842e15b..5266d050 100644
--- a/Hadrons/A2AMatrix.hpp
+++ b/Hadrons/A2AMatrix.hpp
@@ -32,6 +32,10 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Hadrons/Global.hpp>
 #include <Hadrons/TimerArray.hpp>
 #include <Grid/Eigen/unsupported/CXX11/Tensor>
+#ifdef USE_MKL
+#include "mkl.h"
+#include "mkl_cblas.h"
+#endif
 
 #ifndef HADRONS_A2AM_NAME 
 #define HADRONS_A2AM_NAME "a2aMatrix"
@@ -58,6 +62,9 @@ using A2AMatrixSet = Eigen::TensorMap<Eigen::Tensor<T, 5, Eigen::RowMajor>>;
 template <typename T>
 using A2AMatrix = Eigen::Matrix<T, -1, -1, Eigen::RowMajor>;
 
+template <typename T>
+using A2AMatrixTr = Eigen::Matrix<T, -1, -1, Eigen::ColMajor>;
+
 /******************************************************************************
  *                      Abstract class for A2A kernels                        *
  ******************************************************************************/
@@ -150,6 +157,198 @@ private:
     std::vector<IoHelper> nodeIo_;
 };
 
+/******************************************************************************
+ *                       A2A matrix contraction kernels                       *
+ ******************************************************************************/
+class A2AContraction
+{
+public:
+    // accTrMul(acc, a, b): acc += tr(a*b)
+    template <typename C, typename MatLeft, typename MatRight>
+    static inline void accTrMul(C &acc, const MatLeft &a, const MatRight &b)
+    {
+        if ((MatLeft::Options == Eigen::RowMajor) and
+            (MatRight::Options == Eigen::ColMajor))
+        {
+            parallel_for_reduce(ComplexPlus, acc) (unsigned int r = 0; r < a.rows(); ++r)
+            {
+#ifdef USE_MKL
+                ComplexD tmp;
+                dotuRow(tmp, r, a, b);
+                acc += tmp;
+#else
+                acc += a.row(r).conjugate().dot(b.col(r));
+#endif
+            }
+        }
+        else
+        {
+            parallel_for_reduce(ComplexPlus, acc) (unsigned int c = 0; c < a.cols(); ++c)
+            {
+#ifdef USE_MKL
+                ComplexD tmp;
+                dotuCol(tmp, c, a, b);
+                acc += tmp;
+#else
+                acc += a.col(c).conjugate().dot(b.row(c));
+#endif
+            }
+        }
+    }
+
+    // mul(res, a, b): res = a*b
+#ifdef USE_MKL
+    template <template <class, int...> class Mat, int... Opts>
+    static inline void mul(Mat<ComplexD, Opts...> &res, 
+                           const Mat<ComplexD, Opts...> &a, 
+                           const Mat<ComplexD, Opts...> &b)
+    {
+        static const ComplexD one(1., 0.), zero(0., 0.);
+
+        if (Mat<ComplexD, Opts...>::Options == Eigen::RowMajor)
+        {
+            cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),
+                        a.cols(), &one, a.data(), a.cols(), b.data(), b.cols(), &zero,
+                        res.data(), res.cols());
+        }
+        else if (Mat<ComplexD, Opts...>::Options == Eigen::ColMajor)
+        {
+            cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),
+                        a.cols(), &one, a.data(), a.rows(), b.data(), b.rows(), &zero,
+                        res.data(), res.rows());
+        }
+    }
+
+    template <template <class, int...> class Mat, int... Opts>
+    static inline void mul(Mat<ComplexF, Opts...> &res, 
+                           const Mat<ComplexF, Opts...> &a, 
+                           const Mat<ComplexF, Opts...> &b)
+    {
+        static const ComplexF one(1., 0.), zero(0., 0.);
+
+        if (Mat<ComplexF, Opts...>::Options == Eigen::RowMajor)
+        {
+            cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),
+                        a.cols(), &one, a.data(), a.cols(), b.data(), b.cols(), &zero,
+                        res.data(), res.cols());
+        }
+        else if (Mat<ComplexF, Opts...>::Options == Eigen::ColMajor)
+        {
+            cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),
+                        a.cols(), &one, a.data(), a.rows(), b.data(), b.rows(), &zero,
+                        res.data(), res.rows());
+        }
+    }
+#else
+    template <typename Mat>
+    static inline void mul(Mat &res, const Mat &a, const Mat &b)
+    {
+        res = a*b;
+    }
+#endif
+
+private:
+    template <typename C, typename MatLeft, typename MatRight>
+    static inline void makeDotRowPt(C * &aPt, unsigned int &aInc, C * &bPt, 
+                                    unsigned int &bInc, const unsigned int aRow, 
+                                    const MatLeft &a, const MatRight &b)
+    {
+        if (MatLeft::Options == Eigen::RowMajor)
+        {
+            aPt  = a.data() + aRow*a.cols();
+            aInc = 1;
+        }
+        else if (MatLeft::Options == Eigen::ColMajor)
+        {
+            aPt  = a.data() + aRow;
+            aInc = a.rows();
+        }
+        if (MatRight::Options == Eigen::RowMajor)
+        {
+            bPt  = b.data() + aRow;
+            bInc = b.cols();
+        }
+        else if (MatRight::Options == Eigen::ColMajor)
+        {
+            bPt  = b.data() + aRow*b.rows();
+            bInc = 1;
+        }
+    }
+
+#ifdef USE_MKL
+    template <typename C, typename MatLeft, typename MatRight>
+    static inline void makeDotColPt(C * &aPt, unsigned int &aInc, C * &bPt, 
+                                    unsigned int &bInc, const unsigned int aCol, 
+                                    const MatLeft &a, const MatRight &b)
+    {
+        if (MatLeft::Options == Eigen::RowMajor)
+        {
+            aPt  = a.data() + aCol;
+            aInc = a.cols();
+        }
+        else if (MatLeft::Options == Eigen::ColMajor)
+        {
+            aPt  = a.data() + aCol*a.rows();
+            aInc = 1;
+        }
+        if (MatRight::Options == Eigen::RowMajor)
+        {
+            bPt  = b.data() + aCol*b.cols();
+            bInc = 1;
+        }
+        else if (MatRight::Options == Eigen::ColMajor)
+        {
+            bPt  = b.data() + aCol;
+            bInc = b.rows();
+        }
+    }
+
+    template <typename MatLeft, typename MatRight>
+    static inline void dotuRow(ComplexF &res, const unsigned int aRow,
+                               const MatLeft &a, const MatRight &b)
+    {
+        const ComplexF *aPt, *bPt;
+        unsigned int   aInc, bInc;
+
+        makeDotRowPt(aPt, aInc, bPt, bInc, aRow, a, b);
+        cblas_cdotu_sub(a.cols(), aPt, aInc, bPt, bInc, &res);
+    }
+
+    template <typename MatLeft, typename MatRight>
+    static inline void dotuCol(ComplexF &res, const unsigned int aCol,
+                               const MatLeft &a, const MatRight &b)
+    {
+        const ComplexF *aPt, *bPt;
+        unsigned int   aInc, bInc;
+
+        makeDotColPt(aPt, aInc, bPt, bInc, aCol, a, b);
+        cblas_cdotu_sub(a.rows(), aPt, aInc, bPt, bInc, &res);
+    }
+
+    template <typename MatLeft, typename MatRight>
+    static inline void dotuRow(ComplexD &res, const unsigned int aRow,
+                               const MatLeft &a, const MatRight &b)
+    {
+        const ComplexD *aPt, *bPt;
+        unsigned int   aInc, bInc;
+
+        makeDotRowPt(aPt, aInc, bPt, bInc, aRow, a, b);
+        cblas_zdotu_sub(a.cols(), aPt, aInc, bPt, bInc, &res);
+    }
+
+    template <typename MatLeft, typename MatRight>
+    static inline void dotuCol(ComplexD &res, const unsigned int aCol,
+                               const MatLeft &a, const MatRight &b)
+    {
+        const ComplexD *aPt, *bPt;
+        unsigned int   aInc, bInc;
+
+        makeDotColPt(aPt, aInc, bPt, bInc, aCol, a, b);
+        cblas_zdotu_sub(a.rows(), aPt, aInc, bPt, bInc, &res);
+    }
+#endif
+};
+
 /******************************************************************************
  *                     A2AMatrixIo template implementation                    *
  ******************************************************************************/
diff --git a/Hadrons/Utilities/ContractorBenchmark.cc b/Hadrons/Utilities/ContractorBenchmark.cc
index 2a8b033b..83d0e49f 100644
--- a/Hadrons/Utilities/ContractorBenchmark.cc
+++ b/Hadrons/Utilities/ContractorBenchmark.cc
@@ -1,16 +1,12 @@
 #include <Hadrons/Global.hpp>
+#include <Hadrons/A2AMatrix.hpp>
 #ifdef USE_MKL
 #include "mkl.h"
 #include "mkl_cblas.h"
 #endif
-#ifdef AVX2
-#include <immintrin.h>
-#endif
 
 using namespace Grid;
-
-typedef Eigen::Matrix<ComplexD, -1, -1, Eigen::RowMajor> RMCMat;
-typedef Eigen::Matrix<ComplexD, -1, -1, Eigen::ColMajor> CMCMat;
+using namespace Hadrons;
 
 #ifdef GRID_COMMS_MPI3
 #define GET_RANK(rank, nMpi) \
@@ -50,7 +46,7 @@ inline void trBenchmark(const std::string name, const MatLeft &left,
 
     if (rank == 0)
     {
-        std::cout << std::setw(30) << name << ": diff= "
+        std::cout << std::setw(34) << name << ": diff= "
                   << std::setw(12) << std::norm(buf-ref)
                   << std::setw(10) << t/1.0e6 << " sec "
                   << std::setw(10) << flops/t/1.0e3 << " GFlop/s " 
@@ -85,16 +81,17 @@ inline void mulBenchmark(const std::string name, const MatV &left,
 
     if (rank == 0)
     {
-        std::cout << std::setw(30) << name << ": diff= "
-                << std::setw(12) << (buf-ref).squaredNorm()
-                << std::setw(10) << t/1.0e6 << " sec "
-                << std::setw(10) << flops/t/1.0e3 << " GFlop/s " 
-                << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s "
-                << std::endl;
+        std::cout << std::setw(34) << name << ": diff= "
+                  << std::setw(12) << (buf-ref).squaredNorm()
+                  << std::setw(10) << t/1.0e6 << " sec "
+                  << std::setw(10) << flops/t/1.0e3 << " GFlop/s " 
+                  << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s "
+                  << std::endl;
     }
     ::sleep(1);
 }
 
+#ifdef USE_MKL
 template <typename MatLeft, typename MatRight>
 static inline void zdotuRow(ComplexD &res, const unsigned int aRow,
                             const MatLeft &a, const MatRight &b)
@@ -154,6 +151,7 @@ static inline void zdotuCol(ComplexD &res, const unsigned int aCol,
     }
     cblas_zdotu_sub(a.rows(), aPt, aInc, bPt, bInc, &res);
 }
+#endif
 
 template <typename MatLeft, typename MatRight>
 void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigned int nMat)
@@ -192,12 +190,18 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
     }
     BARRIER();
     ref = (left.back()*right.back()).trace();
+    trBenchmark("Hadrons A2AContraction::accTrMul", left, right, ref,
+    [](ComplexD &res, const MatLeft &a, const MatRight &b)
+    { 
+        res = 0.;
+        A2AContraction::accTrMul(res, a, b);
+    });
     trBenchmark("Naive loop rows first", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
     { 
         res = 0.;
         auto nr = a.rows(), nc = a.cols();
-        parallel_for_reduce(ComplexDPlus, res) (unsigned int i = 0; i < nr; ++i)
+        parallel_for_reduce(ComplexPlus, res) (unsigned int i = 0; i < nr; ++i)
         for (unsigned int j = 0; j < nc; ++j)
         {
             res += a(i, j)*b(j, i);
@@ -208,7 +212,7 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
     { 
         res = 0.;
         auto nr = a.rows(), nc = a.cols();
-        parallel_for_reduce(ComplexDPlus, res) (unsigned int j = 0; j < nc; ++j)
+        parallel_for_reduce(ComplexPlus, res) (unsigned int j = 0; j < nc; ++j)
         for (unsigned int i = 0; i < nr; ++i)
         {
             res += a(i, j)*b(j, i);
@@ -224,7 +228,7 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
     { 
         res = 0.;
 
-        parallel_for_reduce(ComplexDPlus, res) (unsigned int r = 0; r < a.rows(); ++r)
+        parallel_for_reduce(ComplexPlus, res) (unsigned int r = 0; r < a.rows(); ++r)
         {
             res += a.row(r).conjugate().dot(b.col(r));
         }
@@ -234,7 +238,7 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
     { 
         res = 0.;
 
-        parallel_for_reduce(ComplexDPlus, res) (unsigned int c = 0; c < a.cols(); ++c)
+        parallel_for_reduce(ComplexPlus, res) (unsigned int c = 0; c < a.cols(); ++c)
         {
             res += a.col(c).conjugate().dot(b.row(c));
         }
@@ -250,7 +254,7 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
     { 
         res = 0.;
 
-        parallel_for_reduce(ComplexDPlus, res) (unsigned int r = 0; r < a.rows(); ++r)
+        parallel_for_reduce(ComplexPlus, res) (unsigned int r = 0; r < a.rows(); ++r)
         {
             ComplexD tmp;
 
@@ -263,7 +267,7 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
     { 
         res = 0.;
 
-        parallel_for_reduce(ComplexDPlus, res) (unsigned int c = 0; c < a.cols(); ++c)
+        parallel_for_reduce(ComplexPlus, res) (unsigned int c = 0; c < a.cols(); ++c)
         {
             ComplexD tmp;
 
@@ -305,6 +309,11 @@ void fullMulBenchmark(const unsigned int ni, const unsigned int nj, const unsign
     }
     BARRIER();
     ref = left.back()*right.back();
+    mulBenchmark("Hadrons A2AContraction::mul", left, right, ref,
+    [](Mat &res, const Mat &a, const Mat &b)
+    { 
+        A2AContraction::mul(res, a, b);
+    });
     mulBenchmark("Eigen A*B", left, right, ref,
     [](Mat &res, const Mat &a, const Mat &b)
     { 
@@ -381,12 +390,12 @@ int main(int argc, char *argv[])
         std::cout << std::endl;
     }
 
-    fullTrBenchmark<RMCMat, RMCMat>(ni, nj, nMat);
-    fullTrBenchmark<RMCMat, CMCMat>(ni, nj, nMat);
-    fullTrBenchmark<CMCMat, RMCMat>(ni, nj, nMat);
-    fullTrBenchmark<CMCMat, CMCMat>(ni, nj, nMat);
-    fullMulBenchmark<RMCMat>(ni, nj, nMat);
-    fullMulBenchmark<CMCMat>(ni, nj, nMat);
+    fullTrBenchmark<A2AMatrix<ComplexD>, A2AMatrix<ComplexD>>(ni, nj, nMat);
+    fullTrBenchmark<A2AMatrix<ComplexD>, A2AMatrixTr<ComplexD>>(ni, nj, nMat);
+    fullTrBenchmark<A2AMatrixTr<ComplexD>, A2AMatrix<ComplexD>>(ni, nj, nMat);
+    fullTrBenchmark<A2AMatrixTr<ComplexD>, A2AMatrixTr<ComplexD>>(ni, nj, nMat);
+    fullMulBenchmark<A2AMatrix<ComplexD>>(ni, nj, nMat);
+    fullMulBenchmark<A2AMatrixTr<ComplexD>>(ni, nj, nMat);
     FINALIZE();
 
     return EXIT_SUCCESS;

From 0c6e581336036dd1323bdf655f14ea73067563d4 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Wed, 7 Nov 2018 19:16:55 +0000
Subject: [PATCH 12/33] Hadrons: first stab at general contraction code, needs
 serious testing

---
 Hadrons/DiskVector.hpp          |   3 +-
 Hadrons/Utilities/Contractor.cc | 124 +++++++++++++++++++++-----------
 2 files changed, 85 insertions(+), 42 deletions(-)

diff --git a/Hadrons/DiskVector.hpp b/Hadrons/DiskVector.hpp
index 1fddbe36..89265326 100644
--- a/Hadrons/DiskVector.hpp
+++ b/Hadrons/DiskVector.hpp
@@ -29,6 +29,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #define Hadrons_DiskVector_hpp_
 
 #include <Hadrons/Global.hpp>
+#include <Hadrons/A2AMatrix.hpp>
 #include <deque>
 #include <sys/stat.h>
 #include <ftw.h>
@@ -143,7 +144,7 @@ private:
  *                      Specialisation for Eigen matrices                     *
  ******************************************************************************/
 template <typename T>
-using EigenDiskVectorMat = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
+using EigenDiskVectorMat = A2AMatrix<T>;
 
 template <typename T>
 class EigenDiskVector: public DiskVectorBase<EigenDiskVectorMat<T>>
diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index a66ad373..5239e3c9 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -32,6 +32,8 @@ using namespace Grid;
 using namespace QCD;
 using namespace Hadrons;
 
+#define TIME_MOD(t) (((t) + par.global.nt) % par.global.nt)
+
 namespace Contractor
 {
     class GlobalPar: Serializable
@@ -58,7 +60,8 @@ namespace Contractor
     public:
         GRID_SERIALIZABLE_CLASS_MEMBERS(ProductPar,
                                         std::string, terms,
-                                        std::vector<std::string>, timeRange);
+                                        std::vector<std::string>, times,
+                                        std::string, translations);
     };
 }
 
@@ -69,6 +72,34 @@ struct ContractorPar
     std::vector<Contractor::ProductPar>    product;
 };
 
+void makeTimeSeq(std::vector<std::vector<unsigned int>> &timeSeq, 
+                 const std::vector<std::set<unsigned int>> &times,
+                 std::vector<unsigned int> &current,
+                 const unsigned int depth)
+{
+    if (depth > 0)
+    {
+        for (auto t: times[times.size() - depth])
+        {
+            current[times.size() - depth] = t;
+            makeTimeSeq(timeSeq, times, current, depth - 1);
+        }
+    }
+    else
+    {
+        timeSeq.push_back(current);
+    }
+}
+
+void makeTimeSeq(std::vector<std::vector<unsigned int>> &timeSeq, 
+                 const std::vector<std::set<unsigned int>> &times)
+{
+    std::vector<unsigned int> current(times.size());
+
+    makeTimeSeq(timeSeq, times, current, times.size());
+}
+
+
 std::set<unsigned int> parseTimeRange(const std::string str, const unsigned int nt)
 {
     std::regex               rex("([0-9]+)|(([0-9]+)\\.\\.([0-9]+))");
@@ -153,7 +184,7 @@ int main(int argc, char* argv[])
         auto   &p = par.a2aMatrix[i];
         double t, size;
 
-        std::cout << "-- Loading '" << p.file << "'..." << std::endl;
+        std::cout << "======== Loading '" << p.file << "'" << std::endl;
 
         A2AMatrixIo<HADRONS_A2AM_IO_TYPE> a2aIo(p.file, p.dataset, par.global.nt);
 
@@ -166,58 +197,69 @@ int main(int argc, char* argv[])
 
     for (auto &p: par.product)
     {
-        std::vector<std::string>            term = strToVec<std::string>(p.terms);
-        std::vector<std::set<unsigned int>> times;
+        std::vector<std::string>               term = strToVec<std::string>(p.terms);
+        std::vector<std::set<unsigned int>>    times;
+        std::vector<std::vector<unsigned int>> timeSeq;
+        std::set<unsigned int>                 translations;
+        std::vector<ComplexD>                  corr(par.global.nt);
+        std::vector<A2AMatrixTr<ComplexD>>     lastTerm(par.global.nt);
+        A2AMatrix<ComplexD>                    prod, buf, tmp;
 
-        if (term.size() != p.timeRange.size())
+        std::cout << "======== Product tr(";
+        for (unsigned int g = 0; g < term.size(); ++g)
+        {
+            std::cout << term[g] << ((g == term.size() - 1) ? ')' : '*');
+        }
+        std::cout << std::endl;
+        if (term.size() != p.times.size() + 1)
         {
             HADRONS_ERROR(Size, "number of terms (" + std::to_string(term.size()) 
-                          + ") different from number of time ranges (" 
-                          + std::to_string(p.timeRange.size()) + ")");
+                          + ") different from number of times (" 
+                          + std::to_string(p.times.size() + 1) + ")");
         }
-        for (auto &s: p.timeRange)
+        for (auto &s: p.times)
         {
             times.push_back(parseTimeRange(s, par.global.nt));
         }
-        
-        // test: just 2-pt function for now
-        if (term.size() != 2)
-        {
-            HADRONS_ERROR(Implementation, "only 2-pt function implemented");
-        }
+        translations = parseTimeRange(p.translations, par.global.nt);
+        makeTimeSeq(timeSeq, times);
+        std::cout << timeSeq.size()*translations.size()*(term.size() - 2) << " A*B, "
+                  << timeSeq.size()*translations.size()*par.global.nt << " tr(A*B)"
+                  << std::endl;
 
-        std::vector<ComplexD> corr(par.global.nt);
-        double                tTrace = 0., flops, bytes;
-
-#define TIME_MOD(t) (((t) + par.global.nt) % par.global.nt)
+        std::cout << "-- caching transposed last term" << std::endl;
         for (unsigned int t = 0; t < par.global.nt; ++t)
         {
-            corr[t] = ComplexD(0., 0.);
+            buf         = a2aMat.at(term.back())[t];
+            lastTerm[t] = buf;
         }
-        for (auto &trans: times[1])
-        for (auto &t0: times[0])
+        for (auto &t: timeSeq)
         {
-            unsigned int tt0 = TIME_MOD(t0 + trans);
-            EigenDiskVector<ComplexD>::Matrix m0 = a2aMat.at(term[0])[tt0];
-            EigenDiskVector<ComplexD>::Matrix m1 = a2aMat.at(term[1])[trans];
-
-            std::cout << "tr(" << term[0] << "[" << tt0 << "]*"
-                      << term[1] << "[" << trans << "])";
-            tTrace  = -usecond();
-            corr[t0] += (m0*m1).trace();
-            tTrace += usecond();
-            flops = 6.*m0.rows()*m0.cols() + 2.*(m0.rows()*m0.cols() - 1.);
-            bytes = 2.*m0.rows()*m0.cols()*sizeof(ComplexD);
-            std::cout << " -- Perf " << flops/tTrace/1.0e3 << " GFlop/s " 
-                      << bytes/tTrace*1.0e6/1024/1024/1024 << " GB/s "<< std::endl;
+            for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
+            {
+                corr[tLast] = 0.;
+            }
+            for (auto &dt: translations)
+            {
+                std::cout << "-- position " << t << ", translation " << dt << std::endl;
+                prod = a2aMat.at(term[0])[TIME_MOD(t[0] + dt)];
+                for (unsigned int i = 1; i < term.size() - 1; ++i)
+                {
+                    const A2AMatrix<ComplexD> &ref = a2aMat.at(term[i])[TIME_MOD(t[i] + dt)];
+                    A2AContraction::mul(tmp, prod, ref);
+                    prod = tmp;
+                }
+                for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
+                {
+                    A2AContraction::accTrMul(corr[tLast], prod, lastTerm[tLast]);
+                }
+            }
+            for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
+            {
+                std::cout << tLast << " " << corr[tLast] << std::endl;
+            }
         }
-        for (unsigned int t = 0; t < par.global.nt; ++t)
-        {
-            std::cout << t << " " << corr[t] << std::endl;
-        }
-#undef TIME_MOD
-   }
-
+    }
     
     return EXIT_SUCCESS;
 }

From 32376f0437c6ed3d05143c7bc4b73f470e1c5940 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Wed, 7 Nov 2018 19:59:11 +0000
Subject: [PATCH 13/33] Hadrons: contractor performances

---
 Hadrons/A2AMatrix.hpp           | 14 ++++++++++++
 Hadrons/Utilities/Contractor.cc | 40 ++++++++++++++++++++++++++++++++-
 2 files changed, 53 insertions(+), 1 deletion(-)

diff --git a/Hadrons/A2AMatrix.hpp b/Hadrons/A2AMatrix.hpp
index 5266d050..a2c2ffcf 100644
--- a/Hadrons/A2AMatrix.hpp
+++ b/Hadrons/A2AMatrix.hpp
@@ -196,6 +196,14 @@ public:
         }
     }
 
+    template <typename C, typename MatLeft, typename MatRight>
+    static inline double accTrMulFlops(C &acc, const MatLeft &a, const MatRight &b)
+    {
+        double n = a.rows()*a.cols();
+
+        return 8.*n;
+    }
+
     // mul(res, a, b): res = a*b
 #ifdef USE_MKL
     template <template <class, int...> class Mat, int... Opts>
@@ -246,7 +254,13 @@ public:
         res = a*b;
     }
 #endif
+    template <typename Mat>
+    static inline double mulFlops(Mat &res, const Mat &a, const Mat &b)
+    {
+        double nr = a.rows(), nc = a.cols();
 
+        return nr*nr*(6.*nc + 2.*(nc - 1.));
+    }
 private:
     template <typename C, typename MatLeft, typename MatRight>
     static inline void makeDotRowPt(C * &aPt, unsigned int &aInc, C * &bPt, 
diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index 5239e3c9..d67d5b9d 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -142,6 +142,13 @@ std::set<unsigned int> parseTimeRange(const std::string str, const unsigned int
     return tSet;
 }
 
+void printPerf(const double flops, const double fusec, const double bytes, const double busec)
+{
+    std::cout << std::setw(10) << flops/fusec/1.0e3 << " GFlop/s " 
+              << std::setw(10) << bytes/busec*1.0e6/1024/1024/1024 << " GB/s"
+              << std::endl;
+}
+
 int main(int argc, char* argv[])
 {
     // parse command line
@@ -204,6 +211,7 @@ int main(int argc, char* argv[])
         std::vector<ComplexD>                  corr(par.global.nt);
         std::vector<A2AMatrixTr<ComplexD>>     lastTerm(par.global.nt);
         A2AMatrix<ComplexD>                    prod, buf, tmp;
+        double                                 fusec, busec, flops, bytes;
 
         std::cout << "======== Product tr(";
         for (unsigned int g = 0; g < term.size(); ++g)
@@ -242,17 +250,47 @@ int main(int argc, char* argv[])
             for (auto &dt: translations)
             {
                 std::cout << "-- position " << t << ", translation " << dt << std::endl;
+                if (term.size() > 2)
+                {
+                    std::cout << "* matrix products" << std::endl;
+                }
+                flops  = 0.;
+                bytes  = 0.;
+                fusec  = 0.;
+                busec  = 0.;
                 prod = a2aMat.at(term[0])[TIME_MOD(t[0] + dt)];
                 for (unsigned int i = 1; i < term.size() - 1; ++i)
                 {
                     const A2AMatrix<ComplexD> &ref = a2aMat.at(term[i])[TIME_MOD(t[i] + dt)];
+                    fusec -= usecond();
+                    busec -= usecond();
                     A2AContraction::mul(tmp, prod, ref);
-                    prod = tmp;
+                    fusec += usecond();
+                    flops += A2AContraction::mulFlops(tmp, prod, ref);
+                    prod   = tmp;
+                    busec += usecond();
+                    bytes += 3.*tmp.rows()*tmp.cols()*sizeof(ComplexD);
                 }
+                if (term.size() > 2)
+                {
+                    printPerf(flops, fusec, bytes, busec);
+                }
+                std::cout << "* traces" << std::endl;
+                flops  = 0.;
+                bytes  = 0.;
+                fusec  = 0.;
+                busec  = 0.;
                 for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
                 {
+                    fusec -= usecond();
+                    busec -= usecond();
                     A2AContraction::accTrMul(corr[tLast], prod, lastTerm[tLast]);
+                    fusec += usecond();
+                    busec += usecond();
+                    flops += A2AContraction::accTrMulFlops(corr[tLast], prod, lastTerm[tLast]);
+                    bytes += 2.*prod.rows()*prod.cols()*sizeof(ComplexD);
                 }
+                printPerf(flops, fusec, bytes, busec);
             }
             for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
             {

From 1ff1422e075b9c2ab5c2f18aa1f3397c3cc0c536 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Wed, 7 Nov 2018 20:02:53 +0000
Subject: [PATCH 14/33] Hadrons: contractor lighter output

---
 Hadrons/Utilities/Contractor.cc | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index d67d5b9d..ba1c139a 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -252,7 +252,7 @@ int main(int argc, char* argv[])
                 std::cout << "-- position " << t << ", translation " << dt << std::endl;
                 if (term.size() > 2)
                 {
-                    std::cout << "* matrix products" << std::endl;
+                    std::cout << "* matrix products ";
                 }
                 flops  = 0.;
                 bytes  = 0.;
@@ -275,7 +275,7 @@ int main(int argc, char* argv[])
                 {
                     printPerf(flops, fusec, bytes, busec);
                 }
-                std::cout << "* traces" << std::endl;
+                std::cout << "* traces ";
                 flops  = 0.;
                 bytes  = 0.;
                 fusec  = 0.;

From 7cd9914f0e9f4efe300fda758e0b9d9d25e18fd6 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Thu, 8 Nov 2018 17:40:57 +0000
Subject: [PATCH 15/33] Hadrons: automatically resize output in MKL A2A matrix
 kernels

---
 Hadrons/A2AMatrix.hpp | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/Hadrons/A2AMatrix.hpp b/Hadrons/A2AMatrix.hpp
index a2c2ffcf..2f958685 100644
--- a/Hadrons/A2AMatrix.hpp
+++ b/Hadrons/A2AMatrix.hpp
@@ -213,6 +213,10 @@ public:
     {
         static const ComplexD one(1., 0.), zero(0., 0.);
 
+        if ((res.rows() != a.rows()) or (res.cols() != b.cols()))
+        {
+            res.resize(a.rows(), b.cols());
+        }
         if (Mat<ComplexD, Opts...>::Options == Eigen::RowMajor)
         {
             cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),
@@ -234,6 +238,10 @@ public:
     {
         static const ComplexF one(1., 0.), zero(0., 0.);
 
+        if ((res.rows() != a.rows()) or (res.cols() != b.cols()))
+        {
+            res.resize(a.rows(), b.cols());
+        }
         if (Mat<ComplexF, Opts...>::Options == Eigen::RowMajor)
         {
             cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, a.rows(), b.cols(),

From 65349b07a76392d3091eb5e28eb70fc0bb122929 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Thu, 8 Nov 2018 18:44:44 +0000
Subject: [PATCH 16/33] Hadrons: simpler A2A perf functions

---
 Hadrons/A2AMatrix.hpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Hadrons/A2AMatrix.hpp b/Hadrons/A2AMatrix.hpp
index 2f958685..831d2a02 100644
--- a/Hadrons/A2AMatrix.hpp
+++ b/Hadrons/A2AMatrix.hpp
@@ -196,8 +196,8 @@ public:
         }
     }
 
-    template <typename C, typename MatLeft, typename MatRight>
-    static inline double accTrMulFlops(C &acc, const MatLeft &a, const MatRight &b)
+    template <typename MatLeft, typename MatRight>
+    static inline double accTrMulFlops(const MatLeft &a, const MatRight &b)
     {
         double n = a.rows()*a.cols();
 
@@ -263,7 +263,7 @@ public:
     }
 #endif
     template <typename Mat>
-    static inline double mulFlops(Mat &res, const Mat &a, const Mat &b)
+    static inline double mulFlops(const Mat &a, const Mat &b)
     {
         double nr = a.rows(), nc = a.cols();
 

From 69b6ba0a733ebbce3efacc7234a003574dd36340 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Thu, 8 Nov 2018 18:46:28 +0000
Subject: [PATCH 17/33] Hadrons: contractor fixes and improvements

---
 Hadrons/Utilities/Contractor.cc | 54 +++++++++++++++++++++++++--------
 1 file changed, 41 insertions(+), 13 deletions(-)

diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index ba1c139a..cb81ba88 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -142,11 +142,38 @@ std::set<unsigned int> parseTimeRange(const std::string str, const unsigned int
     return tSet;
 }
 
-void printPerf(const double flops, const double fusec, const double bytes, const double busec)
+struct Flops
 {
-    std::cout << std::setw(10) << flops/fusec/1.0e3 << " GFlop/s " 
-              << std::setw(10) << bytes/busec*1.0e6/1024/1024/1024 << " GB/s"
-              << std::endl;
+    Flops(const double flops, const double fusec)
+    {
+        gFlopsPerSec = flops/fusec/1.0e3;
+    }
+    
+    double gFlopsPerSec;
+};
+
+inline std::ostream & operator<< (std::ostream& s, const Flops &&f)
+{
+    s << std::setw(10) << f.gFlopsPerSec << " GFlop/s";
+
+    return s;
+}
+
+struct Bytes
+{
+    Bytes(const double bytes, const double busec)
+    {
+        gBytesPerSec = bytes/busec*1.0e6/1024/1024/1024;
+    }
+    
+    double gBytesPerSec;
+};
+
+inline std::ostream & operator<< (std::ostream& s, const Bytes &&b)
+{
+    s << std::setw(10) << b.gBytesPerSec << " GB/s";
+
+    return s;
 }
 
 int main(int argc, char* argv[])
@@ -238,8 +265,9 @@ int main(int argc, char* argv[])
         std::cout << "-- caching transposed last term" << std::endl;
         for (unsigned int t = 0; t < par.global.nt; ++t)
         {
-            buf         = a2aMat.at(term.back())[t];
-            lastTerm[t] = buf;
+            const A2AMatrix<ComplexD> &ref = a2aMat.at(term.back())[t];
+
+            lastTerm[t] = ref;
         }
         for (auto &t: timeSeq)
         {
@@ -252,7 +280,7 @@ int main(int argc, char* argv[])
                 std::cout << "-- position " << t << ", translation " << dt << std::endl;
                 if (term.size() > 2)
                 {
-                    std::cout << "* matrix products ";
+                    std::cout << "*" << std::setw(12) << "products";
                 }
                 flops  = 0.;
                 bytes  = 0.;
@@ -266,16 +294,16 @@ int main(int argc, char* argv[])
                     busec -= usecond();
                     A2AContraction::mul(tmp, prod, ref);
                     fusec += usecond();
-                    flops += A2AContraction::mulFlops(tmp, prod, ref);
+                    flops += A2AContraction::mulFlops(prod, ref);
                     prod   = tmp;
                     busec += usecond();
                     bytes += 3.*tmp.rows()*tmp.cols()*sizeof(ComplexD);
                 }
                 if (term.size() > 2)
                 {
-                    printPerf(flops, fusec, bytes, busec);
+                    std::cout << Flops(flops, fusec) << " " << Bytes(bytes, busec) << std::endl;
                 }
-                std::cout << "* traces ";
+                std::cout << "*" << std::setw(12) << "traces";
                 flops  = 0.;
                 bytes  = 0.;
                 fusec  = 0.;
@@ -284,13 +312,13 @@ int main(int argc, char* argv[])
                 {
                     fusec -= usecond();
                     busec -= usecond();
-                    A2AContraction::accTrMul(corr[tLast], prod, lastTerm[tLast]);
+                    A2AContraction::accTrMul(corr[TIME_MOD(tLast - dt)], prod, lastTerm[tLast]);
                     fusec += usecond();
                     busec += usecond();
-                    flops += A2AContraction::accTrMulFlops(corr[tLast], prod, lastTerm[tLast]);
+                    flops += A2AContraction::accTrMulFlops(prod, lastTerm[tLast]);
                     bytes += 2.*prod.rows()*prod.cols()*sizeof(ComplexD);
                 }
-                printPerf(flops, fusec, bytes, busec);
+                std::cout << Flops(flops, fusec) << " " << Bytes(bytes, busec) << std::endl;
             }
             for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
             {

From 68c76a410d8a0f29e4ea23ade4523cbbafa4f616 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Thu, 8 Nov 2018 19:24:29 +0000
Subject: [PATCH 18/33] Hadrons: more contractor improvements

---
 Hadrons/Utilities/Contractor.cc | 18 ++++++++++++------
 1 file changed, 12 insertions(+), 6 deletions(-)

diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index cb81ba88..40786221 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -269,27 +269,32 @@ int main(int argc, char* argv[])
 
             lastTerm[t] = ref;
         }
-        for (auto &t: timeSeq)
+        for (unsigned int i = 0; i < timeSeq.size(); ++i)
         {
+            unsigned int dti = 0;
+            auto         &t = timeSeq[i];
+
             for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
             {
                 corr[tLast] = 0.;
             }
             for (auto &dt: translations)
             {
-                std::cout << "-- position " << t << ", translation " << dt << std::endl;
+                std::cout << "* Step " << i*translations.size() + dti + 1
+                          << "/" << timeSeq.size()*translations.size()
+                          << " -- positions= " << t << ", dt= " << dt << std::endl;
                 if (term.size() > 2)
                 {
-                    std::cout << "*" << std::setw(12) << "products";
+                    std::cout << std::setw(8) << "products";
                 }
                 flops  = 0.;
                 bytes  = 0.;
                 fusec  = 0.;
                 busec  = 0.;
                 prod = a2aMat.at(term[0])[TIME_MOD(t[0] + dt)];
-                for (unsigned int i = 1; i < term.size() - 1; ++i)
+                for (unsigned int j = 1; j < term.size() - 1; ++j)
                 {
-                    const A2AMatrix<ComplexD> &ref = a2aMat.at(term[i])[TIME_MOD(t[i] + dt)];
+                    const A2AMatrix<ComplexD> &ref = a2aMat.at(term[j])[TIME_MOD(t[j] + dt)];
                     fusec -= usecond();
                     busec -= usecond();
                     A2AContraction::mul(tmp, prod, ref);
@@ -303,7 +308,7 @@ int main(int argc, char* argv[])
                 {
                     std::cout << Flops(flops, fusec) << " " << Bytes(bytes, busec) << std::endl;
                 }
-                std::cout << "*" << std::setw(12) << "traces";
+                std::cout << std::setw(8) << "traces";
                 flops  = 0.;
                 bytes  = 0.;
                 fusec  = 0.;
@@ -319,6 +324,7 @@ int main(int argc, char* argv[])
                     bytes += 2.*prod.rows()*prod.cols()*sizeof(ComplexD);
                 }
                 std::cout << Flops(flops, fusec) << " " << Bytes(bytes, busec) << std::endl;
+                dti++;
             }
             for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
             {

From f5422c7334e4ebd247f20be4bd45c8512531b41f Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Fri, 9 Nov 2018 16:23:53 +0000
Subject: [PATCH 19/33] Hadrons: more contractor instrumentation

---
 Hadrons/Utilities/Contractor.cc | 70 ++++++++++++++++++++++++---------
 1 file changed, 52 insertions(+), 18 deletions(-)

diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index 40786221..f113e954 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -27,6 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Hadrons/Global.hpp>
 #include <Hadrons/A2AMatrix.hpp>
 #include <Hadrons/DiskVector.hpp>
+#include <Hadrons/TimerArray.hpp>
 
 using namespace Grid;
 using namespace QCD;
@@ -142,6 +143,23 @@ std::set<unsigned int> parseTimeRange(const std::string str, const unsigned int
     return tSet;
 }
 
+struct Sec
+{
+    Sec(const double usec)
+    {
+        seconds = usec/1.0e6;
+    }
+    
+    double seconds;
+};
+
+inline std::ostream & operator<< (std::ostream& s, const Sec &&sec)
+{
+    s << std::setw(10) << sec.seconds << " sec";
+
+    return s;
+}
+
 struct Flops
 {
     Flops(const double flops, const double fusec)
@@ -223,7 +241,8 @@ int main(int argc, char* argv[])
         A2AMatrixIo<HADRONS_A2AM_IO_TYPE> a2aIo(p.file, p.dataset, par.global.nt);
 
         a2aIo.load(a2aMat.at(p.name), &t);
-        std::cout << "Read " << a2aIo.getSize() << " bytes in " << t << " usec, " << a2aIo.getSize()/t*1.0e6/1024/1024 << " MB/s" << std::endl;
+        std::cout << "Read " << a2aIo.getSize() << " bytes in " << t/1.0e6 
+                  << " sec, " << a2aIo.getSize()/t*1.0e6/1024/1024 << " MB/s" << std::endl;
     }
 
     // contract
@@ -238,9 +257,10 @@ int main(int argc, char* argv[])
         std::vector<ComplexD>                  corr(par.global.nt);
         std::vector<A2AMatrixTr<ComplexD>>     lastTerm(par.global.nt);
         A2AMatrix<ComplexD>                    prod, buf, tmp;
-        double                                 fusec, busec, flops, bytes;
+        TimerArray                             tAr;
+        double                                 fusec, busec, flops, bytes, tusec;
 
-        std::cout << "======== Product tr(";
+        std::cout << "======== Contraction tr(";
         for (unsigned int g = 0; g < term.size(); ++g)
         {
             std::cout << term[g] << ((g == term.size() - 1) ? ')' : '*');
@@ -262,13 +282,18 @@ int main(int argc, char* argv[])
                   << timeSeq.size()*translations.size()*par.global.nt << " tr(A*B)"
                   << std::endl;
 
-        std::cout << "-- caching transposed last term" << std::endl;
+        std::cout << "* Caching transposed last term" << std::endl;
         for (unsigned int t = 0; t < par.global.nt; ++t)
         {
             const A2AMatrix<ComplexD> &ref = a2aMat.at(term.back())[t];
 
+            tAr.startTimer("Transpose caching");
             lastTerm[t] = ref;
+            tAr.stopTimer("Transpose caching");
         }
+        bytes = par.global.nt*lastTerm[0].rows()*lastTerm[0].cols()*sizeof(ComplexD);
+        std::cout << Sec(tAr.getDTimer("Transpose caching")) << " " 
+                  << Bytes(bytes, tAr.getDTimer("Transpose caching")) << std::endl;
         for (unsigned int i = 0; i < timeSeq.size(); ++i)
         {
             unsigned int dti = 0;
@@ -289,41 +314,50 @@ int main(int argc, char* argv[])
                 }
                 flops  = 0.;
                 bytes  = 0.;
-                fusec  = 0.;
-                busec  = 0.;
+                fusec  = tAr.getDTimer("A*B algebra");
+                busec  = tAr.getDTimer("A*B total");
+                tAr.startTimer("Linear algebra");
+                tAr.startTimer("Disk vector overhead");
                 prod = a2aMat.at(term[0])[TIME_MOD(t[0] + dt)];
+                tAr.stopTimer("Disk vector overhead");
                 for (unsigned int j = 1; j < term.size() - 1; ++j)
                 {
+                    tAr.startTimer("Disk vector overhead");
                     const A2AMatrix<ComplexD> &ref = a2aMat.at(term[j])[TIME_MOD(t[j] + dt)];
-                    fusec -= usecond();
-                    busec -= usecond();
+                    tAr.stopTimer("Disk vector overhead");
+                    
+                    tAr.startTimer("A*B total");
+                    tAr.startTimer("A*B algebra");
                     A2AContraction::mul(tmp, prod, ref);
-                    fusec += usecond();
+                    tAr.stopTimer("A*B algebra");
                     flops += A2AContraction::mulFlops(prod, ref);
                     prod   = tmp;
-                    busec += usecond();
+                    tAr.stopTimer("A*B total");
                     bytes += 3.*tmp.rows()*tmp.cols()*sizeof(ComplexD);
                 }
                 if (term.size() > 2)
                 {
-                    std::cout << Flops(flops, fusec) << " " << Bytes(bytes, busec) << std::endl;
+                    std::cout << Sec(tAr.getDTimer("A*B total") - busec) << " "
+                              << Flops(flops, tAr.getDTimer("A*B algebra") - fusec) << " " 
+                              << Bytes(bytes, tAr.getDTimer("A*B total") - busec) << std::endl;
                 }
                 std::cout << std::setw(8) << "traces";
                 flops  = 0.;
                 bytes  = 0.;
-                fusec  = 0.;
-                busec  = 0.;
+                fusec  = tAr.getDTimer("tr(A*B)");
+                busec  = tAr.getDTimer("tr(A*B)");
                 for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
                 {
-                    fusec -= usecond();
-                    busec -= usecond();
+                    tAr.startTimer("tr(A*B)");
                     A2AContraction::accTrMul(corr[TIME_MOD(tLast - dt)], prod, lastTerm[tLast]);
-                    fusec += usecond();
-                    busec += usecond();
+                    tAr.stopTimer("tr(A*B)");
                     flops += A2AContraction::accTrMulFlops(prod, lastTerm[tLast]);
                     bytes += 2.*prod.rows()*prod.cols()*sizeof(ComplexD);
                 }
-                std::cout << Flops(flops, fusec) << " " << Bytes(bytes, busec) << std::endl;
+                tAr.stopTimer("Linear algebra");
+                std::cout << Sec(tAr.getDTimer("tr(A*B)") - busec) << " "
+                          << Flops(flops, tAr.getDTimer("tr(A*B)") - fusec) << " " 
+                          << Bytes(bytes, tAr.getDTimer("tr(A*B)") - busec) << std::endl;
                 dti++;
             }
             for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)

From cfe281f1a49a4ce7f2457978110a415e23e490b3 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Mon, 12 Nov 2018 15:59:11 +0000
Subject: [PATCH 20/33] Hadrons: diskvectors measure hash performance in debug
 output

---
 Hadrons/DiskVector.hpp | 28 +++++++++++++++++-----------
 1 file changed, 17 insertions(+), 11 deletions(-)

diff --git a/Hadrons/DiskVector.hpp b/Hadrons/DiskVector.hpp
index 89265326..6d96c6d9 100644
--- a/Hadrons/DiskVector.hpp
+++ b/Hadrons/DiskVector.hpp
@@ -165,19 +165,22 @@ private:
         uint32_t      crc, check;
         Eigen::Index  nRow, nCol;
         size_t        matSize;
-        double        t;
+        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);
-        t  = -usecond();
+        tRead  = -usecond();
         f.read(reinterpret_cast<char *>(obj.data()), matSize);
-        t += usecond();
-        DV_DEBUG_MSG(this, "Eigen read " << matSize/t*1.0e6/1024/1024 << " MB/s");
-        check = GridChecksum::crc32(obj.data(), matSize);
-        DV_DEBUG_MSG(this, "Eigen crc32 " << std::hex << check << std::dec);
+        tRead += usecond();
+        tHash  = -usecond();
+        check  = GridChecksum::crc32(obj.data(), matSize);
+        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")
@@ -190,20 +193,23 @@ private:
         uint32_t      crc;
         Eigen::Index  nRow, nCol;
         size_t        matSize;
-        double        t;
+        double        tWrite, tHash;
         
         nRow    = obj.rows();
         nCol    = obj.cols();
         matSize = nRow*nCol*sizeof(T);
+        tHash   = -usecond();
         crc     = GridChecksum::crc32(obj.data(), matSize);
+        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));
-        t  = -usecond();
+        tWrite = -usecond();
         f.write(reinterpret_cast<const char *>(obj.data()), matSize);
-        t += usecond();
-        DV_DEBUG_MSG(this, "Eigen write " << matSize/t*1.0e6/1024/1024 << " MB/s");
-        DV_DEBUG_MSG(this, "Eigen crc32 " << std::hex << crc << std::dec);
+        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");
     }
 };
 

From 344e832a4e41c6d76bafb4c7375a001aaf7ffb18 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Mon, 12 Nov 2018 15:59:54 +0000
Subject: [PATCH 21/33] Hadrons: contractor faster transpose and finer timings

---
 Hadrons/Utilities/Contractor.cc | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index f113e954..3fdbb872 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -260,6 +260,7 @@ int main(int argc, char* argv[])
         TimerArray                             tAr;
         double                                 fusec, busec, flops, bytes, tusec;
 
+        tAr.startTimer("Total");
         std::cout << "======== Contraction tr(";
         for (unsigned int g = 0; g < term.size(); ++g)
         {
@@ -276,6 +277,7 @@ int main(int argc, char* argv[])
         {
             times.push_back(parseTimeRange(s, par.global.nt));
         }
+
         translations = parseTimeRange(p.translations, par.global.nt);
         makeTimeSeq(timeSeq, times);
         std::cout << timeSeq.size()*translations.size()*(term.size() - 2) << " A*B, "
@@ -285,15 +287,23 @@ int main(int argc, char* argv[])
         std::cout << "* Caching transposed last term" << std::endl;
         for (unsigned int t = 0; t < par.global.nt; ++t)
         {
+            tAr.startTimer("Disk vector overhead");
             const A2AMatrix<ComplexD> &ref = a2aMat.at(term.back())[t];
+            tAr.stopTimer("Disk vector overhead");
 
             tAr.startTimer("Transpose caching");
-            lastTerm[t] = ref;
+            lastTerm[t].resize(ref.rows(), ref.cols());
+            parallel_for (unsigned int j = 0; j < ref.cols(); ++j)
+            for (unsigned int i = 0; i < ref.rows(); ++i)
+            {
+                lastTerm[t](i, j) = ref(i, j);
+            }
             tAr.stopTimer("Transpose caching");
         }
         bytes = par.global.nt*lastTerm[0].rows()*lastTerm[0].cols()*sizeof(ComplexD);
         std::cout << Sec(tAr.getDTimer("Transpose caching")) << " " 
                   << Bytes(bytes, tAr.getDTimer("Transpose caching")) << std::endl;
+        std::cout << Sec(tAr.getDTimer("Disk vector overhead")) << std::endl;
         for (unsigned int i = 0; i < timeSeq.size(); ++i)
         {
             unsigned int dti = 0;
@@ -365,6 +375,8 @@ int main(int argc, char* argv[])
                 std::cout << tLast << " " << corr[tLast] << std::endl;
             }
         }
+        tAr.stopTimer("Total");
+        printTimeProfile(tAr.getTimings(), tAr.getTimer("Total"));
     }
     
     return EXIT_SUCCESS;

From 2ac57370f10423bd7e8351d3844afa35890ec3bb Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Mon, 12 Nov 2018 16:04:35 +0000
Subject: [PATCH 22/33] Hadrons: contractor translation average normalisation

---
 Hadrons/Utilities/Contractor.cc | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index 3fdbb872..ffc77cf7 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -371,6 +371,10 @@ int main(int argc, char* argv[])
                 dti++;
             }
             for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
+            {
+                corr[tLast] /= translations.size();
+            }
+            for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
             {
                 std::cout << tLast << " " << corr[tLast] << std::endl;
             }

From 8b007b5c2481882a56bbdbb07e045cc6b5050a62 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Fri, 16 Nov 2018 20:00:29 +0000
Subject: [PATCH 23/33] Hadrons: remove the use of OpenMP reductions

---
 Grid/threads/Threads.h                   |   7 --
 Hadrons/A2AMatrix.hpp                    |  39 +++++--
 Hadrons/Utilities/ContractorBenchmark.cc | 126 ++++++++++++++++++-----
 3 files changed, 128 insertions(+), 44 deletions(-)

diff --git a/Grid/threads/Threads.h b/Grid/threads/Threads.h
index 580278d6..1c13fccb 100644
--- a/Grid/threads/Threads.h
+++ b/Grid/threads/Threads.h
@@ -38,14 +38,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #ifdef GRID_OMP
 #include <omp.h>
 
-// complex reductions
-#pragma omp declare reduction(ComplexPlus:Grid::ComplexD, Grid::vComplexD, Grid::ComplexF, Grid::vComplexF: omp_out += omp_in)
-
 #define PARALLEL_FOR_LOOP        _Pragma("omp parallel for schedule(static)")
 #define PARALLEL_FOR_LOOP_INTERN _Pragma("omp for schedule(static)")
-#define PARALLEL_FOR_REDUCE_HELPER0(x) #x
-#define PARALLEL_FOR_REDUCE_HELPER1(op, var) PARALLEL_FOR_REDUCE_HELPER0(omp parallel for schedule(static) reduction(op:var))
-#define PARALLEL_FOR_LOOP_REDUCE(op, var) _Pragma(PARALLEL_FOR_REDUCE_HELPER1(op, var))
 #define PARALLEL_NESTED_LOOP2 _Pragma("omp parallel for collapse(2)")
 #define PARALLEL_NESTED_LOOP5 _Pragma("omp parallel for collapse(5)")
 #define PARALLEL_REGION       _Pragma("omp parallel")
@@ -63,7 +57,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #define parallel_region    PARALLEL_REGION
 #define parallel_for       PARALLEL_FOR_LOOP for
 #define parallel_for_internal PARALLEL_FOR_LOOP_INTERN for
-#define parallel_for_reduce(op, var) PARALLEL_FOR_LOOP_REDUCE(op, var) for
 #define parallel_for_nest2 PARALLEL_NESTED_LOOP2 for
 #define parallel_for_nest5 PARALLEL_NESTED_LOOP5 for
 
diff --git a/Hadrons/A2AMatrix.hpp b/Hadrons/A2AMatrix.hpp
index 831d2a02..288bb7c5 100644
--- a/Hadrons/A2AMatrix.hpp
+++ b/Hadrons/A2AMatrix.hpp
@@ -167,33 +167,52 @@ public:
     template <typename C, typename MatLeft, typename MatRight>
     static inline void accTrMul(C &acc, const MatLeft &a, const MatRight &b)
     {
+        int            nThreads = GridThread::GetThreads();
+        std::vector<C> tacc(nThreads, 0.);
+
         if ((MatLeft::Options == Eigen::RowMajor) and
             (MatRight::Options == Eigen::ColMajor))
         {
-            parallel_for_reduce(ComplexPlus, acc) (unsigned int r = 0; r < a.rows(); ++r)
+            parallel_for (int thr = 0; thr < nThreads; ++thr)
             {
+                int rt, nr;
+
+                GridThread::GetWork(a.rows(), thr, nr, rt);
+                for (unsigned int r = rt; r < nr + rt; ++r)
+                {
 #ifdef USE_MKL
-                ComplexD tmp;
-                dotuRow(tmp, r, a, b);
-                acc += tmp;
+                    C tmp;
+                    dotuRow(tmp, r, a, b);
+                    tacc[thr] += tmp;
 #else
-                acc += a.row(r).conjugate().dot(b.col(r));
+                    tacc[thr] += a.row(r).conjugate().dot(b.col(r));
 #endif
+                }
             }
         }
         else
         {
-            parallel_for_reduce(ComplexPlus, acc) (unsigned int c = 0; c < a.cols(); ++c)
+            parallel_for (int thr = 0; thr < nThreads; ++thr)
             {
+                int ct, nc;
+
+                GridThread::GetWork(a.cols(), thr, nc, ct);
+                for (unsigned int c = ct; c < nc + ct; ++c)
+                {
 #ifdef USE_MKL
-                ComplexD tmp;
-                dotuCol(tmp, c, a, b);
-                acc += tmp;
+                    C tmp;
+                    dotuCol(tmp, c, a, b);
+                    tacc[thr] += tmp;
 #else
-                acc += a.col(c).conjugate().dot(b.row(c));
+                    tacc[thr] += a.col(c).conjugate().dot(b.row(c));
 #endif
+                }
             }
         }
+        for (int thr = 0; thr < nThreads; ++thr)
+        {
+            acc += tacc[thr];
+        }
     }
 
     template <typename MatLeft, typename MatRight>
diff --git a/Hadrons/Utilities/ContractorBenchmark.cc b/Hadrons/Utilities/ContractorBenchmark.cc
index 83d0e49f..a33f4ab4 100644
--- a/Hadrons/Utilities/ContractorBenchmark.cc
+++ b/Hadrons/Utilities/ContractorBenchmark.cc
@@ -199,23 +199,49 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
     trBenchmark("Naive loop rows first", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
     { 
+        int                   nThreads = GridThread::GetThreads();
+        std::vector<ComplexD> tres(nThreads, 0.);
+
         res = 0.;
-        auto nr = a.rows(), nc = a.cols();
-        parallel_for_reduce(ComplexPlus, res) (unsigned int i = 0; i < nr; ++i)
-        for (unsigned int j = 0; j < nc; ++j)
+        parallel_for (int thr = 0; thr < nThreads; ++thr)
         {
-            res += a(i, j)*b(j, i);
+            int  rt, nr;
+            auto nc = a.cols();
+
+            GridThread::GetWork(a.rows(), thr, nr, rt);
+            for (unsigned int i = rt; i < nr + rt; ++i)
+            for (unsigned int j = 0; j < nc; ++j)
+            {
+                tres[thr] += a(i, j)*b(j, i);
+            }
+        }
+        for (int thr = 0; thr < nThreads; ++thr)
+        {
+            res += tres[thr];
         }
     });
     trBenchmark("Naive loop cols first", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
-    { 
+    {
+        int                   nThreads = GridThread::GetThreads();
+        std::vector<ComplexD> tres(nThreads, 0.);
+
         res = 0.;
-        auto nr = a.rows(), nc = a.cols();
-        parallel_for_reduce(ComplexPlus, res) (unsigned int j = 0; j < nc; ++j)
-        for (unsigned int i = 0; i < nr; ++i)
+        parallel_for (int thr = 0; thr < nThreads; ++thr)
         {
-            res += a(i, j)*b(j, i);
+            int  ct, nc;
+            auto nr = a.rows();
+
+            GridThread::GetWork(a.cols(), thr, nc, ct);
+            for (unsigned int j = ct; j < nc + ct; ++j)
+            for (unsigned int i = 0; i < nr; ++i)
+            {
+                tres[thr] += a(i, j)*b(j, i);
+            }
+        }
+        for (int thr = 0; thr < nThreads; ++thr)
+        {
+            res += tres[thr];
         }
     });
     trBenchmark("Eigen tr(A*B)", left, right, ref,
@@ -225,22 +251,46 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
     });
     trBenchmark("Eigen row-wise dot", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
-    { 
-        res = 0.;
+    {
+        int                   nThreads = GridThread::GetThreads();
+        std::vector<ComplexD> tres(nThreads, 0.);
 
-        parallel_for_reduce(ComplexPlus, res) (unsigned int r = 0; r < a.rows(); ++r)
+        res = 0.;
+        parallel_for (int thr = 0; thr < nThreads; ++thr)
         {
-            res += a.row(r).conjugate().dot(b.col(r));
+            int  rt, nr;
+
+            GridThread::GetWork(a.rows(), thr, nr, rt);
+            for (unsigned int i = rt; i < nr + rt; ++i)
+            {
+                tres[thr] += a.row(i).conjugate().dot(b.col(i));
+            }
+        }
+        for (int thr = 0; thr < nThreads; ++thr)
+        {
+            res += tres[thr];
         }
     });
     trBenchmark("Eigen col-wise dot", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
-    { 
-        res = 0.;
+    {
+        int                   nThreads = GridThread::GetThreads();
+        std::vector<ComplexD> tres(nThreads, 0.);
 
-        parallel_for_reduce(ComplexPlus, res) (unsigned int c = 0; c < a.cols(); ++c)
+        res = 0.;
+        parallel_for (int thr = 0; thr < nThreads; ++thr)
         {
-            res += a.col(c).conjugate().dot(b.row(c));
+            int  ct, nc;
+
+            GridThread::GetWork(a.cols(), thr, nc, ct);
+            for (unsigned int j = ct; j < nc + ct; ++j)
+            {
+                tres[thr] += a.col(j).conjugate().dot(b.row(j));
+            }
+        }
+        for (int thr = 0; thr < nThreads; ++thr)
+        {
+            res += tres[thr];
         }
     });
     trBenchmark("Eigen Hadamard", left, right, ref,
@@ -251,28 +301,50 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
 #ifdef USE_MKL
     trBenchmark("MKL row-wise zdotu", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
-    { 
-        res = 0.;
+    {
+        int                   nThreads = GridThread::GetThreads();
+        std::vector<ComplexD> tres(nThreads, 0.);
 
-        parallel_for_reduce(ComplexPlus, res) (unsigned int r = 0; r < a.rows(); ++r)
+        res = 0.;
+        parallel_for (int thr = 0; thr < nThreads; ++thr)
         {
             ComplexD tmp;
+            int      rt, nr;
 
-            zdotuRow(tmp, r, a, b);
-            res += tmp;
+            GridThread::GetWork(a.rows(), thr, nr, rt);
+            for (unsigned int i = rt; i < nr + rt; ++i)
+            {
+                zdotuRow(tmp, i, a, b);
+                tres[thr] += tmp;
+            }
+        }
+        for (int thr = 0; thr < nThreads; ++thr)
+        {
+            res += tres[thr];
         }
     });
     trBenchmark("MKL col-wise zdotu", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
-    { 
-        res = 0.;
+    {
+        int                   nThreads = GridThread::GetThreads();
+        std::vector<ComplexD> tres(nThreads, 0.);
 
-        parallel_for_reduce(ComplexPlus, res) (unsigned int c = 0; c < a.cols(); ++c)
+        res = 0.;
+        parallel_for (int thr = 0; thr < nThreads; ++thr)
         {
             ComplexD tmp;
+            int      ct, nc;
 
-            zdotuCol(tmp, c, a, b);
-            res += tmp;
+            GridThread::GetWork(a.cols(), thr, nc, ct);
+            for (unsigned int j = ct; j < nc + ct; ++j)
+            {
+                zdotuCol(tmp, j, a, b);
+                tres[thr] += tmp;
+            }
+        }
+        for (int thr = 0; thr < nThreads; ++thr)
+        {
+            res += tres[thr];
         }
     });
 #endif

From f592ec8baac0878cf12874b445e5322bdbf7f467 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Fri, 16 Nov 2018 20:59:49 +0000
Subject: [PATCH 24/33] Hadrons: contractor performance fix

---
 Grid/threads/Threads.h                   |   1 +
 Hadrons/A2AMatrix.hpp                    |  43 +++------
 Hadrons/Utilities/ContractorBenchmark.cc | 112 ++++++++---------------
 3 files changed, 52 insertions(+), 104 deletions(-)

diff --git a/Grid/threads/Threads.h b/Grid/threads/Threads.h
index 1c13fccb..9041d01e 100644
--- a/Grid/threads/Threads.h
+++ b/Grid/threads/Threads.h
@@ -59,6 +59,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #define parallel_for_internal PARALLEL_FOR_LOOP_INTERN for
 #define parallel_for_nest2 PARALLEL_NESTED_LOOP2 for
 #define parallel_for_nest5 PARALLEL_NESTED_LOOP5 for
+#define parallel_critical PARALLEL_CRITICAL
 
 namespace Grid {
 
diff --git a/Hadrons/A2AMatrix.hpp b/Hadrons/A2AMatrix.hpp
index 288bb7c5..7f8791ea 100644
--- a/Hadrons/A2AMatrix.hpp
+++ b/Hadrons/A2AMatrix.hpp
@@ -167,52 +167,39 @@ public:
     template <typename C, typename MatLeft, typename MatRight>
     static inline void accTrMul(C &acc, const MatLeft &a, const MatRight &b)
     {
-        int            nThreads = GridThread::GetThreads();
-        std::vector<C> tacc(nThreads, 0.);
-
         if ((MatLeft::Options == Eigen::RowMajor) and
             (MatRight::Options == Eigen::ColMajor))
         {
-            parallel_for (int thr = 0; thr < nThreads; ++thr)
+            parallel_for (unsigned int r = 0; r < a.rows(); ++r)
             {
-                int rt, nr;
-
-                GridThread::GetWork(a.rows(), thr, nr, rt);
-                for (unsigned int r = rt; r < nr + rt; ++r)
-                {
+                C tmp;
 #ifdef USE_MKL
-                    C tmp;
-                    dotuRow(tmp, r, a, b);
-                    tacc[thr] += tmp;
+                dotuRow(tmp, r, a, b);
 #else
-                    tacc[thr] += a.row(r).conjugate().dot(b.col(r));
+                tmp = a.row(r).conjugate().dot(b.col(r));
 #endif
+                parallel_critical
+                {
+                    acc += tmp;
                 }
             }
         }
         else
         {
-            parallel_for (int thr = 0; thr < nThreads; ++thr)
+            parallel_for (unsigned int c = 0; c < a.cols(); ++c)
             {
-                int ct, nc;
-
-                GridThread::GetWork(a.cols(), thr, nc, ct);
-                for (unsigned int c = ct; c < nc + ct; ++c)
-                {
-#ifdef USE_MKL
-                    C tmp;
-                    dotuCol(tmp, c, a, b);
-                    tacc[thr] += tmp;
+                C tmp;
+#ifdef USE_MKL 
+                dotuCol(tmp, c, a, b);
 #else
-                    tacc[thr] += a.col(c).conjugate().dot(b.row(c));
+                tmp = a.col(c).conjugate().dot(b.row(c));
 #endif
+                parallel_critical
+                {
+                    acc += tmp;
                 }
             }
         }
-        for (int thr = 0; thr < nThreads; ++thr)
-        {
-            acc += tacc[thr];
-        }
     }
 
     template <typename MatLeft, typename MatRight>
diff --git a/Hadrons/Utilities/ContractorBenchmark.cc b/Hadrons/Utilities/ContractorBenchmark.cc
index a33f4ab4..1ec4657b 100644
--- a/Hadrons/Utilities/ContractorBenchmark.cc
+++ b/Hadrons/Utilities/ContractorBenchmark.cc
@@ -199,50 +199,42 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
     trBenchmark("Naive loop rows first", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
     { 
-        int                   nThreads = GridThread::GetThreads();
-        std::vector<ComplexD> tres(nThreads, 0.);
-
+        auto nr = a.rows(), nc = a.cols();
+        
         res = 0.;
-        parallel_for (int thr = 0; thr < nThreads; ++thr)
+        parallel_for (unsigned int i = 0; i < nr; ++i)
         {
-            int  rt, nr;
-            auto nc = a.cols();
+            ComplexD tmp = 0.;
 
-            GridThread::GetWork(a.rows(), thr, nr, rt);
-            for (unsigned int i = rt; i < nr + rt; ++i)
             for (unsigned int j = 0; j < nc; ++j)
             {
-                tres[thr] += a(i, j)*b(j, i);
+                tmp += a(i, j)*b(j, i);
+            }
+            parallel_critical
+            {
+                res += tmp;
             }
-        }
-        for (int thr = 0; thr < nThreads; ++thr)
-        {
-            res += tres[thr];
         }
     });
     trBenchmark("Naive loop cols first", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
     {
-        int                   nThreads = GridThread::GetThreads();
-        std::vector<ComplexD> tres(nThreads, 0.);
-
+        auto nr = a.rows(), nc = a.cols();
+        
         res = 0.;
-        parallel_for (int thr = 0; thr < nThreads; ++thr)
+        parallel_for (unsigned int j = 0; j < nc; ++j)
         {
-            int  ct, nc;
-            auto nr = a.rows();
+            ComplexD tmp = 0.;
 
-            GridThread::GetWork(a.cols(), thr, nc, ct);
-            for (unsigned int j = ct; j < nc + ct; ++j)
             for (unsigned int i = 0; i < nr; ++i)
             {
-                tres[thr] += a(i, j)*b(j, i);
+                tmp += a(i, j)*b(j, i);
+            }        
+            parallel_critical
+            {
+                res += tmp;
             }
         }
-        for (int thr = 0; thr < nThreads; ++thr)
-        {
-            res += tres[thr];
-        }
     });
     trBenchmark("Eigen tr(A*B)", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
@@ -252,46 +244,32 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
     trBenchmark("Eigen row-wise dot", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
     {
-        int                   nThreads = GridThread::GetThreads();
-        std::vector<ComplexD> tres(nThreads, 0.);
-
         res = 0.;
-        parallel_for (int thr = 0; thr < nThreads; ++thr)
+        parallel_for (unsigned int r = 0; r < a.rows(); ++r)
         {
-            int  rt, nr;
+            ComplexD tmp;
 
-            GridThread::GetWork(a.rows(), thr, nr, rt);
-            for (unsigned int i = rt; i < nr + rt; ++i)
+            tmp = a.row(r).conjugate().dot(b.col(r));
+            parallel_critical
             {
-                tres[thr] += a.row(i).conjugate().dot(b.col(i));
+                res += tmp;
             }
         }
-        for (int thr = 0; thr < nThreads; ++thr)
-        {
-            res += tres[thr];
-        }
     });
     trBenchmark("Eigen col-wise dot", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
     {
-        int                   nThreads = GridThread::GetThreads();
-        std::vector<ComplexD> tres(nThreads, 0.);
-
         res = 0.;
-        parallel_for (int thr = 0; thr < nThreads; ++thr)
+        parallel_for (unsigned int c = 0; c < a.cols(); ++c)
         {
-            int  ct, nc;
+            ComplexD tmp;
 
-            GridThread::GetWork(a.cols(), thr, nc, ct);
-            for (unsigned int j = ct; j < nc + ct; ++j)
+            tmp = a.col(c).conjugate().dot(b.row(c));
+            parallel_critical
             {
-                tres[thr] += a.col(j).conjugate().dot(b.row(j));
+                res += tmp;
             }
         }
-        for (int thr = 0; thr < nThreads; ++thr)
-        {
-            res += tres[thr];
-        }
     });
     trBenchmark("Eigen Hadamard", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
@@ -302,50 +280,32 @@ void fullTrBenchmark(const unsigned int ni, const unsigned int nj, const unsigne
     trBenchmark("MKL row-wise zdotu", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
     {
-        int                   nThreads = GridThread::GetThreads();
-        std::vector<ComplexD> tres(nThreads, 0.);
-
         res = 0.;
-        parallel_for (int thr = 0; thr < nThreads; ++thr)
+        parallel_for (unsigned int r = 0; r < a.rows(); ++r)
         {
             ComplexD tmp;
-            int      rt, nr;
 
-            GridThread::GetWork(a.rows(), thr, nr, rt);
-            for (unsigned int i = rt; i < nr + rt; ++i)
+            zdotuRow(tmp, r, a, b);
+            parallel_critical
             {
-                zdotuRow(tmp, i, a, b);
-                tres[thr] += tmp;
+                res += tmp;
             }
         }
-        for (int thr = 0; thr < nThreads; ++thr)
-        {
-            res += tres[thr];
-        }
     });
     trBenchmark("MKL col-wise zdotu", left, right, ref,
     [](ComplexD &res, const MatLeft &a, const MatRight &b)
     {
-        int                   nThreads = GridThread::GetThreads();
-        std::vector<ComplexD> tres(nThreads, 0.);
-
         res = 0.;
-        parallel_for (int thr = 0; thr < nThreads; ++thr)
+        parallel_for (unsigned int c = 0; c < a.cols(); ++c)
         {
             ComplexD tmp;
-            int      ct, nc;
 
-            GridThread::GetWork(a.cols(), thr, nc, ct);
-            for (unsigned int j = ct; j < nc + ct; ++j)
+            zdotuCol(tmp, c, a, b);
+            parallel_critical
             {
-                zdotuCol(tmp, j, a, b);
-                tres[thr] += tmp;
+                res += tmp;
             }
         }
-        for (int thr = 0; thr < nThreads; ++thr)
-        {
-            res += tres[thr];
-        }
     });
 #endif
     BARRIER();

From 33a0bbb17bcab446efbb94edbf0406d1bc819ed5 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Mon, 19 Nov 2018 11:27:57 +0000
Subject: [PATCH 25/33] Const correctness

---
 Grid/lattice/Lattice_transfer.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Grid/lattice/Lattice_transfer.h b/Grid/lattice/Lattice_transfer.h
index 4eca50d6..69c132ed 100644
--- a/Grid/lattice/Lattice_transfer.h
+++ b/Grid/lattice/Lattice_transfer.h
@@ -487,7 +487,7 @@ void InsertSliceLocal(const Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int
 
 
 template<class vobj>
-void ExtractSliceLocal(Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slice_lo,int slice_hi, int orthog)
+void ExtractSliceLocal(Lattice<vobj> &lowDim,const Lattice<vobj> & higherDim,int slice_lo,int slice_hi, int orthog)
 {
   typedef typename vobj::scalar_object sobj;
 

From f22a27d7f92f6d50c9b43bc74f8c54dc1bad0fef Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Mon, 19 Nov 2018 15:45:04 +0000
Subject: [PATCH 26/33] Hadrons: contractor trajectory loop and file output

---
 Hadrons/Global.cc                         |   8 +-
 Hadrons/Global.hpp                        |   6 +-
 Hadrons/Module.hpp                        |   2 +-
 Hadrons/Modules/MScalar/ChargedProp.cc    |   2 +-
 Hadrons/Modules/MScalar/ScalarVP.cc       |   2 +-
 Hadrons/Modules/MScalar/VPCounterTerms.cc |   2 +-
 Hadrons/Utilities/Contractor.cc           | 326 +++++++++++++---------
 7 files changed, 210 insertions(+), 138 deletions(-)

diff --git a/Hadrons/Global.cc b/Hadrons/Global.cc
index a83f422f..d72dbbd6 100644
--- a/Hadrons/Global.cc
+++ b/Hadrons/Global.cc
@@ -166,7 +166,13 @@ std::string Hadrons::dirname(const std::string &s)
 
 void Hadrons::makeFileDir(const std::string filename, GridBase *g)
 {
-    if (g->IsBoss())
+    bool doIt = true;
+
+    if (g)
+    {
+        doIt = g->IsBoss();
+    }
+    if (doIt)
     {
         std::string dir    = dirname(filename);
         int         status = mkdir(dir);
diff --git a/Hadrons/Global.hpp b/Hadrons/Global.hpp
index d8047f09..8a769398 100644
--- a/Hadrons/Global.hpp
+++ b/Hadrons/Global.hpp
@@ -218,15 +218,15 @@ typedef XmlReader ResultReader;
 typedef XmlWriter ResultWriter;
 #endif
 
-#define RESULT_FILE_NAME(name) \
-name + "." + std::to_string(vm().getTrajectory()) + "." + resultFileExt
+#define RESULT_FILE_NAME(name, traj) \
+name + "." + std::to_string(traj) + "." + resultFileExt
 
 // recursive mkdir
 #define MAX_PATH_LENGTH 512u
 int         mkdir(const std::string dirName);
 std::string basename(const std::string &s);
 std::string dirname(const std::string &s);
-void        makeFileDir(const std::string filename, GridBase *g);
+void        makeFileDir(const std::string filename, GridBase *g = nullptr);
 
 // default Schur convention
 #ifndef HADRONS_DEFAULT_SCHUR 
diff --git a/Hadrons/Module.hpp b/Hadrons/Module.hpp
index 752b821e..8d59d06c 100644
--- a/Hadrons/Module.hpp
+++ b/Hadrons/Module.hpp
@@ -144,7 +144,7 @@ if (env().getGrid()->IsBoss() and !ioStem.empty())\
 {\
     makeFileDir(ioStem, env().getGrid());\
     {\
-        ResultWriter _writer(RESULT_FILE_NAME(ioStem));\
+        ResultWriter _writer(RESULT_FILE_NAME(ioStem, vm().getTrajectory()));\
         write(_writer, name, result);\
     }\
 }
diff --git a/Hadrons/Modules/MScalar/ChargedProp.cc b/Hadrons/Modules/MScalar/ChargedProp.cc
index ae5b849f..51cfd585 100644
--- a/Hadrons/Modules/MScalar/ChargedProp.cc
+++ b/Hadrons/Modules/MScalar/ChargedProp.cc
@@ -146,7 +146,7 @@ void TChargedProp::execute(void)
         std::vector<int>    siteCoor;
 
         LOG(Message) << "Saving momentum-projected propagator to '"
-                     << RESULT_FILE_NAME(par().output) << "'..."
+                     << RESULT_FILE_NAME(par().output, vm().getTrajectory()) << "'..."
                      << std::endl;
         result.projection.resize(par().outputMom.size());
         result.lattice_size = env().getGrid()->_fdimensions;
diff --git a/Hadrons/Modules/MScalar/ScalarVP.cc b/Hadrons/Modules/MScalar/ScalarVP.cc
index 8204b7aa..b62ec380 100644
--- a/Hadrons/Modules/MScalar/ScalarVP.cc
+++ b/Hadrons/Modules/MScalar/ScalarVP.cc
@@ -462,7 +462,7 @@ void TScalarVP::execute(void)
     if (!par().output.empty())
     {
         LOG(Message) << "Saving momentum-projected HVP to '"
-                     << RESULT_FILE_NAME(par().output) << "'..."
+                     << RESULT_FILE_NAME(par().output, vm().getTrajectory()) << "'..."
                      << std::endl;
         saveResult(par().output, "HVP", outputData);
     }
diff --git a/Hadrons/Modules/MScalar/VPCounterTerms.cc b/Hadrons/Modules/MScalar/VPCounterTerms.cc
index a5aff51f..974d2171 100644
--- a/Hadrons/Modules/MScalar/VPCounterTerms.cc
+++ b/Hadrons/Modules/MScalar/VPCounterTerms.cc
@@ -239,7 +239,7 @@ void TVPCounterTerms::execute(void)
     if (!par().output.empty())
     {
         LOG(Message) << "Saving momentum-projected correlators to '"
-                     << RESULT_FILE_NAME(par().output) << "'..."
+                     << RESULT_FILE_NAME(par().output, vm().getTrajectory()) << "'..."
                      << std::endl;
         saveResult(par().output, "scalar_loops", outputData);
     }
diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index ffc77cf7..4bfa1e28 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -37,10 +37,20 @@ using namespace Hadrons;
 
 namespace Contractor
 {
+    class TrajRange: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(TrajRange,
+                                        unsigned int, start,
+                                        unsigned int, end,
+                                        unsigned int, step);
+    };
+    
     class GlobalPar: Serializable
     {
     public:
         GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar,
+                                        TrajRange, trajCounter,
                                         unsigned int, nt,
                                         std::string, diskVectorDir,
                                         std::string, output);
@@ -50,7 +60,7 @@ namespace Contractor
     {
     public:
         GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMatrixPar,
-                                        std::string, file,
+                                        std::string, fileStem,
                                         std::string, dataset,
                                         unsigned int, cacheSize,
                                         std::string, name);
@@ -62,7 +72,18 @@ namespace Contractor
         GRID_SERIALIZABLE_CLASS_MEMBERS(ProductPar,
                                         std::string, terms,
                                         std::vector<std::string>, times,
-                                        std::string, translations);
+                                        std::string, translations,
+                                        bool, translationAverage);
+    };
+
+    class CorrelatorResult: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(CorrelatorResult,
+                                        std::vector<Contractor::A2AMatrixPar>,  a2aMatrix,
+                                        ProductPar, contraction,
+                                        std::vector<unsigned int>, times,
+                                        std::vector<ComplexD>, correlator);
     };
 }
 
@@ -100,6 +121,26 @@ void makeTimeSeq(std::vector<std::vector<unsigned int>> &timeSeq,
     makeTimeSeq(timeSeq, times, current, times.size());
 }
 
+void saveCorrelator(const Contractor::CorrelatorResult &result, const std::string dir, 
+                    const unsigned int dt, const unsigned int traj)
+{
+    std::string fileStem = "", us = "_", filename;
+
+    for (unsigned int i = 0; i < result.contraction.terms.size() - 1; i++)
+    {
+        fileStem += result.contraction.terms[i] + us + std::to_string(result.times[i]) + us;
+    }
+    fileStem += result.contraction.terms.back();
+    if (!result.contraction.translationAverage)
+    {
+        fileStem += "_dt_" + std::to_string(dt);
+    }
+    filename = dir + "/" + RESULT_FILE_NAME(fileStem, traj);
+    std::cout << "Saving correlator to '" << filename << "'" << std::endl;
+    makeFileDir(dir);
+    ResultWriter writer(filename);
+    write(writer, fileStem, result);
+}
 
 std::set<unsigned int> parseTimeRange(const std::string str, const unsigned int nt)
 {
@@ -107,7 +148,6 @@ std::set<unsigned int> parseTimeRange(const std::string str, const unsigned int
     std::smatch              sm;
     std::vector<std::string> rstr = strToVec<std::string>(str);
     std::set<unsigned int>   tSet;
-    
 
     for (auto &s: rstr)
     {
@@ -230,158 +270,184 @@ int main(int argc, char* argv[])
         a2aMat.emplace(p.name, EigenDiskVector<ComplexD>(dirName, par.global.nt, p.cacheSize));
     }
 
-    // load data
-    for (unsigned int i = 0; i < a2aMat.size(); ++i)
+    // trajectory loop
+    for (unsigned int traj = par.global.trajCounter.start; 
+         traj < par.global.trajCounter.end; traj += par.global.trajCounter.step)
     {
-        auto   &p = par.a2aMatrix[i];
-        double t, size;
+        std::cout << ":::::::: Trajectory " << traj << std::endl;
 
-        std::cout << "======== Loading '" << p.file << "'" << std::endl;
-
-        A2AMatrixIo<HADRONS_A2AM_IO_TYPE> a2aIo(p.file, p.dataset, par.global.nt);
-
-        a2aIo.load(a2aMat.at(p.name), &t);
-        std::cout << "Read " << a2aIo.getSize() << " bytes in " << t/1.0e6 
-                  << " sec, " << a2aIo.getSize()/t*1.0e6/1024/1024 << " MB/s" << std::endl;
-    }
-
-    // contract
-    EigenDiskVector<ComplexD>::Matrix buf;
-
-    for (auto &p: par.product)
-    {
-        std::vector<std::string>               term = strToVec<std::string>(p.terms);
-        std::vector<std::set<unsigned int>>    times;
-        std::vector<std::vector<unsigned int>> timeSeq;
-        std::set<unsigned int>                 translations;
-        std::vector<ComplexD>                  corr(par.global.nt);
-        std::vector<A2AMatrixTr<ComplexD>>     lastTerm(par.global.nt);
-        A2AMatrix<ComplexD>                    prod, buf, tmp;
-        TimerArray                             tAr;
-        double                                 fusec, busec, flops, bytes, tusec;
-
-        tAr.startTimer("Total");
-        std::cout << "======== Contraction tr(";
-        for (unsigned int g = 0; g < term.size(); ++g)
+        // load data
+        for (auto &p: par.a2aMatrix)
         {
-            std::cout << term[g] << ((g == term.size() - 1) ? ')' : '*');
-        }
-        std::cout << std::endl;
-        if (term.size() != p.times.size() + 1)
-        {
-            HADRONS_ERROR(Size, "number of terms (" + std::to_string(term.size()) 
-                          + ") different from number of times (" 
-                          + std::to_string(p.times.size() + 1) + ")");
-        }
-        for (auto &s: p.times)
-        {
-            times.push_back(parseTimeRange(s, par.global.nt));
+            std::string filename = RESULT_FILE_NAME(p.fileStem, traj);
+            double      t, size;
+
+            std::cout << "======== Loading '" << filename << "'" << std::endl;
+
+            A2AMatrixIo<HADRONS_A2AM_IO_TYPE> a2aIo(filename, p.dataset, par.global.nt);
+
+            a2aIo.load(a2aMat.at(p.name), &t);
+            std::cout << "Read " << a2aIo.getSize() << " bytes in " << t/1.0e6 
+                    << " sec, " << a2aIo.getSize()/t*1.0e6/1024/1024 << " MB/s" << std::endl;
         }
 
-        translations = parseTimeRange(p.translations, par.global.nt);
-        makeTimeSeq(timeSeq, times);
-        std::cout << timeSeq.size()*translations.size()*(term.size() - 2) << " A*B, "
-                  << timeSeq.size()*translations.size()*par.global.nt << " tr(A*B)"
-                  << std::endl;
+        // contract
+        EigenDiskVector<ComplexD>::Matrix buf;
 
-        std::cout << "* Caching transposed last term" << std::endl;
-        for (unsigned int t = 0; t < par.global.nt; ++t)
+        for (auto &p: par.product)
         {
-            tAr.startTimer("Disk vector overhead");
-            const A2AMatrix<ComplexD> &ref = a2aMat.at(term.back())[t];
-            tAr.stopTimer("Disk vector overhead");
+            std::vector<std::string>               term = strToVec<std::string>(p.terms);
+            std::vector<std::set<unsigned int>>    times;
+            std::vector<std::vector<unsigned int>> timeSeq;
+            std::set<unsigned int>                 translations;
+            std::vector<A2AMatrixTr<ComplexD>>     lastTerm(par.global.nt);
+            A2AMatrix<ComplexD>                    prod, buf, tmp;
+            TimerArray                             tAr;
+            double                                 fusec, busec, flops, bytes, tusec;
+            Contractor::CorrelatorResult           result;             
 
-            tAr.startTimer("Transpose caching");
-            lastTerm[t].resize(ref.rows(), ref.cols());
-            parallel_for (unsigned int j = 0; j < ref.cols(); ++j)
-            for (unsigned int i = 0; i < ref.rows(); ++i)
+            tAr.startTimer("Total");
+            std::cout << "======== Contraction tr(";
+            for (unsigned int g = 0; g < term.size(); ++g)
             {
-                lastTerm[t](i, j) = ref(i, j);
+                std::cout << term[g] << ((g == term.size() - 1) ? ')' : '*');
             }
-            tAr.stopTimer("Transpose caching");
-        }
-        bytes = par.global.nt*lastTerm[0].rows()*lastTerm[0].cols()*sizeof(ComplexD);
-        std::cout << Sec(tAr.getDTimer("Transpose caching")) << " " 
-                  << Bytes(bytes, tAr.getDTimer("Transpose caching")) << std::endl;
-        std::cout << Sec(tAr.getDTimer("Disk vector overhead")) << std::endl;
-        for (unsigned int i = 0; i < timeSeq.size(); ++i)
-        {
-            unsigned int dti = 0;
-            auto         &t = timeSeq[i];
-
-            for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
+            std::cout << std::endl;
+            if (term.size() != p.times.size() + 1)
             {
-                corr[tLast] = 0.;
+                HADRONS_ERROR(Size, "number of terms (" + std::to_string(term.size()) 
+                            + ") different from number of times (" 
+                            + std::to_string(p.times.size() + 1) + ")");
             }
-            for (auto &dt: translations)
+            for (auto &s: p.times)
             {
-                std::cout << "* Step " << i*translations.size() + dti + 1
-                          << "/" << timeSeq.size()*translations.size()
-                          << " -- positions= " << t << ", dt= " << dt << std::endl;
-                if (term.size() > 2)
+                times.push_back(parseTimeRange(s, par.global.nt));
+            }
+            for (auto &m: par.a2aMatrix)
+            {
+                if (std::find(result.a2aMatrix.begin(), result.a2aMatrix.end(), m) == result.a2aMatrix.end())
                 {
-                    std::cout << std::setw(8) << "products";
+                    result.a2aMatrix.push_back(m);
                 }
-                flops  = 0.;
-                bytes  = 0.;
-                fusec  = tAr.getDTimer("A*B algebra");
-                busec  = tAr.getDTimer("A*B total");
-                tAr.startTimer("Linear algebra");
+            }
+            result.contraction = p;
+            result.correlator.resize(par.global.nt, 0.);
+
+            translations = parseTimeRange(p.translations, par.global.nt);
+            makeTimeSeq(timeSeq, times);
+            std::cout << timeSeq.size()*translations.size()*(term.size() - 2) << " A*B, "
+                    << timeSeq.size()*translations.size()*par.global.nt << " tr(A*B)"
+                    << std::endl;
+
+            std::cout << "* Caching transposed last term" << std::endl;
+            for (unsigned int t = 0; t < par.global.nt; ++t)
+            {
                 tAr.startTimer("Disk vector overhead");
-                prod = a2aMat.at(term[0])[TIME_MOD(t[0] + dt)];
+                const A2AMatrix<ComplexD> &ref = a2aMat.at(term.back())[t];
                 tAr.stopTimer("Disk vector overhead");
-                for (unsigned int j = 1; j < term.size() - 1; ++j)
+
+                tAr.startTimer("Transpose caching");
+                lastTerm[t].resize(ref.rows(), ref.cols());
+                parallel_for (unsigned int j = 0; j < ref.cols(); ++j)
+                for (unsigned int i = 0; i < ref.rows(); ++i)
                 {
-                    tAr.startTimer("Disk vector overhead");
-                    const A2AMatrix<ComplexD> &ref = a2aMat.at(term[j])[TIME_MOD(t[j] + dt)];
-                    tAr.stopTimer("Disk vector overhead");
-                    
-                    tAr.startTimer("A*B total");
-                    tAr.startTimer("A*B algebra");
-                    A2AContraction::mul(tmp, prod, ref);
-                    tAr.stopTimer("A*B algebra");
-                    flops += A2AContraction::mulFlops(prod, ref);
-                    prod   = tmp;
-                    tAr.stopTimer("A*B total");
-                    bytes += 3.*tmp.rows()*tmp.cols()*sizeof(ComplexD);
+                    lastTerm[t](i, j) = ref(i, j);
                 }
-                if (term.size() > 2)
-                {
-                    std::cout << Sec(tAr.getDTimer("A*B total") - busec) << " "
-                              << Flops(flops, tAr.getDTimer("A*B algebra") - fusec) << " " 
-                              << Bytes(bytes, tAr.getDTimer("A*B total") - busec) << std::endl;
-                }
-                std::cout << std::setw(8) << "traces";
-                flops  = 0.;
-                bytes  = 0.;
-                fusec  = tAr.getDTimer("tr(A*B)");
-                busec  = tAr.getDTimer("tr(A*B)");
+                tAr.stopTimer("Transpose caching");
+            }
+            bytes = par.global.nt*lastTerm[0].rows()*lastTerm[0].cols()*sizeof(ComplexD);
+            std::cout << Sec(tAr.getDTimer("Transpose caching")) << " " 
+                    << Bytes(bytes, tAr.getDTimer("Transpose caching")) << std::endl;
+            std::cout << Sec(tAr.getDTimer("Disk vector overhead")) << std::endl;
+            for (unsigned int i = 0; i < timeSeq.size(); ++i)
+            {
+                unsigned int dti = 0;
+                auto         &t = timeSeq[i];
+
+                result.times = t;
                 for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
                 {
-                    tAr.startTimer("tr(A*B)");
-                    A2AContraction::accTrMul(corr[TIME_MOD(tLast - dt)], prod, lastTerm[tLast]);
-                    tAr.stopTimer("tr(A*B)");
-                    flops += A2AContraction::accTrMulFlops(prod, lastTerm[tLast]);
-                    bytes += 2.*prod.rows()*prod.cols()*sizeof(ComplexD);
+                    result.correlator[tLast] = 0.;
+                }
+                for (auto &dt: translations)
+                {
+                    std::cout << "* Step " << i*translations.size() + dti + 1
+                            << "/" << timeSeq.size()*translations.size()
+                            << " -- positions= " << t << ", dt= " << dt << std::endl;
+                    if (term.size() > 2)
+                    {
+                        std::cout << std::setw(8) << "products";
+                    }
+                    flops  = 0.;
+                    bytes  = 0.;
+                    fusec  = tAr.getDTimer("A*B algebra");
+                    busec  = tAr.getDTimer("A*B total");
+                    tAr.startTimer("Linear algebra");
+                    tAr.startTimer("Disk vector overhead");
+                    prod = a2aMat.at(term[0])[TIME_MOD(t[0] + dt)];
+                    tAr.stopTimer("Disk vector overhead");
+                    for (unsigned int j = 1; j < term.size() - 1; ++j)
+                    {
+                        tAr.startTimer("Disk vector overhead");
+                        const A2AMatrix<ComplexD> &ref = a2aMat.at(term[j])[TIME_MOD(t[j] + dt)];
+                        tAr.stopTimer("Disk vector overhead");
+                        
+                        tAr.startTimer("A*B total");
+                        tAr.startTimer("A*B algebra");
+                        A2AContraction::mul(tmp, prod, ref);
+                        tAr.stopTimer("A*B algebra");
+                        flops += A2AContraction::mulFlops(prod, ref);
+                        prod   = tmp;
+                        tAr.stopTimer("A*B total");
+                        bytes += 3.*tmp.rows()*tmp.cols()*sizeof(ComplexD);
+                    }
+                    if (term.size() > 2)
+                    {
+                        std::cout << Sec(tAr.getDTimer("A*B total") - busec) << " "
+                                << Flops(flops, tAr.getDTimer("A*B algebra") - fusec) << " " 
+                                << Bytes(bytes, tAr.getDTimer("A*B total") - busec) << std::endl;
+                    }
+                    std::cout << std::setw(8) << "traces";
+                    flops  = 0.;
+                    bytes  = 0.;
+                    fusec  = tAr.getDTimer("tr(A*B)");
+                    busec  = tAr.getDTimer("tr(A*B)");
+                    for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
+                    {
+                        tAr.startTimer("tr(A*B)");
+                        A2AContraction::accTrMul(result.correlator[TIME_MOD(tLast - dt)], prod, lastTerm[tLast]);
+                        tAr.stopTimer("tr(A*B)");
+                        flops += A2AContraction::accTrMulFlops(prod, lastTerm[tLast]);
+                        bytes += 2.*prod.rows()*prod.cols()*sizeof(ComplexD);
+                    }
+                    tAr.stopTimer("Linear algebra");
+                    std::cout << Sec(tAr.getDTimer("tr(A*B)") - busec) << " "
+                            << Flops(flops, tAr.getDTimer("tr(A*B)") - fusec) << " " 
+                            << Bytes(bytes, tAr.getDTimer("tr(A*B)") - busec) << std::endl;
+                    if (!p.translationAverage)
+                    {
+                        saveCorrelator(result, par.global.output, dt, traj);
+                        for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
+                        {
+                            result.correlator[tLast] = 0.;
+                        }
+                    }
+                    dti++;
+                }
+                if (p.translationAverage)
+                {
+                    for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
+                    {
+                        result.correlator[tLast] /= translations.size();
+                    }
+                    saveCorrelator(result, par.global.output, 0, traj);
                 }
-                tAr.stopTimer("Linear algebra");
-                std::cout << Sec(tAr.getDTimer("tr(A*B)") - busec) << " "
-                          << Flops(flops, tAr.getDTimer("tr(A*B)") - fusec) << " " 
-                          << Bytes(bytes, tAr.getDTimer("tr(A*B)") - busec) << std::endl;
-                dti++;
-            }
-            for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
-            {
-                corr[tLast] /= translations.size();
-            }
-            for (unsigned int tLast = 0; tLast < par.global.nt; ++tLast)
-            {
-                std::cout << tLast << " " << corr[tLast] << std::endl;
             }
+            tAr.stopTimer("Total");
+            printTimeProfile(tAr.getTimings(), tAr.getTimer("Total"));
         }
-        tAr.stopTimer("Total");
-        printTimeProfile(tAr.getTimings(), tAr.getTimer("Total"));
     }
+
     
     return EXIT_SUCCESS;
 }

From 5d7cc29eaf5484de5931b0511d1e34c1a6a24ada Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Mon, 19 Nov 2018 16:04:01 +0000
Subject: [PATCH 27/33] Hadrons: contractor token @traj@ for trajectory number
 in input file

---
 Hadrons/Global.hpp              | 14 ++++++++++++++
 Hadrons/Utilities/Contractor.cc |  5 +++--
 2 files changed, 17 insertions(+), 2 deletions(-)

diff --git a/Hadrons/Global.hpp b/Hadrons/Global.hpp
index 8a769398..aaf03427 100644
--- a/Hadrons/Global.hpp
+++ b/Hadrons/Global.hpp
@@ -249,6 +249,20 @@ void        makeFileDir(const std::string filename, GridBase *g = nullptr);
 // pretty print time profile
 void printTimeProfile(const std::map<std::string, GridTime> &timing, GridTime total);
 
+// token replacement utility
+template <typename T>
+void tokenReplace(std::string &str, const std::string token,
+                  const T &x, const std::string mark = "@")
+{
+    std::string fullToken = mark + token + mark;
+    
+    auto pos = str.find(fullToken);
+    if (pos != std::string::npos)
+    {
+        str.replace(pos, fullToken.size(), std::to_string(x));
+    }
+}
+
 END_HADRONS_NAMESPACE
 
 #include <Hadrons/Exceptions.hpp>
diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index 4bfa1e28..2bed6bf5 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -60,7 +60,7 @@ namespace Contractor
     {
     public:
         GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMatrixPar,
-                                        std::string, fileStem,
+                                        std::string, file,
                                         std::string, dataset,
                                         unsigned int, cacheSize,
                                         std::string, name);
@@ -279,9 +279,10 @@ int main(int argc, char* argv[])
         // load data
         for (auto &p: par.a2aMatrix)
         {
-            std::string filename = RESULT_FILE_NAME(p.fileStem, traj);
+            std::string filename = p.file;
             double      t, size;
 
+            tokenReplace(filename, "traj", traj);
             std::cout << "======== Loading '" << filename << "'" << std::endl;
 
             A2AMatrixIo<HADRONS_A2AM_IO_TYPE> a2aIo(filename, p.dataset, par.global.nt);

From 2ba19a9e07fb948d473552a424276f4511e4e9b8 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Mon, 19 Nov 2018 16:17:51 +0000
Subject: [PATCH 28/33] Hadrons: contractor IO fix

---
 Hadrons/Utilities/Contractor.cc | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index 2bed6bf5..302eedc7 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -124,13 +124,14 @@ void makeTimeSeq(std::vector<std::vector<unsigned int>> &timeSeq,
 void saveCorrelator(const Contractor::CorrelatorResult &result, const std::string dir, 
                     const unsigned int dt, const unsigned int traj)
 {
-    std::string fileStem = "", us = "_", filename;
+    std::string              fileStem = "", filename;
+    std::vector<std::string> terms = strToVec<std::string>(result.contraction.terms);
 
-    for (unsigned int i = 0; i < result.contraction.terms.size() - 1; i++)
+    for (unsigned int i = 0; i < terms.size() - 1; i++)
     {
-        fileStem += result.contraction.terms[i] + us + std::to_string(result.times[i]) + us;
+        fileStem += terms[i] + "_" + std::to_string(result.times[i]) + "_";
     }
-    fileStem += result.contraction.terms.back();
+    fileStem += terms.back();
     if (!result.contraction.translationAverage)
     {
         fileStem += "_dt_" + std::to_string(dt);

From 494b3c9e572366c2e16d500959336a62e66f09c4 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Mon, 19 Nov 2018 16:26:53 +0000
Subject: [PATCH 29/33] Hadrons: contractor more IO fix

---
 Hadrons/Utilities/Contractor.cc | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index 302eedc7..9ce056f3 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -330,6 +330,7 @@ int main(int argc, char* argv[])
                 if (std::find(result.a2aMatrix.begin(), result.a2aMatrix.end(), m) == result.a2aMatrix.end())
                 {
                     result.a2aMatrix.push_back(m);
+                    tokenReplace(result.a2aMatrix.back().file, "traj", traj);
                 }
             }
             result.contraction = p;

From d77bc8817042d18642883285601561242de717f1 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Mon, 19 Nov 2018 17:21:53 +0000
Subject: [PATCH 30/33] Optional support for faster CRC32C checksum through
 Intel IPP

---
 Grid/util/Sha.h        | 16 +++++++++++++++-
 Hadrons/DiskVector.hpp |  8 ++++++++
 configure.ac           | 24 ++++++++++++++++++++++--
 3 files changed, 45 insertions(+), 3 deletions(-)

diff --git a/Grid/util/Sha.h b/Grid/util/Sha.h
index 15cee338..ee164c34 100644
--- a/Grid/util/Sha.h
+++ b/Grid/util/Sha.h
@@ -28,17 +28,31 @@
 extern "C" {
 #include <openssl/sha.h>
 }
+#ifdef USE_IPP
+#include "ipp.h"
+#endif
 
 #pragma once
 
 class GridChecksum
 {
 public:
-  static inline uint32_t crc32(const void *data,size_t bytes)
+  static inline uint32_t crc32(const void *data, size_t bytes)
   {
     return ::crc32(0L,(unsigned char *)data,bytes);
   }
 
+#ifdef USE_IPP
+  static inline uint32_t crc32c(const void* data, size_t bytes)
+  {
+      uint32_t crc32c = ~(uint32_t)0;
+      ippsCRC32C_8u(reinterpret_cast<const unsigned char *>(data), bytes, &crc32c);
+      ippsSwapBytes_32u_I(&crc32c, 1);
+  
+      return ~crc32c;
+  }
+#endif
+
   template <typename T>
   static inline std::string sha256_string(const std::vector<T> &hash)
   {
diff --git a/Hadrons/DiskVector.hpp b/Hadrons/DiskVector.hpp
index 6d96c6d9..70a00b16 100644
--- a/Hadrons/DiskVector.hpp
+++ b/Hadrons/DiskVector.hpp
@@ -176,7 +176,11 @@ private:
         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 
@@ -199,7 +203,11 @@ private:
         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));
diff --git a/configure.ac b/configure.ac
index d362ed5a..b02b514a 100644
--- a/configure.ac
+++ b/configure.ac
@@ -123,10 +123,13 @@ case ${ac_SFW_FP16} in
       AC_MSG_ERROR(["SFW FP16 option not supported ${ac_SFW_FP16}"]);;
 esac
 
-############### MKL
+############### Intel libraries
 AC_ARG_ENABLE([mkl],
     [AC_HELP_STRING([--enable-mkl=yes|no|prefix], [enable Intel MKL for LAPACK & FFTW])],
     [ac_MKL=${enable_mkl}], [ac_MKL=no])
+AC_ARG_ENABLE([ipp],
+    [AC_HELP_STRING([--enable-ipp=yes|no|prefix], [enable Intel IPP for fast CRC32C])],
+    [ac_IPP=${enable_mkl}], [ac_IPP=no])
 
 case ${ac_MKL} in
     no)
@@ -139,6 +142,17 @@ case ${ac_MKL} in
         AC_DEFINE([USE_MKL], [1], [Define to 1 if you use the Intel MKL]);;
 esac
 
+case ${ac_IPP} in
+    no)
+        ;;
+    yes)
+        AC_DEFINE([USE_IPP], [1], [Define to 1 if you use the Intel IPP]);;
+    *)
+        AM_CXXFLAGS="-I$ac_IPP/include $AM_CXXFLAGS"
+        AM_LDFLAGS="-L$ac_IPP/lib $AM_LDFLAGS"
+        AC_DEFINE([USE_IPP], [1], [Define to 1 if you use the Intel IPP]);;
+esac
+
 ############### HDF5
 AC_ARG_WITH([hdf5],
     [AS_HELP_STRING([--with-hdf5=prefix],
@@ -170,7 +184,13 @@ AC_CHECK_FUNCS([gettimeofday])
 
 if test "${ac_MKL}x" != "nox"; then
     AC_SEARCH_LIBS([mkl_set_interface_layer], [mkl_rt], [],
-                   [AC_MSG_ERROR("MKL enabled but library not found")])
+                   [AC_MSG_ERROR("Intel MKL enabled but library not found")])
+fi
+
+if test "${ac_IPP}x" != "nox"; then
+     AC_SEARCH_LIBS([ippsCRC32C_8u], [ippdc], 
+                   [LIBS="${LIBS} -lippdc -lippvm -lipps -lippcore"],
+                   [AC_MSG_ERROR("Intel IPP enabled but library not found")])
 fi
 
 AC_SEARCH_LIBS([__gmpf_init], [gmp],

From 41032fef44796899a23c451243cb4593a9487da7 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Wed, 21 Nov 2018 18:36:50 +0000
Subject: [PATCH 31/33] Optional RW mode for Hdf5Reader

---
 Grid/serialisation/Hdf5IO.cc | 4 ++--
 Grid/serialisation/Hdf5IO.h  | 2 +-
 Hadrons/A2AMatrix.hpp        | 4 ++--
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/Grid/serialisation/Hdf5IO.cc b/Grid/serialisation/Hdf5IO.cc
index ca2457fa..77396809 100644
--- a/Grid/serialisation/Hdf5IO.cc
+++ b/Grid/serialisation/Hdf5IO.cc
@@ -61,9 +61,9 @@ Group & Hdf5Writer::getGroup(void)
 }
 
 // Reader implementation ///////////////////////////////////////////////////////
-Hdf5Reader::Hdf5Reader(const std::string &fileName)
+Hdf5Reader::Hdf5Reader(const std::string &fileName, const bool readOnly)
 : fileName_(fileName)
-, file_(fileName.c_str(), H5F_ACC_RDWR)
+, file_(fileName.c_str(), readOnly ? H5F_ACC_RDONLY : H5F_ACC_RDWR)
 {
   group_ = file_.openGroup("/");
   readSingleAttribute(dataSetThres_, HDF5_GRID_GUARD "dataset_threshold",
diff --git a/Grid/serialisation/Hdf5IO.h b/Grid/serialisation/Hdf5IO.h
index 1ae2791e..ec26612a 100644
--- a/Grid/serialisation/Hdf5IO.h
+++ b/Grid/serialisation/Hdf5IO.h
@@ -54,7 +54,7 @@ namespace Grid
   class Hdf5Reader: public Reader<Hdf5Reader>
   {
   public:
-    Hdf5Reader(const std::string &fileName);
+    Hdf5Reader(const std::string &fileName, const bool readOnly = true);
     virtual ~Hdf5Reader(void) = default;
     bool push(const std::string &s);
     void pop(void);
diff --git a/Hadrons/A2AMatrix.hpp b/Hadrons/A2AMatrix.hpp
index 7f8791ea..e224a95e 100644
--- a/Hadrons/A2AMatrix.hpp
+++ b/Hadrons/A2AMatrix.hpp
@@ -437,7 +437,7 @@ void A2AMatrixIo<T>::initFile(const MetadataType &d, const unsigned int chunkSiz
     }
 
     // create the dataset
-    Hdf5Reader reader(filename_);
+    Hdf5Reader reader(filename_, false);
 
     push(reader, dataname_);
     auto &group = reader.getGroup();
@@ -457,7 +457,7 @@ void A2AMatrixIo<T>::saveBlock(const T *data,
                                const unsigned int blockSizej)
 {
 #ifdef HAVE_HDF5
-    Hdf5Reader           reader(filename_);
+    Hdf5Reader           reader(filename_, false);
     std::vector<hsize_t> count = {nt_, blockSizei, blockSizej},
                          offset = {0, static_cast<hsize_t>(i),
                                    static_cast<hsize_t>(j)},

From b9e84ecab7471b87ea49e284f84b024c5e748362 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Thu, 22 Nov 2018 15:44:30 +0000
Subject: [PATCH 32/33] Hadrons: minor code cleaning

---
 Hadrons/Utilities/Contractor.cc | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/Hadrons/Utilities/Contractor.cc b/Hadrons/Utilities/Contractor.cc
index 9ce056f3..11ea0bc1 100644
--- a/Hadrons/Utilities/Contractor.cc
+++ b/Hadrons/Utilities/Contractor.cc
@@ -360,8 +360,7 @@ int main(int argc, char* argv[])
             }
             bytes = par.global.nt*lastTerm[0].rows()*lastTerm[0].cols()*sizeof(ComplexD);
             std::cout << Sec(tAr.getDTimer("Transpose caching")) << " " 
-                    << Bytes(bytes, tAr.getDTimer("Transpose caching")) << std::endl;
-            std::cout << Sec(tAr.getDTimer("Disk vector overhead")) << std::endl;
+                      << Bytes(bytes, tAr.getDTimer("Transpose caching")) << std::endl;
             for (unsigned int i = 0; i < timeSeq.size(); ++i)
             {
                 unsigned int dti = 0;
@@ -450,7 +449,6 @@ int main(int argc, char* argv[])
             printTimeProfile(tAr.getTimings(), tAr.getTimer("Total"));
         }
     }
-
     
     return EXIT_SUCCESS;
 }

From e232257cb63c6d5a219d2e452cd2e38584b3a9b2 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Thu, 22 Nov 2018 19:43:49 +0000
Subject: [PATCH 33/33] Hadrons: A2AAslashVector modul cleaning and renaming

---
 Hadrons/Modules.hpp                           |  2 +-
 ...A2AAslashVector.cc => A2AAslashVectors.cc} |  8 +-
 ...AAslashVector.hpp => A2AAslashVectors.hpp} | 99 ++++++++++---------
 Hadrons/modules.inc                           |  4 +-
 4 files changed, 59 insertions(+), 54 deletions(-)
 rename Hadrons/Modules/MSolver/{A2AAslashVector.cc => A2AAslashVectors.cc} (82%)
 rename Hadrons/Modules/MSolver/{A2AAslashVector.hpp => A2AAslashVectors.hpp} (66%)

diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp
index da289307..787fecea 100644
--- a/Hadrons/Modules.hpp
+++ b/Hadrons/Modules.hpp
@@ -24,7 +24,7 @@
 #include <Hadrons/Modules/MSolver/Guesser.hpp>
 #include <Hadrons/Modules/MSolver/RBPrecCG.hpp>
 #include <Hadrons/Modules/MSolver/A2AVectors.hpp>
-#include <Hadrons/Modules/MSolver/A2AAslashVector.hpp>
+#include <Hadrons/Modules/MSolver/A2AAslashVectors.hpp>
 #include <Hadrons/Modules/MGauge/UnitEm.hpp>
 #include <Hadrons/Modules/MGauge/StoutSmearing.hpp>
 #include <Hadrons/Modules/MGauge/Unit.hpp>
diff --git a/Hadrons/Modules/MSolver/A2AAslashVector.cc b/Hadrons/Modules/MSolver/A2AAslashVectors.cc
similarity index 82%
rename from Hadrons/Modules/MSolver/A2AAslashVector.cc
rename to Hadrons/Modules/MSolver/A2AAslashVectors.cc
index 27344397..74e9f59f 100644
--- a/Hadrons/Modules/MSolver/A2AAslashVector.cc
+++ b/Hadrons/Modules/MSolver/A2AAslashVectors.cc
@@ -2,7 +2,7 @@
 
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: Hadrons/Modules/MSolver/A2AAslashVector.cc
+Source file: Hadrons/Modules/MSolver/A2AAslashVectors.cc
 
 Copyright (C) 2015-2018
 
@@ -25,11 +25,11 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-#include <Hadrons/Modules/MSolver/A2AAslashVector.hpp>
+#include <Hadrons/Modules/MSolver/A2AAslashVectors.hpp>
 
 using namespace Grid;
 using namespace Hadrons;
 using namespace MSolver;
 
-template class Grid::Hadrons::MSolver::TA2AAslashVector<FIMPL>;
-template class Grid::Hadrons::MSolver::TA2AAslashVector<ZFIMPL>;
+template class Grid::Hadrons::MSolver::TA2AAslashVectors<FIMPL>;
+template class Grid::Hadrons::MSolver::TA2AAslashVectors<ZFIMPL>;
diff --git a/Hadrons/Modules/MSolver/A2AAslashVector.hpp b/Hadrons/Modules/MSolver/A2AAslashVectors.hpp
similarity index 66%
rename from Hadrons/Modules/MSolver/A2AAslashVector.hpp
rename to Hadrons/Modules/MSolver/A2AAslashVectors.hpp
index 3142d1e7..08c16729 100644
--- a/Hadrons/Modules/MSolver/A2AAslashVector.hpp
+++ b/Hadrons/Modules/MSolver/A2AAslashVectors.hpp
@@ -2,7 +2,7 @@
 
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: Hadrons/Modules/MSolver/A2AAslashVector.hpp
+Source file: Hadrons/Modules/MSolver/A2AAslashVectors.hpp
 
 Copyright (C) 2015-2018
 
@@ -25,13 +25,14 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-#ifndef Hadrons_MSolver_A2AAslashVector_hpp_
-#define Hadrons_MSolver_A2AAslashVector_hpp_
+#ifndef Hadrons_MSolver_A2AAslashVectors_hpp_
+#define Hadrons_MSolver_A2AAslashVectors_hpp_
 
 #include <Hadrons/Global.hpp>
 #include <Hadrons/Module.hpp>
 #include <Hadrons/ModuleFactory.hpp>
 #include <Hadrons/Solver.hpp>
+#include <Hadrons/A2AVectors.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -54,29 +55,30 @@ BEGIN_MODULE_NAMESPACE(MSolver)
 *
 *****************************************************************************/
 
-
-class A2AAslashVectorPar: Serializable
+class A2AAslashVectorsPar: Serializable
 {
 public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(A2AAslashVectorPar,
+  GRID_SERIALIZABLE_CLASS_MEMBERS(A2AAslashVectorsPar,
                                   std::string, vector,
                                   std::string, emField,
-                                  std::string, solver);
+                                  std::string, solver,
+                                  std::string, output,
+                                  bool,        multiFile);
 };
 
 template <typename FImpl>
-class TA2AAslashVector : public Module<A2AAslashVectorPar>
+class TA2AAslashVectors : public Module<A2AAslashVectorsPar>
 {
 public:
     FERM_TYPE_ALIASES(FImpl,);
     SOLVER_TYPE_ALIASES(FImpl,);
 public:
-    typedef PhotonR::GaugeField     EmField;
+    typedef PhotonR::GaugeField EmField;
 public:
     // constructor
-    TA2AAslashVector(const std::string name);
+    TA2AAslashVectors(const std::string name);
     // destructor
-    virtual ~TA2AAslashVector(void) {};
+    virtual ~TA2AAslashVectors(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -88,21 +90,21 @@ private:
     unsigned int Ls_;
 };
 
-MODULE_REGISTER_TMP(A2AAslashVector,TA2AAslashVector<FIMPL>, MSolver);
-MODULE_REGISTER_TMP(ZA2AAslashVector,TA2AAslashVector<ZFIMPL>, MSolver);
+MODULE_REGISTER_TMP(A2AAslashVectors, TA2AAslashVectors<FIMPL>, MSolver);
+MODULE_REGISTER_TMP(ZA2AAslashVectors, TA2AAslashVectors<ZFIMPL>, MSolver);
 
 /******************************************************************************
- *                       TA2AAslashVector implementation                       *
+ *                       TA2AAslashVectors implementation                       *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
 template <typename FImpl>
-TA2AAslashVector<FImpl>::TA2AAslashVector(const std::string name)
-: Module<A2AAslashVectorPar>(name)
+TA2AAslashVectors<FImpl>::TA2AAslashVectors(const std::string name)
+: Module<A2AAslashVectorsPar>(name)
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
 template <typename FImpl>
-std::vector<std::string> TA2AAslashVector<FImpl>::getInput(void)
+std::vector<std::string> TA2AAslashVectors<FImpl>::getInput(void)
 {
     std::vector<std::string> in = {par().vector, par().emField, par().solver};
 
@@ -110,7 +112,7 @@ std::vector<std::string> TA2AAslashVector<FImpl>::getInput(void)
 }
 
 template <typename FImpl>
-std::vector<std::string> TA2AAslashVector<FImpl>::getOutput(void)
+std::vector<std::string> TA2AAslashVectors<FImpl>::getOutput(void)
 {
     std::vector<std::string> out = {getName()};
 
@@ -119,7 +121,7 @@ std::vector<std::string> TA2AAslashVector<FImpl>::getOutput(void)
 
 // setup ///////////////////////////////////////////////////////////////////////
 template <typename FImpl>
-void TA2AAslashVector<FImpl>::setup(void)
+void TA2AAslashVectors<FImpl>::setup(void)
 {
     Ls_  = env().getObjectLs(par().solver);
     auto &vvector = envGet(std::vector<FermionField>, par().vector);
@@ -134,7 +136,7 @@ void TA2AAslashVector<FImpl>::setup(void)
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
-void TA2AAslashVector<FImpl>::execute(void)
+void TA2AAslashVectors<FImpl>::execute(void)
 {
     auto &solver = envGet(Solver, par().solver);
     auto &stoch_photon = envGet(EmField,  par().emField);
@@ -148,42 +150,45 @@ void TA2AAslashVector<FImpl>::execute(void)
 
     Complex ci(0.0,1.0);
 
-
     startTimer("Seq Aslash");
-
-    LOG(Message) << "Calculate Sequential propagator on Aslash * v with the A2A vector " << par().vector
-                  << " and the photon field " << par().emField << std::endl;
-
-
+    LOG(Message) << "Calculate Sequential propagator on Aslash * v with the A2A vector " 
+                 << par().vector << " and the photon field " << par().emField << std::endl;
     for(unsigned int i=0; i<Nmodes; i++)
     {
-	v4dtmp = zero;
-	startTimer("Multiply Aslash");
-	for(unsigned int mu=0;mu<=3;mu++)
-    	{
-	    Gamma gmu(Gamma::gmu[mu]);
-		 v4dtmp +=  ci * PeekIndex<LorentzIndex>(stoch_photon, mu) * (gmu * vvector[i]);
-	}
-	stopTimer("Multiply Aslash");
+        v4dtmp = zero;
+        startTimer("Multiply Aslash");
+        for(unsigned int mu=0;mu<=3;mu++)
+        {
+            Gamma gmu(Gamma::gmu[mu]);
+            v4dtmp +=  ci * PeekIndex<LorentzIndex>(stoch_photon, mu) * (gmu * vvector[i]);
+        }
+        stopTimer("Multiply Aslash");
 
-	if (Ls_ == 1)
-	{
-	    solver(Aslashv[i], v4dtmp);
-	}
-	else
-	{
-	    mat.ImportPhysicalFermionSource(v4dtmp, v5dtmp);
-	    solver(v5dtmp_sol, v5dtmp);
-	    mat.ExportPhysicalFermionSolution(v5dtmp_sol, v4dtmp);
-	    Aslashv[i] = v4dtmp;
-	}
+        startTimer("Inversion");
+        if (Ls_ == 1)
+        {
+            solver(Aslashv[i], v4dtmp);
+        }
+        else
+        {
+            mat.ImportPhysicalFermionSource(v4dtmp, v5dtmp);
+            solver(v5dtmp_sol, v5dtmp);
+            mat.ExportPhysicalFermionSolution(v5dtmp_sol, v4dtmp);
+            Aslashv[i] = v4dtmp;
+        }
+        stopTimer("Inversion");
     }
-
     stopTimer("Seq Aslash");
+    if (!par().output.empty())
+    {
+        startTimer("I/O");
+        A2AVectorsIo::write(par().output, Aslashv, par().multiFile, vm().getTrajectory());
+        stopTimer("I/O");
+    }
 }
 
 END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_MSolver_A2AAslashVector_hpp_
+#endif // Hadrons_MSolver_A2AAslashVectors_hpp_
diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc
index 2add645c..88cfea88 100644
--- a/Hadrons/modules.inc
+++ b/Hadrons/modules.inc
@@ -20,7 +20,7 @@ modules_cc =\
   Modules/MSink/Point.cc \
   Modules/MSink/Smear.cc \
   Modules/MSolver/A2AVectors.cc \
-  Modules/MSolver/A2AAslashVector.cc \
+  Modules/MSolver/A2AAslashVectors.cc \
   Modules/MSolver/RBPrecCG.cc \
   Modules/MSolver/MixedPrecisionRBPrecCG.cc \
   Modules/MSolver/LocalCoherenceLanczos.cc \
@@ -97,7 +97,7 @@ modules_hpp =\
   Modules/MSolver/Guesser.hpp \
   Modules/MSolver/RBPrecCG.hpp \
   Modules/MSolver/A2AVectors.hpp \
-  Modules/MSolver/A2AAslashVector.hpp \
+  Modules/MSolver/A2AAslashVectors.hpp \
   Modules/MGauge/UnitEm.hpp \
   Modules/MGauge/StoutSmearing.hpp \
   Modules/MGauge/Unit.hpp \