Deprecate UVM

Grid/qcd/action/fermion/CayleyFermion5D.h

@@ -90,16 +90,16 @@ public:
   void M5D(const FermionField &psi,
 	   const FermionField &phi,
 	   FermionField &chi,
-	   Vector<Coeff_t> &lower,
-	   Vector<Coeff_t> &diag,
-	   Vector<Coeff_t> &upper);
+	   std::vector<Coeff_t> &lower,
+	   std::vector<Coeff_t> &diag,
+	   std::vector<Coeff_t> &upper);
 
   void M5Ddag(const FermionField &psi,
 	      const FermionField &phi,
 	      FermionField &chi,
-	      Vector<Coeff_t> &lower,
-	      Vector<Coeff_t> &diag,
-	      Vector<Coeff_t> &upper);
+	      std::vector<Coeff_t> &lower,
+	      std::vector<Coeff_t> &diag,
+	      std::vector<Coeff_t> &upper);
 
   virtual void   Instantiatable(void)=0;
 
@@ -119,35 +119,35 @@ public:
   RealD mass_plus, mass_minus;
 
   // Save arguments to SetCoefficientsInternal
-  Vector<Coeff_t> _gamma;
+  std::vector<Coeff_t> _gamma;
   RealD                _zolo_hi;
   RealD                _b;
   RealD                _c;
 
   // Cayley form Moebius (tanh and zolotarev)
-  Vector<Coeff_t> omega;
-  Vector<Coeff_t> bs;    // S dependent coeffs
-  Vector<Coeff_t> cs;
-  Vector<Coeff_t> as;
+  std::vector<Coeff_t> omega;
+  std::vector<Coeff_t> bs;    // S dependent coeffs
+  std::vector<Coeff_t> cs;
+  std::vector<Coeff_t> as;
   // For preconditioning Cayley form
-  Vector<Coeff_t> bee;
-  Vector<Coeff_t> cee;
-  Vector<Coeff_t> aee;
-  Vector<Coeff_t> beo;
-  Vector<Coeff_t> ceo;
-  Vector<Coeff_t> aeo;
+  std::vector<Coeff_t> bee;
+  std::vector<Coeff_t> cee;
+  std::vector<Coeff_t> aee;
+  std::vector<Coeff_t> beo;
+  std::vector<Coeff_t> ceo;
+  std::vector<Coeff_t> aeo;
   // LDU factorisation of the eeoo matrix
-  Vector<Coeff_t> lee;
-  Vector<Coeff_t> leem;
-  Vector<Coeff_t> uee;
-  Vector<Coeff_t> ueem;
-  Vector<Coeff_t> dee;
+  std::vector<Coeff_t> lee;
+  std::vector<Coeff_t> leem;
+  std::vector<Coeff_t> uee;
+  std::vector<Coeff_t> ueem;
+  std::vector<Coeff_t> dee;
 
   // Matrices of 5d ee inverse params
-  Vector<iSinglet<Simd> >  MatpInv;
-  Vector<iSinglet<Simd> >  MatmInv;
-  Vector<iSinglet<Simd> >  MatpInvDag;
-  Vector<iSinglet<Simd> >  MatmInvDag;
+  //  std::vector<iSinglet<Simd> >  MatpInv;
+  //  std::vector<iSinglet<Simd> >  MatmInv;
+  //  std::vector<iSinglet<Simd> >  MatpInvDag;
+  //  std::vector<iSinglet<Simd> >  MatmInvDag;
 
   ///////////////////////////////////////////////////////////////
   // Conserved current utilities
@@ -187,7 +187,7 @@ public:
 protected:
   virtual void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c);
   virtual void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c);
-  virtual void SetCoefficientsInternal(RealD zolo_hi,Vector<Coeff_t> & gamma,RealD b,RealD c);
+  virtual void SetCoefficientsInternal(RealD zolo_hi,std::vector<Coeff_t> & gamma,RealD b,RealD c);
 };
 
 NAMESPACE_END(Grid);

Grid/qcd/action/fermion/ContinuedFractionFermion5D.h

@@ -90,12 +90,12 @@ protected:
   RealD mass;
   RealD R;
   RealD ZoloHiInv;
-  Vector<double> Beta;
-  Vector<double> cc;;
-  Vector<double> cc_d;;
-  Vector<double> sqrt_cc;
-  Vector<double> See;
-  Vector<double> Aee;
+  std::vector<double> Beta;
+  std::vector<double> cc;;
+  std::vector<double> cc_d;;
+  std::vector<double> sqrt_cc;
+  std::vector<double> See;
+  std::vector<double> Aee;
 
 };
 

Grid/qcd/action/fermion/DomainWallEOFAFermion.h

@@ -69,10 +69,10 @@ public:
   // Instantiate different versions depending on Impl
   /////////////////////////////////////////////////////
   void M5D(const FermionField& psi, const FermionField& phi, FermionField& chi,
-	   Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper);
+	   std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
 
   void M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi,
-	      Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper);
+	      std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
 
   virtual void RefreshShiftCoefficients(RealD new_shift);
 
@@ -83,7 +83,7 @@ public:
 			RealD _M5, const ImplParams& p=ImplParams());
 
 protected:
-  void SetCoefficientsInternal(RealD zolo_hi, Vector<Coeff_t>& gamma, RealD b, RealD c);
+  void SetCoefficientsInternal(RealD zolo_hi, std::vector<Coeff_t>& gamma, RealD b, RealD c);
 };
 
 NAMESPACE_END(Grid);

Grid/qcd/action/fermion/ImprovedStaggeredFermion.h

@@ -102,11 +102,11 @@ public:
 		     GaugeField &mat, 
 		     const FermionField &A, const FermionField &B, int dag);
 
-  void DhopInternal(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,DoubledGaugeField &UUU,
+  void DhopInternal(StencilImpl &st, DoubledGaugeField &U,DoubledGaugeField &UUU,
                     const FermionField &in, FermionField &out, int dag);
-  void DhopInternalSerialComms(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,DoubledGaugeField &UUU,
+  void DhopInternalSerialComms(StencilImpl &st, DoubledGaugeField &U,DoubledGaugeField &UUU,
                     const FermionField &in, FermionField &out, int dag);
-  void DhopInternalOverlappedComms(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,DoubledGaugeField &UUU,
+  void DhopInternalOverlappedComms(StencilImpl &st, DoubledGaugeField &U,DoubledGaugeField &UUU,
                     const FermionField &in, FermionField &out, int dag);
 
   //////////////////////////////////////////////////////////////////////////
@@ -164,8 +164,6 @@ public:
   DoubledGaugeField UUUmuEven;
   DoubledGaugeField UUUmuOdd;
 
-  LebesgueOrder Lebesgue;
-  LebesgueOrder LebesgueEvenOdd;
   
   ///////////////////////////////////////////////////////////////
   // Conserved current utilities

Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h

@@ -100,7 +100,6 @@ public:
 		     int dag);
     
   void DhopInternal(StencilImpl & st,
-		    LebesgueOrder &lo,
 		    DoubledGaugeField &U,
 		    DoubledGaugeField &UUU,
 		    const FermionField &in, 
@@ -108,7 +107,6 @@ public:
 		    int dag);
     
     void DhopInternalOverlappedComms(StencilImpl & st,
-		      LebesgueOrder &lo,
 		      DoubledGaugeField &U,
 		      DoubledGaugeField &UUU,
 		      const FermionField &in, 
@@ -116,7 +114,6 @@ public:
 		      int dag);
 
     void DhopInternalSerialComms(StencilImpl & st,
-		      LebesgueOrder &lo,
 		      DoubledGaugeField &U,
 		      DoubledGaugeField &UUU,
 		      const FermionField &in, 
@@ -192,8 +189,6 @@ public:
   DoubledGaugeField UUUmuEven;
   DoubledGaugeField UUUmuOdd;
     
-  LebesgueOrder Lebesgue;
-  LebesgueOrder LebesgueEvenOdd;
     
   // Comms buffer
   //  std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  comm_buf;

Grid/qcd/action/fermion/MobiusEOFAFermion.h

@@ -42,11 +42,11 @@ public:
 
 public:
   // Shift operator coefficients for red-black preconditioned Mobius EOFA
-  Vector<Coeff_t> Mooee_shift;
-  Vector<Coeff_t> MooeeInv_shift_lc;
-  Vector<Coeff_t> MooeeInv_shift_norm;
-  Vector<Coeff_t> MooeeInvDag_shift_lc;
-  Vector<Coeff_t> MooeeInvDag_shift_norm;
+  std::vector<Coeff_t> Mooee_shift;
+  std::vector<Coeff_t> MooeeInv_shift_lc;
+  std::vector<Coeff_t> MooeeInv_shift_norm;
+  std::vector<Coeff_t> MooeeInvDag_shift_lc;
+  std::vector<Coeff_t> MooeeInvDag_shift_norm;
 
   virtual void Instantiatable(void) {};
 
@@ -74,18 +74,18 @@ public:
   // Instantiate different versions depending on Impl
   /////////////////////////////////////////////////////
   void M5D(const FermionField& psi, const FermionField& phi, FermionField& chi,
-	   Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper);
+	   std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
 
   void M5D_shift(const FermionField& psi, const FermionField& phi, FermionField& chi,
-		 Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper,
-		 Vector<Coeff_t>& shift_coeffs);
+		 std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper,
+		 std::vector<Coeff_t>& shift_coeffs);
 
   void M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi,
-	      Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper);
+	      std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
 
   void M5Ddag_shift(const FermionField& psi, const FermionField& phi, FermionField& chi,
-		    Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper,
-		    Vector<Coeff_t>& shift_coeffs);
+		    std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper,
+		    std::vector<Coeff_t>& shift_coeffs);
 
   virtual void RefreshShiftCoefficients(RealD new_shift);
 

Grid/qcd/action/fermion/NaiveStaggeredFermion.h

@@ -102,11 +102,11 @@ public:
 		     GaugeField &mat, 
 		     const FermionField &A, const FermionField &B, int dag);
 
