Merge branch 'develop' of https://github.com/paboyle/Grid into develop

benchmarks/Benchmark_su3.cc

@@ -35,7 +35,8 @@ using namespace Grid::QCD;
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
-#define LMAX (64)
+#define LMAX (40)
+#define LINC (4)
 
   int64_t Nloop=20;
 
@@ -51,7 +52,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
   std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
 
-  for(int lat=2;lat<=LMAX;lat+=2){
+  for(int lat=2;lat<=LMAX;lat+=LINC){
 
       std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
       int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
@@ -83,7 +84,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
   std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
 
-  for(int lat=2;lat<=LMAX;lat+=2){
+  for(int lat=2;lat<=LMAX;lat+=LINC){
 
       std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
       int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
@@ -114,7 +115,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
   std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
 
-  for(int lat=2;lat<=LMAX;lat+=2){
+  for(int lat=2;lat<=LMAX;lat+=LINC){
 
       std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
       int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
@@ -145,7 +146,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
   std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
 
-  for(int lat=2;lat<=LMAX;lat+=2){
+  for(int lat=2;lat<=LMAX;lat+=LINC){
 
       std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
       int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
@@ -165,10 +166,43 @@ int main (int argc, char ** argv)
       double time = (stop-start)/Nloop*1000.0;
       
       double bytes=3*vol*Nc*Nc*sizeof(Complex);
-      double flops=Nc*Nc*(8+8+8)*vol;
+      double flops=Nc*Nc*(6+8+8)*vol;
       std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
 
     }
 
+
+  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
+  std::cout<<GridLogMessage << "= Benchmarking SU3xSU3  CovShiftForward(z,x,y)"<<std::endl;
+  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
+  std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
+  std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
+
+  for(int lat=2;lat<=LMAX;lat+=LINC){
+
+      std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
+      int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
+
+      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
+      GridParallelRNG          pRNG(&Grid);      pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+
+      LatticeColourMatrix z(&Grid); random(pRNG,z);
+      LatticeColourMatrix x(&Grid); random(pRNG,x);
+      LatticeColourMatrix y(&Grid); random(pRNG,y);
+
+      for(int mu=0;mu<=4;mu++){
+	double start=usecond();
+	for(int64_t i=0;i<Nloop;i++){
+	  z = PeriodicBC::CovShiftForward(x,mu,y);
+	}
+	double stop=usecond();
+	double time = (stop-start)/Nloop*1000.0;
+	
+	double bytes=3*vol*Nc*Nc*sizeof(Complex);
+	double flops=Nc*Nc*(6+8+8)*vol;
+	std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
+      }
+    }
+
   Grid_finalize();
 }

lib/lattice/Lattice_transfer.h

@@ -599,6 +599,51 @@ unvectorizeToLexOrdArray(std::vector<sobj> &out, const Lattice<vobj> &in)
     extract1(in_vobj, out_ptrs, 0);
   }
 }
+
+template<typename vobj, typename sobj>
+typename std::enable_if<isSIMDvectorized<vobj>::value && !isSIMDvectorized<sobj>::value, void>::type 
+unvectorizeToRevLexOrdArray(std::vector<sobj> &out, const Lattice<vobj> &in)
+{
+
+  typedef typename vobj::vector_type vtype;
+  
+  GridBase* in_grid = in._grid;
+  out.resize(in_grid->lSites());
+  
+  int ndim = in_grid->Nd();
+  int in_nsimd = vtype::Nsimd();
+
+  std::vector<std::vector<int> > in_icoor(in_nsimd);
+      
+  for(int lane=0; lane < in_nsimd; lane++){
+    in_icoor[lane].resize(ndim);
+    in_grid->iCoorFromIindex(in_icoor[lane], lane);
+  }
+  
+  parallel_for(int in_oidx = 0; in_oidx < in_grid->oSites(); in_oidx++){ //loop over outer index
+    //Assemble vector of pointers to output elements
+    std::vector<sobj*> out_ptrs(in_nsimd);
+
+    std::vector<int> in_ocoor(ndim);
+    in_grid->oCoorFromOindex(in_ocoor, in_oidx);
+
+    std::vector<int> lcoor(in_grid->Nd());
+      
+    for(int lane=0; lane < in_nsimd; lane++){
+      for(int mu=0;mu<ndim;mu++)
+	lcoor[mu] = in_ocoor[mu] + in_grid->_rdimensions[mu]*in_icoor[lane][mu];
+
+      int lex;
+      Lexicographic::IndexFromCoorReversed(lcoor, lex, in_grid->_ldimensions);
+      out_ptrs[lane] = &out[lex];
+    }
+    
+    //Unpack into those ptrs
+    const vobj & in_vobj = in._odata[in_oidx];
+    extract1(in_vobj, out_ptrs, 0);
+  }
+}
+
 //Copy SIMD-vectorized lattice to array of scalar objects in lexicographic order
 template<typename vobj, typename sobj>
 typename std::enable_if<isSIMDvectorized<vobj>::value 
@@ -648,6 +693,54 @@ vectorizeFromLexOrdArray( std::vector<sobj> &in, Lattice<vobj> &out)
   }
 }
 
+template<typename vobj, typename sobj>
+typename std::enable_if<isSIMDvectorized<vobj>::value 
+                    && !isSIMDvectorized<sobj>::value, void>::type 
+vectorizeFromRevLexOrdArray( std::vector<sobj> &in, Lattice<vobj> &out)
+{
+
+  typedef typename vobj::vector_type vtype;
+  
+  GridBase* grid = out._grid;
+  assert(in.size()==grid->lSites());
+  
+  int ndim     = grid->Nd();
+  int nsimd    = vtype::Nsimd();
+
+  std::vector<std::vector<int> > icoor(nsimd);
+      
+  for(int lane=0; lane < nsimd; lane++){
+    icoor[lane].resize(ndim);
+    grid->iCoorFromIindex(icoor[lane],lane);
+  }
+  
+  parallel_for(uint64_t oidx = 0; oidx < grid->oSites(); oidx++){ //loop over outer index
+    //Assemble vector of pointers to output elements
+    std::vector<sobj*> ptrs(nsimd);
+
+    std::vector<int> ocoor(ndim);
+    grid->oCoorFromOindex(ocoor, oidx);
+
+    std::vector<int> lcoor(grid->Nd());
+      
+    for(int lane=0; lane < nsimd; lane++){
+
+      for(int mu=0;mu<ndim;mu++){
+	lcoor[mu] = ocoor[mu] + grid->_rdimensions[mu]*icoor[lane][mu];
+      }
+
+      int lex;
+      Lexicographic::IndexFromCoorReversed(lcoor, lex, grid->_ldimensions);
+      ptrs[lane] = &in[lex];
+    }
+    
+    //pack from those ptrs
+    vobj vecobj;
+    merge1(vecobj, ptrs, 0);
+    out._odata[oidx] = vecobj; 
+  }
+}
+
 //Convert a Lattice from one precision to another
 template<class VobjOut, class VobjIn>
 void precisionChange(Lattice<VobjOut> &out, const Lattice<VobjIn> &in){

lib/qcd/action/fermion/CayleyFermion5D.cc

@@ -52,6 +52,35 @@ namespace QCD {
  { 
  }
 
+///////////////////////////////////////////////////////////////
+// Physical surface field utilities
+///////////////////////////////////////////////////////////////
+template<class Impl>  
+void CayleyFermion5D<Impl>::ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d)
+{
+  int Ls = this->Ls;
+  FermionField tmp(this->FermionGrid());
+  tmp = solution5d;
+  conformable(solution5d._grid,this->FermionGrid());
+  conformable(exported4d._grid,this->GaugeGrid());
+  axpby_ssp_pminus(tmp, 0., solution5d, 1., solution5d, 0, 0);
+  axpby_ssp_pplus (tmp, 1., tmp       , 1., solution5d, 0, Ls-1);
+  ExtractSlice(exported4d, tmp, 0, 0);
+}
+template<class Impl>  
+void CayleyFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d)
+{
+  int Ls = this->Ls;
+  FermionField tmp(this->FermionGrid());
+  conformable(imported5d._grid,this->FermionGrid());
+  conformable(input4d._grid   ,this->GaugeGrid());
+  tmp = zero;
+  InsertSlice(input4d, tmp, 0   , 0);
+  InsertSlice(input4d, tmp, Ls-1, 0);
+  axpby_ssp_pplus (tmp, 0., tmp, 1., tmp, 0, 0);
+  axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1);
+  Dminus(tmp,imported5d);
+}
 template<class Impl>  
 void CayleyFermion5D<Impl>::Dminus(const FermionField &psi, FermionField &chi)
 {

lib/qcd/action/fermion/CayleyFermion5D.h

@@ -83,8 +83,13 @@ namespace Grid {
       virtual void   M5D   (const FermionField &psi, FermionField &chi);
       virtual void   M5Ddag(const FermionField &psi, FermionField &chi);
 
+      ///////////////////////////////////////////////////////////////
+      // Physical surface field utilities
+      ///////////////////////////////////////////////////////////////
       virtual void   Dminus(const FermionField &psi, FermionField &chi);
       virtual void   DminusDag(const FermionField &psi, FermionField &chi);
+      virtual void ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d);
+      virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d);
 
       /////////////////////////////////////////////////////
       // Instantiate different versions depending on Impl

lib/qcd/action/fermion/ContinuedFractionFermion5D.cc

@@ -295,6 +295,27 @@ namespace Grid {
       assert((Ls&0x1)==1); // Odd Ls required
     }
 
+    template<class Impl>
+    void ContinuedFractionFermion5D<Impl>::ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d)
+    {
+      int Ls = this->Ls;
+      conformable(solution5d._grid,this->FermionGrid());
+      conformable(exported4d._grid,this->GaugeGrid());
+      ExtractSlice(exported4d, solution5d, Ls-1, Ls-1);
+    }
+    template<class Impl>
+    void ContinuedFractionFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d)
+    {
+      int Ls = this->Ls;
+      conformable(imported5d._grid,this->FermionGrid());
+      conformable(input4d._grid   ,this->GaugeGrid());
+      FermionField tmp(this->FermionGrid());
+      tmp=zero;
+      InsertSlice(input4d, tmp, Ls-1, Ls-1);
+      tmp=Gamma(Gamma::Algebra::Gamma5)*tmp;
+      this->Dminus(tmp,imported5d);
+    }
+
     FermOpTemplateInstantiate(ContinuedFractionFermion5D);
 
   }

