Updates in tests to make all of Grid compile

Grid/algorithms/iterative/BlockConjugateGradient.h

@@ -157,7 +157,7 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
   Nblock=8;
   std::cout<<GridLogMessage<<" Block Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl;
 
-  X.checkerboard = B.checkerboard;
+  X.Checkerboard() = B.Checkerboard();
   conformable(X, B);
 
   Field tmp(B);
@@ -336,7 +336,7 @@ void CGmultiRHSsolve(LinearOperatorBase<Field> &Linop, const Field &Src, Field &
 
   std::cout<<GridLogMessage<<"MultiRHS Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl;
 
-  Psi.checkerboard = Src.checkerboard;
+  Psi.Checkerboard() = Src.Checkerboard();
   conformable(Psi, Src);
 
   Field P(Src);
@@ -515,7 +515,7 @@ void BlockCGrQsolveVec(LinearOperatorBase<Field> &Linop, const std::vector<Field
   std::cout<<GridLogMessage<<" Block Conjugate Gradient Vec rQ : Nblock "<<Nblock<<std::endl;
 
   for(int b=0;b<Nblock;b++){ 
-    X[b].checkerboard = B[b].checkerboard;
+    X[b].Checkerboard() = B[b].Checkerboard();
     conformable(X[b], B[b]);
     conformable(X[b], X[0]); 
   }

Grid/lattice/Lattice_comparison_utils.h

@@ -1,4 +1,4 @@
-/*************************************************************************************
+    /*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -24,197 +24,198 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     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 */
+    /*  END LEGAL */
-#ifndef GRID_COMPARISON_H
+
-#define GRID_COMPARISON_H
+#pragma once
 
 NAMESPACE_BEGIN(Grid);
 
-/////////////////////////////////////////
+  /////////////////////////////////////////
-// This implementation is a bit poor.
+  // This implementation is a bit poor.
-//
+  //
-// Only support relational logical operations (<, >  etc)
+  // Only support relational logical operations (<, >  etc)
-// on scalar objects. Therefore can strip any tensor structures.
+  // on scalar objects. Therefore can strip any tensor structures.
-//
+  //
-// Should guard this with isGridTensor<> enable if?
+  // Should guard this with isGridTensor<> enable if?
-/////////////////////////////////////////
+  /////////////////////////////////////////
-//
+  //
-// Generic list of functors
+  // Generic list of functors
-//
+  //
-template<class lobj,class robj> class veq {
+  template<class lobj,class robj> class veq {
-public:
+  public:
-  accelerator vInteger operator()(const lobj &lhs, const robj &rhs)
+    vInteger operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) == (rhs);
+      return (lhs) == (rhs);
-  }
+    }
-};
+  };
-template<class lobj,class robj> class vne {
+  template<class lobj,class robj> class vne {
-public:
+  public:
-  accelerator vInteger operator()(const lobj &lhs, const robj &rhs)
+    vInteger operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) != (rhs);
+      return (lhs) != (rhs);
-  }
+    }
-};
+  };
-template<class lobj,class robj> class vlt {
+  template<class lobj,class robj> class vlt {
-public:
+  public:
-  accelerator vInteger operator()(const lobj &lhs, const robj &rhs)
+    vInteger operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) < (rhs);
+      return (lhs) < (rhs);
-  }
+    }
-};
+  };
-template<class lobj,class robj> class vle {
+  template<class lobj,class robj> class vle {
-public:
+  public:
-  accelerator vInteger operator()(const lobj &lhs, const robj &rhs)
+    vInteger operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) <= (rhs);
+      return (lhs) <= (rhs);
-  }
+    }
-};
+  };
-template<class lobj,class robj> class vgt {
+  template<class lobj,class robj> class vgt {
-public:
+  public:
-  accelerator vInteger operator()(const lobj &lhs, const robj &rhs)
+    vInteger operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) > (rhs);
+      return (lhs) > (rhs);
-  }
+    }
-};
+  };
-template<class lobj,class robj> class vge {
+  template<class lobj,class robj> class vge {
-public:
+    public:
-  accelerator vInteger operator()(const lobj &lhs, const robj &rhs)
+    vInteger operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) >= (rhs);
+      return (lhs) >= (rhs);
-  }
+    }
-};
+  };
   
-// Generic list of functors
+  // Generic list of functors
-template<class lobj,class robj> class seq {
+  template<class lobj,class robj> class seq {
-public:
+  public:
-  accelerator Integer operator()(const lobj &lhs, const robj &rhs)
+    Integer operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) == (rhs);
+      return (lhs) == (rhs);
-  }
+    }
-};
+  };
-template<class lobj,class robj> class sne {
+  template<class lobj,class robj> class sne {
-public:
+  public:
-  accelerator Integer operator()(const lobj &lhs, const robj &rhs)
+    Integer operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) != (rhs);
+      return (lhs) != (rhs);
-  }
+    }
-};
+  };
-template<class lobj,class robj> class slt {
+  template<class lobj,class robj> class slt {
-public:
+  public:
-  accelerator Integer operator()(const lobj &lhs, const robj &rhs)
+    Integer operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) < (rhs);
+      return (lhs) < (rhs);
-  }
+    }
-};
+  };
-template<class lobj,class robj> class sle {
+  template<class lobj,class robj> class sle {
-public:
+  public:
-  accelerator Integer operator()(const lobj &lhs, const robj &rhs)
+    Integer operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) <= (rhs);
+      return (lhs) <= (rhs);
-  }
+    }
-};
+  };
-template<class lobj,class robj> class sgt {
+  template<class lobj,class robj> class sgt {
-public:
+  public:
-  accelerator Integer operator()(const lobj &lhs, const robj &rhs)
+    Integer operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) > (rhs);
+      return (lhs) > (rhs);
-  }
+    }
-};
+  };
-template<class lobj,class robj> class sge {
+  template<class lobj,class robj> class sge {
-public:
+  public:
-  accelerator Integer operator()(const lobj &lhs, const robj &rhs)
+    Integer operator()(const lobj &lhs, const robj &rhs)
-  { 
+    { 
-    return (lhs) >= (rhs);
+      return (lhs) >= (rhs);
-  }
+    }
-};
+  };
 
-//////////////////////////////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////////////////////////////
-// Integer and real get extra relational functions.
+  // Integer and real get extra relational functions.
-//////////////////////////////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////////////////////////////////////////////////////////////////////////
-template<class sfunctor, class vsimd,IfNotComplex<vsimd> = 0> 
+  template<class sfunctor, class vsimd,IfNotComplex<vsimd> = 0> 
-accelerator_inline vInteger Comparison(sfunctor sop,const vsimd & lhs, const vsimd & rhs)
+    inline vInteger Comparison(sfunctor sop,const vsimd & lhs, const vsimd & rhs)
-{
+    {
-  typedef typename vsimd::scalar_type scalar;
+      typedef typename vsimd::scalar_type scalar;
-  ExtractBuffer<scalar> vlhs(vsimd::Nsimd());   // Use functors to reduce this to single implementation
+      ExtractBuffer<scalar> vlhs(vsimd::Nsimd());   // Use functors to reduce this to single implementation
-  ExtractBuffer<scalar> vrhs(vsimd::Nsimd());
+      ExtractBuffer<scalar> vrhs(vsimd::Nsimd());
-  ExtractBuffer<Integer> vpred(vsimd::Nsimd());
+      ExtractBuffer<Integer> vpred(vsimd::Nsimd());
-  vInteger ret;
+      vInteger ret;
-  extract<vsimd,scalar>(lhs,vlhs);
+      extract<vsimd,scalar>(lhs,vlhs);
-  extract<vsimd,scalar>(rhs,vrhs);
+      extract<vsimd,scalar>(rhs,vrhs);
-  for(int s=0;s<vsimd::Nsimd();s++){
+      for(int s=0;s<vsimd::Nsimd();s++){
-    vpred[s] = sop(vlhs[s],vrhs[s]);
+	vpred[s] = sop(vlhs[s],vrhs[s]);
-  }
+      }
-  merge<vInteger,Integer>(ret,vpred);
+      merge<vInteger,Integer>(ret,vpred);
-  return ret;
+      return ret;
-}
+    }
 
