Const correctness patch

by Christoph. Differs from the one by Rudy in BFM since it vectorises the twisted
4D solves in pairs.

Grid/algorithms/Algorithms.h

@@ -52,7 +52,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/algorithms/CoarsenedMatrix.h>
 #include <Grid/algorithms/FFT.h>
 
-
 // EigCg
 // Pcg
 // Hdcg

Grid/algorithms/iterative/ConjugateGradient.h

@@ -150,13 +150,13 @@ class ConjugateGradient : public OperatorFunction<Field> {
 	std::cout << GridLogMessage << "\tTrue residual " << true_residual<<std::endl;
 	std::cout << GridLogMessage << "\tTarget " << Tolerance << std::endl;
 
-        std::cout << GridLogMessage << "Time breakdown "<<std::endl;
-	std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tInner      " << InnerTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tAxpyNorm   " << AxpyNormTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
+        std::cout << GridLogPerformance << "Time breakdown "<<std::endl;
+	std::cout << GridLogPerformance << "\tElapsed    " << SolverTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tInner      " << InnerTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tAxpyNorm   " << AxpyNormTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
 
         if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
 

Grid/lattice/Lattice_rng.h

@@ -392,14 +392,10 @@ namespace Grid {
 
     void SeedUniqueString(const std::string &s){
       std::vector<int> seeds;
-      std::stringstream sha;
       seeds = GridChecksum::sha256_seeds(s);
-      for(int i=0;i<seeds.size();i++) { 
-        sha << std::hex << seeds[i];
-      }
       std::cout << GridLogMessage << "Intialising parallel RNG with unique string '" 
                 << s << "'" << std::endl;
-      std::cout << GridLogMessage << "Seed SHA256: " << sha.str() << std::endl;
+      std::cout << GridLogMessage << "Seed SHA256: " << GridChecksum::sha256_string(seeds) << std::endl;
       SeedFixedIntegers(seeds);
     }
     void SeedFixedIntegers(const std::vector<int> &seeds){

Grid/lattice/Lattice_transfer.h

@@ -485,7 +485,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;
 
@@ -520,7 +520,7 @@ void ExtractSliceLocal(Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slic
 
 
 template<class vobj>
-void Replicate(Lattice<vobj> &coarse,Lattice<vobj> & fine)
+void Replicate(const Lattice<vobj> &coarse,Lattice<vobj> & fine)
 {
   typedef typename vobj::scalar_object sobj;
 

Grid/qcd/QCD.h

@@ -90,17 +90,20 @@ namespace QCD {
     // That probably makes for GridRedBlack4dCartesian grid.
 
     // s,sp,c,spc,lc
-    template<typename vtype> using iSinglet                   = iScalar<iScalar<iScalar<vtype> > >;
-    template<typename vtype> using iSpinMatrix                = iScalar<iMatrix<iScalar<vtype>, Ns> >;
-    template<typename vtype> using iColourMatrix              = iScalar<iScalar<iMatrix<vtype, Nc> > > ;
-    template<typename vtype> using iSpinColourMatrix          = iScalar<iMatrix<iMatrix<vtype, Nc>, Ns> >;
-    template<typename vtype> using iLorentzColourMatrix       = iVector<iScalar<iMatrix<vtype, Nc> >, Nd > ;
-    template<typename vtype> using iDoubleStoredColourMatrix  = iVector<iScalar<iMatrix<vtype, Nc> >, Nds > ;
-    template<typename vtype> using iSpinVector                = iScalar<iVector<iScalar<vtype>, Ns> >;
-    template<typename vtype> using iColourVector              = iScalar<iScalar<iVector<vtype, Nc> > >;
-    template<typename vtype> using iSpinColourVector          = iScalar<iVector<iVector<vtype, Nc>, Ns> >;
-    template<typename vtype> using iHalfSpinVector            = iScalar<iVector<iScalar<vtype>, Nhs> >;
-    template<typename vtype> using iHalfSpinColourVector      = iScalar<iVector<iVector<vtype, Nc>, Nhs> >;
+
+    template<typename vtype> using iSinglet                     = iScalar<iScalar<iScalar<vtype> > >;
+    template<typename vtype> using iSpinMatrix                  = iScalar<iMatrix<iScalar<vtype>, Ns> >;
+    template<typename vtype> using iColourMatrix                = iScalar<iScalar<iMatrix<vtype, Nc> > > ;
+    template<typename vtype> using iSpinColourMatrix            = iScalar<iMatrix<iMatrix<vtype, Nc>, Ns> >;
+    template<typename vtype> using iLorentzColourMatrix         = iVector<iScalar<iMatrix<vtype, Nc> >, Nd > ;
+    template<typename vtype> using iDoubleStoredColourMatrix    = iVector<iScalar<iMatrix<vtype, Nc> >, Nds > ;
+    template<typename vtype> using iSpinVector                  = iScalar<iVector<iScalar<vtype>, Ns> >;
+    template<typename vtype> using iColourVector                = iScalar<iScalar<iVector<vtype, Nc> > >;
+    template<typename vtype> using iSpinColourVector            = iScalar<iVector<iVector<vtype, Nc>, Ns> >;
+    template<typename vtype> using iHalfSpinVector              = iScalar<iVector<iScalar<vtype>, Nhs> >;
+    template<typename vtype> using iHalfSpinColourVector        = iScalar<iVector<iVector<vtype, Nc>, Nhs> >;
+    template<typename vtype> using iSpinColourSpinColourMatrix  = iScalar<iMatrix<iMatrix<iMatrix<iMatrix<vtype, Nc>, Ns>, Nc>, Ns> >;
+
 
     template<typename vtype> using iGparitySpinColourVector       = iVector<iVector<iVector<vtype, Nc>, Ns>, Ngp >;
     template<typename vtype> using iGparityHalfSpinColourVector   = iVector<iVector<iVector<vtype, Nc>, Nhs>, Ngp >;
@@ -127,10 +130,28 @@ namespace QCD {
     typedef iSpinColourMatrix<Complex  >    SpinColourMatrix;
     typedef iSpinColourMatrix<ComplexF >    SpinColourMatrixF;
     typedef iSpinColourMatrix<ComplexD >    SpinColourMatrixD;
-
+    
     typedef iSpinColourMatrix<vComplex >    vSpinColourMatrix;
     typedef iSpinColourMatrix<vComplexF>    vSpinColourMatrixF;
     typedef iSpinColourMatrix<vComplexD>    vSpinColourMatrixD;
+    
+    // SpinColourSpinColour matrix
+    typedef iSpinColourSpinColourMatrix<Complex  >    SpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<ComplexF >    SpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<ComplexD >    SpinColourSpinColourMatrixD;
+
+    typedef iSpinColourSpinColourMatrix<vComplex >    vSpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<vComplexF>    vSpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<vComplexD>    vSpinColourSpinColourMatrixD;
+
+    // SpinColourSpinColour matrix
+    typedef iSpinColourSpinColourMatrix<Complex  >    SpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<ComplexF >    SpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<ComplexD >    SpinColourSpinColourMatrixD;
+
+    typedef iSpinColourSpinColourMatrix<vComplex >    vSpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<vComplexF>    vSpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<vComplexD>    vSpinColourSpinColourMatrixD;
 
     // LorentzColour
     typedef iLorentzColourMatrix<Complex  > LorentzColourMatrix;
@@ -229,6 +250,9 @@ namespace QCD {
     typedef Lattice<vSpinColourMatrixF>     LatticeSpinColourMatrixF;
     typedef Lattice<vSpinColourMatrixD>     LatticeSpinColourMatrixD;
 
+    typedef Lattice<vSpinColourSpinColourMatrix>      LatticeSpinColourSpinColourMatrix;
+    typedef Lattice<vSpinColourSpinColourMatrixF>     LatticeSpinColourSpinColourMatrixF;
+    typedef Lattice<vSpinColourSpinColourMatrixD>     LatticeSpinColourSpinColourMatrixD;
 
     typedef Lattice<vLorentzColourMatrix>  LatticeLorentzColourMatrix;
     typedef Lattice<vLorentzColourMatrixF> LatticeLorentzColourMatrixF;

Grid/qcd/action/fermion/CayleyFermion5D.cc

@@ -485,9 +485,13 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,std::vector<Co
     
   double bpc = b+c;
   double bmc = b-c;
+  _b = b;
+  _c = c;
+  _gamma  = gamma; // Save the parameters so we can change mass later.
+  _zolo_hi= zolo_hi;
   for(int i=0; i < Ls; i++){
     as[i] = 1.0;
-    omega[i] = gamma[i]*zolo_hi; //NB reciprocal relative to Chroma NEF code
+    omega[i] = _gamma[i]*_zolo_hi; //NB reciprocal relative to Chroma NEF code
     assert(omega[i]!=Coeff_t(0.0));
     bs[i] = 0.5*(bpc/omega[i] + bmc);
     cs[i] = 0.5*(bpc/omega[i] - bmc);

Grid/qcd/action/fermion/CayleyFermion5D.h

@@ -97,7 +97,10 @@ namespace Grid {
       // Support for MADWF tricks
       ///////////////////////////////////////////////////////////////
       RealD Mass(void) { return mass; };
-      void  SetMass(RealD _mass) { mass=_mass; } ;
+      void  SetMass(RealD _mass) { 
+	mass=_mass; 
+	SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c);  // Reset coeffs
+      } ;
       void  P(const FermionField &psi, FermionField &chi);
       void  Pdag(const FermionField &psi, FermionField &chi);
 
@@ -147,6 +150,12 @@ namespace Grid {
       //    protected:
       RealD mass;
 
+      // Save arguments to SetCoefficientsInternal
+      std::vector<Coeff_t> _gamma;
+      RealD                _zolo_hi;
+      RealD                _b;
+      RealD                _c;
+
       // Cayley form Moebius (tanh and zolotarev)
       std::vector<Coeff_t> omega;
       std::vector<Coeff_t> bs;    // S dependent coeffs

Grid/qcd/action/fermion/Fermion.h

@@ -80,12 +80,24 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
 ///////////////////////////////////////////////////////////////////////////////
 #include <Grid/qcd/action/fermion/g5HermitianLinop.h>
 
+///////////////////////////////////////////////////////////////////////////////
+// Fourier accelerated Pauli Villars inverse support
+///////////////////////////////////////////////////////////////////////////////
+#include <Grid/qcd/action/fermion/WilsonTMFermion5D.h>   
+
+////////////////////////////////////////////////////////////////////////////////
+// Move this group to a DWF specific tools/algorithms subdir? 
+////////////////////////////////////////////////////////////////////////////////
+#include <Grid/qcd/action/fermion/FourierAcceleratedPV.h>
+#include <Grid/qcd/action/fermion/PauliVillarsInverters.h>
+#include <Grid/qcd/action/fermion/Reconstruct5Dprop.h>
+#include <Grid/qcd/action/fermion/MADWF.h>
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // More maintainable to maintain the following typedef list centrally, as more "impl" targets
 // are added, (e.g. extension for gparity, half precision project in comms etc..)
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 
-
 // Cayley 5d
 namespace Grid {
   namespace QCD {

Grid/qcd/action/fermion/FermionOperator.h

@@ -68,7 +68,6 @@ namespace Grid {
       virtual int    ConstEE(void) { return 1; }; // clover returns zero as EE depends on gauge field
       virtual int    isTrivialEE(void) { return 0; };
       virtual RealD  Mass(void) {return 0.0;};
-      virtual void SetMass(RealD _mass) { return; };
 
       // half checkerboard operaions
       virtual void   Meooe       (const FermionField &in, FermionField &out)=0;

Grid/qcd/action/fermion/FermionOperatorImpl.h

@@ -141,6 +141,7 @@ namespace QCD {
   ////////////////////////////////////////////////////////////////////////
   
 #define INHERIT_FIMPL_TYPES(Impl)\
+  typedef Impl Impl_t;							\
   typedef typename Impl::FermionField           FermionField;		\
   typedef typename Impl::PropagatorField     PropagatorField;		\
   typedef typename Impl::DoubledGaugeField DoubledGaugeField;		\

Grid/qcd/action/fermion/FourierAcceleratedPV.h

@@ -0,0 +1,237 @@
+
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/qcd/action/fermion/FourierAcceleratedPV.h
+
+    Copyright (C) 2015
+
+Author: Christoph Lehner (lifted with permission by Peter Boyle, brought back to Grid)
+Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
+
+    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 */
+#pragma once
+namespace Grid {
+namespace QCD {
+
+  template<typename M>
+    void get_real_const_bc(M& m, RealD& _b, RealD& _c) {
+    ComplexD b,c;
+    b=m.bs[0];
+    c=m.cs[0];
+    std::cout << GridLogMessage << "b=" << b << ", c=" << c << std::endl;
+    for (size_t i=1;i<m.bs.size();i++) {
+      assert(m.bs[i] == b);
+      assert(m.cs[i] == c);
+    }
+    assert(b.imag() == 0.0);
+    assert(c.imag() == 0.0);
+    _b = b.real();
+    _c = c.real();
+  }
+
+
+template<typename Vi, typename M, typename G>
+class FourierAcceleratedPV {
+ public:
+
+  ConjugateGradient<Vi> &cg;
+  M& dwfPV;
+  G& Umu;
+  GridCartesian* grid5D;
+  GridRedBlackCartesian* gridRB5D;
+  int group_in_s;
+
+  FourierAcceleratedPV(M& _dwfPV, G& _Umu, ConjugateGradient<Vi> &_cg, int _group_in_s = 2) 
+   : dwfPV(_dwfPV), Umu(_Umu), cg(_cg), group_in_s(_group_in_s) 
+  {
+    assert( dwfPV.FermionGrid()->_fdimensions[0] % (2*group_in_s) == 0);
+    grid5D = QCD::SpaceTimeGrid::makeFiveDimGrid(2*group_in_s, (GridCartesian*)Umu._grid);
+    gridRB5D = QCD::SpaceTimeGrid::makeFiveDimRedBlackGrid(2*group_in_s, (GridCartesian*)Umu._grid);
+  }
+
+  void rotatePV(const Vi& _src, Vi& dst, bool forward) const {
+
+    GridStopWatch gsw1, gsw2;
+
+    typedef typename Vi::scalar_type Coeff_t;
+    int Ls = dst._grid->_fdimensions[0];
+
+    Vi _tmp(dst._grid);
+    double phase = M_PI / (double)Ls;
+    Coeff_t bzero(0.0,0.0);
+
+    FFT theFFT((GridCartesian*)dst._grid);
+
+    if (!forward) {
+      gsw1.Start();
+      for (int s=0;s<Ls;s++) {
+	Coeff_t a(::cos(phase*s),-::sin(phase*s));
+	axpby_ssp(_tmp,a,_src,bzero,_src,s,s);
+      }
+      gsw1.Stop();
+
+      gsw2.Start();
+      theFFT.FFT_dim(dst,_tmp,0,FFT::forward);
+      gsw2.Stop();
+
+    } else {
+
+      gsw2.Start();
+      theFFT.FFT_dim(_tmp,_src,0,FFT::backward);
+      gsw2.Stop();
+
+      gsw1.Start();
+      for (int s=0;s<Ls;s++) {
+	Coeff_t a(::cos(phase*s),::sin(phase*s));
+	axpby_ssp(dst,a,_tmp,bzero,_tmp,s,s);
+      }
+      gsw1.Stop();
+    }
+
+    std::cout << GridLogMessage << "Timing rotatePV: " << gsw1.Elapsed() << ", " << gsw2.Elapsed() << std::endl;
+
+  }
+
+  void pvInv(const Vi& _src, Vi& _dst) const {
+
+    std::cout << GridLogMessage << "Fourier-Accelerated Outer Pauli Villars"<<std::endl;
+
+    typedef typename Vi::scalar_type Coeff_t;
+    int Ls = _dst._grid->_fdimensions[0];
+
+    GridStopWatch gswT;
+    gswT.Start();
+
+    RealD b,c;
+    get_real_const_bc(dwfPV,b,c);
+    RealD M5 = dwfPV.M5;
+    
+    // U(true) Rightinv TMinv U(false) = Minv
+
+    Vi _src_diag(_dst._grid);
+    Vi _src_diag_slice(dwfPV.GaugeGrid());
+    Vi _dst_diag_slice(dwfPV.GaugeGrid());
+    Vi _src_diag_slices(grid5D);
+    Vi _dst_diag_slices(grid5D);
+    Vi _dst_diag(_dst._grid);
+
+    rotatePV(_src,_src_diag,false);
+
+    // now do TM solves
+    Gamma G5(Gamma::Algebra::Gamma5);
+
+    GridStopWatch gswA, gswB;
+
+    gswA.Start();
+
+    typedef typename M::Impl_t Impl;
+    //WilsonTMFermion<Impl> tm(x.Umu,*x.UGridF,*x.UrbGridF,0.0,0.0,solver_outer.parent.par.wparams_f);
+    std::vector<RealD> vmass(grid5D->_fdimensions[0],0.0);
+    std::vector<RealD> vmu(grid5D->_fdimensions[0],0.0);
+
+    WilsonTMFermion5D<Impl> tm(Umu,*grid5D,*gridRB5D,
+			   *(GridCartesian*)dwfPV.GaugeGrid(),
+			   *(GridRedBlackCartesian*)dwfPV.GaugeRedBlackGrid(),
+			   vmass,vmu);
+    
+    //SchurRedBlackDiagTwoSolve<Vi> sol(cg);
+    SchurRedBlackDiagMooeeSolve<Vi> sol(cg); // same performance as DiagTwo
+    gswA.Stop();
+
+    gswB.Start();
+
+    for (int sgroup=0;sgroup<Ls/2/group_in_s;sgroup++) {
+
+      for (int sidx=0;sidx<group_in_s;sidx++) {
+
+	int s = sgroup*group_in_s + sidx;
+	int sprime = Ls-s-1;
+
+	RealD phase = M_PI / (RealD)Ls * (2.0 * s + 1.0);
+	RealD cosp = ::cos(phase);
+	RealD sinp = ::sin(phase);
+	RealD denom = b*b + c*c + 2.0*b*c*cosp;
+	RealD mass = -(b*b*M5 + c*(1.0 - cosp + c*M5) + b*(-1.0 + cosp + 2.0*c*cosp*M5))/denom;
+	RealD mu = (b+c)*sinp/denom;
+
+	vmass[2*sidx + 0] = mass;
+	vmass[2*sidx + 1] = mass;
+	vmu[2*sidx + 0] = mu;
+	vmu[2*sidx + 1] = -mu;
+
+      }
+
+      tm.update(vmass,vmu);
+
+      for (int sidx=0;sidx<group_in_s;sidx++) {
+
+	int s = sgroup*group_in_s + sidx;
+	int sprime = Ls-s-1;
+
+	ExtractSlice(_src_diag_slice,_src_diag,s,0);
+	InsertSlice(_src_diag_slice,_src_diag_slices,2*sidx + 0,0);
+
+	ExtractSlice(_src_diag_slice,_src_diag,sprime,0);
+	InsertSlice(_src_diag_slice,_src_diag_slices,2*sidx + 1,0);
+
+      }
+
+      GridStopWatch gsw;
+      gsw.Start();
+      _dst_diag_slices = zero; // zero guess
+      sol(tm,_src_diag_slices,_dst_diag_slices);
+      gsw.Stop();
+      std::cout << GridLogMessage << "Solve[sgroup=" << sgroup << "] completed in " << gsw.Elapsed() << ", " << gswA.Elapsed() << std::endl;
+
+      for (int sidx=0;sidx<group_in_s;sidx++) {
+
+	int s = sgroup*group_in_s + sidx;
+	int sprime = Ls-s-1;
+
+	RealD phase = M_PI / (RealD)Ls * (2.0 * s + 1.0);
+	RealD cosp = ::cos(phase);
+	RealD sinp = ::sin(phase);
+
+	// now rotate with inverse of
+	Coeff_t pA = b + c*cosp;
+	Coeff_t pB = - Coeff_t(0.0,1.0)*c*sinp;
+	Coeff_t pABden = pA*pA - pB*pB;
+	// (pA + pB * G5) * (pA - pB*G5) = (pA^2 - pB^2)
+      
+	ExtractSlice(_dst_diag_slice,_dst_diag_slices,2*sidx + 0,0);
+	_dst_diag_slice = (pA/pABden) * _dst_diag_slice - (pB/pABden) * (G5 * _dst_diag_slice);
+	InsertSlice(_dst_diag_slice,_dst_diag,s,0);
+	
+	ExtractSlice(_dst_diag_slice,_dst_diag_slices,2*sidx + 1,0);
+	_dst_diag_slice = (pA/pABden) * _dst_diag_slice + (pB/pABden) * (G5 * _dst_diag_slice);
+	InsertSlice(_dst_diag_slice,_dst_diag,sprime,0);
+      }
+    }
+    gswB.Stop();
+
+    rotatePV(_dst_diag,_dst,true);
+
+    gswT.Stop();
+    std::cout << GridLogMessage << "PV completed in " << gswT.Elapsed() << " (Setup: " << gswA.Elapsed() << ", s-loop: " << gswB.Elapsed() << ")" << std::endl;
+  }
+
+};
+}}

Grid/qcd/action/fermion/MADWF.h

@@ -0,0 +1,193 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/MADWF.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    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 */
+#pragma once
+
+namespace Grid {
+namespace QCD {
+
+template <class Fieldi, class Fieldo,IfNotSame<Fieldi,Fieldo> X=0>
+inline void convert(const Fieldi &from,Fieldo &to) 
+{
+  precisionChange(to,from);
+}
+template <class Fieldi, class Fieldo,IfSame<Fieldi,Fieldo> X=0>
+inline void convert(const Fieldi &from,Fieldo &to) 
+{
+  to=from;
+}
+
+template<class Matrixo,class Matrixi,class PVinverter,class SchurSolver, class Guesser> 
+class MADWF 
+{
+ private:
+  typedef typename Matrixo::FermionField FermionFieldo;
+  typedef typename Matrixi::FermionField FermionFieldi;
+
+  PVinverter  & PauliVillarsSolvero;// For the outer field
+  SchurSolver & SchurSolveri;       // For the inner approx field
+  Guesser     & Guesseri;           // To deflate the inner approx solves
+
+  Matrixo & Mato;                   // Action object for outer
+  Matrixi & Mati;                   // Action object for inner
+
+  RealD target_resid;
+  int   maxiter;
+ public:
+
+  MADWF(Matrixo &_Mato,
+	Matrixi &_Mati, 
+	PVinverter &_PauliVillarsSolvero, 
+	SchurSolver &_SchurSolveri,
+	Guesser & _Guesseri,
+	RealD resid,
+	int _maxiter) :
+
+  Mato(_Mato),Mati(_Mati),
+    SchurSolveri(_SchurSolveri),
+    PauliVillarsSolvero(_PauliVillarsSolvero),Guesseri(_Guesseri)
+  {   
+    target_resid=resid;
+    maxiter     =_maxiter; 
+  };
+
+  void operator() (const FermionFieldo &src4,FermionFieldo &sol5)
+  {
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+    std::cout << GridLogMessage<< "  MADWF-like algorithm                           " << std::endl;
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+
+    FermionFieldi    c0i(Mati.GaugeGrid()); // 4d 
+    FermionFieldi    y0i(Mati.GaugeGrid()); // 4d
+    FermionFieldo    c0 (Mato.GaugeGrid()); // 4d 
+    FermionFieldo    y0 (Mato.GaugeGrid()); // 4d
+
+    FermionFieldo    A(Mato.FermionGrid()); // Temporary outer
+    FermionFieldo    B(Mato.FermionGrid()); // Temporary outer
+    FermionFieldo    b(Mato.FermionGrid()); // 5d source
+
+    FermionFieldo    c(Mato.FermionGrid()); // PVinv source; reused so store
+    FermionFieldo    defect(Mato.FermionGrid()); // 5d source
+
+    FermionFieldi   ci(Mati.FermionGrid()); 
+    FermionFieldi   yi(Mati.FermionGrid()); 
+    FermionFieldi   xi(Mati.FermionGrid()); 
+    FermionFieldi srci(Mati.FermionGrid()); 
+    FermionFieldi   Ai(Mati.FermionGrid()); 
+
+    RealD m=Mati.Mass();
+
+    ///////////////////////////////////////
+    //Import source, include Dminus factors
+    ///////////////////////////////////////
+    Mato.ImportPhysicalFermionSource(src4,b); 
+    std::cout << GridLogMessage << " src4 " <<norm2(src4)<<std::endl;
+    std::cout << GridLogMessage << " b    " <<norm2(b)<<std::endl;
+
+    defect = b;
+    sol5=zero;
+    for (int i=0;i<maxiter;i++) {
+
+      ///////////////////////////////////////
+      // Set up c0 from current defect
+      ///////////////////////////////////////
+      PauliVillarsSolvero(Mato,defect,A);
+      Mato.Pdag(A,c);
+      ExtractSlice(c0, c, 0 , 0);
+
+      ////////////////////////////////////////////////
+      // Solve the inner system with surface term c0
+      ////////////////////////////////////////////////
+      ci = zero;  
+      convert(c0,c0i); // Possible precison change
+      InsertSlice(c0i,ci,0, 0);
+
+      // Dwm P y = Dwm x = D(1) P (c0,0,0,0)^T
+      Mati.P(ci,Ai);
+      Mati.SetMass(1.0);      Mati.M(Ai,srci);      Mati.SetMass(m);
+      SchurSolveri(Mati,srci,xi,Guesseri); 
+      Mati.Pdag(xi,yi);
+      ExtractSlice(y0i, yi, 0 , 0);
+      convert(y0i,y0); // Possible precision change
+
+      //////////////////////////////////////
+      // Propagate solution back to outer system
+      // Build Pdag PV^-1 Dm P [-sol4,c2,c3... cL]
+      //////////////////////////////////////
+      c0 = - y0;
+      InsertSlice(c0, c, 0   , 0);
+
+      /////////////////////////////
+      // Reconstruct the bulk solution Pdag PV^-1 Dm P 
+      /////////////////////////////
+      Mato.P(c,B);
+      Mato.M(B,A);
+      PauliVillarsSolvero(Mato,A,B);
+      Mato.Pdag(B,A);
+
+      //////////////////////////////
+      // Reinsert surface prop
+      //////////////////////////////
+      InsertSlice(y0,A,0,0);
+
+      //////////////////////////////
+      // Convert from y back to x 
+      //////////////////////////////
+      Mato.P(A,B);
+
+      //         sol5' = sol5 + M^-1 defect
+      //               = sol5 + M^-1 src - M^-1 M sol5  ...
+      sol5 = sol5 + B;
+      std::cout << GridLogMessage << "***************************************" <<std::endl;
+      std::cout << GridLogMessage << " Sol5 update "<<std::endl;
+      std::cout << GridLogMessage << "***************************************" <<std::endl;
+      std::cout << GridLogMessage << " Sol5 now "<<norm2(sol5)<<std::endl;
+      std::cout << GridLogMessage << " delta    "<<norm2(B)<<std::endl;
+
+       // New defect  = b - M sol5
+       Mato.M(sol5,A);
+       defect = b - A;
+
+       std::cout << GridLogMessage << " defect   "<<norm2(defect)<<std::endl;
+
+       double resid = ::sqrt(norm2(defect) / norm2(b));
+       std::cout << GridLogMessage << "Residual " << i << ": " << resid  << std::endl;
+       std::cout << GridLogMessage << "***************************************" <<std::endl;
+
+       if (resid < target_resid) {
+	 return;
+       }
+    }
+
+    std::cout << GridLogMessage << "MADWF : Exceeded maxiter "<<std::endl;
+    assert(0);
+
+  }
+
+};
+
+}}

Grid/qcd/action/fermion/PauliVillarsInverters.h

@@ -0,0 +1,95 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/SchurRedBlack.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    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 */
+#pragma once
+
+namespace Grid {
+namespace QCD {
+
+template<class Field>
+class PauliVillarsSolverUnprec
+{
+ public:
+  ConjugateGradient<Field> & CG;
+  PauliVillarsSolverUnprec(  ConjugateGradient<Field> &_CG) : CG(_CG){};
+
+  template<class Matrix>
+  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
+  {
+    RealD m = _Matrix.Mass();
+    Field A  (_Matrix.FermionGrid());
+
+    MdagMLinearOperator<Matrix,Field> HermOp(_Matrix);
+
+    _Matrix.SetMass(1.0);
+    _Matrix.Mdag(src,A);
+    CG(HermOp,A,sol);
+    _Matrix.SetMass(m);
+  };
+};
+
+template<class Field,class SchurSolverType>
+class PauliVillarsSolverRBprec
+{
+ public:
+  SchurSolverType & SchurSolver;
+  PauliVillarsSolverRBprec( SchurSolverType &_SchurSolver) : SchurSolver(_SchurSolver){};
+
+  template<class Matrix>
+  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
+  {
+    RealD m = _Matrix.Mass();
+    Field A  (_Matrix.FermionGrid());
+
+    _Matrix.SetMass(1.0);
+    SchurSolver(_Matrix,src,sol);
+    _Matrix.SetMass(m);
+  };
+};
+
+template<class Field,class GaugeField>
+class PauliVillarsSolverFourierAccel
+{
+ public:
+  GaugeField      & Umu;
+  ConjugateGradient<Field> & CG;
+
+  PauliVillarsSolverFourierAccel(GaugeField &_Umu,ConjugateGradient<Field> &_CG) :  Umu(_Umu), CG(_CG)
+  {
+  };
+
+  template<class Matrix>
+  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
+  {
+    FourierAcceleratedPV<Field, Matrix, typename Matrix::GaugeField > faPV(_Matrix,Umu,CG) ;
+    faPV.pvInv(src,sol);
+  };
+};
+
+
+}
+}

Grid/qcd/action/fermion/Reconstruct5Dprop.h

@@ -30,49 +30,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 namespace Grid {
 namespace QCD {
 
-
-template<class Field>
-class PauliVillarsSolverUnprec
-{
- public:
-  ConjugateGradient<Field> & CG;
-  PauliVillarsSolverUnprec(  ConjugateGradient<Field> &_CG) : CG(_CG){};
-
-  template<class Matrix>
-  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
-  {
-    RealD m = _Matrix.Mass();
-    Field A  (_Matrix.FermionGrid());
-
-    MdagMLinearOperator<Matrix,Field> HermOp(_Matrix);
-
-    _Matrix.SetMass(1.0);
-    _Matrix.Mdag(src,A);
-    CG(HermOp,A,sol);
-    _Matrix.SetMass(m);
-  };
-};
-
-template<class Field>
-class PauliVillarsSolverRBprec
-{
- public:
-  ConjugateGradient<Field> & CG;
-  PauliVillarsSolverRBprec(  ConjugateGradient<Field> &_CG) : CG(_CG){};
-
-  template<class Matrix>
-  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
-  {
-    RealD m = _Matrix.Mass();
-    Field A  (_Matrix.FermionGrid());
-
-    _Matrix.SetMass(1.0);
-    SchurRedBlackDiagMooeeSolve<Field> SchurSolver(CG);
-    SchurSolver(_Matrix,src,sol);
-    _Matrix.SetMass(m);
-  };
-};
-
 template<class Field,class PVinverter> class Reconstruct5DfromPhysical {
  private:
   PVinverter & PauliVillarsSolver;
@@ -85,20 +42,12 @@ template<class Field,class PVinverter> class Reconstruct5DfromPhysical {
  // of the Mobius exact AMA corrections.
  //
  // TODO : understand absence of contact term in eqns in Hantao's thesis
- //        sol4 is contact term subtracted.
+ //        sol4 is contact term subtracted, but thesis & Brower's paper suggests not.
  //
- // Options
- // a) Defect correction approach:
- //    1) Compute defect from current soln (initially guess).
- //       This is ...... outerToInner check !!!!
- //    2) Deflated Zmobius solve to get 4d soln
- //       Ensure deflation is working
- //    3) Refine 5d Outer using the inner 4d delta soln
- // 
- // Step 1: localise PV inverse in a routine. [DONE]
+ // Step 1: Localise PV inverse in a routine. [DONE]
  // Step 2: Schur based PV inverse            [DONE]
- // Step 3: Fourier accelerated PV inverse
- // Step 4: 
+ // Step 3: Fourier accelerated PV inverse    [DONE]
+ //
  /////////////////////////////////////////////////////
  
   Reconstruct5DfromPhysical(PVinverter &_PauliVillarsSolver) 

Grid/qcd/action/fermion/WilsonTMFermion5D.h

@@ -0,0 +1,155 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/qcd/action/fermion/WilsonTMFermion5D.h
+
+    Copyright (C) 2015
+
+Author: paboyle <paboyle@ph.ed.ac.uk> ; NB Christoph did similar in GPT
+
+    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 */
+#pragma once 
+
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/WilsonFermion.h>
+
+
+namespace Grid {
+
+  namespace QCD {
+    
+    template<class Impl>
+      class WilsonTMFermion5D : public WilsonFermion5D<Impl>
+      {
+      public:
+	INHERIT_IMPL_TYPES(Impl);
+      public:
+
+	virtual void   Instantiatable(void) {};
+
+	// Constructors
+        WilsonTMFermion5D(GaugeField &_Umu,
+			  GridCartesian         &Fgrid,
+			  GridRedBlackCartesian &Frbgrid, 
+			  GridCartesian         &Ugrid,
+			  GridRedBlackCartesian &Urbgrid, 
+			  const std::vector<RealD> _mass,
+			  const std::vector<RealD> _mu,
+			  const ImplParams &p= ImplParams()
+			  ) :
+	WilsonFermion5D<Impl>(_Umu,
+			      Fgrid,
+			      Frbgrid,
+			      Ugrid,
+			      Urbgrid,
+			      4.0,p)
+	
+	  {
+	    update(_mass,_mu);
+	  }
+
+	virtual void Meooe(const FermionField &in, FermionField &out) {
+	  if (in.checkerboard == Odd) {
+	    this->DhopEO(in, out, DaggerNo);
+	  } else {
+	    this->DhopOE(in, out, DaggerNo);
+	  }
+	}
+
+	virtual void MeooeDag(const FermionField &in, FermionField &out) {
+	  if (in.checkerboard == Odd) {
+	    this->DhopEO(in, out, DaggerYes);
+	  } else {
+	    this->DhopOE(in, out, DaggerYes);
+	  }
+	}	
+	
+	// allow override for twisted mass and clover
+	virtual void Mooee(const FermionField &in, FermionField &out) {
+	  out.checkerboard = in.checkerboard;
+	  //axpibg5x(out,in,a,b); // out = a*in + b*i*G5*in
+	  for (int s=0;s<(int)this->mass.size();s++) {
+	    ComplexD a = 4.0+this->mass[s];
+	    ComplexD b(0.0,this->mu[s]);
+	    axpbg5y_ssp(out,a,in,b,in,s,s);
+	  }
+	}
+
+	virtual void MooeeDag(const FermionField &in, FermionField &out) {
+	  out.checkerboard = in.checkerboard;
+	  for (int s=0;s<(int)this->mass.size();s++) {
+	    ComplexD a = 4.0+this->mass[s];
+	    ComplexD b(0.0,-this->mu[s]);
+	    axpbg5y_ssp(out,a,in,b,in,s,s);
+	  }
+	}
+	virtual void MooeeInv(const FermionField &in, FermionField &out) {
+	  for (int s=0;s<(int)this->mass.size();s++) {
+	    RealD m    = this->mass[s];
+	    RealD tm   = this->mu[s];
+	    RealD mtil = 4.0+this->mass[s];
+	    RealD sq   = mtil*mtil+tm*tm;
+	    ComplexD a    = mtil/sq;
+	    ComplexD b(0.0, -tm /sq);
+	    axpbg5y_ssp(out,a,in,b,in,s,s);
+	  }
+	}
+	virtual void MooeeInvDag(const FermionField &in, FermionField &out) {
+	  for (int s=0;s<(int)this->mass.size();s++) {
+	    RealD m    = this->mass[s];
+	    RealD tm   = this->mu[s];
+	    RealD mtil = 4.0+this->mass[s];
+	    RealD sq   = mtil*mtil+tm*tm;
+	    ComplexD a    = mtil/sq;
+	    ComplexD b(0.0,tm /sq);
+	    axpbg5y_ssp(out,a,in,b,in,s,s);
+	  }
+	}
+
+	virtual RealD M(const FermionField &in, FermionField &out) {
+	  out.checkerboard = in.checkerboard;
+	  this->Dhop(in, out, DaggerNo);
+	  FermionField tmp(out._grid);
+	  for (int s=0;s<(int)this->mass.size();s++) {
+	    ComplexD a = 4.0+this->mass[s];
+	    ComplexD b(0.0,this->mu[s]);
+	    axpbg5y_ssp(tmp,a,in,b,in,s,s);
+	  }
+	  return axpy_norm(out, 1.0, tmp, out);
+	}
+	
+	// needed for fast PV
+	void update(const std::vector<RealD>& _mass, const std::vector<RealD>& _mu) {
+	  assert(_mass.size() == _mu.size());
+	  assert(_mass.size() == this->FermionGrid()->_fdimensions[0]);
+	  this->mass = _mass;
+	  this->mu = _mu;
+	}
+	
+      private:
+	std::vector<RealD> mu;
+	std::vector<RealD> mass;
+	
+      };
+   
+    typedef WilsonTMFermion5D<WilsonImplF> WilsonTMFermion5DF; 
+    typedef WilsonTMFermion5D<WilsonImplD> WilsonTMFermion5DD; 
+
+}}

Grid/qcd/utils/A2Autils.h

@@ -27,12 +27,13 @@ public:
 
   typedef iSpinColourMatrix<vector_type> SpinColourMatrix_v;
 
-  static void MesonField(Eigen::Tensor<ComplexD,5> &mat, 
+  template <typename TensorType> // output: rank 5 tensor, e.g. Eigen::Tensor<ComplexD, 5>
+  static void MesonField(TensorType &mat, 
 			 const FermionField *lhs_wi,
 			 const FermionField *rhs_vj,
 			 std::vector<Gamma::Algebra> gammas,
 			 const std::vector<ComplexField > &mom,
-			 int orthogdim);
+			 int orthogdim, double *t_kernel = nullptr, double *t_gsum = nullptr);
 
   static void PionFieldWVmom(Eigen::Tensor<ComplexD,4> &mat, 
 			     const FermionField *wi,
@@ -59,6 +60,14 @@ public:
 			  const FermionField *vj,
 			  int orthogdim);
 
+  template <typename TensorType> // output: rank 5 tensor, e.g. Eigen::Tensor<ComplexD, 5>
+  static void AslashField(TensorType &mat, 
+        const FermionField *lhs_wi,
+        const FermionField *rhs_vj,
+        const std::vector<ComplexField> &emB0,
+        const std::vector<ComplexField> &emB1,
+        int orthogdim, double *t_kernel = nullptr, double *t_gsum = nullptr);
+
   static void ContractWWVV(std::vector<PropagatorField> &WWVV,
 			   const Eigen::Tensor<ComplexD,3> &WW_sd,
 			   const FermionField *vs,
@@ -92,13 +101,14 @@ public:
 #endif
 };
 
-template<class FImpl>
-void A2Autils<FImpl>::MesonField(Eigen::Tensor<ComplexD,5> &mat, 
+template <class FImpl>
+template <typename TensorType>
+void A2Autils<FImpl>::MesonField(TensorType &mat, 
 				 const FermionField *lhs_wi,
 				 const FermionField *rhs_vj,
 				 std::vector<Gamma::Algebra> gammas,
 				 const std::vector<ComplexField > &mom,
-				 int orthogdim) 
+				 int orthogdim, double *t_kernel, double *t_gsum) 
 {
   typedef typename FImpl::SiteSpinor vobj;
 
@@ -146,6 +156,7 @@ void A2Autils<FImpl>::MesonField(Eigen::Tensor<ComplexD,5> &mat,
   int stride=grid->_slice_stride[orthogdim];
 
   // potentially wasting cores here if local time extent too small
+  if (t_kernel) *t_kernel = -usecond();
   parallel_for(int r=0;r<rd;r++){
 
     int so=r*grid->_ostride[orthogdim]; // base offset for start of plane 
@@ -212,7 +223,7 @@ void A2Autils<FImpl>::MesonField(Eigen::Tensor<ComplexD,5> &mat,
       }
     }}}
   }
-
+  if (t_kernel) *t_kernel += usecond();
   assert(mat.dimension(0) == Nmom);
   assert(mat.dimension(1) == Ngamma);
   assert(mat.dimension(2) == Nt);
@@ -256,9 +267,9 @@ void A2Autils<FImpl>::MesonField(Eigen::Tensor<ComplexD,5> &mat,
   // Vector size is 7 x 16 x 32 x 16 x 16 x sizeof(complex) = 2MB - 60MB depending on volume
   // Healthy size that should suffice
   ////////////////////////////////////////////////////////////////////
-
+  if (t_gsum) *t_gsum = -usecond();
   grid->GlobalSumVector(&mat(0,0,0,0,0),Nmom*Ngamma*Nt*Lblock*Rblock);
-
+  if (t_gsum) *t_gsum += usecond();
 }
 
 
@@ -614,6 +625,189 @@ void A2Autils<FImpl>::PionFieldVV(Eigen::Tensor<ComplexD,3> &mat,
   PionFieldXX(mat,vi,vj,orthogdim,nog5);
 }
 
+// "A-slash" field w_i(x)^dag * i * A_mu * gamma_mu * v_j(x)
+//
+// With:
+//
+// B_0 = A_0 + i A_1
+// B_1 = A_2 + i A_3
+// 
+// then in spin space
+// 
+//                 ( 0          0          -conj(B_1) -B_0 )
+// i * A_mu g_mu = ( 0          0          -conj(B_0)  B_1 )
+//                 ( B_1        B_0        0          0    )
+//                 ( conj(B_0)  -conj(B_1) 0          0    )
+template <class FImpl>
+template <typename TensorType>
+void A2Autils<FImpl>::AslashField(TensorType &mat, 
+          const FermionField *lhs_wi,
+          const FermionField *rhs_vj,
+          const std::vector<ComplexField> &emB0,
+          const std::vector<ComplexField> &emB1,
+          int orthogdim, double *t_kernel, double *t_gsum) 
+{
+    typedef typename FermionField::vector_object vobj;
+    typedef typename vobj::scalar_object         sobj;
+    typedef typename vobj::scalar_type           scalar_type;
+    typedef typename vobj::vector_type           vector_type;
+
+    typedef iSpinMatrix<vector_type> SpinMatrix_v;
+    typedef iSpinMatrix<scalar_type> SpinMatrix_s;
+    typedef iSinglet<vector_type>    Singlet_v;
+    typedef iSinglet<scalar_type>    Singlet_s;
+    
+    int Lblock = mat.dimension(3); 
+    int Rblock = mat.dimension(4);
+
+    GridBase *grid = lhs_wi[0]._grid;
+    
+    const int    Nd = grid->_ndimension;
+    const int Nsimd = grid->Nsimd();
+
+    int Nt  = grid->GlobalDimensions()[orthogdim];
+    int Nem = emB0.size();
+    assert(emB1.size() == Nem);
+
+    int fd=grid->_fdimensions[orthogdim];
+    int ld=grid->_ldimensions[orthogdim];
+    int rd=grid->_rdimensions[orthogdim];
+
+    // will locally sum vectors first
+    // sum across these down to scalars
+    // splitting the SIMD
+    int MFrvol = rd*Lblock*Rblock*Nem;
+    int MFlvol = ld*Lblock*Rblock*Nem;
+
+    Vector<vector_type> lvSum(MFrvol);
+    parallel_for (int r = 0; r < MFrvol; r++)
+    {
+        lvSum[r] = zero;
+    }
+
+    Vector<scalar_type> lsSum(MFlvol);             
+    parallel_for (int r = 0; r < MFlvol; r++)
+    {
+        lsSum[r] = scalar_type(0.0);
+    }
+
+    int e1=    grid->_slice_nblock[orthogdim];
+    int e2=    grid->_slice_block [orthogdim];
+    int stride=grid->_slice_stride[orthogdim];
+
+    // Nested parallelism would be ok
+    // Wasting cores here. Test case r
+    if (t_kernel) *t_kernel = -usecond();
+    parallel_for(int r=0;r<rd;r++)
+    {
+        int so=r*grid->_ostride[orthogdim]; // base offset for start of plane 
+
+        for(int n=0;n<e1;n++)
+        for(int b=0;b<e2;b++)
+        {
+            int ss= so+n*stride+b;
+
+            for(int i=0;i<Lblock;i++)
+            {
+                auto left = conjugate(lhs_wi[i]._odata[ss]);
+
+                for(int j=0;j<Rblock;j++)
+                {
+                    SpinMatrix_v vv;
+                    auto right = rhs_vj[j]._odata[ss];
+
+                    for(int s1=0;s1<Ns;s1++)
+                    for(int s2=0;s2<Ns;s2++)
+                    {
+                        vv()(s1,s2)() = left()(s2)(0) * right()(s1)(0)
+                                        + left()(s2)(1) * right()(s1)(1)
+                                        + left()(s2)(2) * right()(s1)(2);
+                    }
+                    
+                    // After getting the sitewise product do the mom phase loop
+                    int base = Nem*i+Nem*Lblock*j+Nem*Lblock*Rblock*r;
+
+                    for ( int m=0;m<Nem;m++)
+                    {
+                        int idx  = m+base;
+                        auto b0  = emB0[m]._odata[ss];
+                        auto b1  = emB1[m]._odata[ss];
+                        auto cb0 = conjugate(b0);
+                        auto cb1 = conjugate(b1);
+
+                        lvSum[idx] += - vv()(3,0)()*b0()()()  - vv()(2,0)()*cb1()()()
+                                      + vv()(3,1)()*b1()()()  - vv()(2,1)()*cb0()()()
+                                      + vv()(0,2)()*b1()()()  + vv()(1,2)()*b0()()()
+                                      + vv()(0,3)()*cb0()()() - vv()(1,3)()*cb1()()();
+                    }
+                }
+            }
+        }
+    }
+
+    // Sum across simd lanes in the plane, breaking out orthog dir.
+    parallel_for(int rt=0;rt<rd;rt++)
+    {
+        std::vector<int> icoor(Nd);
+        std::vector<scalar_type> extracted(Nsimd);               
+
+        for(int i=0;i<Lblock;i++)
+        for(int j=0;j<Rblock;j++)
+        for(int m=0;m<Nem;m++)
+        {
+
+            int ij_rdx = m+Nem*i+Nem*Lblock*j+Nem*Lblock*Rblock*rt;
+
+            extract<vector_type,scalar_type>(lvSum[ij_rdx],extracted);
+            for(int idx=0;idx<Nsimd;idx++)
+            {
+                grid->iCoorFromIindex(icoor,idx);
+
+                int ldx    = rt+icoor[orthogdim]*rd;
+                int ij_ldx = m+Nem*i+Nem*Lblock*j+Nem*Lblock*Rblock*ldx;
+
+                lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx];
+            }
+        }
+    }
+    if (t_kernel) *t_kernel += usecond();
+
+    // ld loop and local only??
+    int pd = grid->_processors[orthogdim];
+    int pc = grid->_processor_coor[orthogdim];
+    parallel_for_nest2(int lt=0;lt<ld;lt++)
+    {
+        for(int pt=0;pt<pd;pt++)
+        {
+            int t = lt + pt*ld;
+            if (pt == pc)
+            {
+                for(int i=0;i<Lblock;i++)
+                for(int j=0;j<Rblock;j++)
+                for(int m=0;m<Nem;m++)
+                {
+                    int ij_dx = m+Nem*i + Nem*Lblock * j + Nem*Lblock * Rblock * lt;
+
+                    mat(m,0,t,i,j) = lsSum[ij_dx];
+                }
+            } 
+            else 
+            { 
+                const scalar_type zz(0.0);
+
+                for(int i=0;i<Lblock;i++)
+                for(int j=0;j<Rblock;j++)
+                for(int m=0;m<Nem;m++)
+                {
+                    mat(m,0,t,i,j) = zz;
+                }
+            }
+        }
+    }
+    if (t_gsum) *t_gsum = -usecond();
+    grid->GlobalSumVector(&mat(0,0,0,0,0),Nem*Nt*Lblock*Rblock);
+    if (t_gsum) *t_gsum += usecond();
+}
 
 ////////////////////////////////////////////
 // Schematic thoughts about more generalised four quark insertion

Grid/util/Sha.h

@@ -38,7 +38,21 @@ public:
   {
     return ::crc32(0L,(unsigned char *)data,bytes);
   }
-  static inline std::vector<unsigned char> sha256(void *data,size_t bytes)
+  template <typename T>
+  static inline std::string sha256_string(const std::vector<T> &hash)
+  {
+    std::stringstream sha;
+    std::string       s;
+
+    for(unsigned int i = 0; i < hash.size(); i++) 
+    { 
+        sha << std::hex << static_cast<unsigned int>(hash[i]);
+    }
+    s = sha.str();
+
+    return s;
+  }
+  static inline std::vector<unsigned char> sha256(const void *data,size_t bytes)
   {
     std::vector<unsigned char> hash(SHA256_DIGEST_LENGTH);
     SHA256_CTX sha256;

Hadrons/A2AMatrix.hpp

@@ -7,6 +7,7 @@ Source file: Hadrons/A2AMatrix.hpp
 Copyright (C) 2015-2018
 
 Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 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
@@ -29,38 +30,128 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #define A2A_Matrix_hpp_
 
 #include <Hadrons/Global.hpp>
+#include <Hadrons/TimerArray.hpp>
+#include <Grid/Eigen/unsupported/CXX11/Tensor>
+
+#ifndef HADRONS_A2AM_NAME 
+#define HADRONS_A2AM_NAME "a2aMatrix"
+#endif
+
+#define HADRONS_A2AM_PARALLEL_IO
 
 BEGIN_HADRONS_NAMESPACE
 
-template <typename T, typename MetadataType>
+// general A2A matrix set based on Eigen tensors and Grid-allocated memory
+// Dimensions:
+//   0 - ext - external field (momentum, EM field, ...)
+//   1 - str - spin-color structure
+//   2 - t   - timeslice
+//   3 - i   - left  A2A mode index
+//   4 - j   - right A2A mode index
+template <typename T>
+using A2AMatrixSet = Eigen::TensorMap<Eigen::Tensor<T, 5, Eigen::RowMajor>>;
+
+/******************************************************************************
+ *                      Abstract class for A2A kernels                        *
+ ******************************************************************************/
+template <typename T, typename Field>
+class A2AKernel
+{
+public:
+    A2AKernel(void) = default;
+    virtual ~A2AKernel(void) = default;
+    virtual void operator()(A2AMatrixSet<T> &m, const Field *left, const Field *right,
+                          const unsigned int orthogDim, double &time) = 0;
+    virtual double flops(const unsigned int blockSizei, const unsigned int blockSizej) = 0;
+    virtual double bytes(const unsigned int blockSizei, const unsigned int blockSizej) = 0;
+};
+
+/******************************************************************************
+ *                  Class to handle A2A matrix block HDF5 I/O                 *
+ ******************************************************************************/
+template <typename T>
 class A2AMatrixIo
 {
 public:
+    // constructors
     A2AMatrixIo(void) = default;
     A2AMatrixIo(std::string filename, std::string dataname, 
                 const unsigned int nt, const unsigned int ni,
                 const unsigned int nj);
+    // destructor
     ~A2AMatrixIo(void) = default;
+    // file allocation
+    template <typename MetadataType>
     void initFile(const MetadataType &d, const unsigned int chunkSize);
+    // block I/O
     void saveBlock(const T *data, const unsigned int i, const unsigned int j,
                    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);
 private:
     std::string  filename_, dataname_;
     unsigned int nt_, ni_, nj_;
 };
 
-template <typename T, typename MetadataType>
-A2AMatrixIo<T, MetadataType>::A2AMatrixIo(std::string filename, 
-                                          std::string dataname, 
-                                          const unsigned int nt, 
-                                          const unsigned int ni,
-                                          const unsigned int nj)
+/******************************************************************************
+ *                  Wrapper for A2A matrix block computation                  *
+ ******************************************************************************/
+template <typename T, typename Field, typename MetadataType, typename TIo = T>
+class A2AMatrixBlockComputation
+{
+private:
+    struct IoHelper
+    {
+        A2AMatrixIo<TIo> io;
+        MetadataType     md;
+        unsigned int     e, s, i, j;
+    };
+    typedef std::function<std::string(const unsigned int, const unsigned int)>  FilenameFn;
+    typedef std::function<MetadataType(const unsigned int, const unsigned int)> MetadataFn;
+public:
+    // constructor
+    A2AMatrixBlockComputation(GridBase *grid,
+                              const unsigned int orthogDim,
+                              const unsigned int next,
+                              const unsigned int nstr,
+                              const unsigned int blockSize,
+                              const unsigned int cacheBlockSize,
+                              TimerArray *tArray = nullptr);
+    // execution
+    void execute(const std::vector<Field> &left, 
+                 const std::vector<Field> &right,
+                 A2AKernel<T, Field> &kernel,
+                 const FilenameFn &ionameFn,
+                 const FilenameFn &filenameFn,
+                 const MetadataFn &metadataFn);
+private:
+    // I/O handler
+    void saveBlock(const A2AMatrixSet<TIo> &m, IoHelper &h);
+private:
+    TimerArray            *tArray_;
+    GridBase              *grid_;
+    unsigned int          orthogDim_, nt_, next_, nstr_, blockSize_, cacheBlockSize_;
+    Vector<T>             mCache_;
+    Vector<TIo>           mBuf_;
+    std::vector<IoHelper> nodeIo_;
+};
+
+/******************************************************************************
+ *                     A2AMatrixIo template implementation                    *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename T>
+A2AMatrixIo<T>::A2AMatrixIo(std::string filename, std::string dataname, 
+                            const unsigned int nt, const unsigned int ni,
+                            const unsigned int nj)
 : filename_(filename), dataname_(dataname)
 , nt_(nt), ni_(ni), nj_(nj)
 {}
 
-template <typename T, typename MetadataType>
-void A2AMatrixIo<T, MetadataType>::initFile(const MetadataType &d, const unsigned int chunkSize)
+// file allocation /////////////////////////////////////////////////////////////
+template <typename T>
+template <typename MetadataType>
+void A2AMatrixIo<T>::initFile(const MetadataType &d, const unsigned int chunkSize)
 {
 #ifdef HAVE_HDF5
     std::vector<hsize_t>    dim = {static_cast<hsize_t>(nt_), 
@@ -85,18 +176,19 @@ void A2AMatrixIo<T, MetadataType>::initFile(const MetadataType &d, const unsigne
     push(reader, dataname_);
     auto &group = reader.getGroup();
     plist.setChunk(chunk.size(), chunk.data());
-    dataset = group.createDataSet("data", Hdf5Type<T>::type(), dataspace, plist);
+    dataset = group.createDataSet(HADRONS_A2AM_NAME, Hdf5Type<T>::type(), dataspace, plist);
 #else
     HADRONS_ERROR(Implementation, "all-to-all matrix I/O needs HDF5 library");
 #endif
 }
 
-template <typename T, typename MetadataType>
-void A2AMatrixIo<T, MetadataType>::saveBlock(const T *data, 
-                                             const unsigned int i, 
-                                             const unsigned int j,
-                                             const unsigned int blockSizei,
-                                             const unsigned int blockSizej)
+// block I/O ///////////////////////////////////////////////////////////////////
+template <typename T>
+void A2AMatrixIo<T>::saveBlock(const T *data, 
+                               const unsigned int i, 
+                               const unsigned int j,
+                               const unsigned int blockSizei,
+                               const unsigned int blockSizej)
 {
 #ifdef HAVE_HDF5
     Hdf5Reader           reader(filename_);
@@ -111,7 +203,7 @@ void A2AMatrixIo<T, MetadataType>::saveBlock(const T *data,
 
     push(reader, dataname_);
     auto &group = reader.getGroup();
-    dataset     = group.openDataSet("data");
+    dataset     = group.openDataSet(HADRONS_A2AM_NAME);
     dataspace   = dataset.getSpace();
     dataspace.selectHyperslab(H5S_SELECT_SET, count.data(), offset.data(),
                               stride.data(), block.data());
@@ -121,6 +213,193 @@ void A2AMatrixIo<T, MetadataType>::saveBlock(const T *data,
 #endif
 }
 
+template <typename T>
+void A2AMatrixIo<T>::saveBlock(const A2AMatrixSet<T> &m,
+                               const unsigned int ext, const unsigned int str,
+                               const unsigned int i, const unsigned int j)
+{
+    unsigned int blockSizei = m.dimension(3);
+    unsigned int blockSizej = m.dimension(4);
+    unsigned int nstr       = m.dimension(1);
+    size_t       offset     = (ext*nstr + str)*nt_*blockSizei*blockSizej;
+
+    saveBlock(m.data() + offset, i, j, blockSizei, blockSizej);
+}
+
+/******************************************************************************
+ *               A2AMatrixBlockComputation template implementation            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename T, typename Field, typename MetadataType, typename TIo>
+A2AMatrixBlockComputation<T, Field, MetadataType, TIo>
+::A2AMatrixBlockComputation(GridBase *grid,
+                            const unsigned int orthogDim,
+                            const unsigned int next, 
+                            const unsigned int nstr,
+                            const unsigned int blockSize, 
+                            const unsigned int cacheBlockSize,
+                            TimerArray *tArray)
+: grid_(grid), nt_(grid->GlobalDimensions()[orthogDim]), orthogDim_(orthogDim)
+, next_(next), nstr_(nstr), blockSize_(blockSize), cacheBlockSize_(cacheBlockSize)
+, tArray_(tArray)
+{
+    mCache_.resize(nt_*next_*nstr_*cacheBlockSize_*cacheBlockSize_);
+    mBuf_.resize(nt_*next_*nstr_*blockSize_*blockSize_);
+}
+
+#define START_TIMER(name) if (tArray_) tArray_->startTimer(name)
+#define STOP_TIMER(name)  if (tArray_) tArray_->stopTimer(name)
+#define GET_TIMER(name)   ((tArray_ != nullptr) ? tArray_->getDTimer(name) : 0.)
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename T, typename Field, typename MetadataType, typename TIo>
+void A2AMatrixBlockComputation<T, Field, MetadataType, TIo>
+::execute(const std::vector<Field> &left, const std::vector<Field> &right,
+          A2AKernel<T, Field> &kernel, const FilenameFn &ionameFn,
+          const FilenameFn &filenameFn, const MetadataFn &metadataFn)
+{
+    //////////////////////////////////////////////////////////////////////////
+    // i,j   is first  loop over blockSize_ factors
+    // ii,jj is second loop over cacheBlockSize_ factors for high perf contractions
+    // iii,jjj are loops within cacheBlock
+    // Total index is sum of these  i+ii+iii etc...
+    //////////////////////////////////////////////////////////////////////////
+    int    N_i = left.size();
+    int    N_j = right.size();
+    double flops, bytes, t_kernel;
+    double nodes = grid_->NodeCount();
+    
+    int NBlock_i = N_i/blockSize_ + (((N_i % blockSize_) != 0) ? 1 : 0);
+    int NBlock_j = N_j/blockSize_ + (((N_j % blockSize_) != 0) ? 1 : 0);
+
+    for(int i=0;i<N_i;i+=blockSize_)
+    for(int j=0;j<N_j;j+=blockSize_)
+    {
+        // Get the W and V vectors for this block^2 set of terms
+        int N_ii = MIN(N_i-i,blockSize_);
+        int N_jj = MIN(N_j-j,blockSize_);
+        A2AMatrixSet<TIo> mBlock(mBuf_.data(), next_, nstr_, nt_, N_ii, N_jj);
+
+        LOG(Message) << "All-to-all matrix block " 
+                     << j/blockSize_ + NBlock_j*i/blockSize_ + 1 
+                     << "/" << NBlock_i*NBlock_j << " [" << i <<" .. " 
+                     << i+N_ii-1 << ", " << j <<" .. " << j+N_jj-1 << "]" 
+                     << std::endl;
+        // Series of cache blocked chunks of the contractions within this block
+        flops    = 0.0;
+        bytes    = 0.0;
+        t_kernel = 0.0;
+        for(int ii=0;ii<N_ii;ii+=cacheBlockSize_)
+        for(int jj=0;jj<N_jj;jj+=cacheBlockSize_)
+        {
+            double t;
+            int N_iii = MIN(N_ii-ii,cacheBlockSize_);
+            int N_jjj = MIN(N_jj-jj,cacheBlockSize_);
+            A2AMatrixSet<T> mCacheBlock(mCache_.data(), next_, nstr_, nt_, N_iii, N_jjj);
+
+            START_TIMER("kernel");
+            kernel(mCacheBlock, &left[i+ii], &right[j+jj], orthogDim_, t);
+            STOP_TIMER("kernel");
+            t_kernel += t;
+            flops    += kernel.flops(N_iii, N_jjj);
+            bytes    += kernel.bytes(N_iii, N_jjj);
+
+            START_TIMER("cache copy");
+            parallel_for_nest5(int e =0;e<next_;e++)
+            for(int s =0;s< nstr_;s++)
+            for(int t =0;t< nt_;t++)
+            for(int iii=0;iii< N_iii;iii++)
+            for(int jjj=0;jjj< N_jjj;jjj++)
+            {
+                mBlock(e,s,t,ii+iii,jj+jjj) = mCacheBlock(e,s,t,iii,jjj);
+            }
+            STOP_TIMER("cache copy");
+        }
+
+        // perf
+        LOG(Message) << "Kernel perf " << flops/t_kernel/1.0e3/nodes 
+                     << " Gflop/s/node " << std::endl;
+        LOG(Message) << "Kernel perf " << bytes/t_kernel*1.0e6/1024/1024/1024/nodes 
+                     << " GB/s/node "  << std::endl;
+
+        // IO
+        double       blockSize, ioTime;
+        unsigned int myRank = grid_->ThisRank(), nRank  = grid_->RankCount();
+    
+        LOG(Message) << "Writing block to disk" << std::endl;
+        ioTime = -GET_TIMER("IO: write block");
+        START_TIMER("IO: total");
+        makeFileDir(filenameFn(0, 0), grid_);
+#ifdef HADRONS_A2AM_PARALLEL_IO
+        grid_->Barrier();
+        // make task list for current node
+        nodeIo_.clear();
+        for(int f = myRank; f < next_*nstr_; f += nRank)
+        {
+            IoHelper h;
+
+            h.i  = i;
+            h.j  = j;
+            h.e  = f/nstr_;
+            h.s  = f % nstr_;
+            h.io = A2AMatrixIo<TIo>(filenameFn(h.e, h.s), 
+                                    ionameFn(h.e, h.s), nt_, N_i, N_j);
+            h.md = metadataFn(h.e, h.s);
+            nodeIo_.push_back(h);
+        }
+        // parallel IO
+        for (auto &h: nodeIo_)
+        {
+            saveBlock(mBlock, h);
+        }
+        grid_->Barrier();
+#else
+        // serial IO, for testing purposes only
+        for(int e = 0; e < next_; e++)
+        for(int s = 0; s < nstr_; s++)
+        {
+            IoHelper h;
+
+            h.i  = i;
+            h.j  = j;
+            h.e  = e;
+            h.s  = s;
+            h.io = A2AMatrixIo<TIo>(filenameFn(h.e, h.s), 
+                                    ionameFn(h.e, h.s), nt_, N_i, N_j);
+            h.md = metadataFn(h.e, h.s);
+            saveBlock(mfBlock, h);
+        }
+#endif
+        STOP_TIMER("IO: total");
+        blockSize  = static_cast<double>(next_*nstr_*nt_*N_ii*N_jj*sizeof(TIo));
+        ioTime    += GET_TIMER("IO: write block");
+        LOG(Message) << "HDF5 IO done " << sizeString(blockSize) << " in "
+                     << ioTime  << " us (" 
+                     << blockSize/ioTime*1.0e6/1024/1024
+                     << " MB/s)" << std::endl;
+    }
+}
+
+// I/O handler /////////////////////////////////////////////////////////////////
+template <typename T, typename Field, typename MetadataType, typename TIo>
+void A2AMatrixBlockComputation<T, Field, MetadataType, TIo>
+::saveBlock(const A2AMatrixSet<TIo> &m, IoHelper &h)
+{
+    if ((h.i == 0) and (h.j == 0))
+    {
+        START_TIMER("IO: file creation");
+        h.io.initFile(h.md, blockSize_);
+        STOP_TIMER("IO: file creation");
+    }
+    START_TIMER("IO: write block");
+    h.io.saveBlock(m, h.e, h.s, h.i, h.j);
+    STOP_TIMER("IO: write block");
+}
+
+#undef START_TIMER
+#undef STOP_TIMER
+#undef GET_TIMER
+
 END_HADRONS_NAMESPACE
 
 #endif // A2A_Matrix_hpp_

Hadrons/DiskVector.hpp

@@ -34,6 +34,12 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <ftw.h>
 #include <unistd.h>
 
+#ifdef DV_DEBUG
+#define DV_DEBUG_MSG(dv, stream) LOG(Debug) << "diskvector " << (dv) << ": " << stream << std::endl
+#else
+#define DV_DEBUG_MSG(dv, stream)
+#endif
+
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
@@ -56,9 +62,7 @@ public:
         // write to disk and cache
         T &operator=(const T &obj) const
         {
-#ifdef DV_DEBUG
-            LOG(Debug) << "diskvector " << &master_ << ": writing to " << i_ << std::endl;
-#endif
+            DV_DEBUG_MSG(&master_, "writing to " << i_);
             master_.cacheInsert(i_, obj);
             master_.save(master_.filename(i_), obj);
             
@@ -82,6 +86,8 @@ public:
     virtual ~DiskVectorBase(void);
     const T & operator[](const unsigned int i) const;
     RwAccessHelper operator[](const unsigned int i);
+    double hitRatio(void) const;
+    void resetStat(void);
 private:
     virtual void load(T &obj, const std::string filename) const = 0;
     virtual void save(const std::string filename, const T &obj) const = 0;
@@ -93,6 +99,7 @@ private:
 private:
     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
@@ -115,6 +122,7 @@ private:
 
         read(reader, basename(filename), obj);
     }
+
     virtual void save(const std::string filename, const T &obj) const
     {
         Writer writer(filename);
@@ -123,13 +131,76 @@ private:
     }
 };
 
+/******************************************************************************
+ *                      Specialisation for Eigen matrices                     *
+ ******************************************************************************/
+template <typename T>
+using EigenDiskVectorMat = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
+
+template <typename T>
+class EigenDiskVector: public DiskVectorBase<EigenDiskVectorMat<T>>
+{
+public:
+    using DiskVectorBase<EigenDiskVectorMat<T>>::DiskVectorBase;
+    typedef EigenDiskVectorMat<T> Matrix;
+public:
+    T operator()(const unsigned int i, const Eigen::Index j,
+                 const Eigen::Index k) const
+    {
+        return (*this)[i](j, k);
+    }
+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;
+
+        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));
+        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)
+        {
+            HADRONS_ERROR(Io, "checksum failed")
+        }
+    }
+
+    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;
+        
+        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));
+        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");
+    }
+};
+
 /******************************************************************************
  *                       DiskVectorBase implementation                         *
  ******************************************************************************/
-#ifdef DV_DEBUG
-#define DV_DEBUG_MSG(stream) LOG(Debug) << "diskvector " << this << ": " << stream << std::endl
-#endif
-
 template <typename T>
 DiskVectorBase<T>::DiskVectorBase(const std::string dirname, 
                                   const unsigned int size,
@@ -160,28 +231,29 @@ 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 &loads  = *loadsPtr_;
 
-    DV_DEBUG_MSG("accessing " << i << " (RO)");
+    DV_DEBUG_MSG(this, "accessing " << i << " (RO)");
 
     if (i >= size_)
     {
         HADRONS_ERROR(Size, "index out of range");
     }
-
+    const_cast<double &>(access_)++;
     if (cache.find(i) == cache.end())
     {
         // cache miss
-        DV_DEBUG_MSG("cache miss");
+        DV_DEBUG_MSG(this, "cache miss");
         fetch(i);
     }
     else
     {
-        DV_DEBUG_MSG("cache hit");
+        DV_DEBUG_MSG(this, "cache hit");
 
         auto pos = std::find(loads.begin(), loads.end(), i);
 
+        const_cast<double &>(hit_)++;
         loads.erase(pos);
         loads.push_back(i);
     }
@@ -193,7 +265,7 @@ const T & DiskVectorBase<T>::operator[](const unsigned int i) const
     {
         msg += std::to_string(p) + " ";
     }
-    DV_DEBUG_MSG("in cache: " << msg);
+    DV_DEBUG_MSG(this, "in cache: " << msg);
 #endif
 
     return cache.at(i);
@@ -202,7 +274,7 @@ const T & DiskVectorBase<T>::operator[](const unsigned int i) const
 template <typename T>
 typename DiskVectorBase<T>::RwAccessHelper DiskVectorBase<T>::operator[](const unsigned int i)
 {
-    DV_DEBUG_MSG("accessing " << i << " (RW)");
+    DV_DEBUG_MSG(this, "accessing " << i << " (RW)");
 
     if (i >= size_)
     {
@@ -212,6 +284,19 @@ typename DiskVectorBase<T>::RwAccessHelper DiskVectorBase<T>::operator[](const u
     return RwAccessHelper(*this, i);
 }
 
+template <typename T>
+double DiskVectorBase<T>::hitRatio(void) const
+{
+    return hit_/access_;
+}
+
+template <typename T>
+void DiskVectorBase<T>::resetStat(void)
+{
+    access_ = 0.;
+    hit_    = 0.;
+}
+
 template <typename T>
 std::string DiskVectorBase<T>::filename(const unsigned int i) const
 {
@@ -226,7 +311,7 @@ void DiskVectorBase<T>::evict(void) const
 
     if (cache.size() >= cacheSize_)
     {
-        DV_DEBUG_MSG("evicting " << loads.front());
+        DV_DEBUG_MSG(this, "evicting " << loads.front());
         cache.erase(loads.front());
         loads.pop_front();
     }
@@ -239,7 +324,7 @@ void DiskVectorBase<T>::fetch(const unsigned int i) const
     auto &loads = *loadsPtr_;
     struct stat s;
 
-    DV_DEBUG_MSG("loading " << i << " from disk");
+    DV_DEBUG_MSG(this, "loading " << i << " from disk");
 
     evict();
     if(stat(filename(i).c_str(), &s) != 0)
@@ -267,7 +352,7 @@ void DiskVectorBase<T>::cacheInsert(const unsigned int i, const T &obj) const
     {
         msg += std::to_string(p) + " ";
     }
-    DV_DEBUG_MSG("in cache: " << msg);
+    DV_DEBUG_MSG(this, "in cache: " << msg);
 #endif
 }
 

Hadrons/EigenPack.hpp

@@ -39,12 +39,12 @@ BEGIN_HADRONS_NAMESPACE
 #define HADRONS_DEFAULT_LANCZOS_NBASIS 60
 #endif
 
-#define HADRONS_DUMP_EP_METADATA \
+#define HADRONS_DUMP_EP_METADATA(record) \
 LOG(Message) << "Eigenpack metadata:" << std::endl;\
 LOG(Message) << "* operator" << std::endl;\
-LOG(Message) << record.operatorXml << std::endl;\
+LOG(Message) << (record).operatorXml << std::endl;\
 LOG(Message) << "* solver" << std::endl;\
-LOG(Message) << record.solverXml << std::endl;
+LOG(Message) << (record).solverXml << std::endl;
 
 struct PackRecord
 {
@@ -59,66 +59,9 @@ struct VecRecord: Serializable
     VecRecord(void): index(0), eval(0.) {}
 };
 
-template <typename F>
-class EigenPack
+namespace EigenPackIo
 {
-public:
-    typedef F Field;
-public:
-    std::vector<RealD> eval;
-    std::vector<F>     evec;
-    PackRecord         record;
-public:
-    EigenPack(void)          = default;
-    virtual ~EigenPack(void) = default;
-
-    EigenPack(const size_t size, GridBase *grid)
-    {
-        resize(size, grid);
-    }
-
-    void resize(const size_t size, GridBase *grid)
-    {
-        eval.resize(size);
-        evec.resize(size, grid);
-    }
-
-    virtual void read(const std::string fileStem, const bool multiFile, const int traj = -1)
-    {
-        if (multiFile)
-        {
-            for(int k = 0; k < evec.size(); ++k)
-            {
-                basicReadSingle(evec[k], eval[k], evecFilename(fileStem, k, traj), k);
-                if (k == 0)
-                {
-                    HADRONS_DUMP_EP_METADATA;
-                }
-            }
-        }
-        else
-        {
-            basicRead(evec, eval, evecFilename(fileStem, -1, traj), evec.size());
-            HADRONS_DUMP_EP_METADATA;
-        }
-    }
-
-    virtual void write(const std::string fileStem, const bool multiFile, const int traj = -1)
-    {
-        if (multiFile)
-        {
-            for(int k = 0; k < evec.size(); ++k)
-            {
-                basicWriteSingle(evecFilename(fileStem, k, traj), evec[k], eval[k], k);
-            }
-        }
-        else
-        {
-            basicWrite(evecFilename(fileStem, -1, traj), evec, eval, evec.size());
-        }
-    }
-
-    static void readHeader(PackRecord &record, ScidacReader &binReader)
+    inline void readHeader(PackRecord &record, ScidacReader &binReader)
     {
         std::string recordXml;
 
@@ -130,13 +73,75 @@ public:
         xmlReader.readCurrentSubtree(record.solverXml);
     }
 
-    template <typename T>
-    static void readElement(T &evec, VecRecord &vecRecord, ScidacReader &binReader)
+    template <typename T, typename TIo = T>
+    void readElement(T &evec, RealD &eval, const unsigned int index,
+                     ScidacReader &binReader, TIo *ioBuf = nullptr)
     {
-        binReader.readScidacFieldRecord(evec, vecRecord);
+        VecRecord vecRecord;
+
+        LOG(Message) << "Reading eigenvector " << index << std::endl;
+        if (ioBuf == nullptr)
+        {
+            binReader.readScidacFieldRecord(evec, vecRecord);
+        }
+        else
+        {
+            binReader.readScidacFieldRecord(*ioBuf, vecRecord);
+            precisionChange(evec, *ioBuf);
+        }
+        if (vecRecord.index != index)
+        {
+            HADRONS_ERROR(Io, "Eigenvector " + std::to_string(index) + " has a"
+                            + " wrong index (expected " + std::to_string(vecRecord.index) 
+                            + ")");
+        }
+        eval = vecRecord.eval;
     }
 
-    static void writeHeader(ScidacWriter &binWriter, PackRecord &record)
+    template <typename T, typename TIo = T>
+    static void readPack(std::vector<T> &evec, std::vector<RealD> &eval,
+                         PackRecord &record, const std::string filename, 
+                         const unsigned int size, bool multiFile, 
+                         GridBase *gridIo = nullptr)
+    {
+        std::unique_ptr<TIo> ioBuf{nullptr};
+        ScidacReader         binReader;
+
+        if (typeHash<T>() != typeHash<TIo>())
+        {
+            if (gridIo == nullptr)
+            {
+                HADRONS_ERROR(Definition, 
+                              "I/O type different from vector type but null I/O grid passed");
+            }
+            ioBuf.reset(new TIo(gridIo));
+        }
+        if (multiFile)
+        {
+            std::string fullFilename;
+
+            for(int k = 0; k < size; ++k) 
+            {
+                fullFilename = filename + "/v" + std::to_string(k) + ".bin";
+                binReader.open(fullFilename);
+                readHeader(record, binReader);
+                readElement(evec[k], eval[k], k, binReader, ioBuf.get());
+                binReader.close();
+            }
+        }
+        else
+        {
+            binReader.open(filename);
+            readHeader(record, binReader);
+            for(int k = 0; k < size; ++k) 
+            {
+                readElement(evec[k], eval[k], k, binReader, ioBuf.get());
+            }
+            binReader.close();
+        }
+    }
+
+    inline void writeHeader(ScidacWriter &binWriter, PackRecord &record)
     {
         XmlWriter xmlWriter("", "eigenPackPar");
 
@@ -145,165 +150,217 @@ public:
         binWriter.writeLimeObject(1, 1, xmlWriter, "parameters", SCIDAC_FILE_XML);
     }
 
-    template <typename T>
-    static void writeElement(ScidacWriter &binWriter, T &evec, VecRecord &vecRecord)
+    template <typename T, typename TIo = T>
+    void writeElement(ScidacWriter &binWriter, T &evec, RealD &eval, 
+                      const unsigned int index, TIo *ioBuf, 
+                      T *testBuf = nullptr)
     {
-        binWriter.writeScidacFieldRecord(evec, vecRecord, DEFAULT_ASCII_PREC);
-    }
-protected:
-    std::string evecFilename(const std::string stem, const int vec, const int traj)
-    {
-        std::string t = (traj < 0) ? "" : ("." + std::to_string(traj));
+        VecRecord vecRecord;
 
-        if (vec == -1)
+        LOG(Message) << "Writing eigenvector " << index << std::endl;
+        vecRecord.eval  = eval;
+        vecRecord.index = index;
+        if ((ioBuf == nullptr) || (testBuf == nullptr))
         {
-            return stem + t + ".bin";
+            binWriter.writeScidacFieldRecord(evec, vecRecord, DEFAULT_ASCII_PREC);
         }
         else
         {
-            return stem + t + "/v" + std::to_string(vec) + ".bin";
-        }
+            precisionChange(*ioBuf, evec);
+            precisionChange(*testBuf, *ioBuf);
+            *testBuf -= evec;
+            LOG(Message) << "Precision diff norm^2 " << norm2(*testBuf) << std::endl;
+            binWriter.writeScidacFieldRecord(*ioBuf, vecRecord, DEFAULT_ASCII_PREC);
+        }   
     }
-
-    template <typename T>
-    void basicRead(std::vector<T> &evec, std::vector<RealD> &eval,
-                   const std::string filename, const unsigned int size)
+    
+    template <typename T, typename TIo = T>
+    static void writePack(const std::string filename, std::vector<T> &evec, 
+                          std::vector<RealD> &eval, PackRecord &record, 
+                          const unsigned int size, bool multiFile, 
+                          GridBase *gridIo = nullptr)
     {
-        ScidacReader binReader;
+        GridBase             *grid = evec[0]._grid;
+        std::unique_ptr<TIo> ioBuf{nullptr}; 
+        std::unique_ptr<T>   testBuf{nullptr};
+        ScidacWriter         binWriter(grid->IsBoss());
 
-        binReader.open(filename);
-        readHeader(record, binReader);
-        for(int k = 0; k < size; ++k) 
+        if (typeHash<T>() != typeHash<TIo>())
         {
-            VecRecord vecRecord;
-
-            LOG(Message) << "Reading eigenvector " << k << std::endl;
-            readElement(evec[k], vecRecord, binReader);
-            if (vecRecord.index != k)
+            if (gridIo == nullptr)
             {
-                HADRONS_ERROR(Io, "Eigenvector " + std::to_string(k) + " has a"
-                              + " wrong index (expected " + std::to_string(vecRecord.index) 
-                              + ") in file '" + filename + "'");
+                HADRONS_ERROR(Definition, 
+                              "I/O type different from vector type but null I/O grid passed");
             }
-            eval[k] = vecRecord.eval;
+            ioBuf.reset(new TIo(gridIo));
+            testBuf.reset(new T(grid));
         }
-        binReader.close();
-    }
-
-    template <typename T>
-    void basicReadSingle(T &evec, RealD &eval, const std::string filename, 
-                         const unsigned int index)
-    {
-        ScidacReader binReader;
-        VecRecord    vecRecord;
-
-        binReader.open(filename);
-        readHeader(record, binReader);
-        LOG(Message) << "Reading eigenvector " << index << std::endl;
-        readElement(evec, vecRecord, binReader);
-        if (vecRecord.index != index)
+        if (multiFile)
         {
-            HADRONS_ERROR(Io, "Eigenvector " + std::to_string(index) + " has a"
-                          + " wrong index (expected " + std::to_string(vecRecord.index) 
-                          + ") in file '" + filename + "'");
+            std::string fullFilename;
+
+            for(int k = 0; k < size; ++k) 
+            {
+                fullFilename = filename + "/v" + std::to_string(k) + ".bin";
+
+                makeFileDir(fullFilename, grid);
+                binWriter.open(fullFilename);
+                writeHeader(binWriter, record);
+                writeElement(binWriter, evec[k], eval[k], k, ioBuf.get(), testBuf.get());
+                binWriter.close();
+            }
         }
-        eval = vecRecord.eval;
-        binReader.close();
-    }
-
-    template <typename T>
-    void basicWrite(const std::string filename, std::vector<T> &evec, 
-                    const std::vector<RealD> &eval, const unsigned int size)
-    {
-        ScidacWriter binWriter(evec[0]._grid->IsBoss());
-
-        makeFileDir(filename, evec[0]._grid);
-        binWriter.open(filename);
-        writeHeader(binWriter, record);
-        for(int k = 0; k < size; ++k) 
+        else
         {
-            VecRecord vecRecord;
-
-            vecRecord.index = k;
-            vecRecord.eval  = eval[k];
-            LOG(Message) << "Writing eigenvector " << k << std::endl;
-            writeElement(binWriter, evec[k], vecRecord);
+            makeFileDir(filename, grid);
+            binWriter.open(filename);
+            writeHeader(binWriter, record);
+            for(int k = 0; k < size; ++k) 
+            {
+                writeElement(binWriter, evec[k], eval[k], k, ioBuf.get(), testBuf.get());
+            }
+            binWriter.close();
         }
-        binWriter.close();
     }
+}
 
-    template <typename T>
-    void basicWriteSingle(const std::string filename, T &evec, 
-                          const RealD eval, const unsigned int index)
+template <typename F>
+class BaseEigenPack
+{
+public:
+    typedef F Field;
+public:
+    std::vector<RealD> eval;
+    std::vector<F>     evec;
+    PackRecord         record;
+public:
+    BaseEigenPack(void)          = default;
+    BaseEigenPack(const size_t size, GridBase *grid)
     {
-        ScidacWriter binWriter(evec._grid->IsBoss());
-        VecRecord    vecRecord;
-
-        makeFileDir(filename, evec._grid);
-        binWriter.open(filename);
-        writeHeader(binWriter, record);
-        vecRecord.index = index;
-        vecRecord.eval  = eval;
-        LOG(Message) << "Writing eigenvector " << index << std::endl;
-        writeElement(binWriter, evec, vecRecord);
-        binWriter.close();
+        resize(size, grid);
+    }
+    virtual ~BaseEigenPack(void) = default;
+    void resize(const size_t size, GridBase *grid)
+    {
+        eval.resize(size);
+        evec.resize(size, grid);
     }
 };
 
-template <typename FineF, typename CoarseF>
-class CoarseEigenPack: public EigenPack<FineF>
+template <typename F, typename FIo = F>
+class EigenPack: public BaseEigenPack<F>
 {
 public:
-    typedef CoarseF CoarseField;
+    typedef F   Field;
+    typedef FIo FieldIo;
 public:
-    std::vector<RealD>   evalCoarse;
+    EigenPack(void)          = default;
+    virtual ~EigenPack(void) = default;
+
+    EigenPack(const size_t size, GridBase *grid, GridBase *gridIo = nullptr)
+    : BaseEigenPack<F>(size, grid)
+    {
+        if (typeHash<F>() != typeHash<FIo>())
+        {
+            if (gridIo == nullptr)
+            {
+                HADRONS_ERROR(Definition, 
+                              "I/O type different from vector type but null I/O grid passed");
+            }
+        }
+        gridIo_ = gridIo;
+    }
+
+    virtual void read(const std::string fileStem, const bool multiFile, const int traj = -1)
+    {
+        EigenPackIo::readPack<F, FIo>(this->evec, this->eval, this->record, 
+                                      evecFilename(fileStem, traj, multiFile), 
+                                      this->evec.size(), multiFile, gridIo_);
+        HADRONS_DUMP_EP_METADATA(this->record);
+    }
+
+    virtual void write(const std::string fileStem, const bool multiFile, const int traj = -1)
+    {
+        EigenPackIo::writePack<F, FIo>(evecFilename(fileStem, traj, multiFile), 
+                                       this->evec, this->eval, this->record, 
+                                       this->evec.size(), multiFile, gridIo_);
+    }
+protected:
+    std::string evecFilename(const std::string stem, const int traj, const bool multiFile)
+    {
+        std::string t = (traj < 0) ? "" : ("." + std::to_string(traj));
+
+        if (multiFile)
+        {
+            return stem + t;
+        }
+        else
+        {
+            return stem + t + ".bin";
+        }
+    }
+protected:
+    GridBase *gridIo_;
+};
+
+template <typename FineF, typename CoarseF, 
+          typename FineFIo = FineF, typename CoarseFIo = CoarseF>
+class CoarseEigenPack: public EigenPack<FineF, FineFIo>
+{
+public:
+    typedef CoarseF CoarseField;         
     std::vector<CoarseF> evecCoarse;
+    std::vector<RealD>   evalCoarse;
 public:
     CoarseEigenPack(void)          = default;
     virtual ~CoarseEigenPack(void) = default;
 
     CoarseEigenPack(const size_t sizeFine, const size_t sizeCoarse, 
-                    GridBase *gridFine, GridBase *gridCoarse)
+                    GridBase *gridFine, GridBase *gridCoarse,
+                    GridBase *gridFineIo = nullptr, 
+                    GridBase *gridCoarseIo = nullptr)
     {
+        if (typeHash<FineF>() != typeHash<FineFIo>())
+        {
+            if (gridFineIo == nullptr)
+            {
+                HADRONS_ERROR(Definition, 
+                              "Fine I/O type different from vector type but null fine I/O grid passed");
+            }
+        }
+        if (typeHash<CoarseF>() != typeHash<CoarseFIo>())
+        {
+            if (gridCoarseIo == nullptr)
+            {
+                HADRONS_ERROR(Definition, 
+                              "Coarse I/O type different from vector type but null coarse I/O grid passed");
+            }
+        }
+        this->gridIo_ = gridFineIo;
+        gridCoarseIo_ = gridCoarseIo;
         resize(sizeFine, sizeCoarse, gridFine, gridCoarse);
     }
 
     void resize(const size_t sizeFine, const size_t sizeCoarse, 
                 GridBase *gridFine, GridBase *gridCoarse)
     {
-        EigenPack<FineF>::resize(sizeFine, gridFine);
+        EigenPack<FineF, FineFIo>::resize(sizeFine, gridFine);
         evalCoarse.resize(sizeCoarse);
         evecCoarse.resize(sizeCoarse, gridCoarse);
     }
 
     void readFine(const std::string fileStem, const bool multiFile, const int traj = -1)
     {
-        if (multiFile)
-        {
-            for(int k = 0; k < this->evec.size(); ++k)
-            {
-                this->basicReadSingle(this->evec[k], this->eval[k], this->evecFilename(fileStem + "_fine", k, traj), k);
-            }
-        }
-        else
-        {
-            this->basicRead(this->evec, this->eval, this->evecFilename(fileStem + "_fine", -1, traj), this->evec.size());
-        }
+        EigenPack<FineF, FineFIo>::read(fileStem + "_fine", multiFile, traj);
     }
 
     void readCoarse(const std::string fileStem, const bool multiFile, const int traj = -1)
     {
-        if (multiFile)
-        {
-            for(int k = 0; k < evecCoarse.size(); ++k)
-            {
-                this->basicReadSingle(evecCoarse[k], evalCoarse[k], this->evecFilename(fileStem + "_coarse", k, traj), k);
-            }
-        }
-        else
-        {
-            this->basicRead(evecCoarse, evalCoarse, this->evecFilename(fileStem + "_coarse", -1, traj), evecCoarse.size());
-        }
+        PackRecord dummy;
+
+        EigenPackIo::readPack<CoarseF, CoarseFIo>(evecCoarse, evalCoarse, dummy, 
+                              this->evecFilename(fileStem + "_coarse", traj, multiFile), 
+                              evecCoarse.size(), multiFile, gridCoarseIo_);
     }
 
     virtual void read(const std::string fileStem, const bool multiFile, const int traj = -1)
@@ -314,32 +371,14 @@ public:
 
     void writeFine(const std::string fileStem, const bool multiFile, const int traj = -1)
     {
-        if (multiFile)
-        {
-            for(int k = 0; k < this->evec.size(); ++k)
-            {
-                this->basicWriteSingle(this->evecFilename(fileStem + "_fine", k, traj), this->evec[k], this->eval[k], k);
-            }
-        }
-        else
-        {
-            this->basicWrite(this->evecFilename(fileStem + "_fine", -1, traj), this->evec, this->eval, this->evec.size());
-        }
+        EigenPack<FineF, FineFIo>::write(fileStem + "_fine", multiFile, traj);
     }
 
     void writeCoarse(const std::string fileStem, const bool multiFile, const int traj = -1)
     {
-        if (multiFile)
-        {
-            for(int k = 0; k < evecCoarse.size(); ++k)
-            {
-                this->basicWriteSingle(this->evecFilename(fileStem + "_coarse", k, traj), evecCoarse[k], evalCoarse[k], k);
-            }
-        }
-        else
-        {
-            this->basicWrite(this->evecFilename(fileStem + "_coarse", -1, traj), evecCoarse, evalCoarse, evecCoarse.size());
-        }
+        EigenPackIo::writePack<CoarseF, CoarseFIo>(this->evecFilename(fileStem + "_coarse", traj, multiFile), 
+                                                   evecCoarse, evalCoarse, this->record, 
+                                                   evecCoarse.size(), multiFile, gridCoarseIo_);
     }
     
     virtual void write(const std::string fileStem, const bool multiFile, const int traj = -1)
@@ -347,16 +386,25 @@ public:
         writeFine(fileStem, multiFile, traj);
         writeCoarse(fileStem, multiFile, traj);
     }
+private:
+    GridBase *gridCoarseIo_;
 };
 
 template <typename FImpl>
-using FermionEigenPack = EigenPack<typename FImpl::FermionField>;
+using BaseFermionEigenPack = BaseEigenPack<typename FImpl::FermionField>;
 
-template <typename FImpl, int nBasis>
+template <typename FImpl, typename FImplIo = FImpl>
+using FermionEigenPack = EigenPack<typename FImpl::FermionField, typename FImplIo::FermionField>;
+
+template <typename FImpl, int nBasis, typename FImplIo = FImpl>
 using CoarseFermionEigenPack = CoarseEigenPack<
     typename FImpl::FermionField,
     typename LocalCoherenceLanczos<typename FImpl::SiteSpinor, 
                                    typename FImpl::SiteComplex, 
+                                   nBasis>::CoarseField,
+    typename FImplIo::FermionField,
+    typename LocalCoherenceLanczos<typename FImplIo::SiteSpinor, 
+                                   typename FImplIo::SiteComplex, 
                                    nBasis>::CoarseField>;
 
 #undef HADRONS_DUMP_EP_METADATA

Hadrons/Makefile.am

@@ -11,6 +11,7 @@ libHadrons_a_SOURCES = \
 	Exceptions.cc      \
     Global.cc          \
     Module.cc		   \
+	TimerArray.cc      \
 	VirtualMachine.cc
 libHadrons_adir = $(includedir)/Hadrons
 nobase_libHadrons_a_HEADERS = \
@@ -31,4 +32,5 @@ nobase_libHadrons_a_HEADERS = \
 	Modules.hpp               \
 	ModuleFactory.hpp         \
 	Solver.hpp                \
+	TimerArray.hpp            \
 	VirtualMachine.hpp

Hadrons/Module.cc

@@ -66,101 +66,6 @@ void ModuleBase::operator()(void)
     stopAllTimers();
 }
 
-// timers //////////////////////////////////////////////////////////////////////
-void ModuleBase::startTimer(const std::string &name)
-{
-    if (!name.empty())
-    {
-        timer_[name].Start();
-    }
-}
-
-GridTime ModuleBase::getTimer(const std::string &name)
-{
-    GridTime t;
-    
-    if (!name.empty())
-    {
-        try
-        {
-            bool running = timer_.at(name).isRunning();
-
-            if (running) stopTimer(name);
-            t = timer_.at(name).Elapsed();
-            if (running) startTimer(name);
-        }
-        catch (std::out_of_range &)
-        {
-            t = GridTime::zero();
-        }
-    }
-    else
-    {
-        t = GridTime::zero();
-    }
-
-    return t;
-}
-
-double ModuleBase::getDTimer(const std::string &name)
-{
-    return static_cast<double>(getTimer(name).count());
-}
-
-void ModuleBase::startCurrentTimer(const std::string &name)
-{
-    if (!name.empty())
-    {
-        stopCurrentTimer();
-        startTimer(name);
-        currentTimer_ = name;
-    }
-}
-
-void ModuleBase::stopTimer(const std::string &name)
-{
-    if (timer_.at(name).isRunning())
-    {
-        timer_.at(name).Stop();
-    }
-}
-
-void ModuleBase::stopCurrentTimer(void)
-{
-    if (!currentTimer_.empty())
-    {
-        stopTimer(currentTimer_);
-        currentTimer_ = "";
-    }
-}
-
-void ModuleBase::stopAllTimers(void)
-{
-    for (auto &t: timer_)
-    {
-        stopTimer(t.first);
-    }
-    currentTimer_ = "";
-}
-
-void ModuleBase::resetTimers(void)
-{
-    timer_.clear();
-    currentTimer_ = "";
-}
-
-std::map<std::string, GridTime> ModuleBase::getTimings(void)
-{
-    std::map<std::string, GridTime> timing;
-
-    for (auto &t: timer_)
-    {
-        timing[t.first] = t.second.Elapsed();
-    }
-
-    return timing;
-}
-
 std::string ModuleBase::makeSeedString(void)
 {
     std::string seed;

Hadrons/Module.hpp

@@ -30,6 +30,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #define Hadrons_Module_hpp_
 
 #include <Hadrons/Global.hpp>
+#include <Hadrons/TimerArray.hpp>
 #include <Hadrons/VirtualMachine.hpp>
 
 BEGIN_HADRONS_NAMESPACE
@@ -152,7 +153,7 @@ if (env().getGrid()->IsBoss() and !ioStem.empty())\
  *                            Module class                                    *
  ******************************************************************************/
 // base class
-class ModuleBase
+class ModuleBase: public TimerArray
 {
 public:
     // constructor
@@ -180,16 +181,6 @@ public:
     virtual void execute(void) = 0;
     // execution
     void operator()(void);
-    // timers
-    void                            startTimer(const std::string &name);
-    GridTime                        getTimer(const std::string &name);
-    double                          getDTimer(const std::string &name);
-    void                            startCurrentTimer(const std::string &name);
-    void                            stopTimer(const std::string &name);
-    void                            stopCurrentTimer(void);
-    void                            stopAllTimers(void);
-    void                            resetTimers(void);
-    std::map<std::string, GridTime> getTimings(void);
 protected:
     // environment shortcut
     DEFINE_ENV_ALIAS;

Hadrons/Modules.hpp

@@ -1,6 +1,6 @@
 #include <Hadrons/Modules/MContraction/Baryon.hpp>
+#include <Hadrons/Modules/MContraction/A2AAslashField.hpp>
 #include <Hadrons/Modules/MContraction/A2AMesonField.hpp>
-#include <Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp>
 #include <Hadrons/Modules/MContraction/Meson.hpp>
 #include <Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
 #include <Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
@@ -16,6 +16,7 @@
 #include <Hadrons/Modules/MSource/Wall.hpp>
 #include <Hadrons/Modules/MSource/Z2.hpp>
 #include <Hadrons/Modules/MSource/SeqConserved.hpp>
+#include <Hadrons/Modules/MSource/Momentum.hpp>
 #include <Hadrons/Modules/MSink/Smear.hpp>
 #include <Hadrons/Modules/MSink/Point.hpp>
 #include <Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp>
@@ -27,6 +28,7 @@
 #include <Hadrons/Modules/MGauge/StoutSmearing.hpp>
 #include <Hadrons/Modules/MGauge/Unit.hpp>
 #include <Hadrons/Modules/MGauge/Random.hpp>
+#include <Hadrons/Modules/MGauge/GaugeFix.hpp>
 #include <Hadrons/Modules/MGauge/FundtoHirep.hpp>
 #include <Hadrons/Modules/MGauge/StochEm.hpp>
 #include <Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp>
@@ -40,6 +42,9 @@
 #include <Hadrons/Modules/MScalar/ScalarVP.hpp>
 #include <Hadrons/Modules/MScalar/Scalar.hpp>
 #include <Hadrons/Modules/MScalar/ChargedProp.hpp>
+#include <Hadrons/Modules/MNPR/Bilinear.hpp>
+#include <Hadrons/Modules/MNPR/Amputate.hpp>
+#include <Hadrons/Modules/MNPR/FourQuark.hpp>
 #include <Hadrons/Modules/MAction/DWF.hpp>
 #include <Hadrons/Modules/MAction/MobiusDWF.hpp>
 #include <Hadrons/Modules/MAction/Wilson.hpp>
@@ -50,7 +55,6 @@
 #include <Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
 #include <Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
 #include <Hadrons/Modules/MScalarSUN/Div.hpp>
-#include <Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp>
 #include <Hadrons/Modules/MScalarSUN/TrMag.hpp>
 #include <Hadrons/Modules/MScalarSUN/EMT.hpp>
 #include <Hadrons/Modules/MScalarSUN/TwoPoint.hpp>

Hadrons/Modules/MContraction/A2AAslashField.cc

@@ -2,7 +2,7 @@
 
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: Hadrons/Modules/MScalarSUN/TimeMomProbe.cc
+Source file: Hadrons/Modules/MContraction/A2AAslashField.cc
 
 Copyright (C) 2015-2018
 
@@ -25,14 +25,10 @@ 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/MScalarSUN/TimeMomProbe.hpp>
+#include <Hadrons/Modules/MContraction/A2AAslashField.hpp>
 
 using namespace Grid;
 using namespace Hadrons;
-using namespace MScalarSUN;
+using namespace MContraction;
 
-template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<2>>;
-template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<3>>;
-template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<4>>;
-template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<5>>;
-template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<6>>;
+template class Grid::Hadrons::MContraction::TA2AAslashField<FIMPL, PhotonR>;

Hadrons/Modules/MContraction/A2AAslashField.hpp

@@ -0,0 +1,250 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MContraction/A2AAslashField.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+#ifndef Hadrons_MContraction_A2AAslashField_hpp_
+#define Hadrons_MContraction_A2AAslashField_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Hadrons/A2AMatrix.hpp>
+
+#ifndef ASF_IO_TYPE
+#define ASF_IO_TYPE ComplexF
+#endif
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         A2AAslashField                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class A2AAslashFieldPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(A2AAslashFieldPar,
+                                    int, cacheBlock,
+                                    int, block,
+                                    std::string, left,
+                                    std::string, right,
+                                    std::string, output,
+                                    std::vector<std::string>, emField);
+};
+
+class A2AAslashFieldMetadata: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(A2AAslashFieldMetadata,
+                                    std::string, emFieldName);
+};
+
+template <typename T, typename FImpl>
+class AslashFieldKernel: public A2AKernel<T, typename FImpl::FermionField>
+{
+public:
+    typedef typename FImpl::FermionField FermionField;
+public:
+    AslashFieldKernel(const std::vector<LatticeComplex> &emB0,
+                      const std::vector<LatticeComplex> &emB1,
+                      GridBase *grid)
+    : emB0_(emB0), emB1_(emB1), grid_(grid)
+    {
+        vol_ = 1.;
+        for (auto &d: grid_->GlobalDimensions())
+        {
+            vol_ *= d;
+        }
+    }
+
+    virtual ~AslashFieldKernel(void) = default;
+    virtual void operator()(A2AMatrixSet<T> &m, const FermionField *left, 
+                            const FermionField *right,
+                            const unsigned int orthogDim, double &t)
+    {
+        A2Autils<FImpl>::AslashField(m, left, right, emB0_, emB1_, orthogDim, &t);
+    }
+
+    virtual double flops(const unsigned int blockSizei, const unsigned int blockSizej)
+    {
+        return 0.;
+    }
+
+    virtual double bytes(const unsigned int blockSizei, const unsigned int blockSizej)
+    {
+        return 0.;
+    }
+private:
+    const std::vector<LatticeComplex> &emB0_, &emB1_;
+    GridBase                          *grid_;
+    double                            vol_;
+};
+
+template <typename FImpl, typename PhotonImpl>
+class TA2AAslashField: public Module<A2AAslashFieldPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+    typedef typename PhotonImpl::GaugeField EmField;
+    typedef A2AMatrixBlockComputation<Complex, 
+                                      FermionField, 
+                                      A2AAslashFieldMetadata, 
+                                      ASF_IO_TYPE> Computation;
+    typedef AslashFieldKernel<Complex, FImpl> Kernel;
+public:
+    // constructor
+    TA2AAslashField(const std::string name);
+    // destructor
+    virtual ~TA2AAslashField(void) {};
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(A2AAslashField, ARG(TA2AAslashField<FIMPL, PhotonR>), MContraction);
+
+/******************************************************************************
+ *                 TA2AAslashField implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl, typename PhotonImpl>
+TA2AAslashField<FImpl, PhotonImpl>::TA2AAslashField(const std::string name)
+: Module<A2AAslashFieldPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl, typename PhotonImpl>
+std::vector<std::string> TA2AAslashField<FImpl, PhotonImpl>::getInput(void)
+{
+    std::vector<std::string> in = par().emField;
+    
+    in.push_back(par().left);
+    in.push_back(par().right);
+
+    return in;
+}
+
+template <typename FImpl, typename PhotonImpl>
+std::vector<std::string> TA2AAslashField<FImpl, PhotonImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl, typename PhotonImpl>
+void TA2AAslashField<FImpl, PhotonImpl>::setup(void)
+{
+    envTmp(Computation, "computation", 1, envGetGrid(FermionField), 
+           env().getNd() - 1, par().emField.size(), 1, par().block, 
+           par().cacheBlock, this);
+    envTmp(std::vector<ComplexField>, "B0", 1, 
+           par().emField.size(), envGetGrid(ComplexField));
+    envTmp(std::vector<ComplexField>, "B1", 1, 
+           par().emField.size(), envGetGrid(ComplexField));
+    envTmpLat(ComplexField, "Amu");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl, typename PhotonImpl>
+void TA2AAslashField<FImpl, PhotonImpl>::execute(void)
+{
+    auto &left  = envGet(std::vector<FermionField>, par().left);
+    auto &right = envGet(std::vector<FermionField>, par().right);
+
+    int nt         = env().getDim().back();
+    int N_i        = left.size();
+    int N_j        = right.size();
+    int nem        = par().emField.size();
+    int block      = par().block;
+    int cacheBlock = par().cacheBlock;
+
+    LOG(Message) << "Computing all-to-all A-slash fields" << std::endl;
+    LOG(Message) << "Left: '" << par().left << "' Right: '" << par().right << "'" << std::endl;
+    LOG(Message) << "EM fields:" << std::endl;
+    for (auto &name: par().emField)
+    {
+        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)) 
+                 << "/EM field)" << std::endl;
+    
+    // preparing "B" complexified fields
+    startTimer("Complexify EM fields");
+    envGetTmp(std::vector<ComplexField>, B0);
+    envGetTmp(std::vector<ComplexField>, B1);
+    for (unsigned int i = 0; i < par().emField.size(); ++i)
+    {
+        auto &A = envGet(EmField, par().emField[i]);
+        envGetTmp(ComplexField, Amu);
+
+        B0[i]  = peekLorentz(A, 0);
+        B0[i] += timesI(peekLorentz(A, 1));
+        B1[i]  = peekLorentz(A, 2);
+        B1[i] += timesI(peekLorentz(A, 3));
+    }
+    stopTimer("Complexify EM fields");
+
+    // I/O name & metadata lambdas
+    auto ionameFn = [this](const unsigned int em, const unsigned int dummy)
+    {
+        return par().emField[em];
+    };
+
+    auto filenameFn = [this, &ionameFn](const unsigned int em, const unsigned int dummy)
+    {
+        return par().output + "." + std::to_string(vm().getTrajectory()) 
+               + "/" + ionameFn(em, dummy) + ".h5";
+    };
+
+    auto metadataFn = [this](const unsigned int em, const unsigned int dummy)
+    {
+        A2AAslashFieldMetadata md;
+
+        md.emFieldName = par().emField[em];
+        
+        return md;
+    };
+
+    // executing computation
+    Kernel kernel(B0, B1, envGetGrid(FermionField));
+
+    envGetTmp(Computation, computation);
+    computation.execute(left, right, kernel, ionameFn, filenameFn, metadataFn);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MContraction_A2AAslashField_hpp_

Hadrons/Modules/MContraction/A2AMesonField.cc

@@ -33,4 +33,3 @@ using namespace Hadrons;
 using namespace MContraction;
 
 template class Grid::Hadrons::MContraction::TA2AMesonField<FIMPL>;
-template class Grid::Hadrons::MContraction::TA2AMesonField<ZFIMPL>;

Hadrons/Modules/MContraction/A2AMesonField.hpp

@@ -33,12 +33,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Hadrons/Global.hpp>
 #include <Hadrons/Module.hpp>
 #include <Hadrons/ModuleFactory.hpp>
-#include <Hadrons/A2AVectors.hpp>
 #include <Hadrons/A2AMatrix.hpp>
-#include <Hadrons/Modules/MSolver/A2AVectors.hpp>
-#include <Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp>
 
-#define MF_PARALLEL_IO
 #ifndef MF_IO_TYPE
 #define MF_IO_TYPE ComplexF
 #endif
@@ -56,8 +52,8 @@ public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMesonFieldPar,
                                     int, cacheBlock,
                                     int, block,
-                                    std::string, v,
-                                    std::string, w,
+                                    std::string, left,
+                                    std::string, right,
                                     std::string, output,
                                     std::string, gammas,
                                     std::vector<std::string>, mom);
@@ -71,22 +67,59 @@ public:
                                     Gamma::Algebra, gamma);
 };
 
+template <typename T, typename FImpl>
+class MesonFieldKernel: public A2AKernel<T, typename FImpl::FermionField>
+{
+public:
+    typedef typename FImpl::FermionField FermionField;
+public:
+    MesonFieldKernel(const std::vector<Gamma::Algebra> &gamma,
+                     const std::vector<LatticeComplex> &mom,
+                     GridBase *grid)
+    : gamma_(gamma), mom_(mom), grid_(grid)
+    {
+        vol_ = 1.;
+        for (auto &d: grid_->GlobalDimensions())
+        {
+            vol_ *= d;
+        }
+    }
+
+    virtual ~MesonFieldKernel(void) = default;
+    virtual void operator()(A2AMatrixSet<T> &m, const FermionField *left, 
+                            const FermionField *right,
+                            const unsigned int orthogDim, double &t)
+    {
+        A2Autils<FImpl>::MesonField(m, left, right, gamma_, mom_, orthogDim, &t);
+    }
+
+    virtual double flops(const unsigned int blockSizei, const unsigned int blockSizej)
+    {
+        return vol_*(2*8.0+6.0+8.0*mom_.size())*blockSizei*blockSizej*gamma_.size();
+    }
+
+    virtual double bytes(const unsigned int blockSizei, const unsigned int blockSizej)
+    {
+        return vol_*(12.0*sizeof(T))*blockSizei*blockSizej
+               +  vol_*(2.0*sizeof(T)*mom_.size())*blockSizei*blockSizej*gamma_.size();
+    }
+private:
+    const std::vector<Gamma::Algebra> &gamma_;
+    const std::vector<LatticeComplex> &mom_;
+    GridBase                          *grid_;
+    double                            vol_;
+};
+
 template <typename FImpl>
 class TA2AMesonField : public Module<A2AMesonFieldPar>
 {
 public:
     FERM_TYPE_ALIASES(FImpl,);
-    SOLVER_TYPE_ALIASES(FImpl,);
-    typedef Eigen::TensorMap<Eigen::Tensor<Complex, 5, Eigen::RowMajor>>    MesonField;
-    typedef Eigen::TensorMap<Eigen::Tensor<MF_IO_TYPE, 5, Eigen::RowMajor>> MesonFieldIo;
-    typedef A2AMatrixIo<MF_IO_TYPE, A2AMesonFieldMetadata>                  MatrixIo;
-    struct IoHelper
-    {
-        MatrixIo              io;
-        A2AMesonFieldMetadata metadata;
-        size_t                offset;
-        unsigned int          i, j, blockSizei, blockSizej;
-    };
+    typedef A2AMatrixBlockComputation<Complex, 
+                                      FermionField, 
+                                      A2AMesonFieldMetadata, 
+                                      MF_IO_TYPE> Computation;
+    typedef MesonFieldKernel<Complex, FImpl> Kernel;
 public:
     // constructor
     TA2AMesonField(const std::string name);
@@ -100,20 +133,13 @@ public:
     // execution
     virtual void execute(void);
 private:
-    // IO
-    std::string ioname(const unsigned int m, const unsigned int g) const;
-    std::string filename(const unsigned int m, const unsigned int g) const;
-    void saveBlock(const MF_IO_TYPE *data, IoHelper &h);
-private:
-    bool                                               hasPhase_{false};
-    std::string                                        momphName_;
-    std::vector<Gamma::Algebra>                        gamma_;
-    std::vector<std::vector<Real>>                     mom_;
-    std::vector<IoHelper>                              nodeIo_;
+    bool                               hasPhase_{false};
+    std::string                        momphName_;
+    std::vector<Gamma::Algebra>        gamma_;
+    std::vector<std::vector<Real>>     mom_;
 };
 
 MODULE_REGISTER(A2AMesonField, ARG(TA2AMesonField<FIMPL>), MContraction);
-MODULE_REGISTER(ZA2AMesonField, ARG(TA2AMesonField<ZFIMPL>), MContraction);
 
 /******************************************************************************
 *                  TA2AMesonField implementation                             *
@@ -130,7 +156,7 @@ TA2AMesonField<FImpl>::TA2AMesonField(const std::string name)
 template <typename FImpl>
 std::vector<std::string> TA2AMesonField<FImpl>::getInput(void)
 {
-    std::vector<std::string> in = {par().v, par().w};
+    std::vector<std::string> in = {par().left, par().right};
 
     return in;
 }
@@ -186,34 +212,31 @@ void TA2AMesonField<FImpl>::setup(void)
         }
         mom_.push_back(p);
     }
-    
     envCache(std::vector<ComplexField>, momphName_, 1, 
              par().mom.size(), envGetGrid(ComplexField));
     envTmpLat(ComplexField, "coor");
-    // preallocate memory for meson field block
-    auto tgp = env().getDim().back()*gamma_.size()*mom_.size();
-
-    envTmp(Vector<MF_IO_TYPE>, "mfBuf", 1, tgp*par().block*par().block);
-    envTmp(Vector<Complex>, "mfCache", 1, tgp*par().cacheBlock*par().cacheBlock);
+    envTmp(Computation, "computation", 1, envGetGrid(FermionField), 
+           env().getNd() - 1, mom_.size(), gamma_.size(), par().block, 
+           par().cacheBlock, this);
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
 void TA2AMesonField<FImpl>::execute(void)
 {
-    auto &v = envGet(std::vector<FermionField>, par().v);
-    auto &w = envGet(std::vector<FermionField>, par().w);
+    auto &left  = envGet(std::vector<FermionField>, par().left);
+    auto &right = envGet(std::vector<FermionField>, par().right);
 
     int nt         = env().getDim().back();
-    int N_i        = w.size();
-    int N_j        = v.size();
+    int N_i        = left.size();
+    int N_j        = right.size();
     int ngamma     = gamma_.size();
     int nmom       = mom_.size();
     int block      = par().block;
     int cacheBlock = par().cacheBlock;
 
     LOG(Message) << "Computing all-to-all meson fields" << std::endl;
-    LOG(Message) << "W: '" << par().w << "' V: '" << par().v << "'" << std::endl;
+    LOG(Message) << "Left: '" << par().left << "' Right: '" << par().right << "'" << std::endl;
     LOG(Message) << "Momenta:" << std::endl;
     for (auto &p: mom_)
     {
@@ -228,9 +251,6 @@ void TA2AMesonField<FImpl>::execute(void)
                  << " (filesize " << sizeString(nt*N_i*N_j*sizeof(MF_IO_TYPE)) 
                  << "/momentum/bilinear)" << std::endl;
 
-    ///////////////////////////////////////////////
-    // Momentum setup
-    ///////////////////////////////////////////////
     auto &ph = envGet(std::vector<ComplexField>, momphName_);
 
     if (!hasPhase_)
@@ -253,189 +273,43 @@ void TA2AMesonField<FImpl>::execute(void)
         hasPhase_ = true;
         stopTimer("Momentum phases");
     }
-    
-    //////////////////////////////////////////////////////////////////////////
-    // i,j   is first  loop over SchurBlock factors reusing 5D matrices
-    // ii,jj is second loop over cacheBlock factors for high perf contractoin
-    // iii,jjj are loops within cacheBlock
-    // Total index is sum of these  i+ii+iii etc...
-    //////////////////////////////////////////////////////////////////////////
-    
-    double flops;
-    double bytes;
-    double vol      = env().getVolume();
-    double t_kernel = 0.0;
-    double nodes    = env().getGrid()->NodeCount();
-    double tot_kernel;
 
-    envGetTmp(Vector<MF_IO_TYPE>, mfBuf);
-    envGetTmp(Vector<Complex>, mfCache);
-    
-    double t0    = usecond();
-    int NBlock_i = N_i/block + (((N_i % block) != 0) ? 1 : 0);
-    int NBlock_j = N_j/block + (((N_j % block) != 0) ? 1 : 0);
-
-    for(int i=0;i<N_i;i+=block)
-    for(int j=0;j<N_j;j+=block)
+    auto ionameFn = [this](const unsigned int m, const unsigned int g)
     {
-        // Get the W and V vectors for this block^2 set of terms
-        int N_ii = MIN(N_i-i,block);
-        int N_jj = MIN(N_j-j,block);
+        std::stringstream ss;
 
-        LOG(Message) << "Meson field block " 
-                    << j/block + NBlock_j*i/block + 1 
-                    << "/" << NBlock_i*NBlock_j << " [" << i <<" .. " 
-                    << i+N_ii-1 << ", " << j <<" .. " << j+N_jj-1 << "]" 
-                    << std::endl;
-
-        MesonFieldIo mfBlock(mfBuf.data(),nmom,ngamma,nt,N_ii,N_jj);
-
-        // Series of cache blocked chunks of the contractions within this block
-        flops = 0.0;
-        bytes = 0.0;
-        for(int ii=0;ii<N_ii;ii+=cacheBlock)
-        for(int jj=0;jj<N_jj;jj+=cacheBlock)
+        ss << gamma_[g] << "_";
+        for (unsigned int mu = 0; mu < mom_[m].size(); ++mu)
         {
-            int N_iii = MIN(N_ii-ii,cacheBlock);
-            int N_jjj = MIN(N_jj-jj,cacheBlock);
-            MesonField mfCacheBlock(mfCache.data(),nmom,ngamma,nt,N_iii,N_jjj);    
-
-            startTimer("contraction: total");
-            makeMesonFieldBlock(mfCacheBlock, &w[i+ii], &v[j+jj], gamma_, ph, 
-                                env().getNd() - 1, this);
-            stopTimer("contraction: total");
-            
-            // flops for general N_c & N_s
-            flops += vol * ( 2 * 8.0 + 6.0 + 8.0*nmom) * N_iii*N_jjj*ngamma;
-            bytes += vol * (12.0 * sizeof(Complex) ) * N_iii*N_jjj
-                     +  vol * ( 2.0 * sizeof(Complex) *nmom ) * N_iii*N_jjj* ngamma;
-
-            startTimer("cache copy");
-            parallel_for_nest5(int m =0;m< nmom;m++)
-            for(int g =0;g< ngamma;g++)
-            for(int t =0;t< nt;t++)
-            for(int iii=0;iii< N_iii;iii++)
-            for(int jjj=0;jjj< N_jjj;jjj++)
-            {
-                mfBlock(m,g,t,ii+iii,jj+jjj) = mfCacheBlock(m,g,t,iii,jjj);
-            }
-            stopTimer("cache copy");
+            ss << mom_[m][mu] << ((mu == mom_[m].size() - 1) ? "" : "_");
         }
 
-        // perf
-        tot_kernel = getDTimer("contraction: colour trace & mom.")
-                     + getDTimer("contraction: local space sum");
-        t_kernel   = tot_kernel - t_kernel;
-        LOG(Message) << "Kernel perf " << flops/t_kernel/1.0e3/nodes 
-                     << " Gflop/s/node " << std::endl;
-        LOG(Message) << "Kernel perf " << bytes/t_kernel*1.0e6/1024/1024/1024/nodes 
-                     << " GB/s/node "  << std::endl;
-        t_kernel = tot_kernel;
+        return ss.str();
+    };
 
-        // IO
-        if (!par().output.empty())
+    auto filenameFn = [this, &ionameFn](const unsigned int m, const unsigned int g)
+    {
+        return par().output + "." + std::to_string(vm().getTrajectory()) 
+               + "/" + ionameFn(m, g) + ".h5";
+    };
+
+    auto metadataFn = [this](const unsigned int m, const unsigned int g)
+    {
+        A2AMesonFieldMetadata md;
+
+        for (auto pmu: mom_[m])
         {
-            double       blockSize, ioTime;
-            unsigned int myRank = env().getGrid()->ThisRank(),
-                         nRank  = env().getGrid()->RankCount();
+            md.momentum.push_back(pmu);
+        }
+        md.gamma = gamma_[g];
         
-            LOG(Message) << "Writing block to disk" << std::endl;
-            ioTime = -getDTimer("IO: write block");
-            startTimer("IO: total");
-            makeFileDir(filename(0, 0), env().getGrid());
-#ifdef MF_PARALLEL_IO
-            env().getGrid()->Barrier();
-            nodeIo_.clear();
-            for(int f = myRank; f < nmom*ngamma; f += nRank)
-            {
-                const unsigned int    m = f/ngamma, g = f % ngamma;
-                IoHelper              h;
+        return md;
+    };
 
-                h.io = MatrixIo(filename(m, g), ioname(m, g), nt, N_i, N_j);
-                for (auto pmu: mom_[m])
-                {
-                    h.metadata.momentum.push_back(pmu);
-                }
-                h.metadata.gamma = gamma_[g];
-                h.i              = i;
-                h.j              = j;
-                h.blockSizei     = mfBlock.dimension(3);
-                h.blockSizej     = mfBlock.dimension(4);
-                h.offset         = (m*ngamma + g)*nt*h.blockSizei*h.blockSizej;
-                nodeIo_.push_back(h);
-            }
-            // parallel IO
-            for (auto &h: nodeIo_)
-            {
-                saveBlock(mfBlock.data(), h);
-            }
-            env().getGrid()->Barrier();
-#else
-            // serial IO
-            for(int m = 0; m < nmom; m++)
-            for(int g = 0; g < ngamma; g++)
-            {
-                IoHelper h;
+    Kernel      kernel(gamma_, ph, envGetGrid(FermionField));
 
-                h.io = MatrixIo(filename(m, g), ioname(m, g), nt, N_i, N_j);
-                for (auto pmu: mom_[m])
-                {
-                    h.metadata.momentum.push_back(pmu);
-                }
-                h.metadata.gamma = gamma_[g];
-                h.i              = i;
-                h.j              = j;
-                h.blockSizei     = mfBlock.dimension(3);
-                h.blockSizej     = mfBlock.dimension(4);
-                h.offset         = (m*ngamma + g)*nt*h.blockSizei*h.blockSizej;
-                saveBlock(mfBlock.data(), h);
-            }
-#endif
-            stopTimer("IO: total");
-            blockSize  = static_cast<double>(nmom*ngamma*nt*N_ii*N_jj*sizeof(MF_IO_TYPE));
-            ioTime    += getDTimer("IO: write block");
-            LOG(Message) << "HDF5 IO done " << sizeString(blockSize) << " in "
-                         << ioTime  << " us (" 
-                         << blockSize/ioTime*1.0e6/1024/1024
-                         << " MB/s)" << std::endl;
-        }
-    }
-}
-
-// IO
-template <typename FImpl>
-std::string TA2AMesonField<FImpl>::ioname(unsigned int m, unsigned int g) const
-{
-    std::stringstream ss;
-
-    ss << gamma_[g] << "_";
-    for (unsigned int mu = 0; mu < mom_[m].size(); ++mu)
-    {
-        ss << mom_[m][mu] << ((mu == mom_[m].size() - 1) ? "" : "_");
-    }
-
-    return ss.str();
-}
-
-template <typename FImpl>
-std::string TA2AMesonField<FImpl>::filename(unsigned int m, unsigned int g) const
-{
-    return par().output + "." + std::to_string(vm().getTrajectory()) 
-           + "/" + ioname(m, g) + ".h5";
-}
-
-template <typename FImpl>
-void TA2AMesonField<FImpl>::saveBlock(const MF_IO_TYPE *data, IoHelper &h)
-{
-    if ((h.i == 0) and (h.j == 0))
-    {
-        startTimer("IO: file creation");
-        h.io.initFile(h.metadata, par().block);
-        stopTimer("IO: file creation");
-    }
-    startTimer("IO: write block");
-    h.io.saveBlock(data + h.offset, h.i, h.j, h.blockSizei, h.blockSizej);
-    stopTimer("IO: write block");
+    envGetTmp(Computation, computation);
+    computation.execute(left, right, kernel, ionameFn, filenameFn, metadataFn);
 }
 
 END_MODULE_NAMESPACE

Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp

@@ -1,224 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-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 */
-#ifndef Hadrons_MContraction_A2AMesonFieldKernels_hpp_
-#define Hadrons_MContraction_A2AMesonFieldKernels_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Grid/Eigen/unsupported/CXX11/Tensor>
-
-BEGIN_HADRONS_NAMESPACE
-
-BEGIN_MODULE_NAMESPACE(MContraction)
-
-////////////////////////////////////////////////////////////////////////////////
-// Cache blocked arithmetic routine
-// Could move to Grid ???
-////////////////////////////////////////////////////////////////////////////////
-template <typename Field, typename MesonField>
-void makeMesonFieldBlock(MesonField &mat, 
-                         const Field *lhs_wi,
-                         const Field *rhs_vj,
-                         std::vector<Gamma::Algebra> gamma,
-                         const std::vector<LatticeComplex> &mom,
-                         int orthogdim,
-                         ModuleBase *caller = nullptr) 
-{
-    typedef typename Field::vector_object vobj;
-    typedef typename vobj::scalar_object  sobj;
-    typedef typename vobj::scalar_type    scalar_type;
-    typedef typename vobj::vector_type    vector_type;
-
-    typedef iSpinMatrix<vector_type> SpinMatrix_v;
-    typedef iSpinMatrix<scalar_type> SpinMatrix_s;
-    
-    int Lblock = mat.dimension(3); 
-    int Rblock = mat.dimension(4);
-
-    GridBase *grid = lhs_wi[0]._grid;
-    
-    const int    Nd = grid->_ndimension;
-    const int Nsimd = grid->Nsimd();
-
-    int Nt     = grid->GlobalDimensions()[orthogdim];
-    int Ngamma = gamma.size();
-    int Nmom   = mom.size();
-
-    int fd=grid->_fdimensions[orthogdim];
-    int ld=grid->_ldimensions[orthogdim];
-    int rd=grid->_rdimensions[orthogdim];
-
-    // will locally sum vectors first
-    // sum across these down to scalars
-    // splitting the SIMD
-    int MFrvol = rd*Lblock*Rblock*Nmom;
-    int MFlvol = ld*Lblock*Rblock*Nmom;
-
-    Vector<SpinMatrix_v > lvSum(MFrvol);
-    parallel_for (int r = 0; r < MFrvol; r++)
-    {
-        lvSum[r] = zero;
-    }
-
-    Vector<SpinMatrix_s > lsSum(MFlvol);             
-    parallel_for (int r = 0; r < MFlvol; r++)
-    {
-        lsSum[r]=scalar_type(0.0);
-    }
-
-    int e1=    grid->_slice_nblock[orthogdim];
-    int e2=    grid->_slice_block [orthogdim];
-    int stride=grid->_slice_stride[orthogdim];
-
-    if (caller) caller->startTimer("contraction: colour trace & mom.");
-    // Nested parallelism would be ok
-    // Wasting cores here. Test case r
-    parallel_for(int r=0;r<rd;r++)
-    {
-        int so=r*grid->_ostride[orthogdim]; // base offset for start of plane 
-
-        for(int n=0;n<e1;n++)
-        for(int b=0;b<e2;b++)
-        {
-            int ss= so+n*stride+b;
-
-            for(int i=0;i<Lblock;i++)
-            {
-                auto left = conjugate(lhs_wi[i]._odata[ss]);
-
-                for(int j=0;j<Rblock;j++)
-                {
-                    SpinMatrix_v vv;
-                    auto right = rhs_vj[j]._odata[ss];
-
-                    for(int s1=0;s1<Ns;s1++)
-                    for(int s2=0;s2<Ns;s2++)
-                    {
-                        vv()(s1,s2)() = left()(s2)(0) * right()(s1)(0)
-                                        + left()(s2)(1) * right()(s1)(1)
-                                        + left()(s2)(2) * right()(s1)(2);
-                    }
-                    
-                    // After getting the sitewise product do the mom phase loop
-                    int base = Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*r;
-
-                    for ( int m=0;m<Nmom;m++)
-                    {
-                        int idx = m+base;
-                        auto phase = mom[m]._odata[ss];
-                        mac(&lvSum[idx],&vv,&phase);
-                    }
-                }
-            }
-        }
-    }
-    if (caller) caller->stopTimer("contraction: colour trace & mom.");
-
-    // Sum across simd lanes in the plane, breaking out orthog dir.
-    if (caller) caller->startTimer("contraction: local space sum");
-    parallel_for(int rt=0;rt<rd;rt++)
-    {
-        std::vector<int> icoor(Nd);
-        std::vector<SpinMatrix_s> extracted(Nsimd);               
-
-        for(int i=0;i<Lblock;i++)
-        for(int j=0;j<Rblock;j++)
-        for(int m=0;m<Nmom;m++)
-        {
-
-            int ij_rdx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*rt;
-
-            extract(lvSum[ij_rdx],extracted);
-            for(int idx=0;idx<Nsimd;idx++)
-            {
-                grid->iCoorFromIindex(icoor,idx);
-
-                int ldx    = rt+icoor[orthogdim]*rd;
-                int ij_ldx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*ldx;
-
-                lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx];
-            }
-        }
-    }
-    if (caller) caller->stopTimer("contraction: local space sum");
-
-    // ld loop and local only??
-    if (caller) caller->startTimer("contraction: spin trace");
-    int pd = grid->_processors[orthogdim];
-    int pc = grid->_processor_coor[orthogdim];
-    parallel_for_nest2(int lt=0;lt<ld;lt++)
-    {
-        for(int pt=0;pt<pd;pt++)
-        {
-            int t = lt + pt*ld;
-            if (pt == pc)
-            {
-                for(int i=0;i<Lblock;i++)
-                for(int j=0;j<Rblock;j++)
-                for(int m=0;m<Nmom;m++)
-                {
-                    int ij_dx = m+Nmom*i + Nmom*Lblock * j + Nmom*Lblock * Rblock * lt;
-
-                    for(int mu=0;mu<Ngamma;mu++)
-                    {
-                        // this is a bit slow
-                        mat(m,mu,t,i,j) = trace(lsSum[ij_dx]*Gamma(gamma[mu]));
-                    }
-                }
-            } 
-            else 
-            { 
-                const scalar_type zz(0.0);
-
-                for(int i=0;i<Lblock;i++)
-                for(int j=0;j<Rblock;j++)
-                for(int mu=0;mu<Ngamma;mu++)
-                for(int m=0;m<Nmom;m++)
-                {
-                    mat(m,mu,t,i,j) =zz;
-                }
-            }
-        }
-    }
-    if (caller) caller->stopTimer("contraction: spin trace");
-    ////////////////////////////////////////////////////////////////////
-    // This global sum is taking as much as 50% of time on 16 nodes
-    // Vector size is 7 x 16 x 32 x 16 x 16 x sizeof(complex) = 2MB - 60MB depending on volume
-    // Healthy size that should suffice
-    ////////////////////////////////////////////////////////////////////
-    if (caller) caller->startTimer("contraction: global sum");
-    grid->GlobalSumVector(&mat(0,0,0,0,0),Nmom*Ngamma*Nt*Lblock*Rblock);
-    if (caller) caller->stopTimer("contraction: global sum");
-}
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif //Hadrons_MContraction_A2AMesonField_hpp_

Hadrons/Modules/MGauge/GaugeFix.cc

@@ -0,0 +1,36 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MGauge/GaugeFix.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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 */
+
+#include <Hadrons/Modules/MGauge/GaugeFix.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MGauge;
+
+template class Grid::Hadrons::MGauge::TGaugeFix<GIMPL>;

Hadrons/Modules/MGauge/GaugeFix.hpp

@@ -0,0 +1,135 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MGauge/GaugeFix.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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 */
+
+#ifndef Hadrons_MGaugeFix_hpp_
+#define Hadrons_MGaugeFix_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Grid/qcd/utils/GaugeFix.h>
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                              Fix gauge                                    *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MGauge)
+
+class GaugeFixPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(GaugeFixPar,
+                                    std::string, gauge,
+                                    Real,  alpha,
+                                    int, maxiter, 
+                                    Real, Omega_tol, 
+                                    Real, Phi_tol,
+                                    bool, Fourier);
+};
+
+template <typename GImpl>
+class TGaugeFix: public Module<GaugeFixPar>
+{
+public:
+    GAUGE_TYPE_ALIASES(GImpl,);
+public:
+    // constructor
+    TGaugeFix(const std::string name);
+    // destructor
+    virtual ~TGaugeFix(void) {};
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(GaugeFix, TGaugeFix<GIMPL>, MGauge);
+
+/******************************************************************************
+*                            TGaugeFix implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename GImpl>
+TGaugeFix<GImpl>::TGaugeFix(const std::string name)
+: Module<GaugeFixPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename GImpl>
+std::vector<std::string> TGaugeFix<GImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().gauge};
+    return in;
+}
+
+template <typename GImpl>
+std::vector<std::string> TGaugeFix<GImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename GImpl>
+void TGaugeFix<GImpl>::setup(void)
+{
+    envCreateLat(GaugeField, getName());
+}
+
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename GImpl>
+void TGaugeFix<GImpl>::execute(void)
+//Loads the gauge and fixes it
+{
+    std::cout << "executing" << std::endl;
+    LOG(Message) << "Fixing the Gauge" << std::endl;
+    LOG(Message) << par().gauge << std::endl;
+    auto &U     = envGet(GaugeField, par().gauge);
+    auto &Umu   = envGet(GaugeField, getName());
+    LOG(Message) << "Gauge Field fetched" << std::endl;
+    //do we allow maxiter etc to be user set?
+    Real alpha     = par().alpha;
+    int  maxiter   = par().maxiter;
+    Real Omega_tol = par().Omega_tol;
+    Real Phi_tol   = par().Phi_tol;
+    bool Fourier   = par().Fourier;
+    FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(U,alpha,maxiter,Omega_tol,Phi_tol,Fourier);
+    Umu = U;
+    LOG(Message) << "Gauge Fixed" << std::endl;
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MGaugeFix_hpp_

Hadrons/Modules/MIO/LoadEigenPack.cc

@@ -32,4 +32,4 @@ using namespace Hadrons;
 using namespace MIO;
 
 template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL>>;
-
+template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL, FIMPLF>>;

Hadrons/Modules/MIO/LoadEigenPack.hpp

@@ -54,7 +54,9 @@ template <typename Pack>
 class TLoadEigenPack: public Module<LoadEigenPackPar>
 {
 public:
-    typedef EigenPack<typename Pack::Field> BasePack;
+    typedef typename Pack::Field   Field;
+    typedef typename Pack::FieldIo FieldIo;
+    typedef BaseEigenPack<Field>   BasePack;
 public:
     // constructor
     TLoadEigenPack(const std::string name);
@@ -70,6 +72,7 @@ public:
 };
 
 MODULE_REGISTER_TMP(LoadFermionEigenPack, TLoadEigenPack<FermionEigenPack<FIMPL>>, MIO);
+MODULE_REGISTER_TMP(LoadFermionEigenPackIo32, ARG(TLoadEigenPack<FermionEigenPack<FIMPL, FIMPLF>>), MIO);
 
 /******************************************************************************
  *                    TLoadEigenPack implementation                           *
@@ -101,9 +104,14 @@ std::vector<std::string> TLoadEigenPack<Pack>::getOutput(void)
 template <typename Pack>
 void TLoadEigenPack<Pack>::setup(void)
 {
-    env().createGrid(par().Ls);
+    GridBase *gridIo = nullptr;
+
+    if (typeHash<Field>() != typeHash<FieldIo>())
+    {
+        gridIo = envGetRbGrid(FieldIo, par().Ls);
+    }
     envCreateDerived(BasePack, Pack, getName(), par().Ls, par().size, 
-                     env().getRbGrid(par().Ls));
+                     envGetRbGrid(Field, par().Ls), gridIo);
 }
 
 // execution ///////////////////////////////////////////////////////////////////

Hadrons/Modules/MNPR/Amputate.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/Amputate.cc
+
+Copyright (C) 2015-2018
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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 */
+#include <Hadrons/Modules/MNPR/Amputate.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MNPR;
+
+template class Grid::Hadrons::MNPR::TAmputate<FIMPL,FIMPL>;
+

Hadrons/Modules/MNPR/Amputate.hpp

@@ -0,0 +1,200 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/Amputate.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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 */
+
+#ifndef Hadrons_Amputate_hpp_
+#define Hadrons_Amputate_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Grid/Eigen/LU>
+//#include <Grid/qcd/utils/PropagatorUtils.h>
+//#include <Grid/serialisation/Serialisation.h>
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                                TAmputate                                       *
+        Performs bilinear contractions of the type tr[g5*adj(Sout)*g5*G*Sin]
+        Suitable for non exceptional momenta
+******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MNPR)
+
+class AmputatePar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(AmputatePar,
+                                    std::string,    Sin, //need to make this a propogator type?
+                                    std::string,    Sout, //same
+                                    std::string,    vertex,
+                                    std::string,    pin,
+                                    std::string,    pout,
+                                    std::string,    output,
+                                    std::string,    input);
+};
+
+template <typename FImpl1, typename FImpl2>
+class TAmputate: public Module<AmputatePar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl1, 1);
+    FERM_TYPE_ALIASES(FImpl2, 2);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        std::vector<Complex>, Vamp,
+                                        ); 
+    };
+public:
+    // constructor
+    TAmputate(const std::string name);
+    // destructor
+    virtual ~TAmputate(void) {};
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    virtual SpinColourMatrix invertspincolmat(SpinColourMatrix &scmat);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(Amputate, ARG(TAmputate<FIMPL, FIMPL>), MNPR);
+
+/******************************************************************************
+ *                           TAmputate implementation                            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+TAmputate<FImpl1, FImpl2>::TAmputate(const std::string name)
+: Module<AmputatePar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TAmputate<FImpl1, FImpl2>::getInput(void)
+{
+    std::vector<std::string> input = {par().Sin, par().Sout, par().vertex};
+    
+    return input;
+}
+
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TAmputate<FImpl1, FImpl2>::getOutput(void)
+{
+    std::vector<std::string> output = {getName()};
+    
+    
+    return output;
+}
+
+// Invert spin colour matrix using Eigen
+template <typename Fimpl1, typename Fimpl2>
+SpinColourMatrix TAmputate<Fimpl1, Fimpl2>::invertspincolmat(SpinColourMatrix &scmat)
+{
+    Eigen::MatrixXcf scmat_2d(Ns*Nc,Ns*Nc);
+    for(int ic=0; ic<Nc; ic++){
+    for(int jc=0; jc<Nc; jc++){
+        for(int is=0; is<Ns; is++){
+        for(int js=0; js<Ns; js++){
+            scmat_2d(Ns*ic+is,Ns*jc+js) = scmat()(is,js)(ic,jc);
+        }}
+    }}      
+    Eigen::MatrixXcf scmat_2d_inv = scmat_2d.inverse();
+    SpinColourMatrix scmat_inv;
+    for(int ic=0; ic<Nc; ic++){
+    for(int jc=0; jc<Nc; jc++){
+        for(int is=0; is<Ns; is++){
+        for(int js=0; js<Ns; js++){
+            scmat_inv()(is,js)(ic,jc) = scmat_2d_inv(Ns*ic+is,Ns*jc+js);
+        }}
+    }}      
+    return scmat_inv;
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TAmputate<FImpl1, FImpl2>::execute(void)
+{
+    LOG(Message) << "Computing bilinear amputations '" << getName() << "' using"
+                 << " momentum '" << par().Sin << "' and '" << par().Sout << "'"
+                 << std::endl;
+    BinaryWriter                    writer(par().output);
+    PropagatorField1                &Sin = *env().template getObject<PropagatorField1>(par().Sin); //Do these have the phases taken into account?? Don't think so. FIX
+    PropagatorField2                &Sout = *env().template getObject<PropagatorField2>(par().Sout);
+    std::vector<int>                pin  = strToVec<int>(par().pin), pout = strToVec<int>(par().pout);
+    std::vector<Real>               latt_size(pin.begin(), pin.end()); 
+    LatticeComplex                  pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid());
+    LOG(Message) << "Propagators set up " << std::endl;
+    std::vector<SpinColourMatrix>   vertex; // Let's read from file here
+    Gamma                           g5(Gamma::Algebra::Gamma5);
+    Result                          result;
+    LOG(Message) << "reading file - "  << par().input << std::endl;
+    BinaryReader                    reader(par().input); 
+    Complex                         Ci(0.0,1.0);
+
+    std::string svertex;
+    read(reader,"vertex", vertex);
+    LOG(Message) << "vertex read" << std::endl;
+
+    pdotxin=zero;
+    pdotxout=zero;
+    for (unsigned int mu = 0; mu < 4; ++mu)
+    {
+        Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
+        LatticeCoordinate(coor,mu);
+        pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor;
+        pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor;
+    }
+    Sin = Sin*exp(-Ci*pdotxin); //phase corrections
+    Sout = Sout*exp(-Ci*pdotxout);
+
+    SpinColourMatrix Sin_mom = sum(Sin);
+    SpinColourMatrix Sout_mom = sum(Sout);
+    LOG(Message) << "summed over lattice" << std::endl;
+   
+    LOG(Message) << "Lattice -> spincolourmatrix conversion" << std::endl;
+
+    SpinColourMatrix Sin_inv = invertspincolmat(Sin_mom);
+    SpinColourMatrix Sout_inv = invertspincolmat(Sout_mom);
+    LOG(Message) << "Inversions done" << std::endl;
+
+    result.Vamp.resize(Gamma::nGamma/2);
+    for( int mu=0; mu < Gamma::nGamma/2; mu++){
+        Gamma::Algebra gam = mu;
+        result.Vamp[mu] = 1/12.0*trace(adj(Gamma(mu*2+1))*g5*Sout_inv*g5*vertex[mu]*Sin_inv);
+        LOG(Message) << "Vamp[" << mu << "] - " << result.Vamp[mu] << std::endl;
+        }
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Amputate_hpp_

Hadrons/Modules/MNPR/Bilinear.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/Bilinear.cc
+
+Copyright (C) 2015-2018
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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 */
+#include <Hadrons/Modules/MNPR/Bilinear.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MNPR;
+
+template class Grid::Hadrons::MNPR::TBilinear<FIMPL,FIMPL>;
+

Hadrons/Modules/MNPR/Bilinear.hpp

@@ -0,0 +1,225 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/Bilinear.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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 */
+
+#ifndef Hadrons_Bilinear_hpp_
+#define Hadrons_Bilinear_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+//#include <Grid/qcd/utils/PropagatorUtils.h>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                                TBilinear                                       *
+        Performs bilinear contractions of the type tr[g5*adj(Sout)*g5*G*Sin]
+        Suitable for non exceptional momenta in Rome-Southampton NPR
+******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MNPR)
+
+class BilinearPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(BilinearPar,
+                                    std::string,    Sin,
+                                    std::string,    Sout,
+                                    std::string,    pin,
+                                    std::string,    pout,
+                                    std::string,    output);
+};
+
+template <typename FImpl1, typename FImpl2>
+class TBilinear: public Module<BilinearPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl1, 1);
+    FERM_TYPE_ALIASES(FImpl2, 2);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result, 
+                                        std::vector<SpinColourMatrix>, bilinear);
+    };
+public:
+    // constructor
+    TBilinear(const std::string name);
+    // destructor
+    virtual ~TBilinear(void) {};
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    //LatticeSpinColourMatrix PhaseProps(LatticeSpinColourMatrix S, std::vector<Real> p);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(Bilinear, ARG(TBilinear<FIMPL, FIMPL>), MNPR);
+
+/******************************************************************************
+ *                           TBilinear implementation                            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+TBilinear<FImpl1, FImpl2>::TBilinear(const std::string name)
+: Module<BilinearPar>(name)
+{}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TBilinear<FImpl1, FImpl2>::setup(void)
+{
+    //env().template registerLattice<LatticeSpinColourMatrix>(getName());
+    //env().template registerObject<SpinColourMatrix>(getName());
+}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TBilinear<FImpl1, FImpl2>::getInput(void)
+{
+    std::vector<std::string> input = {par().Sin, par().Sout};
+    
+    return input;
+}
+
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TBilinear<FImpl1, FImpl2>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+/*
+/////Phase propagators//////////////////////////
+template <typename FImpl1, typename FImpl2>
+LatticeSpinColourMatrix TBilinear<FImpl1, FImpl2>::PhaseProps(LatticeSpinColourMatrix S, std::vector<Real> p)
+{
+    GridBase *grid = S._grid;
+    LatticeComplex      pdotx(grid),  coor(grid);
+    std::vector<int>   latt_size = grid->_fdimensions; 
+    Complex             Ci(0.0,1.0);
+    pdotx=zero;
+    for (unsigned int mu = 0; mu < 4; ++mu)
+    {
+        Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
+        LatticeCoordinate(coor,mu);
+        pdotx = pdotx +(TwoPiL * p[mu]) * coor;
+    }
+    S = S*exp(-Ci*pdotx);
+    return S;
+}
+*/
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TBilinear<FImpl1, FImpl2>::execute(void)
+{
+/**************************************************************************
+
+Compute the bilinear vertex needed for the NPR.
+V(G) = sum_x  [ g5 * adj(S'(x,p2)) * g5 * G * S'(x,p1) ]_{si,sj,ci,cj}
+G is one of the 16 gamma vertices [I,gmu,g5,g5gmu,sig(mu,nu)]
+
+        * G
+       / \
+    p1/   \p2
+     /     \
+    /       \
+
+Returns a spin-colour matrix, with indices si,sj, ci,cj
+
+Conventions:
+p1 - incoming momenta
+p2 - outgoing momenta
+q = (p1-p2)
+**************************************************************************/
+
+    LOG(Message) << "Computing bilinear contractions '" << getName() << "' using"
+                 << " momentum '" << par().Sin << "' and '" << par().Sout << "'"
+                 << std::endl;
+     
+    BinaryWriter             writer(par().output);
+    
+
+    // Propogators
+    LatticeSpinColourMatrix     &Sin = *env().template getObject<LatticeSpinColourMatrix>(par().Sin);
+    LatticeSpinColourMatrix     &Sout = *env().template getObject<LatticeSpinColourMatrix>(par().Sout);
+    LatticeComplex              pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid());
+    // momentum on legs
+    std::vector<Real>           pin  = strToVec<Real>(par().pin), pout = strToVec<Real>(par().pout);
+    std::vector<Real>           latt_size(pin.begin(), pin.end()); 
+    //bilinears
+    LatticeSpinColourMatrix     bilinear_x(env().getGrid());
+    SpinColourMatrix            bilinear;
+    Gamma                       g5(Gamma::Algebra::Gamma5);
+    Result                      result;
+    Complex                     Ci(0.0,1.0);
+
+    //
+
+    pdotxin=zero;
+    pdotxout=zero;
+    for (unsigned int mu = 0; mu < 4; ++mu)
+    {
+        Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
+        LatticeCoordinate(coor,mu);
+        pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor;
+        pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor;
+    }
+    Sin = Sin*exp(-Ci*pdotxin); //phase corrections
+    Sout = Sout*exp(-Ci*pdotxout);
+    
+    ////Set up gamma vector//////////////////////////
+    std::vector<Gamma> gammavector;
+    for( int i=0; i<Gamma::nGamma; i++){
+        Gamma::Algebra gam = i;
+        gammavector.push_back(Gamma(gam));
+    }
+    result.bilinear.resize(Gamma::nGamma);
+    /////////////////////////////////////////////////
+    //LatticeSpinMatrix temp = g5*Sout;
+    ////////Form Vertex//////////////////////////////
+    for (int i=0; i < Gamma::nGamma; i++){
+        bilinear_x = g5*adj(Sout)*g5*gammavector[i]*Sin; 
+        result.bilinear[i] = sum(bilinear_x); //sum over lattice sites
+    }
+    //////////////////////////////////////////////////
+    write(writer, par().output, result.bilinear);
+    LOG(Message) << "Complete. Writing results to " << par().output << std:: endl;
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Bilinear_hpp_