lib/qcd/action/fermion/ContinuedFractionFermion5D.h

@@ -65,6 +65,14 @@ namespace Grid {
       // Efficient support for multigrid coarsening
       virtual void  Mdir (const FermionField &in, FermionField &out,int dir,int disp);
 
+      ///////////////////////////////////////////////////////////////
+      // Physical surface field utilities
+      ///////////////////////////////////////////////////////////////
+      //      virtual void Dminus(const FermionField &psi, FermionField &chi);     // Inherit trivial case
+      //      virtual void DminusDag(const FermionField &psi, FermionField &chi);  // Inherit trivial case
+      virtual void ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d);
+      virtual void ImportPhysicalFermionSource  (const FermionField &input4d,FermionField &imported5d);
+
       // Constructors
       ContinuedFractionFermion5D(GaugeField &_Umu,
 				 GridCartesian         &FiveDimGrid,

lib/qcd/action/fermion/FermionOperator.h

@@ -128,6 +128,19 @@ namespace Grid {
                                        std::vector<Real> mom,
                                        unsigned int tmin, 
                                        unsigned int tmax)=0;
+      ///////////////////////////////////////////////
+      // Physical field import/export
+      ///////////////////////////////////////////////
+      virtual void Dminus(const FermionField &psi, FermionField &chi)    { chi=psi; }
+      virtual void DminusDag(const FermionField &psi, FermionField &chi) { chi=psi; }
+      virtual void ImportPhysicalFermionSource(const FermionField &input,FermionField &imported)
+      {
+	imported = input;
+      };
+      virtual void ExportPhysicalFermionSolution(const FermionField &solution,FermionField &exported)
+      {
+	exported=solution;
+      };
     };
 
   }