-template<class sfunctor, class vsimd,IfNotComplex<vsimd> = 0> 
+  template<class sfunctor, class vsimd,IfNotComplex<vsimd> = 0> 
-accelerator_inline vInteger Comparison(sfunctor sop,const vsimd & lhs, const typename vsimd::scalar_type & rhs)
+    inline vInteger Comparison(sfunctor sop,const vsimd & lhs, const typename vsimd::scalar_type & rhs)
-{
+    {
-  typedef typename vsimd::scalar_type scalar;
+      typedef typename vsimd::scalar_type scalar;
-  ExtractBuffer<scalar>  vlhs(vsimd::Nsimd());   // Use functors to reduce this to single implementation
+      ExtractBuffer<scalar> vlhs(vsimd::Nsimd());   // Use functors to reduce this to single implementation
-  ExtractBuffer<Integer> vpred(vsimd::Nsimd());
+      ExtractBuffer<Integer> vpred(vsimd::Nsimd());
-  vInteger ret;
+      vInteger ret;
-  extract<vsimd,scalar>(lhs,vlhs);
+      extract<vsimd,scalar>(lhs,vlhs);
-  for(int s=0;s<vsimd::Nsimd();s++){
+      for(int s=0;s<vsimd::Nsimd();s++){
-    vpred[s] = sop(vlhs[s],rhs);
+	vpred[s] = sop(vlhs[s],rhs);
-  }
+      }
-  merge<vInteger,Integer>(ret,vpred);
+      merge<vInteger,Integer>(ret,vpred);
-  return ret;
+      return ret;
-}
+    }
 
-template<class sfunctor, class vsimd,IfNotComplex<vsimd> = 0> 
+  template<class sfunctor, class vsimd,IfNotComplex<vsimd> = 0> 
-accelerator_inline vInteger Comparison(sfunctor sop,const typename vsimd::scalar_type & lhs, const vsimd & rhs)
+    inline vInteger Comparison(sfunctor sop,const typename vsimd::scalar_type & lhs, const vsimd & rhs)
-{
+    {
-  typedef typename vsimd::scalar_type scalar;
+      typedef typename vsimd::scalar_type scalar;
-  ExtractBuffer<scalar> vrhs(vsimd::Nsimd());   // Use functors to reduce this to single implementation
+      ExtractBuffer<scalar> vrhs(vsimd::Nsimd());   // Use functors to reduce this to single implementation
-  ExtractBuffer<Integer> vpred(vsimd::Nsimd());
+      ExtractBuffer<Integer> vpred(vsimd::Nsimd());
-  vInteger ret;
+      vInteger ret;
-  extract<vsimd,scalar>(rhs,vrhs);
+      extract<vsimd,scalar>(rhs,vrhs);
-  for(int s=0;s<vsimd::Nsimd();s++){
+      for(int s=0;s<vsimd::Nsimd();s++){
-    vpred[s] = sop(lhs,vrhs[s]);
+	vpred[s] = sop(lhs,vrhs[s]);
-  }
+      }
-  merge<vInteger,Integer>(ret,vpred);
+      merge<vInteger,Integer>(ret,vpred);
-  return ret;
+      return ret;
-}
+    }
 
 #define DECLARE_RELATIONAL_EQ(op,functor) \
-  template<class vsimd,IfSimd<vsimd> = 0>				\
+  template<class vsimd,IfSimd<vsimd> = 0>\
-  accelerator_inline vInteger operator op (const vsimd & lhs, const vsimd & rhs)	\
+    inline vInteger operator op (const vsimd & lhs, const vsimd & rhs)\
-  {									\
+    {\
-    typedef typename vsimd::scalar_type scalar;				\
+      typedef typename vsimd::scalar_type scalar;\
-    return Comparison(functor<scalar,scalar>(),lhs,rhs);		\
+      return Comparison(functor<scalar,scalar>(),lhs,rhs);\
-  }									\
+    }\
-  template<class vsimd,IfSimd<vsimd> = 0>				\
+  template<class vsimd,IfSimd<vsimd> = 0>\
-  accelerator_inline vInteger operator op (const vsimd & lhs, const typename vsimd::scalar_type & rhs) \
+    inline vInteger operator op (const vsimd & lhs, const typename vsimd::scalar_type & rhs) \
-  {									\
+    {\
-    typedef typename vsimd::scalar_type scalar;				\
+      typedef typename vsimd::scalar_type scalar;\
-    return Comparison(functor<scalar,scalar>(),lhs,rhs);		\
+      return Comparison(functor<scalar,scalar>(),lhs,rhs);\
-  }									\
+    }\
-  template<class vsimd,IfSimd<vsimd> = 0>				\
+  template<class vsimd,IfSimd<vsimd> = 0>\
-  accelerator_inline vInteger operator op (const typename vsimd::scalar_type & lhs, const vsimd & rhs) \
+    inline vInteger operator op (const typename vsimd::scalar_type & lhs, const vsimd & rhs) \
-  {									\
+    {\
-    typedef typename vsimd::scalar_type scalar;				\
+      typedef typename vsimd::scalar_type scalar;\
-    return Comparison(functor<scalar,scalar>(),lhs,rhs);		\
+      return Comparison(functor<scalar,scalar>(),lhs,rhs);\
-  }									\
+    }\
-  template<class vsimd>							\
+  template<class vsimd>\
-  accelerator_inline vInteger operator op(const iScalar<vsimd> &lhs,const iScalar<vsimd> &rhs) \
+    inline vInteger operator op(const iScalar<vsimd> &lhs,const typename vsimd::scalar_type &rhs) \
-  {									\
+    {									\
-    return lhs._internal op rhs._internal;				\
+      return lhs._internal op rhs;					\
-  }									\
+    }									\
-  template<class vsimd>							\
+  template<class vsimd>\
-  accelerator_inline vInteger operator op(const iScalar<vsimd> &lhs,const typename vsimd::scalar_type &rhs) \
+    inline vInteger operator op(const typename vsimd::scalar_type &lhs,const iScalar<vsimd> &rhs) \
-  {									\
+    {									\
-    return lhs._internal op rhs;					\
+      return lhs op rhs._internal;					\
-  }									\
-  template<class vsimd>							\
-  accelerator_inline vInteger operator op(const typename vsimd::scalar_type &lhs,const iScalar<vsimd> &rhs) \
-  {									\
-    return lhs op rhs._internal;					\
     }									\
 
-#define DECLARE_RELATIONAL(op,functor) DECLARE_RELATIONAL_EQ(op,functor)    
+#define DECLARE_RELATIONAL(op,functor) \
+  DECLARE_RELATIONAL_EQ(op,functor)    \
+  template<class vsimd>\
+    inline vInteger operator op(const iScalar<vsimd> &lhs,const iScalar<vsimd> &rhs)\
+    {									\
+      return lhs._internal op rhs._internal;				\
+    }									
 
 DECLARE_RELATIONAL(<,slt);
 DECLARE_RELATIONAL(<=,sle);
@@ -228,4 +229,4 @@ DECLARE_RELATIONAL(!=,sne);
 NAMESPACE_END(Grid);
 
 
-#endif
+

Grid/lattice/Lattice_reduction.h

@@ -477,7 +477,7 @@ static void sliceNorm (std::vector<RealD> &sn,const Lattice<vobj> &rhs,int Ortho
   typedef typename vobj::vector_type vector_type;
   
   int Nblock = rhs.Grid()->GlobalDimensions()[Orthog];
-  Vector<ComplexD> ip(Nblock);
+  std::vector<ComplexD> ip(Nblock);
   sn.resize(Nblock);
   
   sliceInnerProductVector(ip,rhs,rhs,Orthog);
@@ -586,6 +586,10 @@ static void sliceMaddMatrix (Lattice<vobj> &R,Eigen::MatrixXcd &aa,const Lattice
   int block =FullGrid->_slice_block [Orthog];
   int nblock=FullGrid->_slice_nblock[Orthog];
   int ostride=FullGrid->_ostride[Orthog];
+
+  auto X_v=X.View();
+  auto Y_v=Y.View();
+  auto R_v=R.View();
   thread_region
   {
     Vector<vobj> s_x(Nblock);
@@ -595,16 +599,16 @@ static void sliceMaddMatrix (Lattice<vobj> &R,Eigen::MatrixXcd &aa,const Lattice
       int o  = n*stride + b;
 
       for(int i=0;i<Nblock;i++){
-	s_x[i] = X[o+i*ostride];
+	s_x[i] = X_v[o+i*ostride];
       }
 
       vobj dot;
       for(int i=0;i<Nblock;i++){
-	dot = Y[o+i*ostride];
+	dot = Y_v[o+i*ostride];
 	for(int j=0;j<Nblock;j++){
 	  dot = dot + s_x[j]*(scale*aa(j,i));
 	}
-	R[o+i*ostride]=dot;
+	R_v[o+i*ostride]=dot;
       }
     }});
   }
@@ -635,6 +639,8 @@ static void sliceMulMatrix (Lattice<vobj> &R,Eigen::MatrixXcd &aa,const Lattice<
   int block =FullGrid->_slice_block [Orthog];
   int nblock=FullGrid->_slice_nblock[Orthog];
   int ostride=FullGrid->_ostride[Orthog];
+  auto R_v = R.View();
+  auto X_v = X.View();
   thread_region
   {
     std::vector<vobj> s_x(Nblock);
@@ -645,7 +651,7 @@ static void sliceMulMatrix (Lattice<vobj> &R,Eigen::MatrixXcd &aa,const Lattice<
       int o  = n*stride + b;
 
       for(int i=0;i<Nblock;i++){
-	s_x[i] = X[o+i*ostride];
+	s_x[i] = X_v[o+i*ostride];
       }
 
       vobj dot;
@@ -654,7 +660,7 @@ static void sliceMulMatrix (Lattice<vobj> &R,Eigen::MatrixXcd &aa,const Lattice<
 	for(int j=1;j<Nblock;j++){
 	  dot = dot + s_x[j]*(scale*aa(j,i));
 	}
-	R[o+i*ostride]=dot;
+	R_v[o+i*ostride]=dot;
       }
     }});
   }