Hadrons/Modules/MNPR/FourQuark.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/FourQuark.cc
+
+Copyright (C) 2015-2018
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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 */
+#include <Hadrons/Modules/MNPR/FourQuark.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MNPR;
+
+template class Grid::Hadrons::MNPR::TFourQuark<FIMPL,FIMPL>;
+

Hadrons/Modules/MNPR/FourQuark.hpp

@@ -0,0 +1,274 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/FourQuark.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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 */
+
+#ifndef Hadrons_FourQuark_hpp_
+#define Hadrons_FourQuark_hpp_
+
+#include <typeinfo>
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Grid/serialisation/Serialisation.h>
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                                TFourQuark                                       *
+        Performs fourquark contractions of the type tr[g5*adj(Sout)*g5*G*Sin]
+        Suitable for non exceptional momenta
+******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MNPR)
+
+class FourQuarkPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(FourQuarkPar,
+                                    std::string,    Sin, //need to make this a propogator type?
+                                    std::string,    Sout, //same
+                                    std::string,    pin,
+                                    std::string,    pout,
+                                    bool,           fullbasis,
+                                    std::string,    output);
+};
+
+template <typename FImpl1, typename FImpl2>
+class TFourQuark: public Module<FourQuarkPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl1, 1);
+    FERM_TYPE_ALIASES(FImpl2, 2);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        std::vector<SpinColourSpinColourMatrix>, fourquark);
+    };
+public:
+    // constructor
+    TFourQuark(const std::string name);
+    // destructor
+    virtual ~TFourQuark(void) {};
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void tensorprod(LatticeSpinColourSpinColourMatrix &lret, LatticeSpinColourMatrix a, LatticeSpinColourMatrix b);
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(FourQuark, ARG(TFourQuark<FIMPL, FIMPL>), MNPR);
+
+/******************************************************************************
+ *                           TFourQuark implementation                            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+TFourQuark<FImpl1, FImpl2>::TFourQuark(const std::string name)
+: Module<FourQuarkPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TFourQuark<FImpl1, FImpl2>::getInput(void)
+{
+    std::vector<std::string> input = {par().Sin, par().Sout};
+    
+    return input;
+}
+
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TFourQuark<FImpl1, FImpl2>::getOutput(void)
+{
+    std::vector<std::string> output = {getName()};
+    
+    return output;
+}
+
+
+template <typename FImpl1, typename FImpl2>
+void TFourQuark<FImpl1, FImpl2>::tensorprod(LatticeSpinColourSpinColourMatrix &lret, LatticeSpinColourMatrix a, LatticeSpinColourMatrix b)
+{
+#if 0
+            parallel_for(auto site=lret.begin();site<lret.end();site++) {
+                for (int si; si < 4; ++si){
+                for(int sj; sj <4; ++sj){
+                    for (int ci; ci < 3; ++ci){
+                    for (int cj; cj < 3; ++cj){
+                        for (int sk; sk < 4; ++sk){
+                        for(int sl; sl <4; ++sl){
+                            for (int ck; ck < 3; ++ck){
+                            for (int cl; cl < 3; ++cl){
+                        lret[site]()(si,sj)(ci,cj)(sk,sl)(ck,cl)=a[site]()(si,sj)(ci,cj)*b[site]()(sk,sl)(ck,cl);
+                            }}
+                        }}
+                    }}
+                }}
+        }
+#else 
+            // FIXME ; is there a general need for this construct ? In which case we should encapsulate the
+            //         below loops in a helper function.
+            //LOG(Message) << "sp co mat a is - " << a << std::endl;
+            //LOG(Message) << "sp co mat b is - " << b << std::endl;
+            parallel_for(auto site=lret.begin();site<lret.end();site++) {
+            vTComplex left;
+                for(int si=0; si < Ns; ++si){
+                for(int sj=0; sj < Ns; ++sj){
+                    for (int ci=0; ci < Nc; ++ci){
+                    for (int cj=0; cj < Nc; ++cj){
+                      //LOG(Message) << "si, sj, ci, cj -  " << si << ", " << sj  << ", "<< ci  << ", "<< cj << std::endl;
+                      left()()() = a[site]()(si,sj)(ci,cj);
+                      //LOG(Message) << left << std::endl;
+                      lret[site]()(si,sj)(ci,cj)=left()*b[site]();
+                    }}
+                }}
+            }
+#endif      
+}
+
+
+
+
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TFourQuark<FImpl1, FImpl2>::setup(void)
+{
+    envCreateLat(LatticeSpinColourMatrix, getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TFourQuark<FImpl1, FImpl2>::execute(void)
+{
+
+/*********************************************************************************
+
+TFourQuark : Creates the four quark vertex required for the NPR of four-quark ops
+
+V_{Gamma_1,Gamma_2} = sum_x [ ( g5 * adj(S'(x,p2)) * g5 * G1 * S'(x,p1) )_ci,cj;si,sj x ( g5 * adj(S'(x,p2)) * g5 * G2 S'(x,p1) )_ck,cl;sk,cl ]
+
+Create a bilinear vertex for G1 and G2  the spin and colour indices are kept free. Where there are 16 potential Gs.
+We then find the outer product of V1 and V2, keeping the spin and colour indices uncontracted
+Then this is summed over the lattice coordinate
+Result is a SpinColourSpinColourMatrix - with 4 colour and 4 spin indices. 
+We have up to 256 of these including the offdiag (G1 != G2).
+
+        \         /
+         \p1   p1/
+          \     /
+           \   /
+         G1 * * G2
+           /   \
+          /     \
+         /p2   p2\
+        /         \
+
+*********************************************************************************/
+
+
+
+
+    LOG(Message) << "Computing fourquark contractions '" << getName() << "' using"
+                 << " momentum '" << par().Sin << "' and '" << par().Sout << "'"
+                 << std::endl;
+    
+    BinaryWriter             writer(par().output);
+    
+    PropagatorField1                            &Sin = *env().template getObject<PropagatorField1>(par().Sin);
+    PropagatorField2                            &Sout = *env().template getObject<PropagatorField2>(par().Sout);
+    std::vector<Real>                           pin  = strToVec<Real>(par().pin), pout = strToVec<Real>(par().pout);
+    bool                                        fullbasis = par().fullbasis;
+    Gamma                                       g5(Gamma::Algebra::Gamma5);
+    Result                                      result;
+    std::vector<Real>                           latt_size(pin.begin(), pin.end());
+    LatticeComplex                              pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid());
+    LatticeSpinColourMatrix                     bilinear_mu(env().getGrid()), bilinear_nu(env().getGrid());
+    LatticeSpinColourSpinColourMatrix           lret(env().getGrid()); 
+    Complex                         Ci(0.0,1.0);
+
+    //Phase propagators
+    //Sin = Grid::QCD::PropUtils::PhaseProps(Sin,pin);
+    //Sout = Grid::QCD::PropUtils::PhaseProps(Sout,pout);
+    
+    //find p.x for in and out so phase can be accounted for in propagators
+    pdotxin=zero;
+    pdotxout=zero;
+    for (unsigned int mu = 0; mu < 4; ++mu)
+    {
+        Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
+        LatticeCoordinate(coor,mu);
+        pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor;
+        pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor;
+    }
+    Sin = Sin*exp(-Ci*pdotxin); //phase corrections
+    Sout = Sout*exp(-Ci*pdotxout);
+
+
+    //Set up Gammas 
+    std::vector<Gamma> gammavector;
+     for( int i=1; i<Gamma::nGamma; i+=2){
+         Gamma::Algebra gam = i;
+         gammavector.push_back(Gamma(gam));
+       }
+    
+    lret = zero;
+    if (fullbasis == true){ // all combinations of mu and nu
+        result.fourquark.resize(Gamma::nGamma/2*Gamma::nGamma/2);
+        for( int mu=0; mu<Gamma::nGamma/2; mu++){ 
+            bilinear_mu = g5*adj(Sout)*g5*gammavector[mu]*Sin;
+            for ( int nu=0; nu<Gamma::nGamma; nu++){
+                LatticeSpinColourMatrix     bilinear_nu(env().getGrid());
+                bilinear_nu = g5*adj(Sout)*g5*gammavector[nu]*Sin;
+                LOG(Message) << "bilinear_nu for nu = " << nu << " is - " << bilinear_mu << std::endl;
+                result.fourquark[mu*Gamma::nGamma/2 + nu] = zero;
+                tensorprod(lret,bilinear_mu,bilinear_nu);
+                result.fourquark[mu*Gamma::nGamma/2 + nu] = sum(lret);
+            }
+        }
+    } else {
+        result.fourquark.resize(Gamma::nGamma/2);
+        for ( int mu=0; mu<1; mu++){
+        //for( int mu=0; mu<Gamma::nGamma/2; mu++ ){
+            bilinear_mu = g5*adj(Sout)*g5*gammavector[mu]*Sin;
+            //LOG(Message) << "bilinear_mu for mu = " << mu << " is - " << bilinear_mu << std::endl;
+            result.fourquark[mu] = zero;
+            tensorprod(lret,bilinear_mu,bilinear_mu); //tensor outer product
+            result.fourquark[mu] = sum(lret);
+        }
+    }
+    write(writer, "fourquark", result.fourquark);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_FourQuark_hpp_

Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp

@@ -1,268 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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 */
-#ifndef Hadrons_MScalarSUN_TimeMomProbe_hpp_
-#define Hadrons_MScalarSUN_TimeMomProbe_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-#include <Hadrons/Modules/MScalarSUN/Utils.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *          n-point functions O(t,p)*tr(phi(t_1,p_1)*...*phi(t_n,p_n))        *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MScalarSUN)
-
-class TimeMomProbePar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(TimeMomProbePar,
-                                    std::string,              field,
-                                    std::vector<std::string>, op,
-                                    std::vector<std::vector<std::string>>, timeMom,
-                                    std::string,              output);
-};
-
-class TimeMomProbeResult: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(TimeMomProbeResult,
-                                    std::string,                   op,
-                                    std::vector<std::vector<int>>, timeMom,
-                                    std::vector<Complex>,          data);
-};
-
-template <typename SImpl>
-class TTimeMomProbe: public Module<TimeMomProbePar>
-{
-public:
-    typedef typename SImpl::Field                    Field;
-    typedef typename SImpl::SiteField::scalar_object Site;
-    typedef typename SImpl::ComplexField             ComplexField;
-    typedef          std::vector<Complex>            SlicedOp;
-public:
-    // constructor
-    TTimeMomProbe(const std::string name);
-    // destructor
-    virtual ~TTimeMomProbe(void) {};
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-private:
-    void vectorModulo(std::vector<int> &v);
-};
-
-MODULE_REGISTER_TMP(TimeMomProbeSU2, TTimeMomProbe<ScalarNxNAdjImplR<2>>, MScalarSUN);
-MODULE_REGISTER_TMP(TimeMomProbeSU3, TTimeMomProbe<ScalarNxNAdjImplR<3>>, MScalarSUN);
-MODULE_REGISTER_TMP(TimeMomProbeSU4, TTimeMomProbe<ScalarNxNAdjImplR<4>>, MScalarSUN);
-MODULE_REGISTER_TMP(TimeMomProbeSU5, TTimeMomProbe<ScalarNxNAdjImplR<5>>, MScalarSUN);
-MODULE_REGISTER_TMP(TimeMomProbeSU6, TTimeMomProbe<ScalarNxNAdjImplR<6>>, MScalarSUN);
-
-/******************************************************************************
- *                        TTimeMomProbe implementation                        *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-template <typename SImpl>
-TTimeMomProbe<SImpl>::TTimeMomProbe(const std::string name)
-: Module<TimeMomProbePar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-template <typename SImpl>
-std::vector<std::string> TTimeMomProbe<SImpl>::getInput(void)
-{
-    std::vector<std::string> in = par().op;
-    
-    in.push_back(par().field);
-
-    return in;
-}
-
-template <typename SImpl>
-std::vector<std::string> TTimeMomProbe<SImpl>::getOutput(void)
-{
-    std::vector<std::string> out;
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-template <typename SImpl>
-void TTimeMomProbe<SImpl>::setup(void)
-{
-    envTmpLat(ComplexField, "ftBuf");
-    envTmpLat(Field, "ftMatBuf");
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-// NB: time is direction 0
-template <typename SImpl>
-void TTimeMomProbe<SImpl>::vectorModulo(std::vector<int> &v)
-{
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    {
-        auto d = env().getDim(mu);
-        v[mu] = ((v[mu] % d) + d) % d;
-    }
-}
-
-template <typename SImpl>
-void TTimeMomProbe<SImpl>::execute(void)
-{
-    const unsigned int                           nd = env().getNd();
-    const unsigned int                           nt = env().getDim(0);
-    double                                       partVol = 1.;
-    std::set<std::vector<int>>                   timeMomSet;
-    std::vector<std::vector<std::vector<int>>>   timeMom;
-    std::vector<std::vector<int>>                transferMom;
-    FFT                                          fft(envGetGrid(Field));
-    std::vector<int>                             dMask(nd, 1);
-    std::vector<TimeMomProbeResult>              result;
-    std::map<std::string, std::vector<SlicedOp>> slicedOp;
-    std::vector<SlicedOp>                        slicedProbe;
-    auto                                         &phi = envGet(Field, par().field);
-
-    envGetTmp(ComplexField, ftBuf);
-    envGetTmp(Field, ftMatBuf);
-    dMask[0] = 0;
-    for (unsigned int mu = 1; mu < nd; ++mu)
-    {
-        partVol *= env().getDim(mu);
-    }
-    timeMom.resize(par().timeMom.size());
-    for (unsigned int p = 0; p < timeMom.size(); ++p)
-    {
-        for (auto &tms: par().timeMom[p])
-        {
-            std::vector<int> tm = strToVec<int>(tms);
-            
-            timeMom[p].push_back(tm);
-            timeMomSet.insert(tm);
-        }
-        transferMom.push_back(std::vector<int>(nd - 1, 0));
-        for (auto &tm: timeMom[p])
-        {
-            for (unsigned int j = 1; j < nd; ++j)
-            {
-                transferMom[p][j - 1] -= tm[j];
-            }
-        }
-        LOG(Message) << "Probe " << p << " (" << timeMom[p].size() << " points) : " << std::endl;
-        LOG(Message) << "  phi(t_i, p_i) for (t_i, p_i) in " << timeMom[p] << std::endl;
-        LOG(Message) << "  operator with momentum " << transferMom[p] << std::endl;
-    }
-    LOG(Message) << "FFT: field '" << par().field << "'" << std::endl;
-    fft.FFT_dim_mask(ftMatBuf, phi, dMask, FFT::forward);
-    slicedProbe.resize(timeMom.size());
-    for (unsigned int p = 0; p < timeMom.size(); ++p)
-    {
-        std::vector<int> qt;
-
-        LOG(Message) << "Making probe " << p << std::endl;
-        slicedProbe[p].resize(nt);
-        for (unsigned int t = 0; t < nt; ++t)
-        {
-            Site acc;
-            
-            for (unsigned int i = 0; i < timeMom[p].size(); ++i)
-            {
-                Site buf;
-
-                qt     = timeMom[p][i];
-                qt[0] += t;
-                vectorModulo(qt);
-                peekSite(buf, ftMatBuf, qt);
-                if (i == 0)
-                {
-                    acc = buf;
-                }
-                else
-                {
-                    acc *= buf;
-                }
-            }
-            slicedProbe[p][t] = TensorRemove(trace(acc));
-        }
-        //std::cout << slicedProbe[p]<< std::endl;
-    }
-    for (auto &o: par().op)
-    {
-        auto &op = envGet(ComplexField, o);
-
-        slicedOp[o].resize(transferMom.size());
-        LOG(Message) << "FFT: operator '" << o << "'" << std::endl;
-        fft.FFT_dim_mask(ftBuf, op, dMask, FFT::forward);
-        //std::cout << ftBuf << std::endl;
-        for (unsigned int p = 0; p < transferMom.size(); ++p)
-        {
-            std::vector<int> qt(nd, 0);
-
-            for (unsigned int j = 1; j < nd; ++j)
-            {
-                qt[j] = transferMom[p][j - 1];
-            }
-            slicedOp[o][p].resize(nt);
-            for (unsigned int t = 0; t < nt; ++t)
-            {
-                TComplex buf;
-
-                qt[0] = t;
-                vectorModulo(qt);
-                peekSite(buf, ftBuf, qt);
-                slicedOp[o][p][t] = TensorRemove(buf);
-            }
-            //std::cout << ftBuf << std::endl;
-            //std::cout << slicedOp[o][p] << std::endl;
-        }
-    }
-    LOG(Message) << "Making correlators" << std::endl;
-    for (auto &o: par().op)
-    for (unsigned int p = 0; p < timeMom.size(); ++p)
-    {
-        TimeMomProbeResult r;
-
-        LOG(Message) << "  <" << o << " probe_" << p << ">" << std::endl;
-        r.op      = o;
-        r.timeMom = timeMom[p];
-        r.data    = makeTwoPoint(slicedOp[o][p], slicedProbe[p], 1./partVol);
-        result.push_back(r);
-    }
-    saveResult(par().output, "timemomprobe", result);
-}
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MScalarSUN_TimeMomProbe_hpp_

Hadrons/Modules/MScalarSUN/TrMag.hpp

@@ -124,7 +124,8 @@ void TTrMag<SImpl>::execute(void)
     std::vector<TrMagResult> result;
     auto                     &phi = envGet(Field, par().field);
 
-    auto m2 = sum(phi), mn = m2;
+    auto m2 = sum(phi);
+    auto mn = m2;
 
     m2 = -m2*m2;
     mn = 1.;

Hadrons/Modules/MScalarSUN/Utils.hpp

@@ -103,7 +103,7 @@ std::vector<Complex> makeTwoPoint(const std::vector<SinkSite>   &sink,
     {
         for (unsigned int t  = 0; t < nt; ++t)
         {
-            res[dt] += trace(sink[(t+dt)%nt]*source[t]);
+            res[dt] += trace(sink[(t+dt)%nt]*adj(source[t]));
         }
         res[dt] *= factor/static_cast<double>(nt);
     }

Hadrons/Modules/MSolver/A2AVectors.cc

@@ -32,5 +32,5 @@ using namespace Grid;
 using namespace Hadrons;
 using namespace MSolver;
 
-template class Grid::Hadrons::MSolver::TA2AVectors<FIMPL, FermionEigenPack<FIMPL>>;
-template class Grid::Hadrons::MSolver::TA2AVectors<ZFIMPL, FermionEigenPack<ZFIMPL>>;
+template class Grid::Hadrons::MSolver::TA2AVectors<FIMPL, BaseFermionEigenPack<FIMPL>>;
+template class Grid::Hadrons::MSolver::TA2AVectors<ZFIMPL, BaseFermionEigenPack<ZFIMPL>>;

Hadrons/Modules/MSolver/A2AVectors.hpp

@@ -79,9 +79,9 @@ private:
 };
 
 MODULE_REGISTER_TMP(A2AVectors, 
-    ARG(TA2AVectors<FIMPL, FermionEigenPack<FIMPL>>), MSolver);
+    ARG(TA2AVectors<FIMPL, BaseFermionEigenPack<FIMPL>>), MSolver);
 MODULE_REGISTER_TMP(ZA2AVectors, 
-    ARG(TA2AVectors<ZFIMPL, FermionEigenPack<ZFIMPL>>), MSolver);
+    ARG(TA2AVectors<ZFIMPL, BaseFermionEigenPack<ZFIMPL>>), MSolver);
 
 /******************************************************************************
  *                       TA2AVectors implementation                           *

Hadrons/Modules/MSolver/Guesser.hpp

@@ -39,7 +39,7 @@ std::shared_ptr<LinearFunction<typename FImpl::FermionField>>
 makeGuesser(const std::string epackName)
 {
     typedef typename FImpl::FermionField                  FermionField;
-    typedef FermionEigenPack<FImpl>                       EPack;
+    typedef BaseFermionEigenPack<FImpl>                   EPack;
     typedef CoarseFermionEigenPack<FImpl, nBasis>         CoarseEPack;
     typedef DeflatedGuesser<FermionField>                 FineGuesser;
     typedef LocalCoherenceDeflatedGuesser<

Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp

@@ -63,7 +63,7 @@ public:
     typedef LocalCoherenceLanczos<typename FImpl::SiteSpinor, 
                                   typename FImpl::SiteComplex, 
                                   nBasis>                LCL;
-    typedef FermionEigenPack<FImpl>                      BasePack;
+    typedef BaseFermionEigenPack<FImpl>                  BasePack;
     typedef CoarseFermionEigenPack<FImpl, nBasis>        CoarsePack;
     typedef HADRONS_DEFAULT_SCHUR_OP<FMat, FermionField> SchurFMat;
 public:

Hadrons/Modules/MSource/Momentum.cc

@@ -0,0 +1,36 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MSource/Momentum.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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 */
+#include <Hadrons/Modules/MSource/Momentum.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MSource;
+
+template class Grid::Hadrons::MSource::TMomentum<FIMPL>;
+