lib/qcd/action/fermion/PartialFractionFermion5D.cc

@@ -396,6 +396,27 @@ namespace Grid {
       amax=zolo_hi;
     }
 
+    template<class Impl>
+    void PartialFractionFermion5D<Impl>::ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d)
+    {
+      int Ls = this->Ls;
+      conformable(solution5d._grid,this->FermionGrid());
+      conformable(exported4d._grid,this->GaugeGrid());
+      ExtractSlice(exported4d, solution5d, Ls-1, Ls-1);
+    }
+    template<class Impl>
+    void PartialFractionFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d)
+    {
+      int Ls = this->Ls;
+      conformable(imported5d._grid,this->FermionGrid());
+      conformable(input4d._grid   ,this->GaugeGrid());
+      FermionField tmp(this->FermionGrid());
+      tmp=zero;
+      InsertSlice(input4d, tmp, Ls-1, Ls-1);
+      tmp=Gamma(Gamma::Algebra::Gamma5)*tmp;
+      this->Dminus(tmp,imported5d);
+    }
+
       // Constructors
     template<class Impl>
     PartialFractionFermion5D<Impl>::PartialFractionFermion5D(GaugeField &_Umu,

lib/qcd/action/fermion/PartialFractionFermion5D.h

@@ -70,6 +70,12 @@ namespace Grid {
       // Efficient support for multigrid coarsening
       virtual void  Mdir (const FermionField &in, FermionField &out,int dir,int disp);
 
+      ///////////////////////////////////////////////////////////////
+      // Physical surface field utilities
+      ///////////////////////////////////////////////////////////////
+      virtual void ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d);
+      virtual void ImportPhysicalFermionSource  (const FermionField &input4d,FermionField &imported5d);
+
       // Constructors
       PartialFractionFermion5D(GaugeField &_Umu,
 			       GridCartesian         &FiveDimGrid,

lib/qcd/smearing/WilsonFlow.h

@@ -173,8 +173,8 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
         std::cout << "Time to evolve " << diff.count() << " s\n";
         #endif
         std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
-            << step << "  "
+		  << step << "  " << tau(step) << "  " 
-            << energyDensityPlaquette(step,out) << std::endl;
+		  << energyDensityPlaquette(step,out) << std::endl;
          if( step % measure_interval == 0){
          std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "
             << step << "  " 
@@ -193,8 +193,8 @@ void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, Re
         //std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl;
         evolve_step_adaptive(out, maxTau);
         std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
-            << step << "  "
+		  << step << "  " << taus << "  "
-            << energyDensityPlaquette(out) << std::endl;
+		  << energyDensityPlaquette(out) << std::endl;
          if( step % measure_interval == 0){
          std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "
             << step << "  " 

lib/serialisation/Hdf5IO.cc

@@ -55,6 +55,11 @@ void Hdf5Writer::writeDefault(const std::string &s, const char *x)
   writeDefault(s, sx);
 }
 
+Group & Hdf5Writer::getGroup(void)
+{
+  return group_;
+}
+
 // Reader implementation ///////////////////////////////////////////////////////
 Hdf5Reader::Hdf5Reader(const std::string &fileName)
 : fileName_(fileName)
@@ -103,3 +108,8 @@ void Hdf5Reader::readDefault(const std::string &s, std::string &x)
   x.resize(strType.getSize());
   attribute.read(strType, &(x[0]));
 }
+
+Group & Hdf5Reader::getGroup(void)
+{
+  return group_;
+}

lib/serialisation/Hdf5IO.h

@@ -38,6 +38,7 @@ namespace Grid
     template <typename U>
     typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
     writeDefault(const std::string &s, const std::vector<U> &x);
+    H5NS::Group & getGroup(void);
   private:
     template <typename U>
     void writeSingleAttribute(const U &x, const std::string &name,
@@ -65,6 +66,7 @@ namespace Grid
     template <typename U>
     typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
     readDefault(const std::string &s, std::vector<U> &x);
+    H5NS::Group & getGroup(void);
   private:
     template <typename U>
     void readSingleAttribute(U &x, const std::string &name,

lib/serialisation/VectorUtils.h

@@ -30,6 +30,48 @@ namespace Grid {
     typedef typename std::vector<std::vector<typename TensorToVec<T>::type>> type;
   };
 
+  template <typename T>
+  void tensorDim(std::vector<size_t> &dim, const T &t, const bool wipe = true)
+  {
+    if (wipe)
+    {
+      dim.clear();
+    }
+  }
+
+  template <typename T>
+  void tensorDim(std::vector<size_t> &dim, const iScalar<T> &t, const bool wipe = true)
+  {
+    if (wipe)
+    {
+      dim.clear();
+    }
+    tensorDim(dim, t._internal, false);
+  }
+
+  template <typename T, int N>
+  void tensorDim(std::vector<size_t> &dim, const iVector<T, N> &t, const bool wipe = true)
+  {
+    if (wipe)
+    {
+      dim.clear();
+    }
+    dim.push_back(N);
+    tensorDim(dim, t._internal[0], false);
+  }
+
+  template <typename T, int N>
+  void tensorDim(std::vector<size_t> &dim, const iMatrix<T, N> &t, const bool wipe = true)
+  {
+    if (wipe)
+    {
+      dim.clear();
+    }
+    dim.push_back(N);
+    dim.push_back(N);
+    tensorDim(dim, t._internal[0][0], false);
+  }
+
   template <typename T>
   typename TensorToVec<T>::type tensorToVec(const T &t)
   {