@@ -692,6 +698,8 @@ static void sliceInnerProductMatrix(  Eigen::MatrixXcd &mat, const Lattice<vobj>
 
   typedef typename vobj::vector_typeD vector_typeD;
 
+  auto lhs_v=lhs.View();
+  auto rhs_v=rhs.View();
   thread_region
   {
     std::vector<vobj> Left(Nblock);
@@ -704,8 +712,8 @@ static void sliceInnerProductMatrix(  Eigen::MatrixXcd &mat, const Lattice<vobj>
       int o  = n*stride + b;
 
       for(int i=0;i<Nblock;i++){
-	Left [i] = lhs[o+i*ostride];
+	Left [i] = lhs_v[o+i*ostride];
-	Right[i] = rhs[o+i*ostride];
+	Right[i] = rhs_v[o+i*ostride];
       }
 
       for(int i=0;i<Nblock;i++){

Grid/qcd/action/fermion/WilsonTMFermion5D.h

@@ -63,7 +63,7 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
     }
   
   virtual void Meooe(const FermionField &in, FermionField &out) {
-    if (in.checkerboard == Odd) {
+    if (in.Checkerboard() == Odd) {
       this->DhopEO(in, out, DaggerNo);
     } else {
       this->DhopOE(in, out, DaggerNo);
@@ -71,7 +71,7 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
   }
   
   virtual void MeooeDag(const FermionField &in, FermionField &out) {
-    if (in.checkerboard == Odd) {
+    if (in.Checkerboard() == Odd) {
       this->DhopEO(in, out, DaggerYes);
     } else {
       this->DhopOE(in, out, DaggerYes);
@@ -80,7 +80,7 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
   
   // allow override for twisted mass and clover
   virtual void Mooee(const FermionField &in, FermionField &out) {
-    out.checkerboard = in.checkerboard;
+    out.Checkerboard() = in.Checkerboard();
     //axpibg5x(out,in,a,b); // out = a*in + b*i*G5*in
     for (int s=0;s<(int)this->mass.size();s++) {
       ComplexD a = 4.0+this->mass[s];
@@ -90,7 +90,7 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
   }
   
   virtual void MooeeDag(const FermionField &in, FermionField &out) {
-    out.checkerboard = in.checkerboard;
+    out.Checkerboard() = in.Checkerboard();
     for (int s=0;s<(int)this->mass.size();s++) {
       ComplexD a = 4.0+this->mass[s];
       ComplexD b(0.0,-this->mu[s]);
@@ -121,9 +121,9 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
   }
   
   virtual RealD M(const FermionField &in, FermionField &out) {
-    out.checkerboard = in.checkerboard;
+    out.Checkerboard() = in.Checkerboard();
     this->Dhop(in, out, DaggerNo);
-    FermionField tmp(out._grid);
+    FermionField tmp(out.Grid());
     for (int s=0;s<(int)this->mass.size();s++) {
       ComplexD a = 4.0+this->mass[s];
       ComplexD b(0.0,this->mu[s]);

Grid/qcd/action/scalar/ScalarInteractionAction.h

@@ -81,16 +81,20 @@ public:
 
   virtual RealD S(const Field &p)
   {
-    assert(p._grid->Nd() == Ndim);
+    assert(p.Grid()->Nd() == Ndim);
-    static Stencil phiStencil(p._grid, npoint, 0, directions, displacements);
+    static Stencil phiStencil(p.Grid(), npoint, 0, directions, displacements);
     phiStencil.HaloExchange(p, compressor);
-    Field action(p._grid), pshift(p._grid), phisquared(p._grid);
+    Field action(p.Grid()), pshift(p.Grid()), phisquared(p.Grid());
     phisquared = p * p;
     action = (2.0 * Ndim + mass_square) * phisquared - lambda * phisquared * phisquared;
+    
+    
+    auto p_v = p.View();
+    auto action_v = action.View();
     for (int mu = 0; mu < Ndim; mu++)
     {
       //  pshift = Cshift(p, mu, +1);  // not efficient, implement with stencils
-      parallel_for(int i = 0; i < p._grid->oSites(); i++)
+      parallel_for(int i = 0; i < p.Grid()->oSites(); i++)
       {
         int permute_type;
         StencilEntry *SE;
@@ -98,23 +102,20 @@ public:
         const vobj *temp, *t_p;
 
         SE = phiStencil.GetEntry(permute_type, mu, i);
-        t_p = &p._odata[i];
+        t_p = &p_v[i];
         if (SE->_is_local)
         {
-          temp = &p._odata[SE->_offset];
+          temp = &p_v[SE->_offset];
-          if (SE->_permute)
+          if (SE->_permute) {
-          {
             permute(temp2, *temp, permute_type);
-            action._odata[i] -= temp2 * (*t_p) + (*t_p) * temp2;
+            action_v[i] -= temp2 * (*t_p) + (*t_p) * temp2;
-          }
+          } else {
-          else
+            action_v[i] -= (*temp) * (*t_p) + (*t_p) * (*temp);
-          {
-            action._odata[i] -= (*temp) * (*t_p) + (*t_p) * (*temp);
           }
         }
         else
         {
-          action._odata[i] -= phiStencil.CommBuf()[SE->_offset] * (*t_p) + (*t_p) * phiStencil.CommBuf()[SE->_offset];
+          action_v[i] -= phiStencil.CommBuf()[SE->_offset] * (*t_p) + (*t_p) * phiStencil.CommBuf()[SE->_offset];
         }
       }
       //  action -= pshift*p + p*pshift;
@@ -127,12 +128,12 @@ public:
   virtual void deriv(const Field &p, Field &force)
   {
     double t0 = usecond();
-    assert(p._grid->Nd() == Ndim);
+    assert(p.Grid()->Nd() == Ndim);
     force = (2. * Ndim + mass_square) * p - 2. * lambda * p * p * p;
     double interm_t = usecond();
 
     // move this outside
-    static Stencil phiStencil(p._grid, npoint, 0, directions, displacements);
+    static Stencil phiStencil(p.Grid(), npoint, 0, directions, displacements);
 
     phiStencil.HaloExchange(p, compressor);
     double halo_t = usecond();
@@ -145,59 +146,51 @@ public:
     for (int point = 0; point < npoint; point++)
     {
 
-#pragma omp parallel 
+      auto p_v = p.View();
-{
+      auto force_v = force.View();
-        int permute_type;
-        StencilEntry *SE;
-        const vobj *temp;
             
-#pragma omp for schedule(static, chunk)
+      int permute_type;
-      for (int i = 0; i < p._grid->oSites(); i++)
+      StencilEntry *SE;
-      {
+      const vobj *temp;
-        SE = phiStencil.GetEntry(permute_type, point, i);
-        // prefetch next p?
 
-        if (SE->_is_local)
+      parallel_for (int i = 0; i < p.Grid()->oSites(); i++) {
-        {
-          temp = &p._odata[SE->_offset];
 	
-          if (SE->_permute)
+	SE = phiStencil.GetEntry(permute_type, point, i);
-          {
+	// prefetch next p?
+	  
+	if (SE->_is_local) {
+	  temp = &p_v[SE->_offset];
+	    
+          if (SE->_permute) {
             vobj temp2;
             permute(temp2, *temp, permute_type);
-            force._odata[i] -= temp2;
+            force_v[i] -= temp2;
+          } else {
+            force_v[i] -= *temp; // slow part. Dominated by this read/write (BW)
           }
-          else
+        } else {
-          {
+          force_v[i] -= phiStencil.CommBuf()[SE->_offset];
-            force._odata[i] -= *temp; // slow part. Dominated by this read/write (BW)
-          }
-        }
-        else
-        {
-          force._odata[i] -= phiStencil.CommBuf()[SE->_offset];
         }
       }
-
     }
-  }
+    force *= N / g;
-  force *= N / g;
 
-  double t1 = usecond();
+    double t1 = usecond();
-  double total_time = (t1 - t0) / 1e6;
+    double total_time = (t1 - t0) / 1e6;
-  double interm_time = (interm_t - t0) / 1e6;
+    double interm_time = (interm_t - t0) / 1e6;
-  double halo_time = (halo_t - interm_t) / 1e6;
+    double halo_time = (halo_t - interm_t) / 1e6;
-  double stencil_time = (t1 - halo_t) / 1e6;
+    double stencil_time = (t1 - halo_t) / 1e6;
-  std::cout << GridLogIntegrator << "Total time for force computation (s)       : " << total_time << std::endl;
+    std::cout << GridLogIntegrator << "Total time for force computation (s)       : " << total_time << std::endl;
-  std::cout << GridLogIntegrator << "Intermediate time for force computation (s): " << interm_time << std::endl;
+    std::cout << GridLogIntegrator << "Intermediate time for force computation (s): " << interm_time << std::endl;
-  std::cout << GridLogIntegrator << "Halo time in force computation (s)         : " << halo_time << std::endl;
+    std::cout << GridLogIntegrator << "Halo time in force computation (s)         : " << halo_time << std::endl;
-  std::cout << GridLogIntegrator << "Stencil time in force computation (s)      : " << stencil_time << std::endl;
+    std::cout << GridLogIntegrator << "Stencil time in force computation (s)      : " << stencil_time << std::endl;
-  double flops = p._grid->gSites() * (14 * N * N * N + 18 * N * N + 2);
+    double flops = p.Grid()->gSites() * (14 * N * N * N + 18 * N * N + 2);
-  double flops_no_stencil = p._grid->gSites() * (14 * N * N * N + 6 * N * N + 2);
+    double flops_no_stencil = p.Grid()->gSites() * (14 * N * N * N + 6 * N * N + 2);
-  double Gflops = flops / (total_time * 1e9);
+    double Gflops = flops / (total_time * 1e9);
-  double Gflops_no_stencil = flops_no_stencil / (interm_time * 1e9);
+    double Gflops_no_stencil = flops_no_stencil / (interm_time * 1e9);
-  std::cout << GridLogIntegrator << "Flops: " << flops << "  - Gflop/s : " << Gflops << std::endl;
+    std::cout << GridLogIntegrator << "Flops: " << flops << "  - Gflop/s : " << Gflops << std::endl;
-  std::cout << GridLogIntegrator << "Flops NS: " << flops_no_stencil << "  - Gflop/s NS: " << Gflops_no_stencil << std::endl;
+    std::cout << GridLogIntegrator << "Flops NS: " << flops_no_stencil << "  - Gflop/s NS: " << Gflops_no_stencil << std::endl;
-}
+  }
 };
 
 NAMESPACE_END(Grid);

Grid/qcd/hmc/checkpointers/ILDGCheckpointer.h

@@ -73,7 +73,7 @@ public:
     if ((traj % Params.saveInterval) == 0) {
       std::string config, rng;
       this->build_filenames(traj, Params, config, rng);
-      GridBase *grid = U._grid;
+      GridBase *grid = U.Grid();
       uint32_t nersc_csum,scidac_csuma,scidac_csumb;
       BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
       IldgWriter _IldgWriter(grid->IsBoss());

Grid/qcd/hmc/checkpointers/ScidacCheckpointer.h

@@ -75,7 +75,7 @@ class ScidacHmcCheckpointer : public BaseHmcCheckpointer<Implementation> {
     if ((traj % Params.saveInterval) == 0) {
       std::string config, rng;
       this->build_filenames(traj, Params, config, rng);
-      GridBase *grid = U._grid;
+      GridBase *grid = U.Grid();
       uint32_t nersc_csum,scidac_csuma,scidac_csumb;
       BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
       ScidacWriter _ScidacWriter(grid->IsBoss());

Grid/qcd/utils/WilsonLoops.h

@@ -128,12 +128,12 @@ public:
   // average over all x,y,z the temporal loop
   //////////////////////////////////////////////////
   static ComplexD avgPolyakovLoop(const GaugeField &Umu) {  //assume Nd=4
-    GaugeMat Ut(Umu._grid), P(Umu._grid);
+    GaugeMat Ut(Umu.Grid()), P(Umu.Grid());
     ComplexD out;
-    int T = Umu._grid->GlobalDimensions()[3];
+    int T = Umu.Grid()->GlobalDimensions()[3];
-    int X = Umu._grid->GlobalDimensions()[0];
+    int X = Umu.Grid()->GlobalDimensions()[0];
-    int Y = Umu._grid->GlobalDimensions()[1];
+    int Y = Umu.Grid()->GlobalDimensions()[1];
-    int Z = Umu._grid->GlobalDimensions()[2];
+    int Z = Umu.Grid()->GlobalDimensions()[2];
 
     Ut = peekLorentz(Umu,3); //Select temporal direction
     P = Ut;

Grid/tensors/Tensor_logical.h

@@ -55,13 +55,13 @@ LOGICAL_BINOP(||);
 LOGICAL_BINOP(&&);
 
 template <class T>
-strong_inline bool operator==(const iScalar<T> &t1, const iScalar<T> &t2)
+accelerator_inline bool operator==(const iScalar<T> &t1, const iScalar<T> &t2)
 {
   return (t1._internal == t2._internal);
 }
 
 template <class T, int N>
-strong_inline bool operator==(const iVector<T, N> &t1, const iVector<T, N> &t2)
+accelerator_inline bool operator==(const iVector<T, N> &t1, const iVector<T, N> &t2)
 {
   bool res = true;
 
@@ -74,7 +74,7 @@ strong_inline bool operator==(const iVector<T, N> &t1, const iVector<T, N> &t2)
 }
 
 template <class T, int N>
-strong_inline bool operator==(const iMatrix<T, N> &t1, const iMatrix<T, N> &t2)
+accelerator_inline bool operator==(const iMatrix<T, N> &t1, const iMatrix<T, N> &t2)
 {
   bool res = true;
 

tests/IO/Test_serialisation.cc

@@ -45,8 +45,8 @@ public:
                           bool , b,
                           std::vector<double>, array,
                           std::vector<std::vector<double> >, twodimarray,
-                          std::vector<std::vector<std::vector<Complex> > >, cmplx3darray,
+			  std::vector<std::vector<std::vector<Complex> > >, cmplx3darray,
-                          SpinColourMatrix, scm
+			  SpinColourMatrix, scm
                           );
   myclass() {}
   myclass(int i)

tests/core/Test_staggered5DvecF.cc

@@ -30,15 +30,14 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 using namespace std;
 using namespace Grid;
-using namespace Grid::QCD;
 
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
+  auto latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  auto mpi_layout  = GridDefaultMpi();
 
   const int Ls=16;
   GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
@@ -76,9 +75,9 @@ int main (int argc, char ** argv)
   src = zero;
   pokeSite(cv,src,site);
   */
-  FermionField result(FGrid); result=zero;
+  FermionField result(FGrid); result=Zero();
-  FermionField    tmp(FGrid);    tmp=zero;
+  FermionField    tmp(FGrid);    tmp=Zero();
-  FermionField    err(FGrid);    tmp=zero;
+  FermionField    err(FGrid);    tmp=Zero();
   FermionField phi   (FGrid); random(pRNG5,phi);
   FermionField chi   (FGrid); random(pRNG5,chi);
 
@@ -88,7 +87,7 @@ int main (int argc, char ** argv)
   /*
   for(int mu=1;mu<4;mu++){
     auto tmp = PeekIndex<LorentzIndex>(Umu,mu);
-        tmp = zero;
+        tmp = Zero();
     PokeIndex<LorentzIndex>(Umu,tmp,mu);
   }
   */
@@ -129,9 +128,9 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "Calling vectorised staggered operator"<<std::endl;
 
 #ifdef AVX512
-  QCD::StaggeredKernelsStatic::Opt=QCD::StaggeredKernelsStatic::OptInlineAsm;
+  StaggeredKernelsStatic::Opt=StaggeredKernelsStatic::OptInlineAsm;
 #else
-  QCD::StaggeredKernelsStatic::Opt=QCD::StaggeredKernelsStatic::OptGeneric;
+  StaggeredKernelsStatic::Opt=StaggeredKernelsStatic::OptGeneric;
 #endif
 
   t0=usecond();
@@ -150,9 +149,9 @@ int main (int argc, char ** argv)
 
   
   FermionField ssrc  (sFGrid);  localConvert(src,ssrc);
-  FermionField sresult(sFGrid); sresult=zero;
+  FermionField sresult(sFGrid); sresult=Zero();
 
-  QCD::StaggeredKernelsStatic::Opt=QCD::StaggeredKernelsStatic::OptHandUnroll;
+  StaggeredKernelsStatic::Opt=StaggeredKernelsStatic::OptHandUnroll;
   t0=usecond();
   for(int i=0;i<ncall1;i++){
     sDs.Dhop(ssrc,sresult,0);
@@ -167,9 +166,9 @@ int main (int argc, char ** argv)
 
 
 #ifdef AVX512
-  QCD::StaggeredKernelsStatic::Opt=QCD::StaggeredKernelsStatic::OptInlineAsm;
+  StaggeredKernelsStatic::Opt=StaggeredKernelsStatic::OptInlineAsm;
 #else
-  QCD::StaggeredKernelsStatic::Opt=QCD::StaggeredKernelsStatic::OptGeneric;
+  StaggeredKernelsStatic::Opt=StaggeredKernelsStatic::OptGeneric;
 #endif
 
   err = tmp-result; 

tests/core/Test_wilson_clover.cc

@@ -35,9 +35,9 @@ int main(int argc, char **argv)
 {
   Grid_init(&argc, &argv);
 
-  std::vector<int> latt_size = GridDefaultLatt();
+  auto latt_size = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
-  std::vector<int> mpi_layout = GridDefaultMpi();
+  auto mpi_layout = GridDefaultMpi();
   GridCartesian Grid(latt_size, simd_layout, mpi_layout);
   GridRedBlackCartesian RBGrid(&Grid);
 
@@ -59,17 +59,17 @@ int main(int argc, char **argv)
   FermionField src(&Grid);
   random(pRNG, src);
   FermionField result(&Grid);
-  result = zero;
+  result = Zero();
   FermionField result2(&Grid);
-  result2 = zero;
+  result2 = Zero();
   FermionField ref(&Grid);
-  ref = zero;
+  ref = Zero();
   FermionField tmp(&Grid);
-  tmp = zero;
+  tmp = Zero();
   FermionField err(&Grid);
-  err = zero;
+  err = Zero();
   FermionField err2(&Grid);
-  err2 = zero;
+  err2 = Zero();
   FermionField phi(&Grid);
   random(pRNG, phi);
   FermionField chi(&Grid);
@@ -214,9 +214,9 @@ int main(int argc, char **argv)
   std::cout << GridLogMessage << "= Testing gauge covariance Clover term with EO preconditioning  " << std::endl;
   std::cout << GridLogMessage << "================================================================" << std::endl;
 
-  chi = zero;
+  chi = Zero();
-  phi = zero;
+  phi = Zero();
-  tmp = zero;
+  tmp = Zero();
   pickCheckerboard(Even, chi_e, chi);
   pickCheckerboard(Odd, chi_o, chi);
   pickCheckerboard(Even, phi_e, phi);
@@ -236,9 +236,9 @@ int main(int argc, char **argv)
   LatticeColourMatrix Omega(&Grid);
   LatticeColourMatrix ShiftedOmega(&Grid);
   LatticeGaugeField U_prime(&Grid);
-  U_prime = zero;
+  U_prime = Zero();
   LatticeColourMatrix U_prime_mu(&Grid);
-  U_prime_mu = zero;
+  U_prime_mu = Zero();
   SU<Nc>::LieRandomize(pRNG2, Omega, 1.0);
   for (int mu = 0; mu < Nd; mu++)
   {
@@ -269,8 +269,8 @@ int main(int argc, char **argv)
   std::cout << GridLogMessage << "= Testing gauge covariance Clover term w/o EO preconditioning  " << std::endl;
   std::cout << GridLogMessage << "================================================================" << std::endl;
 
-  chi = zero;
+  chi = Zero();
-  phi = zero;
+  phi = Zero();
 
   WilsonFermionR Dw(Umu, Grid, RBGrid, mass, params);
   Dw.ImportGauge(Umu);
@@ -293,9 +293,9 @@ int main(int argc, char **argv)
   std::cout << GridLogMessage << "= Testing Mooee(csw=0) Clover to reproduce Mooee Wilson   " << std::endl;
   std::cout << GridLogMessage << "==========================================================" << std::endl;
 
-  chi = zero;
+  chi = Zero();
-  phi = zero;
+  phi = Zero();
-  err = zero;
+  err = Zero();
   WilsonCloverFermionR Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); //  <-- Notice: csw=0
   Dwc_csw0.ImportGauge(Umu);
 
@@ -323,9 +323,9 @@ int main(int argc, char **argv)
   std::cout << GridLogMessage << "= Testing EO operator is equal to the unprec              " << std::endl;
   std::cout << GridLogMessage << "==========================================================" << std::endl;
 
-  chi = zero;
+  chi = Zero();
-  phi = zero;
+  phi = Zero();
-  err = zero;
+  err = Zero();
 
   pickCheckerboard(Even, phi_e, phi);
   pickCheckerboard(Odd, phi_o, phi);

tests/debug/Test_cayley_cg.cc

@@ -30,19 +30,9 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 using namespace std;
 using namespace Grid;
- ;
 
-template<class d>
-struct scal {
-  d internal;
-};
 
-  Gamma::Algebra Gmu [] = {
+
-    Gamma::Algebra::GammaX,
-    Gamma::Algebra::GammaY,
-    Gamma::Algebra::GammaZ,
-    Gamma::Algebra::GammaT
-  };
 
 template<class What> 
 void  TestCGinversions(What & Ddwf, 
@@ -261,7 +251,7 @@ void  TestReconstruct5D(What & Ddwf,
 			GridParallelRNG *RNG5)
 {
   LatticeFermion src4   (UGrid); random(*RNG4,src4);
-  LatticeFermion res4   (UGrid); res4 = zero;
+  LatticeFermion res4   (UGrid); res4 = Zero();
 
   LatticeFermion src   (FGrid);
   LatticeFermion src_NE(FGrid);
@@ -332,7 +322,7 @@ void  TestReconstruct5DFA(What & Ddwf,
 			  GridParallelRNG *RNG5)
 {
   LatticeFermion src4   (UGrid); random(*RNG4,src4);
-  LatticeFermion res4   (UGrid); res4 = zero;
+  LatticeFermion res4   (UGrid); res4 = Zero();
 
   LatticeFermion src   (FGrid);
   LatticeFermion src_NE(FGrid);

tests/debug/Test_split_laplacian.cc

@@ -38,13 +38,13 @@ int main (int argc, char ** argv)
 
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
+  Coordinate latt_size   = GridDefaultLatt();
   int nd   = latt_size.size();
   int ndm1 = nd-1;
 
-  std::vector<int> simd_layout = GridDefaultSimd(nd,vComplex::Nsimd());
+  Coordinate simd_layout = GridDefaultSimd(nd,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  Coordinate mpi_layout  = GridDefaultMpi();
-  std::vector<int> mpi_split (mpi_layout.size(),1);
+  Coordinate mpi_split (mpi_layout.size(),1);
 
   std::cout << " Full " << GridCmdVectorIntToString(latt_size)  << " subgrid"         <<std::endl;
   std::cout << " Full " << GridCmdVectorIntToString(mpi_layout) << " sub communicator"<<std::endl;
@@ -54,9 +54,9 @@ int main (int argc, char ** argv)
 						    simd_layout,
 						    mpi_layout);
 
-  std::vector<int> latt_m  = latt_size;   latt_m[nd-1] = 1;
+  Coordinate latt_m  = latt_size;   latt_m[nd-1] = 1;
-  std::vector<int> mpi_m   = mpi_layout;  mpi_m [nd-1] = 1;
+  Coordinate mpi_m   = mpi_layout;  mpi_m [nd-1] = 1;
-  std::vector<int> simd_m  = GridDefaultSimd(ndm1,vComplex::Nsimd()); simd_m.push_back(1);
+  Coordinate simd_m  = GridDefaultSimd(ndm1,vComplex::Nsimd()); simd_m.push_back(1);
 
 
   std::cout << " Requesting " << GridCmdVectorIntToString(latt_m)<< " subgrid"         <<std::endl;
@@ -88,7 +88,7 @@ int main (int argc, char ** argv)
 
   auto local_latt = GridN->LocalDimensions();
 
-  Full_cpy = zero;
+  Full_cpy = Zero();
   std::vector<int> seeds({1,2,3,4});
   GridParallelRNG          RNG(GridN);  RNG.SeedFixedIntegers(seeds);
 

tests/forces/Test_gp_plaq_force.cc

@@ -29,15 +29,14 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 using namespace std;
 using namespace Grid;
-using namespace Grid::QCD;
 
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
+  auto latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  auto mpi_layout  = GridDefaultMpi();
 
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
   GridRedBlackCartesian     RBGrid(&Grid);
@@ -85,8 +84,11 @@ int main (int argc, char ** argv)
     PokeIndex<LorentzIndex>(mom,mommu,mu);
 
     // fourth order exponential approx
-    parallel_for(auto i=mom.begin();i<mom.end();i++){ // exp(pmu dt) * Umu
+    auto Uprime_v = Uprime.View();
-      Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt ;
+    auto U_v = U.View();
+    auto mom_v = mom.View();
+    parallel_for(auto i=mom_v.begin();i<mom_v.end();i++){ // exp(pmu dt) * Umu
+      Uprime_v[i](mu) = U_v[i](mu) + mom_v[i](mu)*U_v[i](mu)*dt ;
     }
   }
 
@@ -96,7 +98,7 @@ int main (int argc, char ** argv)
   // Use derivative to estimate dS
   //////////////////////////////////////////////
 
-  LatticeComplex dS(&Grid); dS = zero;
+  LatticeComplex dS(&Grid); dS = Zero();
 
   for(int mu=0;mu<Nd;mu++){
 

tests/forces/Test_wilsonclover_force.cc

@@ -35,9 +35,9 @@ int main(int argc, char **argv)
 {
   Grid_init(&argc, &argv);
 
-  std::vector<int> latt_size = GridDefaultLatt();
+  auto latt_size = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
-  std::vector<int> mpi_layout = GridDefaultMpi();
+  auto mpi_layout = GridDefaultMpi();
 
   GridCartesian Grid(latt_size, simd_layout, mpi_layout);
   GridRedBlackCartesian RBGrid(&Grid);
@@ -63,9 +63,8 @@ int main(int argc, char **argv)
 
   LatticeGaugeField U(&Grid);
 
-  std::vector<int> site = {0, 0, 0, 0};
   SU3::HotConfiguration(pRNG, U);
-  //SU3::ColdConfiguration(pRNG, U);// Clover term zero
+  //SU3::ColdConfiguration(pRNG, U);// Clover term Zero()
 
   ////////////////////////////////////
   // Unmodified matrix element
@@ -107,9 +106,12 @@ int main(int argc, char **argv)
     Hmom -= real(sum(trace(mommu * mommu)));
     PokeIndex<LorentzIndex>(mom, mommu, mu);
 
-    parallel_for(int ss = 0; ss < mom._grid->oSites(); ss++)
+    auto Uprime_v = Uprime.View();
+    auto U_v      = U.View();
+    auto mom_v    = mom.View();
+    parallel_for(int ss = 0; ss < mom.Grid()->oSites(); ss++)
     {
-      Uprime[ss]._internal[mu] = ProjectOnGroup(Exponentiate(mom[ss]._internal[mu], dt, 12) * U[ss]._internal[mu]);
+      Uprime_v[ss]._internal[mu] = ProjectOnGroup(Exponentiate(mom_v[ss]._internal[mu], dt, 12) * U_v[ss]._internal[mu]);
     }
   }
 
@@ -125,11 +127,11 @@ int main(int argc, char **argv)
   //////////////////////////////////////////////
 
   LatticeComplex dS(&Grid);
-  dS = zero;
+  dS = Zero();
   LatticeComplex dSmom(&Grid);
-  dSmom = zero;
+  dSmom = Zero();
   LatticeComplex dSmom2(&Grid);
-  dSmom2 = zero;
+  dSmom2 = Zero();
 
   for (int mu = 0; mu < Nd; mu++)
   {

tests/lanczos/Test_WCMultiRep_lanczos.cc

@@ -30,7 +30,6 @@ directory
 
 using namespace std;
 using namespace Grid;
-using namespace Grid::QCD;
 
 
 //typedef WilsonCloverFermionR FermionOp;
@@ -154,11 +153,11 @@ NerscHmcCheckpointer<PeriodicGimplR> Checkpoint(CPparams);
   std::vector<ASymmFermionField>   ASevec(Nm, FGrid);
   
   for (int i = 0; i < 1; i++) {
-    std::cout << i << " / " << Nm << "Fund: grid pointer " << Fundevec[i]._grid
+    std::cout << i << " / " << Nm << "Fund: grid pointer " << Fundevec[i].Grid()
               << std::endl;
   };
   for (int i = 0; i < 1; i++) {
-    std::cout << i << " / " << Nm << "AS: grid pointer " << ASevec[i]._grid
+    std::cout << i << " / " << Nm << "AS: grid pointer " << ASevec[i].Grid()
               << std::endl;
   };
   

tests/lanczos/Test_compressed_lanczos.cc

@@ -36,7 +36,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 using namespace std;
 using namespace Grid;
-using namespace Grid::QCD;
 
 template<class Fobj,class CComplex,int nbasis>
 class LocalCoherenceLanczosScidac : public LocalCoherenceLanczos<Fobj,CComplex,nbasis>
@@ -58,7 +57,7 @@ public:
   {
     assert(this->subspace.size()==nbasis);
     emptyUserRecord record;
-    Grid::QCD::ScidacWriter WR(this->_FineGrid->IsBoss());
+    Grid::ScidacWriter WR(this->_FineGrid->IsBoss());
     WR.open(evecs_file);
     for(int k=0;k<nbasis;k++) {
       WR.writeScidacFieldRecord(this->subspace[k],record);
@@ -82,10 +81,10 @@ public:
     
     std::cout << GridLogIRL<< "checkpointFineRestore:  Reading evecs from "<<evecs_file<<std::endl;
     emptyUserRecord record;
-    Grid::QCD::ScidacReader RD ;
+    Grid::ScidacReader RD ;
     RD.open(evecs_file);
     for(int k=0;k<nbasis;k++) {
-      this->subspace[k].checkerboard=this->_checkerboard;
+      this->subspace[k].Checkerboard()=this->_checkerboard;
       RD.readScidacFieldRecord(this->subspace[k],record);
       
     }
@@ -96,7 +95,7 @@ public:
   {
     int n = this->evec_coarse.size();
     emptyUserRecord record;
-    Grid::QCD::ScidacWriter WR(this->_CoarseGrid->IsBoss());
+    Grid::ScidacWriter WR(this->_CoarseGrid->IsBoss());
     WR.open(evecs_file);
     for(int k=0;k<n;k++) {
       WR.writeScidacFieldRecord(this->evec_coarse[k],record);
@@ -119,7 +118,7 @@ public:
     assert(this->evals_coarse.size()==nvec);
     emptyUserRecord record;
     std::cout << GridLogIRL<< "checkpointCoarseRestore:  Reading evecs from "<<evecs_file<<std::endl;
-    Grid::QCD::ScidacReader RD ;
+    Grid::ScidacReader RD ;
     RD.open(evecs_file);
     for(int k=0;k<nvec;k++) {
       RD.readScidacFieldRecord(this->evec_coarse[k],record);
@@ -150,7 +149,7 @@ int main (int argc, char ** argv) {
   int     Ls = (int)Params.omega.size();
   RealD mass = Params.mass;
   RealD M5   = Params.M5;
-  std::vector<int> blockSize = Params.blockSize;
+  auto blockSize = Params.blockSize;
 
   // Grids
   GridCartesian         * UGrid     = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
@@ -160,11 +159,10 @@ int main (int argc, char ** argv) {
   GridCartesian         * FGrid     = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
   GridRedBlackCartesian * FrbGrid   = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
-  std::vector<int> fineLatt     = GridDefaultLatt();
+  auto fineLatt     = GridDefaultLatt();
   int dims=fineLatt.size();
   assert(blockSize.size()==dims+1);
-  std::vector<int> coarseLatt(dims);
+  Coordinate coarseLatt(dims);
-  std::vector<int> coarseLatt5d ;
 
   for (int d=0;d<coarseLatt.size();d++){
     coarseLatt[d] = fineLatt[d]/blockSize[d];    assert(coarseLatt[d]*blockSize[d]==fineLatt[d]);

tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc

@@ -58,7 +58,7 @@ public:
   {
     assert(this->subspace.size()==nbasis);
     emptyUserRecord record;
-    Grid::QCD::ScidacWriter WR(this->_FineGrid->IsBoss());
+    Grid::ScidacWriter WR(this->_FineGrid->IsBoss());
     WR.open(evecs_file);
     for(int k=0;k<nbasis;k++) {
       WR.writeScidacFieldRecord(this->subspace[k],record);
@@ -85,7 +85,7 @@ public:
     Grid::ScidacReader RD ;
     RD.open(evecs_file);
     for(int k=0;k<nbasis;k++) {
-      this->subspace[k].checkerboard=this->_checkerboard;
+      this->subspace[k].Checkerboard()=this->_checkerboard;
       RD.readScidacFieldRecord(this->subspace[k],record);
     }
     RD.close();
@@ -95,7 +95,7 @@ public:
   {
     int n = this->evec_coarse.size();
     emptyUserRecord record;
-    Grid::QCD::ScidacWriter WR(this->_CoarseGrid->IsBoss());
+    Grid::ScidacWriter WR(this->_CoarseGrid->IsBoss());
     WR.open(evecs_file);
     for(int k=0;k<n;k++) {
       WR.writeScidacFieldRecord(this->evec_coarse[k],record);

tests/solver/Test_dwf_mrhs_cg_mpi.cc

@@ -228,7 +228,7 @@ int main (int argc, char ** argv)
     std::cout << GridLogMessage<<" resid["<<n<<"]  "<< norm2(tmp)/norm2(src[n])<<std::endl;
   }
 
-  for(int s=0;s<nrhs;s++) result[s]=zero;
+  for(int s=0;s<nrhs;s++) result[s]=Zero();
   int blockDim = 0;//not used for BlockCGVec
   BlockConjugateGradient<FermionField>    BCGV  (BlockCGVec,blockDim,stp,10000);
   BCGV.PrintInterval=10;

tests/solver/Test_mobius_bcg.cc

@@ -30,7 +30,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 using namespace std;
 using namespace Grid;
-using namespace Grid::QCD;
 
 int main (int argc, char ** argv)
 {
@@ -42,12 +41,12 @@ int main (int argc, char ** argv)
 
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
+  auto latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  auto mpi_layout  = GridDefaultMpi();
-  std::vector<int> mpi_split (mpi_layout.size(),1);
+  Coordinate mpi_split (mpi_layout.size(),1);
-  std::vector<int> split_coor (mpi_layout.size(),1);
+  Coordinate split_coor (mpi_layout.size(),1);
-  std::vector<int> split_dim (mpi_layout.size(),1);
+  Coordinate split_dim (mpi_layout.size(),1);
 
   std::vector<ComplexD> boundary_phases(Nd,1.);
   boundary_phases[Nd-1]=-1.;
@@ -108,7 +107,7 @@ int main (int argc, char ** argv)
   FermionField tmp(FGrid);
   std::cout << GridLogMessage << "Made the Fermion Fields"<<std::endl;
 
-  for(int s=0;s<nrhs;s++) result[s]=zero;
+  for(int s=0;s<nrhs;s++) result[s]=Zero();
   GridParallelRNG pRNG5(FGrid);  pRNG5.SeedFixedIntegers(seeds);
   for(int s=0;s<nrhs;s++) {
     random(pRNG5,src[s]);
@@ -131,7 +130,6 @@ int main (int argc, char ** argv)
     std::cout << GridLogMessage << "Intialised 4D RNG "<<std::endl;
     SU3::HotConfiguration(pRNG,Umu);
     std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
-    std::cout << " Site zero "<< Umu._odata[0]   <<std::endl;
   } 
 
   /////////////////
@@ -170,7 +168,7 @@ int main (int argc, char ** argv)
   MdagMLinearOperator<MobiusFermionR,FermionField> HermOp(Ddwf);
   MdagMLinearOperator<MobiusFermionR,FermionField> HermOpCk(Dchk);
   ConjugateGradient<FermionField> CG((stp),100000);
-  s_res = zero;
+  s_res = Zero();
 
   CG(HermOp,s_src,s_res);
 
@@ -202,7 +200,7 @@ int main (int argc, char ** argv)
 
 
   for(int s=0;s<nrhs;s++){
-    result[s]=zero;
+    result[s]=Zero();
   }
 
   /////////////////////////////////////////////////////////////

tests/solver/Test_mobius_bcg_nosplit.cc

@@ -43,9 +43,9 @@ int main (int argc, char ** argv)
 
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
+  auto latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  auto mpi_layout  = GridDefaultMpi();
 
   std::vector<ComplexD> boundary_phases(Nd,1.);
   boundary_phases[Nd-1]=-1.;
@@ -72,7 +72,7 @@ int main (int argc, char ** argv)
   FermionField tmp(FGrid);
   std::cout << GridLogMessage << "Made the Fermion Fields"<<std::endl;
 
-  for(int s=0;s<nrhs;s++) result[s]=zero;
+  for(int s=0;s<nrhs;s++) result[s]=Zero();
   GridParallelRNG pRNG5(FGrid);  pRNG5.SeedFixedIntegers(seeds);
   for(int s=0;s<nrhs;s++) {
     random(pRNG5,src[s]);
@@ -116,7 +116,7 @@ int main (int argc, char ** argv)
   ConjugateGradient<FermionField> CG((stp),100000);
 
   for(int rhs=0;rhs<1;rhs++){
-    result[rhs] = zero;
+    result[rhs] = Zero();
     CG(HermOp,src[rhs],result[rhs]);
   }
 
@@ -129,7 +129,7 @@ int main (int argc, char ** argv)
   /////////////////////////////////////////////////////////////
   int blockDim = 0;//not used for BlockCGVec
   for(int s=0;s<nrhs;s++){
-      result[s]=zero;
+      result[s]=Zero();
    }
   BlockConjugateGradient<FermionField>    BCGV  (BlockCGrQVec,blockDim,stp,100000);
   {

tests/solver/Test_mobius_bcg_phys_nosplit.cc

@@ -43,9 +43,9 @@ int main (int argc, char ** argv)
 
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
+  auto latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  auto mpi_layout  = GridDefaultMpi();
 
   std::vector<ComplexD> boundary_phases(Nd,1.);
   boundary_phases[Nd-1]=-1.;
@@ -73,7 +73,7 @@ int main (int argc, char ** argv)
   FermionField tmp(FGrid);
   std::cout << GridLogMessage << "Made the Fermion Fields"<<std::endl;
 
-  for(int s=0;s<nrhs;s++) result[s]=zero;
+  for(int s=0;s<nrhs;s++) result[s]=Zero();
   GridParallelRNG pRNG4(UGrid);  pRNG4.SeedFixedIntegers(seeds);
   for(int s=0;s<nrhs;s++) {
     random(pRNG4,src4[s]);
@@ -120,7 +120,7 @@ int main (int argc, char ** argv)
   ConjugateGradient<FermionField> CG((stp),100000);
 
   for(int rhs=0;rhs<1;rhs++){
-    result[rhs] = zero;
+    result[rhs] = Zero();
     //    CG(HermOp,src[rhs],result[rhs]);
   }
 
@@ -133,7 +133,7 @@ int main (int argc, char ** argv)
   /////////////////////////////////////////////////////////////
   int blockDim = 0;//not used for BlockCGVec
   for(int s=0;s<nrhs;s++){
-    result[s]=zero;
+    result[s]=Zero();
   }
   BlockConjugateGradient<FermionField>    BCGV  (BlockCGrQVec,blockDim,stp,100000);
   {

tests/solver/Test_mobius_bcg_prec_nosplit.cc

@@ -43,9 +43,9 @@ int main (int argc, char ** argv)
 
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
+  auto latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  auto mpi_layout  = GridDefaultMpi();
 
   std::vector<ComplexD> boundary_phases(Nd,1.);
   boundary_phases[Nd-1]=-1.;
@@ -72,7 +72,7 @@ int main (int argc, char ** argv)
   FermionField tmp(FGrid);
   std::cout << GridLogMessage << "Made the Fermion Fields"<<std::endl;
 
-  for(int s=0;s<nrhs;s++) result[s]=zero;
+  for(int s=0;s<nrhs;s++) result[s]=Zero();
   GridParallelRNG pRNG5(FGrid);  pRNG5.SeedFixedIntegers(seeds);
   for(int s=0;s<nrhs;s++) {
     random(pRNG5,src[s]);
@@ -116,7 +116,7 @@ int main (int argc, char ** argv)
   ConjugateGradient<FermionField> CG((stp),100000);
 
   for(int rhs=0;rhs<1;rhs++){
-    result[rhs] = zero;
+    result[rhs] = Zero();
     CG(HermOp,src[rhs],result[rhs]);
   }
 
@@ -129,7 +129,7 @@ int main (int argc, char ** argv)
   /////////////////////////////////////////////////////////////
   int blockDim = 0;//not used for BlockCGVec
   for(int s=0;s<nrhs;s++){
-    result[s]=zero;
+    result[s]=Zero();
   }
 
 

tests/solver/Test_staggered_block_cg_prec.cc

@@ -72,7 +72,7 @@ int main (int argc, char ** argv)
 #if 1
   random(pRNG5,src);
 #else
-  src=zero;
+  src=Zero();
   ComplexField coor(FGrid);
   LatticeCoordinate(coor,0);
   for(int ss=0;ss<FGrid->oSites();ss++){
@@ -196,7 +196,7 @@ int main (int argc, char ** argv)
   std::cout << GridLogMessage << " Calling Block CG for "<<Ls <<" right hand sides" <<std::endl;
   std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
   Ds.ZeroCounters();
-  result_o=zero;
+  result_o=Zero();
   {
     double t1=usecond();
     BCG(HermOp,src_o,result_o);
@@ -216,7 +216,7 @@ int main (int argc, char ** argv)
   std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
   std::vector<FermionField> src_v   (Ls,UrbGrid);
   std::vector<FermionField> result_v(Ls,UrbGrid);
-  for(int s=0;s<Ls;s++) result_v[s] = zero;
+  for(int s=0;s<Ls;s++) result_v[s] = Zero();
   for(int s=0;s<Ls;s++) {
     FermionField src4(UGrid);
     ExtractSlice(src4,src,s,0);

tests/solver/Test_staggered_multishift.cc

@@ -31,7 +31,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 using namespace std;
 using namespace Grid;
-using namespace Grid::QCD;
 
 template<class d>
 struct scal {
@@ -52,9 +51,9 @@ int main (int argc, char ** argv)
 
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size   = GridDefaultLatt();
+  auto latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
+  auto mpi_layout  = GridDefaultMpi();
 
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
   GridRedBlackCartesian     RBGrid(&Grid);