Hadrons/Modules/MSource/Momentum.hpp

@@ -0,0 +1,149 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MSource/Momentum.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+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 */
+#ifndef Hadrons_Momentum_hpp_
+#define Hadrons_Momentum_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/* 
+Plane Wave source
+-----------------
+src_x = e^i2pi/L * p *position
+*/
+
+/******************************************************************************
+ *                          Plane Wave source                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class MomentumPar: Serializable
+{
+public:
+//What is meant by serializable in this context
+GRID_SERIALIZABLE_CLASS_MEMBERS(MomentumPar,
+std::string, mom);
+};
+
+
+template <typename FImpl>
+class TMomentum: public Module<MomentumPar>
+{
+public:
+FERM_TYPE_ALIASES(FImpl,);
+public:
+// constructor
+TMomentum(const std::string name);
+// destructor
+virtual ~TMomentum(void) {};
+// dependency relation
+virtual std::vector<std::string> getInput(void);
+virtual std::vector<std::string> getOutput(void);
+// setup
+virtual void setup(void);
+// execution
+virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(Momentum, TMomentum<FIMPL>, MSource);
+//MODULE_REGISTER_NS(Momentum, TMomentum, MSource);
+
+/******************************************************************************
+*                       TMomentum template implementation                     *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TMomentum<FImpl>::TMomentum(const std::string name)
+: Module<MomentumPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TMomentum<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TMomentum<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    return out;
+}
+
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TMomentum<FImpl>::setup(void)
+{
+    envCreateLat(PropagatorField, getName());
+}
+
+
+//execution//////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TMomentum<FImpl>::execute(void)
+{
+    LOG(Message) << "Generating planewave momentum source with momentum " << par().mom << std::endl;
+    //what does this env do?
+    PropagatorField &src = envGet(PropagatorField, getName());
+    Lattice<iScalar<vInteger>> t(env().getGrid());
+    LatticeComplex             C(env().getGrid()), coor(env().getGrid());
+    std::vector<Real>          p;
+    std::vector<Real> latt_size(GridDefaultLatt().begin(), GridDefaultLatt().end()); 
+    Complex                    i(0.0,1.0);
+
+    LOG(Message) << " " << std::endl;
+    //get the momentum from parameters
+    p  = strToVec<Real>(par().mom);
+    C = zero;
+    LOG(Message) << "momentum converted from string - " << std::to_string(p[0]) <<std::to_string(p[1]) <<std::to_string(p[2]) <<   std::to_string(p[3]) << std::endl;
+    for(int mu=0;mu<4;mu++){
+    Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
+    LatticeCoordinate(coor,mu);
+    C = C +(TwoPiL * p[mu]) * coor;
+    }
+    C = exp(C*i);
+    LOG(Message) << "exponential of pdotx taken " << std::endl;
+    src = src + C;
+    LOG(Message) << "source created" << std::endl;
+
+}
+
+
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Momentum_hpp_

Hadrons/Modules/MSource/Point.hpp

@@ -126,6 +126,11 @@ void TPoint<FImpl>::execute(void)
     auto             &src     = envGet(PropagatorField, getName());
     SitePropagator   id;
     
+    if (position.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "position has " + std::to_string(position.size())
+                      + " components (must have " + std::to_string(env().getNd()) + ")");
+    }
     id  = 1.;
     src = zero;
     pokeSite(id, src, position);

Hadrons/TimerArray.cc

@@ -0,0 +1,126 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/TimerArray.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+#include <Hadrons/TimerArray.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+void TimerArray::startTimer(const std::string &name)
+{
+    if (!name.empty())
+    {
+        timer_[name].Start();
+    }
+}
+
+GridTime TimerArray::getTimer(const std::string &name)
+{
+    GridTime t;
+    
+    if (!name.empty())
+    {
+        try
+        {
+            bool running = timer_.at(name).isRunning();
+
+            if (running) stopTimer(name);
+            t = timer_.at(name).Elapsed();
+            if (running) startTimer(name);
+        }
+        catch (std::out_of_range &)
+        {
+            t = GridTime::zero();
+        }
+    }
+    else
+    {
+        t = GridTime::zero();
+    }
+
+    return t;
+}
+
+double TimerArray::getDTimer(const std::string &name)
+{
+    return static_cast<double>(getTimer(name).count());
+}
+
+void TimerArray::startCurrentTimer(const std::string &name)
+{
+    if (!name.empty())
+    {
+        stopCurrentTimer();
+        startTimer(name);
+        currentTimer_ = name;
+    }
+}
+
+void TimerArray::stopTimer(const std::string &name)
+{
+    if (timer_.at(name).isRunning())
+    {
+        timer_.at(name).Stop();
+    }
+}
+
+void TimerArray::stopCurrentTimer(void)
+{
+    if (!currentTimer_.empty())
+    {
+        stopTimer(currentTimer_);
+        currentTimer_ = "";
+    }
+}
+
+void TimerArray::stopAllTimers(void)
+{
+    for (auto &t: timer_)
+    {
+        stopTimer(t.first);
+    }
+    currentTimer_ = "";
+}
+
+void TimerArray::resetTimers(void)
+{
+    timer_.clear();
+    currentTimer_ = "";
+}
+
+std::map<std::string, GridTime> TimerArray::getTimings(void)
+{
+    std::map<std::string, GridTime> timing;
+
+    for (auto &t: timer_)
+    {
+        timing[t.first] = t.second.Elapsed();
+    }
+
+    return timing;
+}

Hadrons/TimerArray.hpp

@@ -0,0 +1,56 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/TimerArray.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+#ifndef Hadrons_TimerArray_hpp_
+#define Hadrons_TimerArray_hpp_
+
+#include <Hadrons/Global.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+class TimerArray
+{
+public:
+    TimerArray(void) = default;
+    virtual ~TimerArray(void) = default;
+    void                            startTimer(const std::string &name);
+    GridTime                        getTimer(const std::string &name);
+    double                          getDTimer(const std::string &name);
+    void                            startCurrentTimer(const std::string &name);
+    void                            stopTimer(const std::string &name);
+    void                            stopCurrentTimer(void);
+    void                            stopAllTimers(void);
+    void                            resetTimers(void);
+    std::map<std::string, GridTime> getTimings(void);
+private:
+    std::string                          currentTimer_;
+    std::map<std::string, GridStopWatch> timer_; 
+};
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_TimerArray_hpp_

Hadrons/Utilities/EigenPackCast.cc

@@ -1,3 +1,30 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Utilities/EigenPackCast.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
 #include <Hadrons/EigenPack.hpp>
 #include <Hadrons/Environment.hpp>
 
@@ -59,8 +86,12 @@ void convert(const std::string outFilename, const std::string inFilename,
         }
     }
 
-    FOut bufOut(gOut);
-    FIn  bufIn(gIn), testIn(gIn);
+    FOut         bufOut(gOut);
+    FIn          bufIn(gIn), testIn(gIn);
+    ScidacWriter binWriter(gOut->IsBoss());
+    ScidacReader binReader;
+    PackRecord   record;
+    RealD        eval;
 
     LOG(Message) << "==== EIGENPACK CONVERSION" << std::endl;
     LOG(Message) << "Lattice       : " << gIn->GlobalDimensions() << std::endl;
@@ -75,10 +106,6 @@ void convert(const std::string outFilename, const std::string inFilename,
     {
         for(unsigned int k = 0; k < size; ++k)
         {
-            ScidacWriter binWriter(gOut->IsBoss());
-            ScidacReader binReader;
-            PackRecord   record;
-            VecRecord    vecRecord;
             std::string  outV = outFilename + "/v" + std::to_string(k) + ".bin";
             std::string  inV  = inFilename + "/v" + std::to_string(k) + ".bin";
 
@@ -88,40 +115,25 @@ void convert(const std::string outFilename, const std::string inFilename,
             makeFileDir(outV, gOut);
             binWriter.open(outV);
             binReader.open(inV);
-            EPIn::readHeader(record, binReader);
-            EPOut::writeHeader(binWriter, record);
-            EPIn::readElement(bufIn, vecRecord, binReader);
-            precisionChange(bufOut, bufIn);
-            precisionChange(testIn, bufOut);
-            testIn -= bufIn;
-            LOG(Message) << "Diff norm^2: " << norm2(testIn) << std::endl;
-            EPOut::writeElement(binWriter, bufOut, vecRecord);
+            EigenPackIo::readHeader(record, binReader);
+            EigenPackIo::writeHeader(binWriter, record);
+            EigenPackIo::readElement<FIn>(bufIn, eval, k, binReader);
+            EigenPackIo::writeElement<FIn, FOut>(binWriter, bufIn, eval, k, &bufOut, &testIn);
             binWriter.close();
             binReader.close();
         }
     }
     else
     {
-        ScidacWriter binWriter(gOut->IsBoss());
-        ScidacReader binReader;
-        PackRecord   record;
-
         makeFileDir(outFilename, gOut);
         binWriter.open(outFilename);
         binReader.open(inFilename);
-        EPIn::readHeader(record, binReader);
-        EPOut::writeHeader(binWriter, record);
+        EigenPackIo::readHeader(record, binReader);
+        EigenPackIo::writeHeader(binWriter, record);
         for(unsigned int k = 0; k < size; ++k)
         {
-            VecRecord vecRecord;
-
-            LOG(Message) << "==== Converting vector " << k << std::endl;
-            EPIn::readElement(bufIn, vecRecord, binReader);
-            precisionChange(bufOut, bufIn);
-            precisionChange(testIn, bufOut);
-            testIn -= bufIn;
-            LOG(Message) << "Diff norm^2: " << norm2(testIn) << std::endl;
-            EPOut::writeElement(binWriter, bufOut, vecRecord);
+            EigenPackIo::readElement<FIn>(bufIn, eval, k, binReader);
+            EigenPackIo::writeElement<FIn, FOut>(binWriter, bufIn, eval, k, &bufOut, &testIn);
         }
         binWriter.close();
         binReader.close();

Hadrons/Utilities/HadronsXmlRun.cc

@@ -2,7 +2,7 @@
 
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: Hadrons/HadronsXmlRun.cc
+Source file: Hadrons/Utilities/HadronsXmlRun.cc
 
 Copyright (C) 2015-2018
 

Hadrons/Utilities/Makefile.am

@@ -1,10 +1,8 @@
-AM_LDFLAGS += -L../../Hadrons
-
 bin_PROGRAMS = HadronsXmlRun HadronsFermionEP64To32
 
 HadronsXmlRun_SOURCES = HadronsXmlRun.cc
-HadronsXmlRun_LDADD   = -lHadrons -lGrid
+HadronsXmlRun_LDADD   = ../libHadrons.a ../../Grid/libGrid.a
 
 HadronsFermionEP64To32_SOURCES  = EigenPackCast.cc
 HadronsFermionEP64To32_CXXFLAGS = $(AM_CXXFLAGS) -DFIN=WilsonImplD::FermionField -DFOUT=WilsonImplF::FermionField
-HadronsFermionEP64To32_LDADD    = -lHadrons -lGrid
+HadronsFermionEP64To32_LDADD    = ../libHadrons.a ../../Grid/libGrid.a

Hadrons/modules.inc

@@ -4,12 +4,14 @@ modules_cc =\
   Modules/MContraction/Meson.cc \
   Modules/MContraction/WeakNeutral4ptDisc.cc \
   Modules/MContraction/WeakHamiltonianNonEye.cc \
+  Modules/MContraction/A2AAslashField.cc \
   Modules/MContraction/WardIdentity.cc \
   Modules/MContraction/A2AMesonField.cc \
   Modules/MContraction/DiscLoop.cc \
   Modules/MContraction/Gamma3pt.cc \
   Modules/MFermion/FreeProp.cc \
   Modules/MFermion/GaugeProp.cc \
+  Modules/MSource/Momentum.cc \
   Modules/MSource/Point.cc \
   Modules/MSource/Wall.cc \
   Modules/MSource/SeqConserved.cc \
@@ -27,6 +29,7 @@ modules_cc =\
   Modules/MGauge/StochEm.cc \
   Modules/MGauge/Random.cc \
   Modules/MGauge/FundtoHirep.cc \
+  Modules/MGauge/GaugeFix.cc \
   Modules/MNoise/TimeDilutedSpinColorDiagonal.cc \
   Modules/MUtilities/RandomVectors.cc \
   Modules/MUtilities/TestSeqGamma.cc \
@@ -37,6 +40,9 @@ modules_cc =\
   Modules/MScalar/VPCounterTerms.cc \
   Modules/MScalar/ChargedProp.cc \
   Modules/MScalar/ScalarVP.cc \
+  Modules/MNPR/Amputate.cc \
+  Modules/MNPR/Bilinear.cc \
+  Modules/MNPR/FourQuark.cc \
   Modules/MAction/Wilson.cc \
   Modules/MAction/MobiusDWF.cc \
   Modules/MAction/ZMobiusDWF.cc \
@@ -45,7 +51,6 @@ modules_cc =\
   Modules/MAction/ScaledDWF.cc \
   Modules/MScalarSUN/TrPhi.cc \
   Modules/MScalarSUN/Grad.cc \
-  Modules/MScalarSUN/TimeMomProbe.cc \
   Modules/MScalarSUN/TrMag.cc \
   Modules/MScalarSUN/TrKinetic.cc \
   Modules/MScalarSUN/EMT.cc \
@@ -63,8 +68,8 @@ modules_cc =\
 
 modules_hpp =\
   Modules/MContraction/Baryon.hpp \
+  Modules/MContraction/A2AAslashField.hpp \
   Modules/MContraction/A2AMesonField.hpp \
-  Modules/MContraction/A2AMesonFieldKernels.hpp \
   Modules/MContraction/Meson.hpp \
   Modules/MContraction/WeakHamiltonian.hpp \
   Modules/MContraction/WeakHamiltonianNonEye.hpp \
@@ -80,6 +85,7 @@ modules_hpp =\
   Modules/MSource/Wall.hpp \
   Modules/MSource/Z2.hpp \
   Modules/MSource/SeqConserved.hpp \
+  Modules/MSource/Momentum.hpp \
   Modules/MSink/Smear.hpp \
   Modules/MSink/Point.hpp \
   Modules/MSolver/MixedPrecisionRBPrecCG.hpp \
@@ -91,6 +97,7 @@ modules_hpp =\
   Modules/MGauge/StoutSmearing.hpp \
   Modules/MGauge/Unit.hpp \
   Modules/MGauge/Random.hpp \
+  Modules/MGauge/GaugeFix.hpp \
   Modules/MGauge/FundtoHirep.hpp \
   Modules/MGauge/StochEm.hpp \
   Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp \
@@ -104,6 +111,9 @@ modules_hpp =\
   Modules/MScalar/ScalarVP.hpp \
   Modules/MScalar/Scalar.hpp \
   Modules/MScalar/ChargedProp.hpp \
+  Modules/MNPR/Bilinear.hpp \
+  Modules/MNPR/Amputate.hpp \
+  Modules/MNPR/FourQuark.hpp \
   Modules/MAction/DWF.hpp \
   Modules/MAction/MobiusDWF.hpp \
   Modules/MAction/Wilson.hpp \
@@ -114,7 +124,6 @@ modules_hpp =\
   Modules/MScalarSUN/TwoPointNPR.hpp \
   Modules/MScalarSUN/ShiftProbe.hpp \
   Modules/MScalarSUN/Div.hpp \
-  Modules/MScalarSUN/TimeMomProbe.hpp \
   Modules/MScalarSUN/TrMag.hpp \
   Modules/MScalarSUN/EMT.hpp \
   Modules/MScalarSUN/TwoPoint.hpp \

benchmarks/Benchmark_IO.cc

@@ -1,108 +1,48 @@
-#include <Grid/Grid.h>
-#ifdef HAVE_LIME
 
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
+#include "Benchmark_IO.hpp"
 
-#define MSG cout << GridLogMessage
-#define SEP \
-"============================================================================="
 #ifndef BENCH_IO_LMAX
 #define BENCH_IO_LMAX 40
 #endif
 
-typedef function<void(const string, LatticeFermion &)> WriterFn;
-typedef function<void(LatticeFermion &, const string)> ReaderFn;
+using namespace Grid;
+using namespace QCD;
 
-string filestem(const int l)
+std::string filestem(const int l)
 {
-  return "iobench_l" + to_string(l);
-}
-
-void limeWrite(const string filestem, LatticeFermion &vec)
-{
-  emptyUserRecord record;
-  ScidacWriter    binWriter(vec._grid->IsBoss());
-
-  binWriter.open(filestem + ".bin");
-  binWriter.writeScidacFieldRecord(vec, record);
-  binWriter.close();
-}
-
-void limeRead(LatticeFermion &vec, const string filestem)
-{
-  emptyUserRecord record;
-  ScidacReader    binReader;
-
-  binReader.open(filestem + ".bin");
-  binReader.readScidacFieldRecord(vec, record);
-  binReader.close();
-}
-
-void writeBenchmark(const int l, const WriterFn &write)
-{
-  auto                      mpi  = GridDefaultMpi();
-  auto                      simd = GridDefaultSimd(Nd, vComplex::Nsimd());
-  vector<int>               latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]};
-  unique_ptr<GridCartesian> gPt(SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi));
-  GridCartesian             *g = gPt.get();
-  GridParallelRNG           rng(g);
-  LatticeFermion            vec(g);
-  emptyUserRecord           record;
-  ScidacWriter              binWriter(g->IsBoss());
-
-  cout << "-- Local volume " << l << "^4" << endl;
-  random(rng, vec);
-  write(filestem(l), vec);
-}
-
-void readBenchmark(const int l, const ReaderFn &read)
-{
-  auto                      mpi  = GridDefaultMpi();
-  auto                      simd = GridDefaultSimd(Nd, vComplex::Nsimd());
-  vector<int>               latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]};
-  unique_ptr<GridCartesian> gPt(SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi));
-  GridCartesian             *g = gPt.get();
-  LatticeFermion            vec(g);
-  emptyUserRecord           record;
-  ScidacReader              binReader;
-
-  cout << "-- Local volume " << l << "^4" << endl;
-  read(vec, filestem(l));
+  return "iobench_l" + std::to_string(l);
 }
 
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  auto simd = GridDefaultSimd(Nd,vComplex::Nsimd());
-  auto mpi  = GridDefaultMpi();
-
   int64_t threads = GridThread::GetThreads();
-  MSG << "Grid is setup to use " << threads << " threads" << endl;
-  MSG << SEP << endl;
-  MSG << "Benchmark Lime write" << endl;
-  MSG << SEP << endl;
+  MSG << "Grid is setup to use " << threads << " threads" << std::endl;
+  MSG << SEP << std::endl;
+  MSG << "Benchmark Lime write" << std::endl;
+  MSG << SEP << std::endl;
   for (int l = 4; l <= BENCH_IO_LMAX; l += 2)
   {
-    writeBenchmark(l, limeWrite);
+    auto             mpi  = GridDefaultMpi();
+    std::vector<int> latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]};
+
+    std::cout << "-- Local volume " << l << "^4" << std::endl;
+    writeBenchmark<LatticeFermion>(latt, filestem(l), limeWrite<LatticeFermion>);
   }
 
-  MSG << "Benchmark Lime read" << endl;
-  MSG << SEP << endl;
+  MSG << "Benchmark Lime read" << std::endl;
+  MSG << SEP << std::endl;
   for (int l = 4; l <= BENCH_IO_LMAX; l += 2)
   {
-    readBenchmark(l, limeRead);
+    auto             mpi  = GridDefaultMpi();
+    std::vector<int> latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]};
+
+    std::cout << "-- Local volume " << l << "^4" << std::endl;
+    readBenchmark<LatticeFermion>(latt, filestem(l), limeRead<LatticeFermion>);
   }
 
   Grid_finalize();
 
   return EXIT_SUCCESS;
 }
-#else
-int main (int argc, char ** argv)
-{
-  return EXIT_SUCCESS;
-}
-#endif

benchmarks/Benchmark_IO.hpp

@@ -0,0 +1,107 @@
+#ifndef Benchmark_IO_hpp_
+#define Benchmark_IO_hpp_
+
+#include <Grid/Grid.h>
+
+#define MSG std::cout << GridLogMessage
+#define SEP \
+"============================================================================="
+
+namespace Grid {
+
+template <typename Field>
+using WriterFn = std::function<void(const std::string, Field &)> ;
+template <typename Field>
+using ReaderFn = std::function<void(Field &, const std::string)>;
+
+template <typename Field>
+void limeWrite(const std::string filestem, Field &vec)
+{
+  emptyUserRecord   record;
+  QCD::ScidacWriter binWriter(vec._grid->IsBoss());
+
+  binWriter.open(filestem + ".bin");
+  binWriter.writeScidacFieldRecord(vec, record);
+  binWriter.close();
+}
+
+template <typename Field>
+void limeRead(Field &vec, const std::string filestem)
+{
+  emptyUserRecord   record;
+  QCD::ScidacReader binReader;
+
+  binReader.open(filestem + ".bin");
+  binReader.readScidacFieldRecord(vec, record);
+  binReader.close();
+}
+
+inline void makeGrid(std::shared_ptr<GridBase> &gPt, 
+                     const std::shared_ptr<GridCartesian> &gBasePt,
+                     const unsigned int Ls = 1, const bool rb = false)
+{
+  if (rb)
+  {
+    if (Ls > 1)
+    {
+      gPt.reset(QCD::SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, gBasePt.get()));
+    }
+    else
+    {
+      gPt.reset(QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(gBasePt.get()));
+    }
+  }
+  else
+  {
+    if (Ls > 1)
+    {
+        gPt.reset(QCD::SpaceTimeGrid::makeFiveDimGrid(Ls, gBasePt.get()));
+    }
+    else
+    {
+        gPt = gBasePt;
+    }
+  }
+}
+
+template <typename Field>
+void writeBenchmark(const std::vector<int> &latt, const std::string filename,
+                    const WriterFn<Field> &write, 
+                    const unsigned int Ls = 1, const bool rb = false)
+{
+  auto                           mpi  = GridDefaultMpi();
+  auto                           simd = GridDefaultSimd(latt.size(), Field::vector_type::Nsimd());
+  std::shared_ptr<GridCartesian> gBasePt(QCD::SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi));
+  std::shared_ptr<GridBase>      gPt;
+
+  makeGrid(gPt, gBasePt, Ls, rb);
+
+  GridBase                       *g = gPt.get();
+  GridParallelRNG                rng(g);
+  Field                          vec(g);
+
+  random(rng, vec);
+  write(filename, vec);
+}
+
+template <typename Field>
+void readBenchmark(const std::vector<int> &latt, const std::string filename,
+                   const ReaderFn<Field> &read, 
+                   const unsigned int Ls = 1, const bool rb = false)
+{
+  auto                           mpi  = GridDefaultMpi();
+  auto                           simd = GridDefaultSimd(latt.size(), Field::vector_type::Nsimd());
+  std::shared_ptr<GridCartesian> gBasePt(QCD::SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi));
+  std::shared_ptr<GridBase>      gPt;
+
+  makeGrid(gPt, gBasePt, Ls, rb);
+
+  GridBase                       *g = gPt.get();
+  Field                          vec(g);
+
+  read(vec, filename);
+}
+
+}
+
+#endif // Benchmark_IO_hpp_

benchmarks/Benchmark_IO_vs_dir.cc

@@ -0,0 +1,63 @@
+#include "Benchmark_IO.hpp"
+
+#define MSG std::cout << GridLogMessage
+#define SEP \
+"============================================================================="
+#ifndef BENCH_IO_LMAX
+#define BENCH_IO_LMAX 40
+#endif
+
+using namespace Grid;
+using namespace QCD;
+
+int main (int argc, char ** argv)
+{
+  std::vector<std::string> dir;
+  unsigned int             Ls;
+  bool                     rb;
+  if (argc < 4)
+  {
+    std::cerr << "usage: " << argv[0] << " <Ls> <RB {0|1}> <dir1> [<dir2> ... <dirn>] [Grid options]";
+    std::cerr << std::endl;
+  }
+  Ls = std::stoi(argv[1]);
+  rb = (std::string(argv[2]) == "1");
+  for (unsigned int i = 3; i < argc; ++i)
+  {
+    std::string a = argv[i];
+
+    if (a[0] != '-')
+    {
+      dir.push_back(std::string(argv[i]));
+    }
+    else
+    {
+      break;
+    }
+  }
+  Grid_init(&argc,&argv);
+
+
+  int64_t threads = GridThread::GetThreads();
+  MSG << "Grid is setup to use " << threads << " threads" << std::endl;
+  MSG << SEP << std::endl;
+  MSG << "Benchmark Lime write" << std::endl;
+  MSG << SEP << std::endl;
+  for (auto &d: dir)
+  {
+    MSG << "-- Directory " << d << std::endl;
+    writeBenchmark<LatticeFermion>(GridDefaultLatt(), d + "/ioBench", limeWrite<LatticeFermion>, Ls, rb);
+  }
+
+  MSG << "Benchmark Lime read" << std::endl;
+  MSG << SEP << std::endl;
+  for (auto &d: dir)
+  {
+    MSG << "-- Directory " << d << std::endl;
+    readBenchmark<LatticeFermion>(GridDefaultLatt(), d + "/ioBench", limeRead<LatticeFermion>, Ls, rb);
+  }
+
+  Grid_finalize();
+
+  return EXIT_SUCCESS;
+}

tests/debug/Test_cayley_cg.cc

@@ -1,5 +1,4 @@
     /*************************************************************************************
-
     Grid physics library, www.github.com/paboyle/Grid 
 
     Source file: ./tests/Test_cayley_cg.cc
@@ -27,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     *************************************************************************************/
     /*  END LEGAL */
 #include <Grid/Grid.h>