-  void DhopInternal(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+  void DhopInternal(StencilImpl &st, DoubledGaugeField &U,
                     const FermionField &in, FermionField &out, int dag);
-  void DhopInternalSerialComms(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+  void DhopInternalSerialComms(StencilImpl &st, DoubledGaugeField &U,
 			       const FermionField &in, FermionField &out, int dag);
-  void DhopInternalOverlappedComms(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+  void DhopInternalOverlappedComms(StencilImpl &st, DoubledGaugeField &U,
 				   const FermionField &in, FermionField &out, int dag);
 
   //////////////////////////////////////////////////////////////////////////
@@ -152,9 +152,6 @@ public:
   DoubledGaugeField UmuEven;
   DoubledGaugeField UmuOdd;
 
-  LebesgueOrder Lebesgue;
-  LebesgueOrder LebesgueEvenOdd;
-  
   ///////////////////////////////////////////////////////////////
   // Conserved current utilities
   ///////////////////////////////////////////////////////////////

Grid/qcd/action/fermion/PartialFractionFermion5D.h

@@ -94,8 +94,8 @@ protected:
   RealD R;
   RealD amax;
   RealD scale;
-  Vector<double> p; 
-  Vector<double> q;
+  std::vector<double> p; 
+  std::vector<double> q;
 
 };
 

Grid/qcd/action/fermion/SchurDiagTwoKappa.h

@@ -35,7 +35,7 @@ template<class Matrix, class Field>
 class KappaSimilarityTransform {
 public:
   INHERIT_IMPL_TYPES(Matrix);
-  Vector<Coeff_t> kappa, kappaDag, kappaInv, kappaInvDag;
+  std::vector<Coeff_t> kappa, kappaDag, kappaInv, kappaInvDag;
 
   KappaSimilarityTransform (Matrix &zmob) {
     for (int i=0;i<(int)zmob.bs.size();i++) {

Grid/qcd/action/fermion/StaggeredKernels.h

@@ -49,10 +49,10 @@ template<class Impl> class StaggeredKernels : public FermionOperator<Impl> , pub
    
  public:
 
-  void DhopImproved(StencilImpl &st, LebesgueOrder &lo, 
+  void DhopImproved(StencilImpl &st,
 		    DoubledGaugeField &U, DoubledGaugeField &UUU, 
 		    const FermionField &in, FermionField &out, int dag, int interior,int exterior);
-  void DhopNaive(StencilImpl &st, LebesgueOrder &lo, 
+  void DhopNaive(StencilImpl &st,
 		 DoubledGaugeField &U,
 		 const FermionField &in, FermionField &out, int dag, int interior,int exterior);
   

Grid/qcd/action/fermion/WilsonCompressor.h

@@ -47,7 +47,7 @@ public:
   static int PartialCompressionFactor(GridBase *grid) { return 1;}
 #endif
   template<class vobj,class cobj,class compressor>
-  static void Gather_plane_simple (commVector<std::pair<int,int> >& table,
+  static void Gather_plane_simple (deviceVector<std::pair<int,int> >& table,
 				   const Lattice<vobj> &rhs,
 				   cobj *buffer,
 				   compressor &compress,
@@ -109,7 +109,7 @@ public:
   // Reorder the fifth dim to be s=Ls-1 , s=0, s=1,...,Ls-2.
   ////////////////////////////////////////////////////////////////////////////////////////////
   template<class vobj,class cobj,class compressor>
-  static void Gather_plane_exchange(commVector<std::pair<int,int> >& table,const Lattice<vobj> &rhs,
+  static void Gather_plane_exchange(deviceVector<std::pair<int,int> >& table,const Lattice<vobj> &rhs,
 				    std::vector<cobj *> pointers,int dimension,int plane,int cbmask,
 				    compressor &compress,int type,int partial)
   {
@@ -197,7 +197,7 @@ public:
 #endif
   
   template<class vobj,class cobj,class compressor>
-  static void Gather_plane_simple (commVector<std::pair<int,int> >& table,
+  static void Gather_plane_simple (deviceVector<std::pair<int,int> >& table,
 					 const Lattice<vobj> &rhs,
 					 cobj *buffer,
 					 compressor &compress,
@@ -208,7 +208,7 @@ public:
     else        FaceGatherSimple::Gather_plane_simple(table,rhs,buffer,compress,off,so,partial);
   }
   template<class vobj,class cobj,class compressor>
-  static void Gather_plane_exchange(commVector<std::pair<int,int> >& table,const Lattice<vobj> &rhs,
+  static void Gather_plane_exchange(deviceVector<std::pair<int,int> >& table,const Lattice<vobj> &rhs,
 				    std::vector<cobj *> pointers,int dimension,int plane,int cbmask,
 				    compressor &compress,int type,int partial)
   {
@@ -402,7 +402,6 @@ public:
 
   typedef CartesianStencil<vobj,cobj,Parameters> Base;
   typedef typename Base::View_type View_type;
-  typedef typename Base::StencilVector StencilVector;
 
   //  Vector<int> surface_list;
   WilsonStencil(GridBase *grid,

Grid/qcd/action/fermion/WilsonFermion.h

@@ -126,14 +126,17 @@ public:
   void DerivInternal(StencilImpl &st, DoubledGaugeField &U, GaugeField &mat,
                      const FermionField &A, const FermionField &B, int dag);
 
-  void DhopInternal(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+  void DhopInternal(StencilImpl &st,
+		    DoubledGaugeField &U,
                     const FermionField &in, FermionField &out, int dag);
 
-  void DhopInternalSerial(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
-                    const FermionField &in, FermionField &out, int dag);
+  void DhopInternalSerial(StencilImpl &st,
+			  DoubledGaugeField &U,
+			  const FermionField &in, FermionField &out, int dag);
 
-  void DhopInternalOverlappedComms(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
-                    const FermionField &in, FermionField &out, int dag);
+  void DhopInternalOverlappedComms(StencilImpl &st,
+				   DoubledGaugeField &U,
+				   const FermionField &in, FermionField &out, int dag);
 
   // Constructor
   WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
@@ -168,9 +171,6 @@ public:
   DoubledGaugeField UmuEven;
   DoubledGaugeField UmuOdd;
 
-  LebesgueOrder Lebesgue;
-  LebesgueOrder LebesgueEvenOdd;
-
   WilsonAnisotropyCoefficients anisotropyCoeff;
 
   ///////////////////////////////////////////////////////////////

Grid/qcd/action/fermion/WilsonFermion5D.h

@@ -135,21 +135,18 @@ public:
 		     int dag);
     
   void DhopInternal(StencilImpl & st,
-		    LebesgueOrder &lo,
 		    DoubledGaugeField &U,
 		    const FermionField &in, 
 		    FermionField &out,
 		    int dag);
 
   void DhopInternalOverlappedComms(StencilImpl & st,
-				   LebesgueOrder &lo,
 				   DoubledGaugeField &U,
 				   const FermionField &in, 
 				   FermionField &out,
 				   int dag);
 
   void DhopInternalSerialComms(StencilImpl & st,
-			       LebesgueOrder &lo,
 			       DoubledGaugeField &U,
 			       const FermionField &in, 
 			       FermionField &out,
@@ -203,9 +200,6 @@ public:
   DoubledGaugeField UmuEven;
   DoubledGaugeField UmuOdd;
     
-  LebesgueOrder Lebesgue;
-  LebesgueOrder LebesgueEvenOdd;
-    
   // Comms buffer
   //  std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  comm_buf;
 

Grid/qcd/action/fermion/ZMobiusFermion.h

@@ -58,7 +58,7 @@ public:
   {
     //    RealD eps = 1.0;
     std::cout<<GridLogMessage << "ZMobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<this->Ls<<" gamma passed in"<<std::endl;
-    Vector<Coeff_t> zgamma(this->Ls);
+    std::vector<Coeff_t> zgamma(this->Ls);
     for(int s=0;s<this->Ls;s++){
       zgamma[s] = gamma[s];
     }

Grid/qcd/action/fermion/deprecated/CayleyFermion5Dvec.h

@@ -1,3 +1,5 @@
+#if 0
+
 /*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
@@ -818,3 +820,5 @@ CayleyFermion5D<Impl>::MooeeInternal(const FermionField &psi, FermionField &chi,
 }
 
 NAMESPACE_END(Grid);
+
+#endif

Grid/qcd/action/fermion/deprecated/Lebesgue.cc

@@ -0,0 +1,245 @@
+#if 0
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/stencil/Lebesgue.cc
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: paboyle <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 <Grid/GridCore.h>
+#include <algorithm>
+
+NAMESPACE_BEGIN(Grid);
+
+int LebesgueOrder::UseLebesgueOrder;
+#ifdef KNL
+std::vector<int> LebesgueOrder::Block({8,2,2,2});
+#else
+std::vector<int> LebesgueOrder::Block({2,2,2,2});
+#endif
+LebesgueOrder::IndexInteger LebesgueOrder::alignup(IndexInteger n){
+  n--;           // 1000 0011 --> 1000 0010
+  n |= n >> 1;   // 1000 0010 | 0100 0001 = 1100 0011
+  n |= n >> 2;   // 1100 0011 | 0011 0000 = 1111 0011
+  n |= n >> 4;   // 1111 0011 | 0000 1111 = 1111 1111
+  n |= n >> 8;   // ... (At this point all bits are 1, so further bitwise-or
+  n |= n >> 16;  //      operations produce no effect.)
+  n++;           // 1111 1111 --> 1 0000 0000
+  return n;
+};
+
+LebesgueOrder::LebesgueOrder(GridBase *_grid) 
+{
+  grid = _grid;
+  if ( Block[0]==0) ZGraph();
+  else if ( Block[1]==0) NoBlocking();
+  else CartesianBlocking();
+
+  if (0) {
+    std::cout << "Thread Interleaving"<<std::endl;
+    ThreadInterleave();
+  } 
+}
+void LebesgueOrder::ThreadInterleave(void)
+{
+  Vector<IndexInteger> reorder = _LebesgueReorder;
+  Vector<IndexInteger> throrder;
+  int vol = _LebesgueReorder.size();
+  int threads = GridThread::GetThreads();
+  int blockbits=3;
+  //  int blocklen = 8;
+  //  int msk      = 0x7;
+
+  for(int t=0;t<threads;t++){
+    for(int ss=0;ss<vol;ss++){
+      if ( ( ss >> blockbits) % threads == t ) { 
+	throrder.push_back(reorder[ss]);
+      }
+    }
+  }
+  _LebesgueReorder = throrder;
+}
+void LebesgueOrder::NoBlocking(void) 
+{
+  std::cout<<GridLogDebug<<"Lexicographic : no cache blocking"<<std::endl;
+  _LebesgueReorder.resize(0);
+  for ( int s = 0 ; s!= grid->oSites();s++){
+    _LebesgueReorder.push_back(s);
+  }
+}
+void LebesgueOrder::CartesianBlocking(void) 
+{
+  _LebesgueReorder.resize(0);
+
+  //    std::cout << GridLogDebug << " CartesianBlocking ";
+  //    for(int d=0;d<Block.size();d++) std::cout <<Block[d]<<" ";
+  //    std::cout<<std::endl; 
+
+  IndexInteger ND = grid->_ndimension;
+
+  assert(ND==4);
+  assert(ND==Block.size());
+
+  Coordinate dims(ND);
+  Coordinate xo(ND,0);
+  Coordinate xi(ND,0);
+
+  for(IndexInteger mu=0;mu<ND;mu++){
+    dims[mu] = grid->_rdimensions[mu];
+  }
+
+  IterateO(ND,ND-1,xo,xi,dims);
+};
+
+void LebesgueOrder::IterateO(int ND,int dim,
+			     Coordinate & xo,
+			     Coordinate & xi,
+			     Coordinate &dims)
+{
+  for(xo[dim]=0;xo[dim]<dims[dim];xo[dim]+=Block[dim]){
+    if ( dim > 0 ) {
+      IterateO(ND,dim-1,xo,xi,dims);
+    } else {
+      IterateI(ND,ND-1,xo,xi,dims);
+    }
+  }
+};
+
+void LebesgueOrder::IterateI(int ND,
+			     int dim,
+			     Coordinate & xo,
+			     Coordinate & xi,
+			     Coordinate &dims)
+{
+  Coordinate x(ND);
+  for(xi[dim]=0;xi[dim]<std::min(dims[dim]-xo[dim],Block[dim]);xi[dim]++){
+    if ( dim > 0 ) {
+      IterateI(ND,dim-1,xo,xi,dims);
+    } else {
+      for(int d=0;d<ND;d++){
+	x[d]=xi[d]+xo[d];
+	//	std::cout << x[d]<<" ";
+      }
+      //      std::cout << "\n";
+      IndexInteger index;
+      Lexicographic::IndexFromCoor(x,index,grid->_rdimensions);
+      _LebesgueReorder.push_back(index);
+    }
+  }
+}
+
+void LebesgueOrder::ZGraph(void) 
+{
+  _LebesgueReorder.resize(0);
+
+  std::cout << GridLogDebug << " Lebesgue order "<<std::endl;
+  // Align up dimensions to power of two.
+  const IndexInteger one=1;
+
+  IndexInteger ND = grid->_ndimension;
+  Coordinate dims(ND);
+  Coordinate adims(ND);
+  std::vector<std::vector<IndexInteger> > bitlist(ND);
+  
+  for(IndexInteger mu=0;mu<ND;mu++){
+    dims[mu] = grid->_rdimensions[mu];
+    assert ( dims[mu] != 0 );
+    adims[mu] = alignup(dims[mu]);
+  }
+  
+  // List which bits of padded volume coordinate contribute; this strategy 
+  // i) avoids recursion 
+  // ii) has loop lengths at most the width of a 32 bit word.
+  int sitebit=0;
+  for(int bit=0;bit<32;bit++){
+    IndexInteger mask = one<<bit;
+    for(int mu=0;mu<ND;mu++){   // mu 0 takes bit 0; mu 1 takes bit 1 etc...
+      if ( mask&(adims[mu]-1) ){
+	bitlist[mu].push_back(sitebit);
+	sitebit++;
+      }
+    }
+  }
+  
+  // Work out padded and unpadded volumes
+  IndexInteger avol = 1;
+  for(int mu=0;mu<ND;mu++) avol = avol * adims[mu];
+  
+  IndexInteger vol = 1;
+  for(int mu=0;mu<ND;mu++) vol = vol * dims[mu];
+  
+  // Loop over padded volume, following Lebesgue curve
+  // We interleave the bits from sequential "mu".
+  std::vector<IndexInteger> ax(ND);
+  
+  for(IndexInteger asite=0;asite<avol;asite++){
+    
+    // Start with zero and collect bits
+    for(int mu=0;mu<ND;mu++) ax[mu] = 0;
+    
+    int contained = 1;
+    for(int mu=0;mu<ND;mu++){
+      
+      // Build the coordinate on the aligned volume
+      for(int bit=0;bit<bitlist[mu].size();bit++){
+	int sbit=bitlist[mu][bit];
+	
+	if(asite&(one<<sbit)){
+	  ax[mu]|=one<<bit;
+	}
+      }
+      
+      // Is it contained in original box
+      if ( ax[mu]>dims[mu]-1 ) contained = 0;
+      
+    }
+    
+    if ( contained ) {
+      int site = ax[0]
+	+        dims[0]*ax[1]
+	+dims[0]*dims[1]*ax[2]
+	+dims[0]*dims[1]*dims[2]*ax[3];
+
+      assert(site < vol);
+      _LebesgueReorder.push_back(site);
+    }
+  }
+  assert( _LebesgueReorder.size() == vol );
+
+  /*
+    std::vector<int> coor(4);
+    for(IndexInteger asite=0;asite<vol;asite++){
+    grid->oCoorFromOindex (coor,_LebesgueReorder[asite]);
+    std::cout << " site "<<asite << "->" << _LebesgueReorder[asite]<< " = ["
+    << coor[0]<<","
+    << coor[1]<<","
+    << coor[2]<<","
+    << coor[3]<<"]"
+    <<std::endl;
+    }
+  */
+}
+NAMESPACE_END(Grid);
+
+#endif

Grid/qcd/action/fermion/deprecated/Lebesgue.h

@@ -0,0 +1,81 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/stencil/Lebesgue.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: paboyle <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 GRID_LEBESGUE_H
+#define GRID_LEBESGUE_H
+
+#include<vector>
+
+// Lebesgue, Morton, Z-graph ordering assistance
+NAMESPACE_BEGIN(Grid);
+  
+class LebesgueOrder { 
+public:
+
+  typedef int32_t IndexInteger;
+  static int UseLebesgueOrder;
+  GridBase *grid;
+
+public:
+  LebesgueOrder(GridBase *_grid);
+
+  inline IndexInteger Reorder(IndexInteger ss) { 
+    return _LebesgueReorder[ss] ;
+  };
+
+  ////////////////////////////
+  // Space filling fractal for cache oblivious
+  ////////////////////////////
+  void ZGraph(void);
+  IndexInteger alignup(IndexInteger n);
+
+  /////////////////////////////////
+  // Cartesian stencil blocking strategy
+  /////////////////////////////////
+  static std::vector<int> Block;
+  void NoBlocking(void);
+  void CartesianBlocking(void);
+  void IterateO(int ND,int dim,
+		Coordinate & xo,
+		Coordinate & xi,
+		Coordinate &dims);
+  void IterateI(int ND,int dim,
+		Coordinate & xo,
+		Coordinate & xi,
+		Coordinate &dims);
+
+  void ThreadInterleave(void);
+
+private:
+  deviceVector<IndexInteger> _LebesgueReorder;
+
+};    
+
+NAMESPACE_END(Grid);
+
+#endif

Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h

@@ -156,18 +156,18 @@ template<class Impl>
 void CayleyFermion5D<Impl>::M5D   (const FermionField &psi, FermionField &chi)
 {
   int Ls=this->Ls;
-  Vector<Coeff_t> diag (Ls,1.0);
-  Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1]=mass_minus;
-  Vector<Coeff_t> lower(Ls,-1.0); lower[0]   =mass_plus;
+  std::vector<Coeff_t> diag (Ls,1.0);
+  std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1]=mass_minus;
+  std::vector<Coeff_t> lower(Ls,-1.0); lower[0]   =mass_plus;
   M5D(psi,chi,chi,lower,diag,upper);
 }
 template<class Impl>
 void CayleyFermion5D<Impl>::Meooe5D    (const FermionField &psi, FermionField &Din)
 {
   int Ls=this->Ls;
-  Vector<Coeff_t> diag = bs;
-  Vector<Coeff_t> upper= cs;
-  Vector<Coeff_t> lower= cs; 
+  std::vector<Coeff_t> diag = bs;
+  std::vector<Coeff_t> upper= cs;
+  std::vector<Coeff_t> lower= cs; 
   upper[Ls-1]=-mass_minus*upper[Ls-1];
   lower[0]   =-mass_plus*lower[0];
   M5D(psi,psi,Din,lower,diag,upper);
@@ -176,9 +176,9 @@ void CayleyFermion5D<Impl>::Meooe5D    (const FermionField &psi, FermionField &D
 template<class Impl> void CayleyFermion5D<Impl>::Meo5D     (const FermionField &psi, FermionField &chi)
 {
   int Ls=this->Ls;
-  Vector<Coeff_t> diag = beo;
-  Vector<Coeff_t> upper(Ls);
-  Vector<Coeff_t> lower(Ls);
+  std::vector<Coeff_t> diag = beo;
+  std::vector<Coeff_t> upper(Ls);
+  std::vector<Coeff_t> lower(Ls);
   for(int i=0;i<Ls;i++) {
     upper[i]=-ceo[i];
     lower[i]=-ceo[i];
@@ -191,9 +191,9 @@ template<class Impl>
 void CayleyFermion5D<Impl>::Mooee       (const FermionField &psi, FermionField &chi)
 {
   int Ls=this->Ls;
-  Vector<Coeff_t> diag = bee;
-  Vector<Coeff_t> upper(Ls);
-  Vector<Coeff_t> lower(Ls);
+  std::vector<Coeff_t> diag = bee;
+  std::vector<Coeff_t> upper(Ls);
+  std::vector<Coeff_t> lower(Ls);
   for(int i=0;i<Ls;i++) {
     upper[i]=-cee[i];
     lower[i]=-cee[i];
@@ -206,9 +206,9 @@ template<class Impl>
 void CayleyFermion5D<Impl>::MooeeDag    (const FermionField &psi, FermionField &chi)
 {
   int Ls=this->Ls;
-  Vector<Coeff_t> diag = bee;
-  Vector<Coeff_t> upper(Ls);
-  Vector<Coeff_t> lower(Ls);
+  std::vector<Coeff_t> diag = bee;
+  std::vector<Coeff_t> upper(Ls);
+  std::vector<Coeff_t> lower(Ls);
 
   for (int s=0;s<Ls;s++){
     // Assemble the 5d matrix
@@ -236,9 +236,9 @@ template<class Impl>
 void CayleyFermion5D<Impl>::M5Ddag (const FermionField &psi, FermionField &chi)
 {
   int Ls=this->Ls;
-  Vector<Coeff_t> diag(Ls,1.0);
-  Vector<Coeff_t> upper(Ls,-1.0);
-  Vector<Coeff_t> lower(Ls,-1.0);
+  std::vector<Coeff_t> diag(Ls,1.0);
+  std::vector<Coeff_t> upper(Ls,-1.0);
+  std::vector<Coeff_t> lower(Ls,-1.0);
   upper[Ls-1]=-mass_plus*upper[Ls-1];
   lower[0]   =-mass_minus*lower[0];
   M5Ddag(psi,chi,chi,lower,diag,upper);
@@ -248,9 +248,9 @@ template<class Impl>
 void CayleyFermion5D<Impl>::MeooeDag5D    (const FermionField &psi, FermionField &Din)
 {
   int Ls=this->Ls;
-  Vector<Coeff_t> diag =bs;
-  Vector<Coeff_t> upper=cs;
-  Vector<Coeff_t> lower=cs; 
+  std::vector<Coeff_t> diag =bs;
+  std::vector<Coeff_t> upper=cs;
+  std::vector<Coeff_t> lower=cs; 
 
   for (int s=0;s<Ls;s++){
     if ( s== 0 ) {
@@ -394,7 +394,7 @@ void CayleyFermion5D<Impl>::MeoDeriv(GaugeField &mat,const FermionField &U,const
 template<class Impl>
 void CayleyFermion5D<Impl>::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c)
 {
-  Vector<Coeff_t> gamma(this->Ls);
+  std::vector<Coeff_t> gamma(this->Ls);
   for(int s=0;s<this->Ls;s++) gamma[s] = zdata->gamma[s];
   SetCoefficientsInternal(1.0,gamma,b,c);
 }
@@ -402,13 +402,13 @@ void CayleyFermion5D<Impl>::SetCoefficientsTanh(Approx::zolotarev_data *zdata,Re
 template<class Impl>
 void CayleyFermion5D<Impl>::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c)
 {
-  Vector<Coeff_t> gamma(this->Ls);
+  std::vector<Coeff_t> gamma(this->Ls);
   for(int s=0;s<this->Ls;s++) gamma[s] = zdata->gamma[s];
   SetCoefficientsInternal(zolo_hi,gamma,b,c);
 }
 //Zolo
 template<class Impl>
-void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,Vector<Coeff_t> & gamma,RealD b,RealD c)
+void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,std::vector<Coeff_t> & gamma,RealD b,RealD c)
 {
   int Ls=this->Ls;
 

Grid/qcd/action/fermion/implementation/CayleyFermion5Dcache.h

@@ -43,9 +43,9 @@ void
 CayleyFermion5D<Impl>::M5D(const FermionField &psi_i,
 			       const FermionField &phi_i, 
 			       FermionField &chi_i,
-			       Vector<Coeff_t> &lower,
-			       Vector<Coeff_t> &diag,
-			       Vector<Coeff_t> &upper)
+			       std::vector<Coeff_t> &lower,
+			       std::vector<Coeff_t> &diag,
+			       std::vector<Coeff_t> &upper)
 {
   
   chi_i.Checkerboard()=psi_i.Checkerboard();
@@ -55,12 +55,16 @@ CayleyFermion5D<Impl>::M5D(const FermionField &psi_i,
   autoView(chi , chi_i,AcceleratorWrite);
   assert(phi.Checkerboard() == psi.Checkerboard());
 
-  auto pdiag = &diag[0];
-  auto pupper = &upper[0];
-  auto plower = &lower[0];
-
   int Ls =this->Ls;
 
+  static deviceVector<Coeff_t> d_diag(Ls) ; acceleratorCopyToDevice(&diag[0] ,&d_diag[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_upper(Ls); acceleratorCopyToDevice(&upper[0],&d_upper[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_lower(Ls); acceleratorCopyToDevice(&lower[0],&d_lower[0],Ls*sizeof(Coeff_t));
+  
+  auto pdiag = &d_diag[0];
+  auto pupper = &d_upper[0];
+  auto plower = &d_lower[0];
+
   // 10 = 3 complex mult + 2 complex add
   // Flops = 10.0*(Nc*Ns) *Ls*vol (/2 for red black counting)
   uint64_t nloop = grid->oSites();
@@ -82,9 +86,9 @@ void
 CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi_i,
 			      const FermionField &phi_i, 
 			      FermionField &chi_i,
-			      Vector<Coeff_t> &lower,
-			      Vector<Coeff_t> &diag,
-			      Vector<Coeff_t> &upper)
+			      std::vector<Coeff_t> &lower,
+			      std::vector<Coeff_t> &diag,
+			      std::vector<Coeff_t> &upper)
 {
   chi_i.Checkerboard()=psi_i.Checkerboard();
   GridBase *grid=psi_i.Grid();
@@ -93,12 +97,16 @@ CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi_i,
   autoView(chi , chi_i,AcceleratorWrite);
   assert(phi.Checkerboard() == psi.Checkerboard());
 
-  auto pdiag = &diag[0];
-  auto pupper = &upper[0];
-  auto plower = &lower[0];
-
   int Ls=this->Ls;
 
+  static deviceVector<Coeff_t> d_diag(Ls) ; acceleratorCopyToDevice(&diag[0] ,&d_diag[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_upper(Ls); acceleratorCopyToDevice(&upper[0],&d_upper[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_lower(Ls); acceleratorCopyToDevice(&lower[0],&d_lower[0],Ls*sizeof(Coeff_t));
+  
+  auto pdiag = &d_diag[0];
+  auto pupper = &d_upper[0];
+  auto plower = &d_lower[0];
+
   // Flops = 6.0*(Nc*Ns) *Ls*vol
   uint64_t nloop = grid->oSites();
   accelerator_for(sss,nloop,Simd::Nsimd(),{
@@ -126,11 +134,17 @@ CayleyFermion5D<Impl>::MooeeInv    (const FermionField &psi_i, FermionField &chi
 
   int Ls=this->Ls;
 
-  auto plee  = & lee [0];
-  auto pdee  = & dee [0];
-  auto puee  = & uee [0];
-  auto pleem = & leem[0];
-  auto pueem = & ueem[0];
+  static deviceVector<Coeff_t> d_lee(Ls); acceleratorCopyToDevice(&lee[0],&d_lee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_dee(Ls); acceleratorCopyToDevice(&dee[0],&d_dee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_uee(Ls); acceleratorCopyToDevice(&uee[0],&d_uee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_leem(Ls); acceleratorCopyToDevice(&leem[0],&d_leem[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_ueem(Ls); acceleratorCopyToDevice(&ueem[0],&d_ueem[0],Ls*sizeof(Coeff_t));
+
+  auto plee  = & d_lee [0];
+  auto pdee  = & d_dee [0];
+  auto puee  = & d_uee [0];
+  auto pleem = & d_leem[0];
+  auto pueem = & d_ueem[0];
 
   uint64_t nloop = grid->oSites()/Ls;
   accelerator_for(sss,nloop,Simd::Nsimd(),{
@@ -182,11 +196,17 @@ CayleyFermion5D<Impl>::MooeeInvDag (const FermionField &psi_i, FermionField &chi
   autoView(psi , psi_i,AcceleratorRead);
   autoView(chi , chi_i,AcceleratorWrite);
 
-  auto plee  = & lee [0];
-  auto pdee  = & dee [0];
-  auto puee  = & uee [0];
-  auto pleem = & leem[0];
-  auto pueem = & ueem[0];
+  static deviceVector<Coeff_t> d_lee(Ls); acceleratorCopyToDevice(&lee[0],&d_lee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_dee(Ls); acceleratorCopyToDevice(&dee[0],&d_dee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_uee(Ls); acceleratorCopyToDevice(&uee[0],&d_uee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_leem(Ls); acceleratorCopyToDevice(&leem[0],&d_leem[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_ueem(Ls); acceleratorCopyToDevice(&ueem[0],&d_ueem[0],Ls*sizeof(Coeff_t));
+
+  auto plee  = & d_lee [0];
+  auto pdee  = & d_dee [0];
+  auto puee  = & d_uee [0];
+  auto pleem = & d_leem[0];
+  auto pueem = & d_ueem[0];
 
   assert(psi.Checkerboard() == psi.Checkerboard());
 

Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionCache.h

@@ -41,7 +41,7 @@ NAMESPACE_BEGIN(Grid);
 // Pplus  backwards..
 template<class Impl>
 void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi_i, const FermionField& phi_i,FermionField& chi_i, 
-				      Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
+				      std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
 {
   chi_i.Checkerboard() = psi_i.Checkerboard();
   int Ls = this->Ls;
@@ -50,9 +50,15 @@ void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi_i, const FermionFi
   autoView( psi , psi_i, AcceleratorRead);
   autoView( chi , chi_i, AcceleratorWrite);
   assert(phi.Checkerboard() == psi.Checkerboard());
-  auto pdiag = &diag[0];
-  auto pupper = &upper[0];
-  auto plower = &lower[0];
+
+  static deviceVector<Coeff_t> d_diag(Ls); acceleratorCopyToDevice(&diag[0],&d_diag[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_upper(Ls);acceleratorCopyToDevice(&upper[0],&d_upper[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_lower(Ls);acceleratorCopyToDevice(&lower[0],&d_lower[0],Ls*sizeof(Coeff_t));
+  
+  auto pdiag = &d_diag[0];
+  auto pupper = &d_upper[0];
+  auto plower = &d_lower[0];
+
   // Flops = 6.0*(Nc*Ns) *Ls*vol
   
   auto nloop=grid->oSites()/Ls;
@@ -73,7 +79,7 @@ void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi_i, const FermionFi
 
 template<class Impl>
 void DomainWallEOFAFermion<Impl>::M5Ddag(const FermionField& psi_i, const FermionField& phi_i, FermionField& chi_i, 
-					 Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
+					 std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
 {
   chi_i.Checkerboard() = psi_i.Checkerboard();
   GridBase* grid = psi_i.Grid();
@@ -83,9 +89,14 @@ void DomainWallEOFAFermion<Impl>::M5Ddag(const FermionField& psi_i, const Fermio
   autoView( phi , phi_i, AcceleratorRead);
   autoView( chi , chi_i, AcceleratorWrite);
   assert(phi.Checkerboard() == psi.Checkerboard());
-  auto pdiag = &diag[0];
-  auto pupper = &upper[0];
-  auto plower = &lower[0];
+
+  static deviceVector<Coeff_t> d_diag(Ls); acceleratorCopyToDevice(&diag[0],&d_diag[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_upper(Ls);acceleratorCopyToDevice(&upper[0],&d_upper[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_lower(Ls);acceleratorCopyToDevice(&lower[0],&d_lower[0],Ls*sizeof(Coeff_t));
+  
+  auto pdiag = &d_diag[0];
+  auto pupper = &d_upper[0];
+  auto plower = &d_lower[0];
 
   // Flops = 6.0*(Nc*Ns) *Ls*vol
 
@@ -114,13 +125,18 @@ void DomainWallEOFAFermion<Impl>::MooeeInv(const FermionField& psi_i, FermionFie
   autoView( chi, chi_i, AcceleratorWrite);
   int Ls = this->Ls;
 
-  auto plee  = & this->lee[0];
-  auto pdee  = & this->dee[0];
-  auto puee  = & this->uee[0];
-
-  auto pleem = & this->leem[0];
-  auto pueem = & this->ueem[0];
+  static deviceVector<Coeff_t> d_lee(Ls); acceleratorCopyToDevice(&this->lee[0],&d_lee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_dee(Ls); acceleratorCopyToDevice(&this->dee[0],&d_dee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_uee(Ls); acceleratorCopyToDevice(&this->uee[0],&d_uee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_leem(Ls); acceleratorCopyToDevice(&this->leem[0],&d_leem[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_ueem(Ls); acceleratorCopyToDevice(&this->ueem[0],&d_ueem[0],Ls*sizeof(Coeff_t));
 
+  auto plee  = & d_lee [0];
+  auto pdee  = & d_dee [0];
+  auto puee  = & d_uee [0];
+  auto pleem = & d_leem[0];
+  auto pueem = & d_ueem[0];
+  
   uint64_t nloop=grid->oSites()/Ls;
   accelerator_for(sss,nloop,Simd::Nsimd(),{
     uint64_t ss=sss*Ls;

Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h

@@ -131,9 +131,9 @@ void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi, FermionField& chi
     else{ shiftm = -shift*(mq3-mq2); }
   }
 
-  Vector<Coeff_t> diag(Ls,1.0);
-  Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftm;
-  Vector<Coeff_t> lower(Ls,-1.0); lower[0]    = mq1 + shiftp;
+  std::vector<Coeff_t> diag(Ls,1.0);
+  std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftm;
+  std::vector<Coeff_t> lower(Ls,-1.0); lower[0]    = mq1 + shiftp;
 
 #if(0)
   std::cout << GridLogMessage << "DomainWallEOFAFermion::M5D(FF&,FF&):" << std::endl;
@@ -168,9 +168,9 @@ void DomainWallEOFAFermion<Impl>::M5Ddag(const FermionField& psi, FermionField&
     else{ shiftm = -shift*(mq3-mq2); }
   }
 
-  Vector<Coeff_t> diag(Ls,1.0);
-  Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftp;
-  Vector<Coeff_t> lower(Ls,-1.0); lower[0]    = mq1 + shiftm;
+  std::vector<Coeff_t> diag(Ls,1.0);
+  std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftp;
+  std::vector<Coeff_t> lower(Ls,-1.0); lower[0]    = mq1 + shiftm;
 
   this->M5Ddag(psi, chi, chi, lower, diag, upper);
 }
@@ -181,9 +181,9 @@ void DomainWallEOFAFermion<Impl>::Mooee(const FermionField& psi, FermionField& c
 {
   int Ls = this->Ls;
 
-  Vector<Coeff_t> diag = this->bee;
-  Vector<Coeff_t> upper(Ls);
-  Vector<Coeff_t> lower(Ls);
+  std::vector<Coeff_t> diag = this->bee;
+  std::vector<Coeff_t> upper(Ls);
+  std::vector<Coeff_t> lower(Ls);
 
   for(int s=0; s<Ls; s++){
     upper[s] = -this->cee[s];
@@ -200,9 +200,9 @@ void DomainWallEOFAFermion<Impl>::MooeeDag(const FermionField& psi, FermionField
 {
   int Ls = this->Ls;
 
-  Vector<Coeff_t> diag = this->bee;
-  Vector<Coeff_t> upper(Ls);
-  Vector<Coeff_t> lower(Ls);
+  std::vector<Coeff_t> diag = this->bee;
+  std::vector<Coeff_t> upper(Ls);
+  std::vector<Coeff_t> lower(Ls);
 
   for(int s=0; s<Ls; s++){
     upper[s] = -this->cee[s];
@@ -218,7 +218,7 @@ void DomainWallEOFAFermion<Impl>::MooeeDag(const FermionField& psi, FermionField
 
 //Zolo
 template<class Impl>
-void DomainWallEOFAFermion<Impl>::SetCoefficientsInternal(RealD zolo_hi, Vector<Coeff_t>& gamma, RealD b, RealD c)
+void DomainWallEOFAFermion<Impl>::SetCoefficientsInternal(RealD zolo_hi, std::vector<Coeff_t>& gamma, RealD b, RealD c)
 {
   int   Ls    = this->Ls;
   int   pm    = this->pm;

Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermion5DImplementation.h

@@ -61,8 +61,6 @@ ImprovedStaggeredFermion5D<Impl>::ImprovedStaggeredFermion5D(GridCartesian
   UUUmu(&FourDimGrid),
   UUUmuEven(&FourDimRedBlackGrid),
   UUUmuOdd(&FourDimRedBlackGrid),
-  Lebesgue(&FourDimGrid),
-  LebesgueEvenOdd(&FourDimRedBlackGrid),
   _tmp(&FiveDimRedBlackGrid)
 {
 
@@ -277,18 +275,18 @@ void ImprovedStaggeredFermion5D<Impl>::DhopDerivOE(GaugeField &mat,
 
 /*CHANGE */
 template<class Impl>
-void ImprovedStaggeredFermion5D<Impl>::DhopInternal(StencilImpl & st, LebesgueOrder &lo,
+void ImprovedStaggeredFermion5D<Impl>::DhopInternal(StencilImpl & st, 
 						    DoubledGaugeField & U,DoubledGaugeField & UUU,
 						    const FermionField &in, FermionField &out,int dag)
 {
   if ( StaggeredKernelsStatic::Comms == StaggeredKernelsStatic::CommsAndCompute )
-    DhopInternalOverlappedComms(st,lo,U,UUU,in,out,dag);
+    DhopInternalOverlappedComms(st,U,UUU,in,out,dag);
   else
-    DhopInternalSerialComms(st,lo,U,UUU,in,out,dag);
+    DhopInternalSerialComms(st,U,UUU,in,out,dag);
 }
 
 template<class Impl>
-void ImprovedStaggeredFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl & st, LebesgueOrder &lo,
+void ImprovedStaggeredFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl & st, 
 								   DoubledGaugeField & U,DoubledGaugeField & UUU,
 								   const FermionField &in, FermionField &out,int dag)
 {
@@ -313,7 +311,7 @@ void ImprovedStaggeredFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl &
   {
     int interior=1;
     int exterior=0;
-    Kernels::DhopImproved(st,lo,U,UUU,in,out,dag,interior,exterior);
+    Kernels::DhopImproved(st,U,UUU,in,out,dag,interior,exterior);
   }
 
   st.CommsMerge(compressor);
@@ -323,12 +321,12 @@ void ImprovedStaggeredFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl &
   {
     int interior=0;
     int exterior=1;
-    Kernels::DhopImproved(st,lo,U,UUU,in,out,dag,interior,exterior);
+    Kernels::DhopImproved(st,U,UUU,in,out,dag,interior,exterior);
   }
 }
 
 template<class Impl>
-void ImprovedStaggeredFermion5D<Impl>::DhopInternalSerialComms(StencilImpl & st, LebesgueOrder &lo,
+void ImprovedStaggeredFermion5D<Impl>::DhopInternalSerialComms(StencilImpl & st, 
 						    DoubledGaugeField & U,DoubledGaugeField & UUU,
 						    const FermionField &in, FermionField &out,int dag)
 {
@@ -341,7 +339,7 @@ void ImprovedStaggeredFermion5D<Impl>::DhopInternalSerialComms(StencilImpl & st,
   {
     int interior=1;
     int exterior=1;
-    Kernels::DhopImproved(st,lo,U,UUU,in,out,dag,interior,exterior);
+    Kernels::DhopImproved(st,U,UUU,in,out,dag,interior,exterior);
   }
 }
 /*CHANGE END*/
@@ -357,7 +355,7 @@ void ImprovedStaggeredFermion5D<Impl>::DhopOE(const FermionField &in, FermionFie
   assert(in.Checkerboard()==Even);
   out.Checkerboard() = Odd;
 
-  DhopInternal(StencilEven,LebesgueEvenOdd,UmuOdd,UUUmuOdd,in,out,dag);
+  DhopInternal(StencilEven,UmuOdd,UUUmuOdd,in,out,dag);
 }
 template<class Impl>
 void ImprovedStaggeredFermion5D<Impl>::DhopEO(const FermionField &in, FermionField &out,int dag)
@@ -368,7 +366,7 @@ void ImprovedStaggeredFermion5D<Impl>::DhopEO(const FermionField &in, FermionFie
   assert(in.Checkerboard()==Odd);
   out.Checkerboard() = Even;
 
-  DhopInternal(StencilOdd,LebesgueEvenOdd,UmuEven,UUUmuEven,in,out,dag);
+  DhopInternal(StencilOdd,UmuEven,UUUmuEven,in,out,dag);
 }
 template<class Impl>
 void ImprovedStaggeredFermion5D<Impl>::Dhop(const FermionField &in, FermionField &out,int dag)
@@ -378,7 +376,7 @@ void ImprovedStaggeredFermion5D<Impl>::Dhop(const FermionField &in, FermionField
 
   out.Checkerboard() = in.Checkerboard();
 
-  DhopInternal(Stencil,Lebesgue,Umu,UUUmu,in,out,dag);
+  DhopInternal(Stencil,Umu,UUUmu,in,out,dag);
 }
 
 /////////////////////////////////////////////////////////////////////////

Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermionImplementation.h

@@ -48,8 +48,6 @@ ImprovedStaggeredFermion<Impl>::ImprovedStaggeredFermion(GridCartesian &Fgrid, G
     StencilEven(&Hgrid, npoint, Even, directions, displacements,p),  // source is Even
     StencilOdd(&Hgrid, npoint, Odd, directions, displacements,p),  // source is Odd
     mass(_mass),
-    Lebesgue(_grid),
-    LebesgueEvenOdd(_cbgrid),
     Umu(&Fgrid),
     UmuEven(&Hgrid),
     UmuOdd(&Hgrid),
@@ -339,7 +337,7 @@ void ImprovedStaggeredFermion<Impl>::Dhop(const FermionField &in, FermionField &
 
   out.Checkerboard() = in.Checkerboard();
 
-  DhopInternal(Stencil, Lebesgue, Umu, UUUmu, in, out, dag);
+  DhopInternal(Stencil, Umu, UUUmu, in, out, dag);
 }
 
 template <class Impl>
@@ -351,7 +349,7 @@ void ImprovedStaggeredFermion<Impl>::DhopOE(const FermionField &in, FermionField
   assert(in.Checkerboard() == Even);
   out.Checkerboard() = Odd;
 
-  DhopInternal(StencilEven, LebesgueEvenOdd, UmuOdd, UUUmuOdd, in, out, dag);
+  DhopInternal(StencilEven, UmuOdd, UUUmuOdd, in, out, dag);
 }
 
 template <class Impl>
@@ -363,7 +361,7 @@ void ImprovedStaggeredFermion<Impl>::DhopEO(const FermionField &in, FermionField
   assert(in.Checkerboard() == Odd);
   out.Checkerboard() = Even;
 
-  DhopInternal(StencilOdd, LebesgueEvenOdd, UmuEven, UUUmuEven, in, out, dag);
+  DhopInternal(StencilOdd, UmuEven, UUUmuEven, in, out, dag);
 }
 
 template <class Impl>
@@ -394,19 +392,19 @@ void ImprovedStaggeredFermion<Impl>::DhopDir(const FermionField &in, FermionFiel
 
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder &lo,
+void ImprovedStaggeredFermion<Impl>::DhopInternal(StencilImpl &st, 
 						  DoubledGaugeField &U,
 						  DoubledGaugeField &UUU,
 						  const FermionField &in,
 						  FermionField &out, int dag) 
 {
   if ( StaggeredKernelsStatic::Comms == StaggeredKernelsStatic::CommsAndCompute )
-    DhopInternalOverlappedComms(st,lo,U,UUU,in,out,dag);
+    DhopInternalOverlappedComms(st,U,UUU,in,out,dag);
   else
-    DhopInternalSerialComms(st,lo,U,UUU,in,out,dag);
+    DhopInternalSerialComms(st,U,UUU,in,out,dag);
 }
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, LebesgueOrder &lo,
+void ImprovedStaggeredFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, 
 								 DoubledGaugeField &U,
 								 DoubledGaugeField &UUU,
 								 const FermionField &in,
@@ -429,7 +427,7 @@ void ImprovedStaggeredFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st
   {
     int interior=1;
     int exterior=0;
-    Kernels::DhopImproved(st,lo,U,UUU,in,out,dag,interior,exterior);
+    Kernels::DhopImproved(st,U,UUU,in,out,dag,interior,exterior);
   }
 
   st.CommunicateComplete(requests);
@@ -440,13 +438,13 @@ void ImprovedStaggeredFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st
   {
     int interior=0;
     int exterior=1;
-    Kernels::DhopImproved(st,lo,U,UUU,in,out,dag,interior,exterior);
+    Kernels::DhopImproved(st,U,UUU,in,out,dag,interior,exterior);
   }
 }
 
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::DhopInternalSerialComms(StencilImpl &st, LebesgueOrder &lo,
+void ImprovedStaggeredFermion<Impl>::DhopInternalSerialComms(StencilImpl &st, 
 							     DoubledGaugeField &U,
 							     DoubledGaugeField &UUU,
 							     const FermionField &in,
@@ -460,7 +458,7 @@ void ImprovedStaggeredFermion<Impl>::DhopInternalSerialComms(StencilImpl &st, Le
   {
     int interior=1;
     int exterior=1;
-    Kernels::DhopImproved(st,lo,U,UUU,in,out,dag,interior,exterior);
+    Kernels::DhopImproved(st,U,UUU,in,out,dag,interior,exterior);
   }
 };
 

Grid/qcd/action/fermion/implementation/MobiusEOFAFermionCache.h

@@ -39,7 +39,7 @@ NAMESPACE_BEGIN(Grid);
  
 template<class Impl>
 void MobiusEOFAFermion<Impl>::M5D(const FermionField &psi_i, const FermionField &phi_i, FermionField &chi_i,
-				  Vector<Coeff_t> &lower, Vector<Coeff_t> &diag, Vector<Coeff_t> &upper)
+				  std::vector<Coeff_t> &lower, std::vector<Coeff_t> &diag, std::vector<Coeff_t> &upper)
 {
   chi_i.Checkerboard() = psi_i.Checkerboard();
   GridBase *grid = psi_i.Grid();
@@ -50,9 +50,13 @@ void MobiusEOFAFermion<Impl>::M5D(const FermionField &psi_i, const FermionField
 
   assert(phi.Checkerboard() == psi.Checkerboard());
 
-  auto pdiag = &diag[0];
-  auto pupper = &upper[0];
-  auto plower = &lower[0];
+  static deviceVector<Coeff_t> d_diag(Ls); acceleratorCopyToDevice(&diag[0],&d_diag[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_upper(Ls);acceleratorCopyToDevice(&upper[0],&d_upper[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_lower(Ls);acceleratorCopyToDevice(&lower[0],&d_lower[0],Ls*sizeof(Coeff_t));
+  
+  auto pdiag = &d_diag[0];
+  auto pupper = &d_upper[0];
+  auto plower = &d_lower[0];
 
   // Flops = 6.0*(Nc*Ns) *Ls*vol
   int nloop = grid->oSites()/Ls;
@@ -74,8 +78,8 @@ void MobiusEOFAFermion<Impl>::M5D(const FermionField &psi_i, const FermionField
 
 template<class Impl>
 void MobiusEOFAFermion<Impl>::M5D_shift(const FermionField &psi_i, const FermionField &phi_i, FermionField &chi_i,
-					Vector<Coeff_t> &lower, Vector<Coeff_t> &diag, Vector<Coeff_t> &upper,
-					Vector<Coeff_t> &shift_coeffs)
+					std::vector<Coeff_t> &lower, std::vector<Coeff_t> &diag, std::vector<Coeff_t> &upper,
+					std::vector<Coeff_t> &shift_coeffs)
 {
   chi_i.Checkerboard() = psi_i.Checkerboard();
   GridBase *grid = psi_i.Grid();
@@ -86,13 +90,18 @@ void MobiusEOFAFermion<Impl>::M5D_shift(const FermionField &psi_i, const Fermion
 
   auto pm  = this->pm;
   int shift_s = (pm == 1) ? (Ls-1) : 0; // s-component modified by shift operator
-
+  
   assert(phi.Checkerboard() == psi.Checkerboard());
 
-  auto pdiag = &diag[0];
-  auto pupper = &upper[0];
-  auto plower = &lower[0];
-  auto pshift_coeffs = &shift_coeffs[0];
+  static deviceVector<Coeff_t> d_diag(Ls); acceleratorCopyToDevice(&diag[0],&d_diag[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_upper(Ls);acceleratorCopyToDevice(&upper[0],&d_upper[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_lower(Ls);acceleratorCopyToDevice(&lower[0],&d_lower[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_shift_coeffs(Ls);acceleratorCopyToDevice(&shift_coeffs[0],&d_shift_coeffs[0],Ls*sizeof(Coeff_t));
+  
+  auto pdiag = &d_diag[0];
+  auto pupper = &d_upper[0];
+  auto plower = &d_lower[0];
+  auto pshift_coeffs = &d_shift_coeffs[0];
 
   // Flops = 6.0*(Nc*Ns) *Ls*vol
   int nloop = grid->oSites()/Ls;
@@ -119,7 +128,7 @@ void MobiusEOFAFermion<Impl>::M5D_shift(const FermionField &psi_i, const Fermion
 
 template<class Impl>
 void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField &psi_i, const FermionField &phi_i, FermionField &chi_i,
-				     Vector<Coeff_t> &lower, Vector<Coeff_t> &diag, Vector<Coeff_t> &upper)
+				     std::vector<Coeff_t> &lower, std::vector<Coeff_t> &diag, std::vector<Coeff_t> &upper)
 {
   chi_i.Checkerboard() = psi_i.Checkerboard();
   GridBase *grid = psi_i.Grid();
@@ -130,9 +139,13 @@ void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField &psi_i, const FermionFie
 
   assert(phi.Checkerboard() == psi.Checkerboard());
 
-  auto pdiag = &diag[0];
-  auto pupper = &upper[0];
-  auto plower = &lower[0];
+  static deviceVector<Coeff_t> d_diag(Ls); acceleratorCopyToDevice(&diag[0],&d_diag[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_upper(Ls);acceleratorCopyToDevice(&upper[0],&d_upper[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_lower(Ls);acceleratorCopyToDevice(&lower[0],&d_lower[0],Ls*sizeof(Coeff_t));
+  
+  auto pdiag = &d_diag[0];
+  auto pupper = &d_upper[0];
+  auto plower = &d_lower[0];
 
   // Flops = 6.0*(Nc*Ns) *Ls*vol
   int nloop = grid->oSites()/Ls;
@@ -154,8 +167,8 @@ void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField &psi_i, const FermionFie
 
 template<class Impl>
 void MobiusEOFAFermion<Impl>::M5Ddag_shift(const FermionField &psi_i, const FermionField &phi_i, FermionField &chi_i,
-					   Vector<Coeff_t> &lower, Vector<Coeff_t> &diag, Vector<Coeff_t> &upper,
-					   Vector<Coeff_t> &shift_coeffs)
+					   std::vector<Coeff_t> &lower, std::vector<Coeff_t> &diag, std::vector<Coeff_t> &upper,
+					   std::vector<Coeff_t> &shift_coeffs)
 {
   chi_i.Checkerboard() = psi_i.Checkerboard();
   GridBase *grid = psi_i.Grid();
@@ -167,10 +180,15 @@ void MobiusEOFAFermion<Impl>::M5Ddag_shift(const FermionField &psi_i, const Ferm
 
   assert(phi.Checkerboard() == psi.Checkerboard());
 
-  auto pdiag = &diag[0];
-  auto pupper = &upper[0];
-  auto plower = &lower[0];
-  auto pshift_coeffs = &shift_coeffs[0];
+  static deviceVector<Coeff_t> d_diag(Ls); acceleratorCopyToDevice(&diag[0],&d_diag[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_upper(Ls);acceleratorCopyToDevice(&upper[0],&d_upper[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_lower(Ls);acceleratorCopyToDevice(&lower[0],&d_lower[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_shift_coeffs(Ls);acceleratorCopyToDevice(&shift_coeffs[0],&d_shift_coeffs[0],Ls*sizeof(Coeff_t));
+  
+  auto pdiag = &d_diag[0];
+  auto pupper = &d_upper[0];
+  auto plower = &d_lower[0];
+  auto pshift_coeffs = &d_shift_coeffs[0];
 
   // Flops = 6.0*(Nc*Ns) *Ls*vol
   auto pm = this->pm;
@@ -212,11 +230,17 @@ void MobiusEOFAFermion<Impl>::MooeeInv(const FermionField &psi_i, FermionField &
   autoView(psi , psi_i, AcceleratorRead);
   autoView(chi , chi_i, AcceleratorWrite);
 
-  auto plee = & this->lee [0];
-  auto pdee = & this->dee [0];
-  auto puee = & this->uee [0];
-  auto pleem= & this->leem[0];
-  auto pueem= & this->ueem[0];
+  static deviceVector<Coeff_t> d_lee(Ls); acceleratorCopyToDevice(&this->lee[0],&d_lee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_dee(Ls); acceleratorCopyToDevice(&this->dee[0],&d_dee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_uee(Ls); acceleratorCopyToDevice(&this->uee[0],&d_uee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_leem(Ls); acceleratorCopyToDevice(&this->leem[0],&d_leem[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_ueem(Ls); acceleratorCopyToDevice(&this->ueem[0],&d_ueem[0],Ls*sizeof(Coeff_t));
+
+  auto plee  = & d_lee [0];
+  auto pdee  = & d_dee [0];
+  auto puee  = & d_uee [0];
+  auto pleem = & d_leem[0];
+  auto pueem = & d_ueem[0];
 
   if(this->shift != 0.0){ MooeeInv_shift(psi_i,chi_i); return; }
 
@@ -268,14 +292,24 @@ void MobiusEOFAFermion<Impl>::MooeeInv_shift(const FermionField &psi_i, FermionF
   autoView(psi , psi_i, AcceleratorRead);
   autoView(chi , chi_i, AcceleratorWrite);
 
+  // Move into object and constructor
+  static deviceVector<Coeff_t> d_lee(Ls); acceleratorCopyToDevice(&this->lee[0],&d_lee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_dee(Ls); acceleratorCopyToDevice(&this->dee[0],&d_dee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_uee(Ls); acceleratorCopyToDevice(&this->uee[0],&d_uee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_leem(Ls); acceleratorCopyToDevice(&this->leem[0],&d_leem[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_ueem(Ls); acceleratorCopyToDevice(&this->ueem[0],&d_ueem[0],Ls*sizeof(Coeff_t));
+
   auto pm = this->pm;
-  auto plee = & this->lee [0];
-  auto pdee = & this->dee [0];
-  auto puee = & this->uee [0];
-  auto pleem= & this->leem[0];
-  auto pueem= & this->ueem[0];
-  auto pMooeeInv_shift_lc   = &MooeeInv_shift_lc[0];
-  auto pMooeeInv_shift_norm = &MooeeInv_shift_norm[0];
+  auto plee  = & d_lee [0];
+  auto pdee  = & d_dee [0];
+  auto puee  = & d_uee [0];
+  auto pleem = & d_leem[0];
+  auto pueem = & d_ueem[0];
+
+  static deviceVector<Coeff_t> d_MooeeInv_shift_lc(Ls); acceleratorCopyToDevice(&MooeeInv_shift_lc[0],&d_MooeeInv_shift_lc[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_MooeeInv_shift_norm(Ls); acceleratorCopyToDevice(&MooeeInv_shift_norm[0],&d_MooeeInv_shift_norm[0],Ls*sizeof(Coeff_t));
+  auto pMooeeInv_shift_lc   = &d_MooeeInv_shift_lc[0];
+  auto pMooeeInv_shift_norm = &d_MooeeInv_shift_norm[0];
 
   int nloop = grid->oSites()/Ls;
   accelerator_for(sss,nloop,Simd::Nsimd(),{
@@ -333,11 +367,17 @@ void MobiusEOFAFermion<Impl>::MooeeInvDag(const FermionField &psi_i, FermionFiel
   autoView(psi , psi_i, AcceleratorRead);
   autoView(chi , chi_i, AcceleratorWrite);
 
-  auto plee = & this->lee [0];
-  auto pdee = & this->dee [0];
-  auto puee = & this->uee [0];
-  auto pleem= & this->leem[0];
-  auto pueem= & this->ueem[0];
+  static deviceVector<Coeff_t> d_lee(Ls); acceleratorCopyToDevice(&this->lee[0],&d_lee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_dee(Ls); acceleratorCopyToDevice(&this->dee[0],&d_dee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_uee(Ls); acceleratorCopyToDevice(&this->uee[0],&d_uee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_leem(Ls); acceleratorCopyToDevice(&this->leem[0],&d_leem[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_ueem(Ls); acceleratorCopyToDevice(&this->ueem[0],&d_ueem[0],Ls*sizeof(Coeff_t));
+
+  auto plee  = & d_lee [0];
+  auto pdee  = & d_dee [0];
+  auto puee  = & d_uee [0];
+  auto pleem = & d_leem[0];
+  auto pueem = & d_ueem[0];
 
   int nloop = grid->oSites()/Ls;
   accelerator_for(sss,nloop,Simd::Nsimd(),{
@@ -386,14 +426,28 @@ void MobiusEOFAFermion<Impl>::MooeeInvDag_shift(const FermionField &psi_i, Fermi
   autoView(chi , chi_i, AcceleratorWrite);
   int Ls = this->Ls;
 
+  static deviceVector<Coeff_t> d_lee(Ls); acceleratorCopyToDevice(&this->lee[0],&d_lee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_dee(Ls); acceleratorCopyToDevice(&this->dee[0],&d_dee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_uee(Ls); acceleratorCopyToDevice(&this->uee[0],&d_uee[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_leem(Ls); acceleratorCopyToDevice(&this->leem[0],&d_leem[0],Ls*sizeof(Coeff_t));
+  static deviceVector<Coeff_t> d_ueem(Ls); acceleratorCopyToDevice(&this->ueem[0],&d_ueem[0],Ls*sizeof(Coeff_t));
+
   auto pm = this->pm;
-  auto plee = & this->lee [0];
-  auto pdee = & this->dee [0];
-  auto puee = & this->uee [0];
-  auto pleem= & this->leem[0];
-  auto pueem= & this->ueem[0];
-  auto pMooeeInvDag_shift_lc   = &MooeeInvDag_shift_lc[0];
-  auto pMooeeInvDag_shift_norm = &MooeeInvDag_shift_norm[0];
+  auto plee  = & d_lee [0];
+  auto pdee  = & d_dee [0];
+  auto puee  = & d_uee [0];
+  auto pleem = & d_leem[0];
+  auto pueem = & d_ueem[0];
+
+  static deviceVector<Coeff_t> d_MooeeInvDag_shift_lc(Ls);
+  static deviceVector<Coeff_t> d_MooeeInvDag_shift_norm(Ls);
+  acceleratorCopyToDevice(&MooeeInvDag_shift_lc[0],&d_MooeeInvDag_shift_lc[0],Ls*sizeof(Coeff_t));
+  acceleratorCopyToDevice(&MooeeInvDag_shift_norm[0],&d_MooeeInvDag_shift_norm[0],Ls*sizeof(Coeff_t));
+  auto pMooeeInvDag_shift_lc   = &d_MooeeInvDag_shift_lc[0];
+  auto pMooeeInvDag_shift_norm = &d_MooeeInvDag_shift_norm[0];
+
+  //  auto pMooeeInvDag_shift_lc   = &MooeeInvDag_shift_lc[0];
+  //  auto pMooeeInvDag_shift_norm = &MooeeInvDag_shift_norm[0];
 
   int nloop = grid->oSites()/Ls;
   accelerator_for(sss,nloop,Simd::Nsimd(),{

Grid/qcd/action/fermion/implementation/MobiusEOFAFermionImplementation.h

@@ -196,9 +196,9 @@ void MobiusEOFAFermion<Impl>::M5D(const FermionField& psi, FermionField& chi)
 {
   int Ls = this->Ls;
 
-  Vector<Coeff_t> diag(Ls,1.0);
-  Vector<Coeff_t> upper(Ls,-1.0);  upper[Ls-1] = this->mq1;
-  Vector<Coeff_t> lower(Ls,-1.0);  lower[0]    = this->mq1;
+  std::vector<Coeff_t> diag(Ls,1.0);
+  std::vector<Coeff_t> upper(Ls,-1.0);  upper[Ls-1] = this->mq1;
+  std::vector<Coeff_t> lower(Ls,-1.0);  lower[0]    = this->mq1;
 
   // no shift term
   if(this->shift == 0.0){ this->M5D(psi, chi, chi, lower, diag, upper); }
@@ -212,9 +212,9 @@ void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField& psi, FermionField& chi)
 {
   int Ls = this->Ls;
 
-  Vector<Coeff_t> diag(Ls,1.0);
-  Vector<Coeff_t> upper(Ls,-1.0);  upper[Ls-1] = this->mq1;
-  Vector<Coeff_t> lower(Ls,-1.0);  lower[0]    = this->mq1;
+  std::vector<Coeff_t> diag(Ls,1.0);
+  std::vector<Coeff_t> upper(Ls,-1.0);  upper[Ls-1] = this->mq1;
+  std::vector<Coeff_t> lower(Ls,-1.0);  lower[0]    = this->mq1;
 
   // no shift term
   if(this->shift == 0.0){ this->M5Ddag(psi, chi, chi, lower, diag, upper); }
@@ -230,9 +230,9 @@ void MobiusEOFAFermion<Impl>::Mooee(const FermionField& psi, FermionField& chi)
   int Ls = this->Ls;
 
   // coefficients of Mooee
-  Vector<Coeff_t> diag = this->bee;
-  Vector<Coeff_t> upper(Ls);
-  Vector<Coeff_t> lower(Ls);
+  std::vector<Coeff_t> diag = this->bee;
+  std::vector<Coeff_t> upper(Ls);
+  std::vector<Coeff_t> lower(Ls);
   for(int s=0; s<Ls; s++){
     upper[s] = -this->cee[s];
     lower[s] = -this->cee[s];
@@ -253,9 +253,9 @@ void MobiusEOFAFermion<Impl>::MooeeDag(const FermionField& psi, FermionField& ch
   int Ls = this->Ls;
 
   // coefficients of MooeeDag
-  Vector<Coeff_t> diag = this->bee;
-  Vector<Coeff_t> upper(Ls);
-  Vector<Coeff_t> lower(Ls);
+  std::vector<Coeff_t> diag = this->bee;
+  std::vector<Coeff_t> upper(Ls);
+  std::vector<Coeff_t> lower(Ls);
   for(int s=0; s<Ls; s++){
     if(s==0) {
       upper[s] = -this->cee[s+1];
@@ -314,10 +314,10 @@ void MobiusEOFAFermion<Impl>::SetCoefficientsPrecondShiftOps()
   // Tridiagonal solve for MooeeInvDag_shift_lc
   {
     Coeff_t m(0.0);
-    Vector<Coeff_t> d = Mooee_shift;
-    Vector<Coeff_t> u(Ls,0.0);
-    Vector<Coeff_t> y(Ls,0.0);
-    Vector<Coeff_t> q(Ls,0.0);
+    std::vector<Coeff_t> d = Mooee_shift;
+    std::vector<Coeff_t> u(Ls,0.0);
+    std::vector<Coeff_t> y(Ls,0.0);
+    std::vector<Coeff_t> q(Ls,0.0);
     if(pm == 1){ u[0] = 1.0; }
     else{ u[Ls-1] = 1.0; }
 

Grid/qcd/action/fermion/implementation/NaiveStaggeredFermionImplementation.h

@@ -48,8 +48,6 @@ NaiveStaggeredFermion<Impl>::NaiveStaggeredFermion(GridCartesian &Fgrid, GridRed
     StencilEven(&Hgrid, npoint, Even, directions, displacements,p),  // source is Even
     StencilOdd(&Hgrid, npoint, Odd, directions, displacements,p),  // source is Odd
     mass(_mass),
-    Lebesgue(_grid),
-    LebesgueEvenOdd(_cbgrid),
     Umu(&Fgrid),
     UmuEven(&Hgrid),
     UmuOdd(&Hgrid),
@@ -268,7 +266,7 @@ void NaiveStaggeredFermion<Impl>::Dhop(const FermionField &in, FermionField &out
 
   out.Checkerboard() = in.Checkerboard();
 
-  DhopInternal(Stencil, Lebesgue, Umu, in, out, dag);
+  DhopInternal(Stencil, Umu, in, out, dag);
 }
 
 template <class Impl>
@@ -280,7 +278,7 @@ void NaiveStaggeredFermion<Impl>::DhopOE(const FermionField &in, FermionField &o
   assert(in.Checkerboard() == Even);
   out.Checkerboard() = Odd;
 
-  DhopInternal(StencilEven, LebesgueEvenOdd, UmuOdd, in, out, dag);
+  DhopInternal(StencilEven, UmuOdd, in, out, dag);
 }
 
 template <class Impl>
@@ -292,7 +290,7 @@ void NaiveStaggeredFermion<Impl>::DhopEO(const FermionField &in, FermionField &o
   assert(in.Checkerboard() == Odd);
   out.Checkerboard() = Even;
 
-  DhopInternal(StencilOdd, LebesgueEvenOdd, UmuEven, in, out, dag);
+  DhopInternal(StencilOdd, UmuEven, in, out, dag);
 }
 
 template <class Impl>
@@ -323,18 +321,18 @@ void NaiveStaggeredFermion<Impl>::DhopDir(const FermionField &in, FermionField &
 
 
 template <class Impl>
-void NaiveStaggeredFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder &lo,
+void NaiveStaggeredFermion<Impl>::DhopInternal(StencilImpl &st,
 					       DoubledGaugeField &U,
 					       const FermionField &in,
 					       FermionField &out, int dag) 
 {
   if ( StaggeredKernelsStatic::Comms == StaggeredKernelsStatic::CommsAndCompute )
-    DhopInternalOverlappedComms(st,lo,U,in,out,dag);
+    DhopInternalOverlappedComms(st,U,in,out,dag);
   else
-    DhopInternalSerialComms(st,lo,U,in,out,dag);
+    DhopInternalSerialComms(st,U,in,out,dag);
 }
 template <class Impl>
-void NaiveStaggeredFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, LebesgueOrder &lo,
+void NaiveStaggeredFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st,
 							      DoubledGaugeField &U,
 							      const FermionField &in,
 							      FermionField &out, int dag) 
@@ -356,7 +354,7 @@ void NaiveStaggeredFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, L
   {
     int interior=1;
     int exterior=0;
-    Kernels::DhopNaive(st,lo,U,in,out,dag,interior,exterior);
+    Kernels::DhopNaive(st,U,in,out,dag,interior,exterior);
   }
 
   st.CommunicateComplete(requests);
@@ -367,12 +365,12 @@ void NaiveStaggeredFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, L
   {
     int interior=0;
     int exterior=1;
-    Kernels::DhopNaive(st,lo,U,in,out,dag,interior,exterior);
+    Kernels::DhopNaive(st,U,in,out,dag,interior,exterior);
   }
 }
 
 template <class Impl>
-void NaiveStaggeredFermion<Impl>::DhopInternalSerialComms(StencilImpl &st, LebesgueOrder &lo,
+void NaiveStaggeredFermion<Impl>::DhopInternalSerialComms(StencilImpl &st,
 							  DoubledGaugeField &U,
 							  const FermionField &in,
 							  FermionField &out, int dag) 
@@ -385,7 +383,7 @@ void NaiveStaggeredFermion<Impl>::DhopInternalSerialComms(StencilImpl &st, Lebes
   {
     int interior=1;
     int exterior=1;
-    Kernels::DhopNaive(st,lo,U,in,out,dag,interior,exterior);
+    Kernels::DhopNaive(st,U,in,out,dag,interior,exterior);
   }
 };
 

Grid/qcd/action/fermion/implementation/StaggeredKernelsHand.h

@@ -375,23 +375,6 @@ void StaggeredKernels<Impl>::DhopSiteHandExt(StencilView &st,
   }
 }
 
-/*
-#define DHOP_SITE_HAND_INSTANTIATE(IMPL)				\
-  template void StaggeredKernels<IMPL>::DhopSiteHand(StencilImpl &st, LebesgueOrder &lo, \
-						     DoubledGaugeFieldView &U,DoubledGaugeFieldView &UUU, \
-						     SiteSpinor *buf, int LLs, int sU, \
-						     const FermionFieldView &in, FermionFieldView &out, int dag); \
-									\
-  template void StaggeredKernels<IMPL>::DhopSiteHandInt(StencilImpl &st, LebesgueOrder &lo, \
-						     DoubledGaugeFieldView &U,DoubledGaugeFieldView &UUU, \
-						     SiteSpinor *buf, int LLs, int sU, \
-						     const FermionFieldView &in, FermionFieldView &out, int dag); \
-									\
-  template void StaggeredKernels<IMPL>::DhopSiteHandExt(StencilImpl &st, LebesgueOrder &lo, \
-						     DoubledGaugeFieldView &U,DoubledGaugeFieldView &UUU, \
-						     SiteSpinor *buf, int LLs, int sU, \
-						     const FermionFieldView &in, FermionFieldView &out, int dag); \
-*/
 #undef LOAD_CHI
 #undef HAND_DECLARATIONS
 

Grid/qcd/action/fermion/implementation/StaggeredKernelsImplementation.h

@@ -256,7 +256,7 @@ void StaggeredKernels<Impl>::DhopDirKernel(StencilImpl &st, DoubledGaugeFieldVie
   });
 
 template <class Impl> 
-void StaggeredKernels<Impl>::DhopImproved(StencilImpl &st, LebesgueOrder &lo, 
+void StaggeredKernels<Impl>::DhopImproved(StencilImpl &st, 
 					  DoubledGaugeField &U, DoubledGaugeField &UUU, 
 					  const FermionField &in, FermionField &out, int dag, int interior,int exterior)
 {
@@ -294,7 +294,7 @@ void StaggeredKernels<Impl>::DhopImproved(StencilImpl &st, LebesgueOrder &lo,
   assert(0 && " Kernel optimisation case not covered ");
 }
 template <class Impl> 
-void StaggeredKernels<Impl>::DhopNaive(StencilImpl &st, LebesgueOrder &lo, 
+void StaggeredKernels<Impl>::DhopNaive(StencilImpl &st, 
 				       DoubledGaugeField &U,
 				       const FermionField &in, FermionField &out, int dag, int interior,int exterior)
 {

Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h

@@ -58,15 +58,9 @@ WilsonFermion5D<Impl>::WilsonFermion5D(GaugeField &_Umu,
   Umu(_FourDimGrid),
   UmuEven(_FourDimRedBlackGrid),
   UmuOdd (_FourDimRedBlackGrid),
-  Lebesgue(_FourDimGrid),
-  LebesgueEvenOdd(_FourDimRedBlackGrid),
   _tmp(&FiveDimRedBlackGrid),
   Dirichlet(0)
 {
-  Stencil.lo     = &Lebesgue;
-  StencilEven.lo = &LebesgueEvenOdd;
-  StencilOdd.lo  = &LebesgueEvenOdd;
-  
   // some assertions
   assert(FiveDimGrid._ndimension==5);
   assert(FourDimGrid._ndimension==4);
@@ -305,19 +299,19 @@ void WilsonFermion5D<Impl>::DhopDerivOE(GaugeField &mat,
 }
 
 template<class Impl>
-void WilsonFermion5D<Impl>::DhopInternal(StencilImpl & st, LebesgueOrder &lo,
+void WilsonFermion5D<Impl>::DhopInternal(StencilImpl & st,
                                          DoubledGaugeField & U,
                                          const FermionField &in, FermionField &out,int dag)
 {
   if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsAndCompute )
-    DhopInternalOverlappedComms(st,lo,U,in,out,dag);
+    DhopInternalOverlappedComms(st,U,in,out,dag);
   else 
-    DhopInternalSerialComms(st,lo,U,in,out,dag);
+    DhopInternalSerialComms(st,U,in,out,dag);
 }
 
 
 template<class Impl>
-void WilsonFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl & st, LebesgueOrder &lo,
+void WilsonFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl & st,
 							DoubledGaugeField & U,
 							const FermionField &in, FermionField &out,int dag)
 {
@@ -331,10 +325,12 @@ void WilsonFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl & st, Lebesg
   // Start comms  // Gather intranode and extra node differentiated??
   /////////////////////////////
   {
+    std::cout << " WilsonFermion5D gather " <<std::endl;
     GRID_TRACE("Gather");
     st.HaloExchangeOptGather(in,compressor); // Put the barrier in the routine
   }
   
+  std::cout << " WilsonFermion5D Communicate Begin " <<std::endl;
   std::vector<std::vector<CommsRequest_t> > requests;
   auto id=traceStart("Communicate overlapped");
   st.CommunicateBegin(requests);
@@ -343,6 +339,7 @@ void WilsonFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl & st, Lebesg
   // Overlap with comms
   /////////////////////////////
   {
+  std::cout << " WilsonFermion5D Comms merge " <<std::endl;
     GRID_TRACE("MergeSHM");
     st.CommsMergeSHM(compressor);// Could do this inside parallel region overlapped with comms
   }
@@ -350,6 +347,7 @@ void WilsonFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl & st, Lebesg
   /////////////////////////////
   // do the compute interior
   /////////////////////////////
+  std::cout << " WilsonFermion5D Interior " <<std::endl;
   int Opt = WilsonKernelsStatic::Opt; // Why pass this. Kernels should know
   if (dag == DaggerYes) {
     GRID_TRACE("DhopDagInterior");
@@ -362,6 +360,7 @@ void WilsonFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl & st, Lebesg
   /////////////////////////////
   // Complete comms
   /////////////////////////////
+  std::cout << " WilsonFermion5D Comms Complete " <<std::endl;
   st.CommunicateComplete(requests);
   traceStop(id);
 
@@ -369,11 +368,13 @@ void WilsonFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl & st, Lebesg
   // do the compute exterior
   /////////////////////////////
   {
+    std::cout << " WilsonFermion5D Comms Merge " <<std::endl;
     GRID_TRACE("Merge");
     st.CommsMerge(compressor);
   }
   
 
+  std::cout << " WilsonFermion5D Exterior " <<std::endl;
   if (dag == DaggerYes) {
     GRID_TRACE("DhopDagExterior");
     Kernels::DhopDagKernel(Opt,st,U,st.CommBuf(),LLs,U.oSites(),in,out,0,1);
@@ -381,11 +382,12 @@ void WilsonFermion5D<Impl>::DhopInternalOverlappedComms(StencilImpl & st, Lebesg
     GRID_TRACE("DhopExterior");
     Kernels::DhopKernel   (Opt,st,U,st.CommBuf(),LLs,U.oSites(),in,out,0,1);
   }
+  std::cout << " WilsonFermion5D Done " <<std::endl;
 }
 
 
 template<class Impl>
-void WilsonFermion5D<Impl>::DhopInternalSerialComms(StencilImpl & st, LebesgueOrder &lo,
+void WilsonFermion5D<Impl>::DhopInternalSerialComms(StencilImpl & st, 
 						    DoubledGaugeField & U,
 						    const FermionField &in, 
 						    FermionField &out,int dag)
@@ -395,11 +397,13 @@ void WilsonFermion5D<Impl>::DhopInternalSerialComms(StencilImpl & st, LebesgueOr
 
   int LLs = in.Grid()->_rdimensions[0];
 
+  std::cout << " WilsonFermion5D Halo exch " <<std::endl;
   {
     GRID_TRACE("HaloExchange");
     st.HaloExchangeOpt(in,compressor);
   }
   
+  std::cout << " WilsonFermion5D Dhop " <<std::endl;
   int Opt = WilsonKernelsStatic::Opt;
   if (dag == DaggerYes) {
     GRID_TRACE("DhopDag");
@@ -408,6 +412,7 @@ void WilsonFermion5D<Impl>::DhopInternalSerialComms(StencilImpl & st, LebesgueOr
     GRID_TRACE("Dhop");
     Kernels::DhopKernel(Opt,st,U,st.CommBuf(),LLs,U.oSites(),in,out);
   }
+  std::cout << " WilsonFermion5D Done " <<std::endl;
 }
 
 
@@ -420,7 +425,7 @@ void WilsonFermion5D<Impl>::DhopOE(const FermionField &in, FermionField &out,int
   assert(in.Checkerboard()==Even);
   out.Checkerboard() = Odd;
 
-  DhopInternal(StencilEven,LebesgueEvenOdd,UmuOdd,in,out,dag);
+  DhopInternal(StencilEven,UmuOdd,in,out,dag);
 }
 template<class Impl>
 void WilsonFermion5D<Impl>::DhopEO(const FermionField &in, FermionField &out,int dag)
@@ -431,7 +436,7 @@ void WilsonFermion5D<Impl>::DhopEO(const FermionField &in, FermionField &out,int
   assert(in.Checkerboard()==Odd);
   out.Checkerboard() = Even;
 
-  DhopInternal(StencilOdd,LebesgueEvenOdd,UmuEven,in,out,dag);
+  DhopInternal(StencilOdd,UmuEven,in,out,dag);
 }
 template<class Impl>
 void WilsonFermion5D<Impl>::Dhop(const FermionField &in, FermionField &out,int dag)
@@ -441,7 +446,7 @@ void WilsonFermion5D<Impl>::Dhop(const FermionField &in, FermionField &out,int d
 
   out.Checkerboard() = in.Checkerboard();
 
-  DhopInternal(Stencil,Lebesgue,Umu,in,out,dag);
+  DhopInternal(Stencil,Umu,in,out,dag);
 }
 template<class Impl>
 void WilsonFermion5D<Impl>::DW(const FermionField &in, FermionField &out,int dag)

Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h

@@ -52,17 +52,12 @@ WilsonFermion<Impl>::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
     StencilEven(&Hgrid, npoint, Even, directions,displacements,p),  // source is Even
     StencilOdd(&Hgrid, npoint, Odd, directions,displacements,p),  // source is Odd
     mass(_mass),
-    Lebesgue(_grid),
-    LebesgueEvenOdd(_cbgrid),
     Umu(&Fgrid),
     UmuEven(&Hgrid),
     UmuOdd(&Hgrid),
       _tmp(&Hgrid),
       anisotropyCoeff(anis)
 {
-  Stencil.lo     = &Lebesgue;
-  StencilEven.lo = &LebesgueEvenOdd;
-  StencilOdd.lo  = &LebesgueEvenOdd;
   // Allocate the required comms buffer
   ImportGauge(_Umu);
   if  (anisotropyCoeff.isAnisotropic){
@@ -314,7 +309,7 @@ void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out, int da
 
   out.Checkerboard() = in.Checkerboard();
 
-  DhopInternal(Stencil, Lebesgue, Umu, in, out, dag);
+  DhopInternal(Stencil, Umu, in, out, dag);
 }
 
 template <class Impl>
@@ -326,7 +321,7 @@ void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out, int
   assert(in.Checkerboard() == Even);
   out.Checkerboard() = Odd;
 
-  DhopInternal(StencilEven, LebesgueEvenOdd, UmuOdd, in, out, dag);
+  DhopInternal(StencilEven, UmuOdd, in, out, dag);
 }
 
 template <class Impl>
@@ -338,7 +333,7 @@ void WilsonFermion<Impl>::DhopEO(const FermionField &in, FermionField &out,int d
   assert(in.Checkerboard() == Odd);
   out.Checkerboard() = Even;
 
-  DhopInternal(StencilOdd, LebesgueEvenOdd, UmuEven, in, out, dag);
+  DhopInternal(StencilOdd, UmuEven, in, out, dag);
 }
 
 template <class Impl>
@@ -391,21 +386,21 @@ void WilsonFermion<Impl>::DhopDirCalc(const FermionField &in, FermionField &out,
 };
 
 template <class Impl>
-void WilsonFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder &lo,
+void WilsonFermion<Impl>::DhopInternal(StencilImpl &st, 
                                        DoubledGaugeField &U,
                                        const FermionField &in,
                                        FermionField &out, int dag)
 {
 #ifdef GRID_OMP
   if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsAndCompute )
-    DhopInternalOverlappedComms(st,lo,U,in,out,dag);
+    DhopInternalOverlappedComms(st,U,in,out,dag);
   else
 #endif
-    DhopInternalSerial(st,lo,U,in,out,dag);
+    DhopInternalSerial(st,U,in,out,dag);
 }
 
 template <class Impl>
-void WilsonFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, LebesgueOrder &lo,
+void WilsonFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, 
 						      DoubledGaugeField &U,
 						      const FermionField &in,
 						      FermionField &out, int dag)
@@ -474,10 +469,10 @@ void WilsonFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, LebesgueO
 
 
 template <class Impl>
-void WilsonFermion<Impl>::DhopInternalSerial(StencilImpl &st, LebesgueOrder &lo,
-                                       DoubledGaugeField &U,
-                                       const FermionField &in,
-                                       FermionField &out, int dag)
+void WilsonFermion<Impl>::DhopInternalSerial(StencilImpl &st, 
+					     DoubledGaugeField &U,
+					     const FermionField &in,
+					     FermionField &out, int dag)
 {
   GRID_TRACE("DhopSerial");
   assert((dag == DaggerNo) || (dag == DaggerYes));

Grid/qcd/action/fermion/implementation/WilsonKernelsAsmAvx512.h

@@ -40,11 +40,11 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 /// Switch off the 5d vectorised code optimisations
 #undef DWFVEC5D
 
-static Vector<vComplexF> signsF;
+static std::vector<vComplexF> signsF;
 
   template<typename vtype>    
-  int setupSigns(Vector<vtype>& signs ){
-    Vector<vtype> bother(2);
+  int setupSigns(std::vector<vtype>& signs ){
+    std::vector<vtype> bother(2);
     signs = bother;
     vrsign(signs[0]);
     visign(signs[1]);
@@ -364,7 +364,7 @@ WilsonKernels<ZDomainWallVec5dImplF>::AsmDhopSiteDagExt(StencilView &st, Doubled
 
 #include <simd/Intel512double.h>
     
-static Vector<vComplexD> signsD;
+static std::vector<vComplexD> signsD;
 static int signInitD = setupSigns(signsD);
     
 #define MAYBEPERM(A,perm) if (perm) { A ; }

Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h

@@ -434,7 +434,7 @@ void WilsonKernels<Impl>::DhopDirKernel( StencilImpl &st, DoubledGaugeField &U,S
 
 #define ASM_CALL(A)							\
   thread_for( sss, Nsite, {						\
-    int ss = st.lo->Reorder(sss);					\
+    int ss = sss; /*st.lo->Reorder(sss);*/			\
     int sU = ss;							\
     int sF = ss*Ls;							\
     WilsonKernels<Impl>::A(st_v,U_v,buf,sF,sU,Ls,1,in_v,out_v);		\