-#include <Grid/algorithms/iterative/Reconstruct5Dprop.h>
+#include <Grid/qcd/action/fermion/Reconstruct5Dprop.h>
 
 using namespace std;
 using namespace Grid;
@@ -47,6 +46,7 @@ struct scal {
 
 template<class What> 
 void  TestCGinversions(What & Ddwf, 
+		       LatticeGaugeField &Umu,
 		       GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
 		       GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
 		       RealD mass, RealD M5,
@@ -78,12 +78,23 @@ void  TestCGprec(What & Ddwf,
 
 template<class What> 
 void  TestReconstruct5D(What & Ddwf, 
+			LatticeGaugeField &Umu,
 			GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
 			GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
 			RealD mass, RealD M5,
 			GridParallelRNG *RNG4,
 			GridParallelRNG *RNG5);
 
+template<class What,class WhatF> 
+void  TestReconstruct5DFA(What & Ddwf, 
+			  WhatF & DdwfF, 
+			  LatticeGaugeField &Umu,
+			  GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+			  GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+			  RealD mass, RealD M5,
+			  GridParallelRNG *RNG4,
+			  GridParallelRNG *RNG5);
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
@@ -92,63 +103,104 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
 
   const int Ls=8;
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), 
+								   GridDefaultSimd(Nd,vComplex::Nsimd()),
+								   GridDefaultMpi());
   GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
   GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
 
+  GridCartesian         * UGridF   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), 
+								    GridDefaultSimd(Nd,vComplexF::Nsimd()),
+								    GridDefaultMpi());
+  GridRedBlackCartesian * UrbGridF = SpaceTimeGrid::makeFourDimRedBlackGrid(UGridF);
+  GridCartesian         * FGridF   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGridF);
+  GridRedBlackCartesian * FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGridF);
+
+
   std::vector<int> seeds4({1,2,3,4});
   std::vector<int> seeds5({5,6,7,8});
   GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
   GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
 
   LatticeGaugeField Umu(UGrid);
+  LatticeGaugeFieldF UmuF(UGridF);
   SU3::HotConfiguration(RNG4,Umu);
+  precisionChange(UmuF,Umu);
   std::vector<LatticeColourMatrix> U(4,UGrid);
 
   RealD mass=0.1;
   RealD M5  =1.8;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"DomainWallFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  TestCGinversions<DomainWallFermionR>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  DomainWallFermionF DdwfF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5);
+  TestCGinversions<DomainWallFermionR>(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5DFA<DomainWallFermionR,DomainWallFermionF>(Ddwf,DdwfF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   RealD b=1.5;// Scale factor b+c=2, b-c=1
   RealD c=0.5;
   std::vector<ComplexD> gamma(Ls,ComplexD(1.0,0.0));
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"MobiusFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   MobiusFermionR Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
-  TestCGinversions<MobiusFermionR>(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  MobiusFermionF DmobF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5,b,c);
+  TestCGinversions<MobiusFermionR>(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5DFA<MobiusFermionR,MobiusFermionF>(Dmob,DmobF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"ZMobiusFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   ZMobiusFermionR ZDmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,gamma,b,c);
-  TestCGinversions<ZMobiusFermionR>(ZDmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestCGinversions<ZMobiusFermionR>(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5D<ZMobiusFermionR>(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"MobiusZolotarevFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   MobiusZolotarevFermionR Dzolo(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,0.1,2.0);
-  TestCGinversions<MobiusZolotarevFermionR>(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestCGinversions<MobiusZolotarevFermionR>(Dzolo,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5D<MobiusZolotarevFermionR>(Dzolo,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"ScaledShamirFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   ScaledShamirFermionR Dsham(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,2.0);
-  TestCGinversions<ScaledShamirFermionR>(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  ScaledShamirFermionF DshamF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5,2.0);
+  TestCGinversions<ScaledShamirFermionR>(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5DFA<ScaledShamirFermionR,ScaledShamirFermionF>(Dsham,DshamF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"ShamirZolotarevFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   ShamirZolotarevFermionR Dshamz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0);
-  TestCGinversions<ShamirZolotarevFermionR>(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestCGinversions<ShamirZolotarevFermionR>(Dshamz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5D<ShamirZolotarevFermionR>(Dshamz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"OverlapWilsonCayleyTanhFermion test"<<std::endl;
-  OverlapWilsonCayleyTanhFermionR Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
-  TestCGinversions<OverlapWilsonCayleyTanhFermionR>(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
+  OverlapWilsonCayleyTanhFermionR Dov (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
+  OverlapWilsonCayleyTanhFermionF DovF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5,1.0);
+  TestCGinversions<OverlapWilsonCayleyTanhFermionR>(Dov,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5DFA<OverlapWilsonCayleyTanhFermionR,OverlapWilsonCayleyTanhFermionF>(Dov,DovF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"OverlapWilsonCayleyZolotarevFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   OverlapWilsonCayleyZolotarevFermionR Dovz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0);
-  TestCGinversions<OverlapWilsonCayleyZolotarevFermionR>(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestCGinversions<OverlapWilsonCayleyZolotarevFermionR>(Dovz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5D<OverlapWilsonCayleyZolotarevFermionR>(Dovz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   Grid_finalize();
 }
 template<class What> 
 void  TestCGinversions(What & Ddwf, 
+		       LatticeGaugeField &Umu,
 		       GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
 		       GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
 		       RealD mass, RealD M5,
@@ -161,11 +213,9 @@ void  TestCGinversions(What & Ddwf,
   TestCGprec<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
   std::cout<<GridLogMessage << "Testing red black Schur inverter"<<std::endl;
   TestCGschur<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
-
-  std::cout<<GridLogMessage << "Testing 5D PV reconstruction"<<std::endl;
-  TestReconstruct5D<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
 }
 
+
 template<class What> 
 void  TestCGunprec(What & Ddwf, 
 		   GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
@@ -203,6 +253,7 @@ void  TestCGprec(What & Ddwf,
 
 template<class What> 
 void  TestReconstruct5D(What & Ddwf, 
+			LatticeGaugeField & Umu,
 			GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
 			GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
 			RealD mass, RealD M5,
@@ -216,9 +267,13 @@ void  TestReconstruct5D(What & Ddwf,
   LatticeFermion src_NE(FGrid);
   LatticeFermion result(FGrid);
   LatticeFermion result_rec(FGrid);
+  LatticeFermion result_madwf(FGrid);
 
   MdagMLinearOperator<What,LatticeFermion> HermOp(Ddwf);
-  ConjugateGradient<LatticeFermion> CG(1.0e-12,10000);
+  double Resid = 1.0e-12;
+  double Residi = 1.0e-6;
+  ConjugateGradient<LatticeFermion> CG(Resid,10000);
+  ConjugateGradient<LatticeFermion> CGi(Residi,10000);
 
   Ddwf.ImportPhysicalFermionSource(src4,src);
   Ddwf.Mdag(src,src_NE);
@@ -232,9 +287,16 @@ void  TestReconstruct5D(What & Ddwf,
 
   std::cout <<GridLogMessage<< " Reconstructing " <<std::endl;
 
-  //  typedef PauliVillarsSolverUnprec<LatticeFermion> PVinverter;
-  typedef PauliVillarsSolverRBprec<LatticeFermion> PVinverter;
-  PVinverter PVinverse(CG);
+  ////////////////////////////
+  // RBprec PV inverse
+  ////////////////////////////
+  typedef LatticeFermion Field;
+  typedef SchurRedBlackDiagTwoSolve<Field> SchurSolverType; 
+  typedef SchurRedBlackDiagTwoSolve<Field> SchurSolverTypei; 
+  typedef PauliVillarsSolverRBprec<Field,SchurSolverType> PVinverter;
+  SchurSolverType SchurSolver(CG);
+  PVinverter      PVinverse(SchurSolver);
+
   Reconstruct5DfromPhysical<LatticeFermion,PVinverter> reconstructor(PVinverse);
 
   reconstructor(Ddwf,res4,src4,result_rec);
@@ -245,6 +307,88 @@ void  TestReconstruct5D(What & Ddwf,
   result_rec = result_rec - result;
   std::cout <<GridLogMessage << "Difference "<<norm2(result_rec)<<std::endl;
 
+  //////////////////////////////
+  // Now try MADWF
+  //////////////////////////////
+  SchurSolverTypei SchurSolveri(CGi);
+  ZeroGuesser<LatticeFermion> Guess;
+  MADWF<What,What,PVinverter,SchurSolverTypei,ZeroGuesser<LatticeFermion> > 
+    madwf(Ddwf,Ddwf,PVinverse,SchurSolveri,Guess,Resid,10);
+  
+  madwf(src4,result_madwf);
+  result_madwf = result_madwf - result;
+  std::cout <<GridLogMessage << "Difference "<<norm2(result_madwf)<<std::endl;
+
+
+}
+template<class What,class WhatF> 
+void  TestReconstruct5DFA(What & Ddwf, 
+			  WhatF & DdwfF, 
+			  LatticeGaugeField & Umu,
+			  GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+			  GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+			  RealD mass, RealD M5,
+			  GridParallelRNG *RNG4,
+			  GridParallelRNG *RNG5)
+{
+  LatticeFermion src4   (UGrid); random(*RNG4,src4);
+  LatticeFermion res4   (UGrid); res4 = zero;
+
+  LatticeFermion src   (FGrid);
+  LatticeFermion src_NE(FGrid);
+  LatticeFermion result(FGrid);
+  LatticeFermion result_rec(FGrid);
+  LatticeFermion result_madwf(FGrid);
+
+  MdagMLinearOperator<What,LatticeFermion> HermOp(Ddwf);
+  double Resid = 1.0e-12;
+  double Residi = 1.0e-5;
+  ConjugateGradient<LatticeFermion> CG(Resid,10000);
+  ConjugateGradient<LatticeFermionF> CGi(Residi,10000);
+
+  Ddwf.ImportPhysicalFermionSource(src4,src);
+  Ddwf.Mdag(src,src_NE);
+  CG(HermOp,src_NE,result);
+
+  Ddwf.ExportPhysicalFermionSolution(result, res4);
+
+  Ddwf.M(result,src_NE);
+  src_NE = src_NE - src;
+  std::cout <<GridLogMessage<< " True residual is " << norm2(src_NE)<<std::endl;
+
+  std::cout <<GridLogMessage<< " Reconstructing " <<std::endl;
+
+  ////////////////////////////
+  // Fourier accel PV inverse
+  ////////////////////////////
+  typedef LatticeFermion Field;
+  typedef LatticeFermionF FieldF;
+  typedef SchurRedBlackDiagTwoSolve<FieldF> SchurSolverTypei; 
+  typedef PauliVillarsSolverFourierAccel<LatticeFermion,LatticeGaugeField> PVinverter;
+  PVinverter PVinverse(Umu,CG);
+
+  Reconstruct5DfromPhysical<LatticeFermion,PVinverter> reconstructor(PVinverse);
+
+  reconstructor(Ddwf,res4,src4,result_rec);
+
+  std::cout <<GridLogMessage << "Result     "<<norm2(result)<<std::endl;
+  std::cout <<GridLogMessage << "Result_rec "<<norm2(result_rec)<<std::endl;
+
+  result_rec = result_rec - result;
+  std::cout <<GridLogMessage << "Difference "<<norm2(result_rec)<<std::endl;
+
+  //////////////////////////////
+  // Now try MADWF
+  //////////////////////////////
+  SchurSolverTypei SchurSolver(CGi);
+  ZeroGuesser<LatticeFermionF> Guess;
+  MADWF<What,WhatF,PVinverter,SchurSolverTypei,ZeroGuesser<LatticeFermionF> > 
+    madwf(Ddwf,DdwfF,PVinverse,SchurSolver,Guess,Resid,10);
+  
+  madwf(src4,result_madwf);
+  result_madwf = result_madwf - result;
+  std::cout <<GridLogMessage << "Difference "<<norm2(result_madwf)<<std::endl;
+
 }
 
 

tests/hadrons/Test_diskvector.cc

@@ -91,6 +91,22 @@ int main(int argc, char *argv[])
     v13r = v[13];
     LOG(Message) << "v[13] correct? " 
                  << ((v13r == v13w) ? "yes" : "no" ) << std::endl;
+    LOG(Message) << "hit ratio " << v.hitRatio() << std::endl;
+
+    EigenDiskVector<ComplexD>         w("eigendiskvector_test", 1000, 4);
+    EigenDiskVector<ComplexD>::Matrix m,n;
+
+    w[2] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
+    m    = w[2];
+    w[3] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
+    w[4] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
+    w[5] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
+    w[6] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
+    w[7] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
+    n    = w[2];
+    LOG(Message) << "w[2] correct? " 
+                 << ((m == n) ? "yes" : "no" ) << std::endl;
+    LOG(Message) << "hit ratio " << w.hitRatio() << std::endl;
 
     Grid_finalize();