Got the NERSC IO working and fixed a bug in cshift.

configure

@@ -4999,6 +4999,28 @@ _ACEOF
 fi
 done
 
+for ac_func in ntohll
+do :
+  ac_fn_c_check_func "$LINENO" "ntohll" "ac_cv_func_ntohll"
+if test "x$ac_cv_func_ntohll" = xyes; then :
+  cat >>confdefs.h <<_ACEOF
+#define HAVE_NTOHLL 1
+_ACEOF
+
+fi
+done
+
+for ac_func in be64toh
+do :
+  ac_fn_c_check_func "$LINENO" "be64toh" "ac_cv_func_be64toh"
+if test "x$ac_cv_func_be64toh" = xyes; then :
+  cat >>confdefs.h <<_ACEOF
+#define HAVE_BE64TOH 1
+_ACEOF
+
+fi
+done
+
 
 # Check whether --enable-simd was given.
 if test "${enable_simd+set}" = set; then :

configure.ac

@@ -25,6 +25,8 @@ AC_TYPE_UINT64_T
 
 # Checks for library functions.
 AC_CHECK_FUNCS([gettimeofday])
+AC_CHECK_FUNCS([ntohll])
+AC_CHECK_FUNCS([be64toh])
 
 AC_ARG_ENABLE([simd],[AC_HELP_STRING([--enable-simd=SSE|AVX|AVX2|AVX512],[Select instructions])],[ac_SIMD=${enable_simd}],[ac_SIMD=AVX2])
 

lib/Grid.h

@@ -51,6 +51,7 @@
 #include <Grid_where.h>
 #include <Grid_stencil.h>
 #include <qcd/Grid_QCD.h>
+#include <parallelIO/GridNerscIO.h>
 
 namespace Grid {
 

lib/Grid_communicator.h

@@ -35,6 +35,7 @@ class CartesianCommunicator {
     // Grid information queries
     /////////////////////////////////
     int                      IsBoss(void)            { return _processor==0; };
+    int                      BossRank(void)          { return 0; };
     int                      ThisRank(void)          { return _processor; };
     const std::vector<int> & ThisProcessorCoor(void) { return _processor_coor; };
     const std::vector<int> & ProcessorGrid(void)     { return _processors; };
@@ -49,6 +50,8 @@ class CartesianCommunicator {
     void GlobalSum(RealD &);
     void GlobalSumVector(RealD *,int N);
 
+    void GlobalSum(uint32_t &);
+
     void GlobalSum(ComplexF &c)
     {
       GlobalSumVector((float *)&c,2);
@@ -68,12 +71,10 @@ class CartesianCommunicator {
     }
     
     template<class obj> void GlobalSum(obj &o){
-
       typedef typename obj::scalar_type scalar_type;
       int words = sizeof(obj)/sizeof(scalar_type);
-
       scalar_type * ptr = (scalar_type *)& o;
-      GlobalSum(ptr,words);
+      GlobalSumVector(ptr,words);
     }
     ////////////////////////////////////////////////////////////
     // Face exchange

lib/Grid_comparison.h

@@ -8,42 +8,42 @@ namespace Grid {
     public:
       vInteger operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs == rhs;
+	  return TensorRemove(lhs) == TensorRemove(rhs);
 	}
     };
     template<class lobj,class robj> class vne {
     public:
       vInteger operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs != rhs;
+	  return TensorRemove(lhs) != TensorRemove(rhs);
 	}
     };
     template<class lobj,class robj> class vlt {
     public:
       vInteger operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs < rhs;
+	  return TensorRemove(lhs) < TensorRemove(rhs);
 	}
     };
     template<class lobj,class robj> class vle {
     public:
       vInteger operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs <= rhs;
+	  return TensorRemove(lhs) <= TensorRemove(rhs);
 	}
     };
     template<class lobj,class robj> class vgt {
     public:
       vInteger operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs > rhs;
+	  return TensorRemove(lhs) > TensorRemove(rhs);
 	}
     };
     template<class lobj,class robj> class vge {
     public:
       vInteger operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs >= rhs;
+	  return TensorRemove(lhs) >= TensorRemove(rhs);
 	}
     };
 
@@ -52,42 +52,42 @@ namespace Grid {
     public:
       Integer operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs == rhs;
+	  return TensorRemove(lhs) == TensorRemove(rhs);
 	}
     };
     template<class lobj,class robj> class sne {
     public:
       Integer operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs != rhs;
+	  return TensorRemove(lhs) != TensorRemove(rhs);
 	}
     };
     template<class lobj,class robj> class slt {
     public:
       Integer operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs < rhs;
+	  return TensorRemove(lhs) < TensorRemove(rhs);
 	}
     };
     template<class lobj,class robj> class sle {
     public:
       Integer operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs <= rhs;
+	  return TensorRemove(lhs) <= TensorRemove(rhs);
 	}
     };
     template<class lobj,class robj> class sgt {
     public:
       Integer operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs > rhs;
+	  return TensorRemove(lhs) > TensorRemove(rhs);
 	}
     };
     template<class lobj,class robj> class sge {
     public:
       Integer operator()(const lobj &lhs, const robj &rhs)
 	{ 
-	  return lhs >= rhs;
+	  return TensorRemove(lhs) >= TensorRemove(rhs);
 	}
     };
 

lib/Grid_config.h

@@ -16,6 +16,9 @@
 /* GRID_COMMS_NONE */
 /* #undef GRID_COMMS_NONE */
 
+/* Define to 1 if you have the `be64toh' function. */
+/* #undef HAVE_BE64TOH */
+
 /* Define to 1 if you have the `gettimeofday' function. */
 #define HAVE_GETTIMEOFDAY 1
 
@@ -31,6 +34,9 @@
 /* Define to 1 if you have the <memory.h> header file. */
 #define HAVE_MEMORY_H 1
 
+/* Define to 1 if you have the `ntohll' function. */
+/* #undef HAVE_NTOHLL */
+
 /* Define to 1 if you have the <stdint.h> header file. */
 #define HAVE_STDINT_H 1
 

lib/Grid_config.h.in

@@ -15,6 +15,9 @@
 /* GRID_COMMS_NONE */
 #undef GRID_COMMS_NONE
 
+/* Define to 1 if you have the `be64toh' function. */
+#undef HAVE_BE64TOH
+
 /* Define to 1 if you have the `gettimeofday' function. */
 #undef HAVE_GETTIMEOFDAY
 
@@ -30,6 +33,9 @@
 /* Define to 1 if you have the <memory.h> header file. */
 #undef HAVE_MEMORY_H
 
+/* Define to 1 if you have the `ntohll' function. */
+#undef HAVE_NTOHLL
+
 /* Define to 1 if you have the <stdint.h> header file. */
 #undef HAVE_STDINT_H
 

lib/Grid_lattice.h

@@ -37,7 +37,8 @@ public:
     // what about a default grid?
     //////////////////////////////////////////////////////////////////
     Lattice(GridBase *grid) : _grid(grid) {
-        _odata.reserve(_grid->oSites());
+      //        _odata.reserve(_grid->oSites());
+        _odata.resize(_grid->oSites());
         assert((((uint64_t)&_odata[0])&0xF) ==0);
         checkerboard=0;
     }
@@ -49,6 +50,15 @@ public:
         }
         return *this;
     }
+    template<class robj> inline Lattice<vobj> & operator = (const Lattice<robj> & r){
+      conformable(*this,r);
+      std::cout<<"Lattice operator ="<<std::endl;
+
+        for(int ss=0;ss<_grid->oSites();ss++){
+            this->_odata[ss]=r._odata[ss];
+        }
+        return *this;
+    }
 
     // *=,+=,-= operators inherit behvour from correspond */+/- operation
     template<class T> inline Lattice<vobj> &operator *=(const T &r) {

lib/cartesian/Grid_cartesian_base.h

@@ -6,34 +6,23 @@
 
 namespace Grid{
 
-  class GridRNG ; // Forward declaration;
+  //////////////////////////////////////////////////////////////////////
-
+  // Commicator provides information on the processor grid
+  //////////////////////////////////////////////////////////////////////
+  //    unsigned long _ndimension;
+  //    std::vector<int> _processors; // processor grid
+  //    int              _processor;  // linear processor rank
+  //    std::vector<int> _processor_coor;  // linear processor rank
+  //////////////////////////////////////////////////////////////////////
   class GridBase : public CartesianCommunicator {
+
 public:
 
     // Give Lattice access
     template<class object> friend class Lattice;
 
-    GridRNG *_rng;
-
     GridBase(std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};
             
-    //FIXME 
-    // protected:
-    // Lattice wide random support. not yet fully implemented. Need seed strategy
-    // and one generator per site.
-    // std::default_random_engine generator;
-    //    static std::mt19937  generator( 9 );
-    
-    //////////////////////////////////////////////////////////////////////
-    // Commicator provides information on the processor grid
-    //////////////////////////////////////////////////////////////////////
-    //    unsigned long _ndimension;
-    //    std::vector<int> _processors; // processor grid
-    //    int              _processor;  // linear processor rank
-    //    std::vector<int> _processor_coor;  // linear processor rank
-    //////////////////////////////////////////////////////////////////////
-
     // Physics Grid information.
     std::vector<int> _simd_layout;// Which dimensions get relayed out over simd lanes.
     std::vector<int> _fdimensions;// Global dimensions of array prior to cb removal
@@ -101,6 +90,13 @@ public:
       for(int d=0;d<_ndimension;d++) idx+=_ostride[d]*ocoor[d];
       return idx;
     }
+    inline void CoorFromIndex (std::vector<int>& coor,int index,std::vector<int> &dims){
+      coor.resize(_ndimension);
+      for(int d=0;d<_ndimension;d++){
+	coor[d] = index % dims[d];
+	index   = index / dims[d];
+      }
+    }
     inline void oCoorFromOindex (std::vector<int>& coor,int Oindex){
       coor.resize(_ndimension);
       for(int d=0;d<_ndimension;d++){
@@ -146,6 +142,7 @@ public:
     inline int lSites(void) { return _isites*_osites; }; 
     inline int gSites(void) { return _isites*_osites*_Nprocessors; }; 
     inline int Nd    (void) { return _ndimension;};
+
     inline const std::vector<int> &FullDimensions(void)         { return _fdimensions;};
     inline const std::vector<int> &GlobalDimensions(void)       { return _gdimensions;};
     inline const std::vector<int> &LocalDimensions(void)        { return _ldimensions;};
@@ -175,10 +172,10 @@ public:
       std::vector<int> coor(_ndimension);
 
       ProcessorCoorFromRank(rank,coor);
-      for(int mu=0;mu<_ndimension;mu++) gcoor[mu] = _ldimensions[mu]&coor[mu];
+      for(int mu=0;mu<_ndimension;mu++) gcoor[mu] = _ldimensions[mu]*coor[mu];
 
       iCoorFromIindex(coor,i_idx);
-      for(int mu=0;mu<_ndimension;mu++) gcoor[mu] += _rdimensions[mu]&coor[mu];
+      for(int mu=0;mu<_ndimension;mu++) gcoor[mu] += _rdimensions[mu]*coor[mu];
 
       oCoorFromOindex (coor,o_idx);
       for(int mu=0;mu<_ndimension;mu++) gcoor[mu] += coor[mu];

lib/communicator/Grid_communicator_fake.cc

@@ -17,6 +17,7 @@ CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
 void CartesianCommunicator::GlobalSum(float &){}
 void CartesianCommunicator::GlobalSumVector(float *,int N){}
 void CartesianCommunicator::GlobalSum(double &){}
+void CartesianCommunicator::GlobalSum(uint32_t &){}
 void CartesianCommunicator::GlobalSumVector(double *,int N){}
 
 // Basic Halo comms primitive

lib/communicator/Grid_communicator_mpi.cc

@@ -28,6 +28,9 @@ CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
   assert(Size==_Nprocessors);
 }
 
+void CartesianCommunicator::GlobalSum(uint32_t &u){
+  MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_SUM,communicator);
+}
 void CartesianCommunicator::GlobalSum(float &f){
   MPI_Allreduce(MPI_IN_PLACE,&f,1,MPI_FLOAT,MPI_SUM,communicator);
 }

lib/cshift/Grid_cshift_mpi.h

@@ -188,8 +188,8 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,Lattice<vobj> &r
 
     if ( comm_any ) {
 
-      for(int i=0;i<Nsimd;i++){
+      for(int i=0;i<Nsimd;i++){         // there is a reversal in extract merge
-	pointers[i] = (scalar_type *)&send_buf_extract[i][0];
+	pointers[i] = (scalar_type *)&send_buf_extract[Nsimd-1-i][0];
       }
       Gather_plane_extract(rhs,pointers,dimension,sx,cbmask);
 
@@ -238,9 +238,9 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,Lattice<vobj> &r
       for(int i=0;i<Nsimd;i++){
 	if ( permute_slice ) {
 	  PermuteMap=i^toggle_bit;
-	  pointers[i] = rpointers[PermuteMap];
+	  pointers[Nsimd-1-i] = rpointers[PermuteMap];
 	} else {
-	  pointers[i] = rpointers[i];
+	  pointers[Nsimd-1-i] = rpointers[i];
 	}
       }
 

lib/lattice/Grid_lattice_coordinate.h

@@ -5,16 +5,17 @@ namespace Grid {
 
     template<class iobj> inline void LatticeCoordinate(Lattice<iobj> &l,int mu)
     {
+      typedef typename iobj::scalar_type scalar_type;
+      typedef typename iobj::vector_type vector_type;
       GridBase *grid = l._grid;
       int Nsimd = grid->iSites();
       std::vector<int> gcoor;
-      std::vector<Integer> mergebuf(Nsimd);
+      std::vector<scalar_type> mergebuf(Nsimd);
-      std::vector<Integer *> mergeptr(Nsimd);
+      std::vector<scalar_type *> mergeptr(Nsimd);
-      vInteger vI;
+      vector_type vI;
       for(int o=0;o<grid->oSites();o++){
 	for(int i=0;i<grid->iSites();i++){
 	  grid->RankIndexToGlobalCoor(grid->ThisRank(),o,i,gcoor);
-	  //	  grid->RankIndexToGlobalCoor(0,o,i,gcoor);
 	  mergebuf[i]=gcoor[mu];
 	  mergeptr[i]=&mergebuf[i];
 	}

lib/lattice/Grid_lattice_peekpoke.h

@@ -43,6 +43,7 @@ namespace Grid {
         }
         return ret;
     };
+
     ////////////////////////////////////////////////////////////////////////////////////////////////////
     // Poke internal indices of a Lattice object
     ////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -88,26 +89,26 @@ namespace Grid {
       assert( sizeof(sobj)*Nsimd == sizeof(vobj));
 
       int rank,odx,idx;
-      grid->GlobalCoorToRankIndex(rank,odx,idx,site);
-
       // Optional to broadcast from node 0.
-      grid->Broadcast(0,s);
+      grid->GlobalCoorToRankIndex(rank,odx,idx,site);
+      grid->Broadcast(grid->BossRank(),s);
 
       std::vector<sobj> buf(Nsimd);
       std::vector<scalar_type *> pointers(Nsimd);  
-      for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
 
       // extract-modify-merge cycle is easiest way and this is not perf critical
-      extract(l._odata[odx],pointers);
+      if ( rank == grid->ThisRank() ) {
-      
+	for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
-      buf[idx] = s;
+	extract(l._odata[odx],pointers);
-
+	buf[idx] = s;
-      for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
+	for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
-      merge(l._odata[odx],pointers);
+	merge(l._odata[odx],pointers);
+      }
 
       return;
     };
 
+
     //////////////////////////////////////////////////////////
     // Peek a scalar object from the SIMD array
     //////////////////////////////////////////////////////////
@@ -138,6 +139,69 @@ namespace Grid {
       return;
     };
 
+    //////////////////////////////////////////////////////////
+    // Peek a scalar object from the SIMD array
+    //////////////////////////////////////////////////////////
+    template<class vobj,class sobj>
+    void peekLocalSite(sobj &s,Lattice<vobj> &l,std::vector<int> &site){
+        
+      GridBase *grid=l._grid;
+
+      typedef typename vobj::scalar_type scalar_type;
+      typedef typename vobj::vector_type vector_type;
+
+      int Nsimd = grid->Nsimd();
+
+      assert( l.checkerboard== l._grid->CheckerBoard(site));
+      assert( sizeof(sobj)*Nsimd == sizeof(vobj));
+
+      int odx,idx;
+      idx= grid->iIndex(site);
+      odx= grid->oIndex(site);
+
+      std::vector<sobj> buf(Nsimd);
+      std::vector<scalar_type *> pointers(Nsimd);  
+      for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
+
+      extract(l._odata[odx],pointers);
+      
+      s = buf[idx];
+
+      return;
+    };
+
+    template<class vobj,class sobj>
+    void pokeLocalSite(const sobj &s,Lattice<vobj> &l,std::vector<int> &site){
+
+      GridBase *grid=l._grid;
+
+      typedef typename vobj::scalar_type scalar_type;
+      typedef typename vobj::vector_type vector_type;
+
+      int Nsimd = grid->Nsimd();
+
+      assert( l.checkerboard== l._grid->CheckerBoard(site));
+      assert( sizeof(sobj)*Nsimd == sizeof(vobj));
+
+      int odx,idx;
+      idx= grid->iIndex(site);
+      odx= grid->oIndex(site);
+
+      std::vector<sobj> buf(Nsimd);
+      std::vector<scalar_type *> pointers(Nsimd);  
+      for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
+
+      // extract-modify-merge cycle is easiest way and this is not perf critical
+      extract(l._odata[odx],pointers);
+      
+      buf[idx] = s;
+
+      for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
+      merge(l._odata[odx],pointers);
+
+      return;
+    };
+
 }
 #endif
 

lib/lattice/Grid_lattice_reduction.h

@@ -40,6 +40,36 @@ namespace Grid {
       return nrm;
     }
 
+    template<class vobj>
+      inline typename vobj::scalar_object sum(const Lattice<vobj> &arg){
+
+      GridBase *grid=arg._grid;
+      int Nsimd = grid->Nsimd();
+
+      typedef typename vobj::scalar_object sobj;
+      typedef typename vobj::scalar_type   scalar_type;
+
+      vobj vsum;
+      sobj ssum;
+
+      vsum=zero;
+      ssum=zero;
+      for(int ss=0;ss<arg._grid->oSites(); ss++){
+	vsum = vsum + arg._odata[ss];
+      }
+      
+      std::vector<sobj>               buf(Nsimd);
+      std::vector<scalar_type *> pointers(Nsimd);  
+      for(int i=0;i<Nsimd;i++) pointers[i] = (scalar_type *)&buf[i];
+      extract(vsum,pointers);
+
+      for(int i=0;i<Nsimd;i++) ssum = ssum + buf[i];
+
+      arg._grid->GlobalSum(ssum);
+
+      return ssum;
+    }
+
 }
 #endif
 

lib/math/Grid_math_tensors.h

@@ -13,10 +13,12 @@ template<class vtype> class iScalar
 public:
   vtype _internal;
 
-  typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
+  typedef typename GridTypeMapper<vtype>::scalar_type   scalar_type;
   typedef typename GridTypeMapper<vtype>::vector_type vector_type;
   typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
   typedef iScalar<tensor_reduced_v> tensor_reduced;
+  typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
+  typedef iScalar<recurse_scalar_object> scalar_object;
 
   enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
 
@@ -27,7 +29,7 @@ public:
     
     iScalar(scalar_type s) : _internal(s) {};// recurse down and hit the constructor for vector_type
 
-    iScalar(Zero &z){ *this = zero; };
+    iScalar(const Zero &z){ *this = zero; };
 
     iScalar<vtype> & operator= (const Zero &hero){
         zeroit(*this);
@@ -73,11 +75,18 @@ public:
     operator ComplexD () const { return(TensorRemove(_internal)); };
     operator RealD () const { return(real(TensorRemove(_internal))); }
 
+
+    template<class T,typename std::enable_if<isGridTensor<T>::notvalue, T>::type* = nullptr > inline auto operator = (T arg) -> iScalar<vtype>
+    { 
+      _internal = arg;
+      return *this;
+    }
+
 };
 ///////////////////////////////////////////////////////////
 // Allows to turn scalar<scalar<scalar<double>>>> back to double.
 ///////////////////////////////////////////////////////////
-template<class T> inline typename std::enable_if<isGridTensor<T>::notvalue, T>::type TensorRemove(T arg) { return arg;}
+template<class T>     inline typename std::enable_if<isGridTensor<T>::notvalue, T>::type TensorRemove(T arg) { return arg;}
 template<class vtype> inline auto TensorRemove(iScalar<vtype> arg) -> decltype(TensorRemove(arg._internal))
 {
   return TensorRemove(arg._internal);
@@ -91,14 +100,17 @@ public:
   typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
   typedef typename GridTypeMapper<vtype>::vector_type vector_type;
   typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
+  typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
+  typedef iScalar<tensor_reduced_v> tensor_reduced;
+  typedef iVector<recurse_scalar_object,N> scalar_object;
+
 
   enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
-  typedef iScalar<tensor_reduced_v> tensor_reduced;
 
-    iVector(Zero &z){ *this = zero; };
+    iVector(const Zero &z){ *this = zero; };
     iVector() {};// Empty constructure
 
-    iVector<vtype,N> & operator= (Zero &hero){
+    iVector<vtype,N> & operator= (const Zero &hero){
         zeroit(*this);
         return *this;
     }
@@ -153,18 +165,27 @@ public:
 
   typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
   typedef typename GridTypeMapper<vtype>::vector_type vector_type;
-
   typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
+  typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
+  typedef iScalar<tensor_reduced_v> tensor_reduced;
+  typedef iMatrix<recurse_scalar_object,N> scalar_object;
 
   enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
-  typedef iScalar<tensor_reduced_v> tensor_reduced;
 
-  iMatrix(Zero &z){ *this = zero; };
+  iMatrix(const Zero &z){ *this = zero; };
   iMatrix() {};
-  iMatrix<vtype,N> & operator= (Zero &hero){
+  iMatrix<vtype,N> & operator= (const Zero &hero){
     zeroit(*this);
     return *this;
   }
+  template<class T,typename std::enable_if<isGridTensor<T>::notvalue, T>::type* = nullptr > inline auto operator = (T arg) -> iMatrix<vtype,N>
+    { 
+      zeroit(*this);
+      for(int i=0;i<N;i++)
+	_internal[i][i] = arg;
+      return *this;
+    }
+
   friend void zeroit(iMatrix<vtype,N> &that){
     for(int i=0;i<N;i++){
       for(int j=0;j<N;j++){

lib/math/Grid_math_traits.h

@@ -25,6 +25,7 @@ namespace Grid {
     typedef typename T::scalar_type scalar_type;
     typedef typename T::vector_type vector_type;
     typedef typename T::tensor_reduced tensor_reduced;
+    typedef typename T::scalar_object scalar_object;
     enum { TensorLevel = T::TensorLevel };
   };
 
@@ -36,6 +37,7 @@ namespace Grid {
     typedef RealF scalar_type;
     typedef RealF vector_type;
     typedef RealF tensor_reduced ;
+    typedef RealF scalar_object;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<RealD> {
@@ -43,6 +45,7 @@ namespace Grid {
     typedef RealD scalar_type;
     typedef RealD vector_type;
     typedef RealD tensor_reduced;
+    typedef RealD scalar_object;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<ComplexF> {
@@ -50,6 +53,7 @@ namespace Grid {
     typedef ComplexF scalar_type;
     typedef ComplexF vector_type;
     typedef ComplexF tensor_reduced;
+    typedef ComplexF scalar_object;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<ComplexD> {
@@ -57,6 +61,7 @@ namespace Grid {
     typedef ComplexD scalar_type;
     typedef ComplexD vector_type;
     typedef ComplexD tensor_reduced;
+    typedef ComplexD scalar_object;
     enum { TensorLevel = 0 };
   };
 
@@ -65,6 +70,7 @@ namespace Grid {
     typedef RealF  scalar_type;
     typedef vRealF vector_type;
     typedef vRealF tensor_reduced;
+    typedef RealF  scalar_object;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<vRealD> {
@@ -72,6 +78,7 @@ namespace Grid {
     typedef RealD  scalar_type;
     typedef vRealD vector_type;
     typedef vRealD tensor_reduced;
+    typedef RealD  scalar_object;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<vComplexF> {
@@ -79,6 +86,7 @@ namespace Grid {
     typedef ComplexF  scalar_type;
     typedef vComplexF vector_type;
     typedef vComplexF tensor_reduced;
+    typedef ComplexF  scalar_object;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<vComplexD> {
@@ -86,6 +94,7 @@ namespace Grid {
     typedef ComplexD  scalar_type;
     typedef vComplexD vector_type;
     typedef vComplexD tensor_reduced;
+    typedef ComplexD  scalar_object;
     enum { TensorLevel = 0 };
   };
   template<> class GridTypeMapper<vInteger> {
@@ -93,6 +102,7 @@ namespace Grid {
     typedef Integer  scalar_type;
     typedef vInteger vector_type;
     typedef vInteger tensor_reduced;
+    typedef Integer  scalar_object;
     enum { TensorLevel = 0 };
   };
 
@@ -102,6 +112,10 @@ namespace Grid {
     static const bool notvalue = false;
   };
 
+  template<> struct isGridTensor<int > {
+    static const bool value = false;
+    static const bool notvalue = true;
+  };
   template<> struct isGridTensor<RealD > {
     static const bool value = false;
     static const bool notvalue = true;

lib/parallelIO/GridNerscIO.h

@@ -0,0 +1,524 @@
+#ifndef GRID_NERSC_IO_H
+#define GRID_NERSC_IO_H
+
+#include <algorithm>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+
+#include <arpa/inet.h>
+
+namespace Grid {
+
+  using namespace QCD;
+
+////////////////////////////////////////////////////////////////////////////////
+// Some data types for intermediate storage
+////////////////////////////////////////////////////////////////////////////////
+  template<typename vtype> using iLorentzColour2x3 = iVector<iVector<iVector<vtype, Nc>, 2>, 4 >;
+typedef iLorentzColour2x3<Complex>  LorentzColour2x3;
+typedef iLorentzColour2x3<ComplexF> LorentzColour2x3F;
+typedef iLorentzColour2x3<ComplexD> LorentzColour2x3D;
+
+////////////////////////////////////////////////////////////////////////////////
+// header specification/interpretation
+////////////////////////////////////////////////////////////////////////////////
+class NerscField {
+ public:
+    // header strings (not in order)
+    int dimension[4];
+    std::string boundary[4]; 
+    int data_start;
+    std::string hdr_version;
+    std::string storage_format;
+    // Checks on data
+    double link_trace;
+    double plaquette;
+    uint32_t checksum;
+    unsigned int sequence_number;
+    std::string data_type;
+    std::string ensemble_id ;
+    std::string ensemble_label ;
+    std::string creator ;
+    std::string creator_hardware ;
+    std::string creation_date ;
+    std::string archive_date ;
+    std::string floating_point;
+};
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Write and read from fstream; comput header offset for payload
+////////////////////////////////////////////////////////////////////////////////
+inline unsigned int writeNerscHeader(NerscField &field,std::string file)
+{
+  std::ofstream fout(file,std::ios::out);
+  
+  fout.seekp(0,std::ios::beg);
+  fout << "BEGIN_HEADER"      << std::endl;
+  fout << "HDR_VERSION = "    << field.hdr_version    << std::endl;
+  fout << "DATATYPE = "       << field.data_type      << std::endl;
+  fout << "STORAGE_FORMAT = " << field.storage_format << std::endl;
+
+  for(int i=0;i<4;i++){
+    fout << "DIMENSION_" << i+1 << " = " << field.dimension[i] << std::endl ;
+  }
+  // just to keep the space and write it later
+  fout << "LINK_TRACE = " << std::setprecision(10) << field.link_trace << std::endl;
+  fout << "PLAQUETTE  = " << std::setprecision(10) << field.plaquette  << std::endl;
+  for(int i=0;i<4;i++){
+    fout << "BOUNDARY_"<<i+1<<" = " << field.boundary[i] << std::endl;
+  }
+  fout << "CHECKSUM = "<< std::hex << std::setw(16) << 0 << field.checksum << std::endl;
+
+  fout << "ENSEMBLE_ID = "     << field.ensemble_id      << std::endl;
+  fout << "ENSEMBLE_LABEL = "  << field.ensemble_label   << std::endl;
+  fout << "SEQUENCE_NUMBER = " << field.sequence_number  << std::endl;
+  fout << "CREATOR = "         << field.creator          << std::endl;
+  fout << "CREATOR_HARDWARE = "<< field.creator_hardware << std::endl;
+  fout << "CREATION_DATE = "   << field.creation_date    << std::endl;
+  fout << "ARCHIVE_DATE = "    << field.archive_date     << std::endl;
+  fout << "FLOATING_POINT = "  << field.floating_point   << std::endl;
+  fout << "END_HEADER"         << std::endl;
+  field.data_start = fout.tellp();
+  return field.data_start;
+}
+
+// A little helper
+inline void removeWhitespace(std::string &key)
+{
+  key.erase(std::remove_if(key.begin(), key.end(), ::isspace),key.end());
+}
+// for the header-reader
+inline int readNerscHeader(std::string file,GridBase *grid,  NerscField &field)
+{
+  int offset=0;
+  std::map<std::string,std::string> header;
+  std::string line;
+
+  //////////////////////////////////////////////////
+  // read the header
+  //////////////////////////////////////////////////
+  std::ifstream fin(file);
+
+  getline(fin,line); // read one line and insist is 
+
+  removeWhitespace(line);
+  assert(line==std::string("BEGIN_HEADER"));
+
+  do {
+    getline(fin,line); // read one line
+    int eq = line.find("=");
+    if(eq >0) {
+      std::string key=line.substr(0,eq);
+      std::string val=line.substr(eq+1);
+      removeWhitespace(key);
+      removeWhitespace(val);
+      
+      header[key] = val;
+    }
+  } while( line.find("END_HEADER") == std::string::npos );
+
+  field.data_start = fin.tellg();
+
+  //////////////////////////////////////////////////
+  // chomp the values
+  //////////////////////////////////////////////////
+
+  field.hdr_version    = header[std::string("HDR_VERSION")];
+  field.data_type      = header[std::string("DATATYPE")];
+  field.storage_format = header[std::string("STORAGE_FORMAT")];
+  
+  field.dimension[0] = std::stol(header[std::string("DIMENSION_1")]);
+  field.dimension[1] = std::stol(header[std::string("DIMENSION_2")]);
+  field.dimension[2] = std::stol(header[std::string("DIMENSION_3")]);
+  field.dimension[3] = std::stol(header[std::string("DIMENSION_4")]);
+
+  assert(grid->_ndimension == 4);
+  for(int d=0;d<4;d++){
+    assert(grid->_fdimensions[d]==field.dimension[d]);
+  }
+
+  field.link_trace = std::stod(header[std::string("LINK_TRACE")]);
+  field.plaquette  = std::stod(header[std::string("PLAQUETTE")]);
+
+  field.boundary[0] = header[std::string("BOUNDARY_1")];
+  field.boundary[1] = header[std::string("BOUNDARY_2")];
+  field.boundary[2] = header[std::string("BOUNDARY_3")];
+  field.boundary[3] = header[std::string("BOUNDARY_4")];
+
+  field.checksum = std::stoul(header[std::string("CHECKSUM")],0,16);
+  field.ensemble_id      = header[std::string("ENSEMBLE_ID")];
+  field.ensemble_label   = header[std::string("ENSEMBLE_LABEL")];
+  field.sequence_number  = std::stol(header[std::string("SEQUENCE_NUMBER")]);
+  field.creator          = header[std::string("CREATOR")];
+  field.creator_hardware = header[std::string("CREATOR_HARDWARE")];
+  field.creation_date    = header[std::string("CREATION_DATE")];
+  field.archive_date     = header[std::string("ARCHIVE_DATE")];
+  field.floating_point   = header[std::string("FLOATING_POINT")];
+
+  return field.data_start;
+}
+
+
+//////////////////////////////////////////////////////////////////////
+// Utilities
+//////////////////////////////////////////////////////////////////////
+inline void reconstruct3(LorentzColourMatrix & cm)
+{
+  const int x=0;
+  const int y=1;
+  const int z=2;
+  for(int mu=0;mu<4;mu++){
+    cm(mu)()(2,x) = adj(cm(mu)()(0,y)*cm(mu)()(1,z)-cm(mu)()(0,z)*cm(mu)()(1,y)); //x= yz-zy
+    cm(mu)()(2,y) = adj(cm(mu)()(0,z)*cm(mu)()(1,x)-cm(mu)()(0,x)*cm(mu)()(1,z)); //y= zx-xz
+    cm(mu)()(2,z) = adj(cm(mu)()(0,x)*cm(mu)()(1,y)-cm(mu)()(0,y)*cm(mu)()(1,x)); //z= xy-yx
+  }
+}
+
+
+ void inline be32toh(void *file_object,uint32_t bytes)
+ {
+   uint32_t * f = (uint32_t *)file_object;
+   for(int i=0;i*sizeof(uint32_t)<bytes;i++){  
+     f[i] = ntohl(f[i]);
+   }
+ }
+ void inline le32toh(void *file_object,uint32_t bytes)
+ {
+   uint32_t *fp = (uint32_t *)file_object;
+
+   uint32_t f;
+
+   for(int i=0;i*sizeof(uint32_t)<bytes;i++){  
+     f = fp[i];
+     // got network order and the network to host
+     f = ((f&0xFF)<<24) | ((f&0xFF00)<<8) | ((f&0xFF0000)>>8) | ((f&0xFF000000UL)>>24) ; 
+     fp[i] = ntohl(f);
+   }
+ }
+ void inline be64toh(void *file_object,uint32_t bytes)
+ {
+   uint64_t * f = (uint64_t *)file_object;
+   for(int i=0;i*sizeof(uint64_t)<bytes;i++){  
+     f[i] = ntohll(f[i]);
+   }
+ }
+ void inline le64toh(void *file_object,uint32_t bytes)
+ {
+   uint64_t *fp = (uint64_t *)file_object;
+   uint64_t f,g;
+
+   for(int i=0;i*sizeof(uint64_t)<bytes;i++){  
+     f = fp[i];
+     // got network order and the network to host
+     g = ((f&0xFF)<<24) | ((f&0xFF00)<<8) | ((f&0xFF0000)>>8) | ((f&0xFF000000UL)>>24) ; 
+     g = g << 32;
+     f = f >> 32;
+     g|= ((f&0xFF)<<24) | ((f&0xFF00)<<8) | ((f&0xFF0000)>>8) | ((f&0xFF000000UL)>>24) ; 
+     fp[i] = ntohl(g);
+   }
+ }
+
+inline void NerscChecksum(uint32_t *buf,uint32_t buf_size,uint32_t &csum)
+{
+  for(int i=0;i*sizeof(uint32_t)<buf_size;i++){
+    csum=csum+buf[i];
+  }
+}
+
+  template<class fobj,class sobj>
+  struct NerscSimpleMunger{
+    void operator() (fobj &in,sobj &out,uint32_t &csum){
+
+      for(int mu=0;mu<4;mu++){
+      for(int i=0;i<3;i++){
+      for(int j=0;j<3;j++){
+	out(mu)()(i,j) = in(mu)()(i,j);
+      }}}
+
+      NerscChecksum((uint32_t *)&in,sizeof(in),csum); 
+    };
+  };
+  template<class fobj,class sobj>
+  struct NerscSimpleUnmunger{
+    void operator() (sobj &in,fobj &out,uint32_t &csum){
+      for(int mu=0;mu<4;mu++){
+      for(int i=0;i<3;i++){
+      for(int j=0;j<3;j++){
+	out(mu)()(i,j) = in(mu)()(i,j);
+      }}}
+      NerscChecksum((uint32_t *)&out,sizeof(out),csum); 
+    };
+  };
+ 
+ template<class fobj,class sobj>
+ struct Nersc3x2munger{
+   void operator() (fobj &in,sobj &out,uint32_t &csum){
+     
+     NerscChecksum((uint32_t *)&in,sizeof(in),csum); 
+
+     for(int mu=0;mu<4;mu++){
+       for(int i=0;i<2;i++){
+       for(int j=0;j<3;j++){
+	   out(mu)()(i,j) = in(mu)(i)(j);
+       }}
+     }
+     reconstruct3(out);
+   }
+ };
+
+ template<class fobj,class sobj>
+ struct Nersc3x2unmunger{
+
+  void operator() (sobj &in,fobj &out,uint32_t &csum){
+
+    NerscChecksum((uint32_t *)&out,sizeof(out),csum); 
+
+    for(int mu=0;mu<4;mu++){
+      for(int i=0;i<2;i++){
+      for(int j=0;j<3;j++){
+	out(mu)(i)(j) = in(mu)()(i,j);
+      }}
+    }
+  }
+};
+
+////////////////////////////////////////////////////////////////////////////
+// Template wizardry to map types to strings for NERSC in an extensible way
+////////////////////////////////////////////////////////////////////////////
+ template<class vobj> struct NerscDataType {
+   static void DataType     (std::string &str) { str = std::string("4D_BINARY_UNKNOWN"); };
+   static void FloatingPoint(std::string &str) { str = std::string("IEEE64BIG"); };
+ };
+
+ template<> struct NerscDataType<iColourMatrix<ComplexD> > {
+   static void DataType     (std::string &str) { str = std::string("4D_SU3_GAUGE_3X3"); };
+   static void FloatingPoint(std::string &str) { str = std::string("IEEE64BIG");};
+ };
+
+ template<> struct NerscDataType<iColourMatrix<ComplexF> > {
+   static void DataType     (std::string &str) { str = std::string("4D_SU3_GAUGE_3X3"); };
+   static void FloatingPoint(std::string &str) { str = std::string("IEEE32BIG");};
+ };
+
+//////////////////////////////////////////////////////////////////////
+// Bit and Physical Checksumming and QA of data
+//////////////////////////////////////////////////////////////////////
+/*
+template<class vobj> inline uint32_t NerscChecksum(Lattice<vobj> & data)
+{
+  uint32_t sum;
+  for(int ss=0;ss<data._grid->Osites();ss++){
+    uint32_t *iptr = (uint32_t *)& data._odata[0] ;
+    for(int i=0;i<sizeof(vobj);i+=sizeof(uint32_t)){
+      sum=sum+iptr[i];
+    }
+  }
+  data._grid->globalSum(sum);
+  return sum;
+}
+*/
+template<class vobj> inline void NerscPhysicalCharacteristics(Lattice<vobj> & data,NerscField &header)
+{
+  header.data_type      = NerscDataType<vobj>::DataType;
+  header.floating_point = NerscDataType<vobj>::FloatingPoint;
+  return;
+}
+
+ template<> inline void NerscPhysicalCharacteristics(LatticeGaugeField & data,NerscField &header)
+{
+  NerscDataType<decltype(data._odata[0])>::DataType(header.data_type);
+  NerscDataType<decltype(data._odata[0])>::FloatingPoint(header.floating_point);
+  header.link_trace=1.0;
+  header.plaquette =1.0;
+}
+
+template<class vobj> inline void NerscStatisics(Lattice<vobj> & data,NerscField &header)
+{
+  assert(data._grid->_ndimension==4);
+
+  for(int d=0;d<4;d++)
+    header.dimension[d] = data._grid->_fdimensions[d];
+
+  // compute checksum and any physical properties contained for this type
+  //  header.checksum = NerscChecksum(data);
+
+  NerscPhysicalCharacteristics(data,header);
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Now the meat: the object readers
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+template<class vobj,class sobj,class fobj,class munger>
+inline void readNerscObject(Lattice<vobj> &Umu,std::string file,munger munge,int offset,std::string &format)
+{
+  GridBase *grid = Umu._grid;
+
+  int ieee32big = (format == std::string("IEEE32BIG"));
+  int ieee32    = (format == std::string("IEEE32"));
+  int ieee64big = (format == std::string("IEEE64BIG"));
+  int ieee64    = (format == std::string("IEEE64"));
+
+  // Find the location of each site and send to primary node
+  //  for(int site=0; site < Layout::vol(); ++site){
+  //     multi1d<int> coord = crtesn(site, Layout::lattSize());
+  //     for(int dd=0; dd<Nd; dd++){        /* dir */
+  //        cfg_in.readArray(su3_buffer, float_size, mat_size);
+  //
+  // Above from Chroma; defines loop order now that NERSC doc no longer
+  // available (how short sighted is that?)
+  {
+    std::ifstream fin(file,std::ios::binary|std::ios::in);
+    fin.seekg(offset);
+
+    Umu = zero;
+    uint32_t csum=0;
+    fobj file_object;
+    sobj munged;
+    
+    for(int t=0;t<grid->_fdimensions[3];t++){
+    for(int z=0;z<grid->_fdimensions[2];z++){
+    for(int y=0;y<grid->_fdimensions[1];y++){
+    for(int x=0;x<grid->_fdimensions[0];x++){
+
+      std::vector<int> site({x,y,z,t});
+
+      if ( grid->IsBoss() ) {
+	fin.read((char *)&file_object,sizeof(file_object));
+
+	if(ieee32big) be32toh((void *)&file_object,sizeof(file_object));
+	if(ieee32)    le32toh((void *)&file_object,sizeof(file_object));
+	if(ieee64big) be64toh((void *)&file_object,sizeof(file_object));
+	if(ieee64)    le64toh((void *)&file_object,sizeof(file_object));
+
+	munge(file_object,munged,csum);
+      }
+      // The boss who read the file has their value poked
+      pokeSite(munged,Umu,site);
+    }}}}
+  }
+}
+
+template<class vobj,class sobj,class fobj,class munger>
+inline void writeNerscObject(Lattice<vobj> &Umu,std::string file,munger munge,int offset,
+                       int sequence,double lt,double pl)
+{
+  GridBase *grid = Umu._grid;
+  NerscField header;
+  
+  //////////////////////////////////////////////////
+  // First write the header; this is in wrong place
+  //////////////////////////////////////////////////
+  assert(grid->_ndimension == 4);
+  for(int d=0;d<4;d++){
+    header.dimension[d]=grid->_fdimensions[d];
+    header.boundary [d]=std::string("PERIODIC");; 
+  }
+  header.hdr_version=std::string("WHATDAHECK");
+  //  header.storage_format=storage_format<vobj>::string; // use template specialisation
+  //  header.data_type=data_type<vobj>::string;
+  header.storage_format=std::string("debug");
+  header.data_type     =std::string("debug");
+
+  //FIXME; use template specialisation to fill these out
+  header.link_trace   =lt;
+  header.plaquette    =pl;
+  header.checksum     =0;
+
+  //
+  header.sequence_number =sequence;
+  header.ensemble_id     =std::string("UKQCD");
+  header.ensemble_label  =std::string("UKQCD");
+  header.creator         =std::string("Tadahito");
+  header.creator_hardware=std::string("BlueGene/Q");
+  header.creation_date   =std::string("AnnoDomini");
+  header.archive_date    =std::string("AnnoDomini");
+  header.floating_point  =std::string("IEEE64BIG");
+  //  header.data_start=;
+  //  unsigned int checksum;
+
+  //////////////////////////////////////////////////
+  // Now write the body
+  //////////////////////////////////////////////////
+  {
+    std::ofstream fout(file,std::ios::binary|std::ios::in);
+    fout.seekp(offset);
+
+    Umu = zero;
+    uint32_t csum=0;
+    fobj file_object;
+    sobj unmunged;
+    for(int t=0;t<grid->_fdimensions[3];t++){
+    for(int z=0;z<grid->_fdimensions[2];z++){
+    for(int y=0;y<grid->_fdimensions[1];y++){
+    for(int x=0;x<grid->_fdimensions[0];x++){
+      std::vector<int> site({x,y,z,t});
+      peekSite(unmunged,Umu,site);
+      munge(unmunged,file_object,csum);
+      // broadcast & insert
+      fout.write((char *)&file_object,sizeof(file_object));
+    }}}}
+  }
+}
+
+
+
+inline void readNerscConfiguration(LatticeGaugeField &Umu,NerscField& header,std::string file)
+{
+  GridBase *grid = Umu._grid;
+
+  int offset = readNerscHeader(file,Umu._grid,header);
+
+  std::string format(header.floating_point);
+
+  int ieee32big = (format == std::string("IEEE32BIG"));
+  int ieee32    = (format == std::string("IEEE32"));
+  int ieee64big = (format == std::string("IEEE64BIG"));
+  int ieee64    = (format == std::string("IEEE64"));
+
+  // depending on datatype, set up munger;
+  // munger is a function of <floating point, Real, data_type>
+  if ( header.data_type == std::string("4D_SU3_GAUGE") ) {
+    if ( ieee32 || ieee32big ) {
+      readNerscObject<vLorentzColourMatrix, LorentzColourMatrix, LorentzColour2x3F> 
+	(Umu,file,
+	 Nersc3x2munger<LorentzColour2x3F,LorentzColourMatrix>(),
+	 offset,format);
+    }
+    if ( ieee64 || ieee64big ) {
+      readNerscObject<vLorentzColourMatrix, LorentzColourMatrix, LorentzColour2x3D> 
+	(Umu,file,
+	 Nersc3x2munger<LorentzColour2x3D,LorentzColourMatrix>(),
+	 offset,format);
+    }
+  } else if ( header.data_type == std::string("4D_SU3_GAUGE_3X3") ) {
+    if ( ieee32 || ieee32big ) {
+      readNerscObject<vLorentzColourMatrix,LorentzColourMatrix,LorentzColourMatrixF>
+	(Umu,file,NerscSimpleMunger<LorentzColourMatrixF,LorentzColourMatrix>(),offset,format);
+    }
+    if ( ieee64 || ieee64big ) {
+      readNerscObject<vLorentzColourMatrix,LorentzColourMatrix,LorentzColourMatrixD>
+	(Umu,file,NerscSimpleMunger<LorentzColourMatrixD,LorentzColourMatrix>(),offset,format);
+    }
+  } else {
+    assert(0);
+  }
+
+}
+
+template<class vobj>
+inline void writeNerscConfiguration(Lattice<vobj> &Umu,NerscField &header,std::string file)
+{
+  GridBase &grid = Umu._grid;
+  
+  NerscStatisics(Umu,header);
+
+  int offset = writeNerscHeader(header,file);
+
+  writeNerscObject(Umu,NerscSimpleMunger<vobj,vobj>(),offset);
+}
+
+
+}
+#endif

lib/qcd/Grid_QCD.h

@@ -42,6 +42,8 @@ namespace QCD {
     typedef iColourMatrix<Complex >        ColourMatrix;
     typedef iSpinColourMatrix<Complex >    SpinColourMatrix;
     typedef iLorentzColourMatrix<Complex > LorentzColourMatrix;
+    typedef iLorentzColourMatrix<ComplexF > LorentzColourMatrixF;
+    typedef iLorentzColourMatrix<ComplexD > LorentzColourMatrixD;
 
     typedef iSpinVector<Complex >       SpinVector;
     typedef iColourVector<Complex >     ColourVector;
@@ -66,7 +68,7 @@ namespace QCD {
 
     typedef Lattice<vTReal>              LatticeReal;
     typedef Lattice<vTComplex>           LatticeComplex;
-    typedef Lattice<vInteger>            LatticeInteger; // Predicates for "where"
+    typedef Lattice<vTInteger>            LatticeInteger; // Predicates for "where"
     
     typedef Lattice<vColourMatrix>     LatticeColourMatrix;
     typedef Lattice<vSpinMatrix>       LatticeSpinMatrix;

tests/Grid_cshift.cc

@@ -0,0 +1,84 @@
+#include <Grid.h>
+#include <parallelIO/GridNerscIO.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> simd_layout({1,1,2,2});
+  std::vector<int> mpi_layout ({2,2,2,2});
+  std::vector<int> latt_size  ({8,8,8,16});
+    
+  GridCartesian     Fine(latt_size,simd_layout,mpi_layout);
+  GridRNG           FineRNG(&Fine);
+
+  LatticeComplex U(&Fine);
+  LatticeComplex ShiftU(&Fine);
+
+  LatticeComplex lex(&Fine);
+  lex=zero;
+  Integer stride =1;
+  {
+    double nrm;
+    LatticeComplex coor(&Fine);
+
+    for(int d=0;d<4;d++){
+      LatticeCoordinate(coor,d);
+      lex = lex + coor*stride;
+      stride=stride*latt_size[d];
+    }
+    U=lex;
+  }
+
+  TComplex cm;
+  
+  for(int dir=0;dir<4;dir++){
+    for(int shift=0;shift<latt_size[dir];shift++){
+      if ( Fine.IsBoss() ) 
+	std::cout<<"Shifting by "<<shift<<" in direction"<<dir<<std::endl;
+
+	ShiftU  = Cshift(U,dir,shift);    // Shift everything
+
+	std::vector<int> coor(4);
+
+	for(coor[3]=0;coor[3]<latt_size[3];coor[3]++){
+	for(coor[2]=0;coor[2]<latt_size[2];coor[2]++){
+	for(coor[1]=0;coor[1]<latt_size[1];coor[1]++){
+	for(coor[0]=0;coor[0]<latt_size[0];coor[0]++){
+	  
+	  peekSite(cm,ShiftU,coor);
+
+	  double nrm=norm2(U);
+
+	  std::vector<int> scoor(coor);
+	  scoor[dir] = (scoor[dir]+shift)%latt_size[dir];
+	  
+	  Integer slex = scoor[0]
+	    + latt_size[0]*scoor[1]
+	    + latt_size[0]*latt_size[1]*scoor[2]
+	    + latt_size[0]*latt_size[1]*latt_size[2]*scoor[3];
+
+	  Complex scm(slex);
+	  
+	  nrm = abs(scm-cm()()());
+	  std::vector<int> peer(4);
+	  int index=real(cm);
+	  Fine.CoorFromIndex(peer,index,latt_size);
+
+	  if (nrm > 0){
+	    cout<<"FAIL shift "<< shift<<" in dir "<< dir<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] = "<< cm()()()<<" expect "<<scm<<"  "<<nrm<<endl;
+	    cout<<"Got    "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
+	    index=real(scm);
+	    Fine.CoorFromIndex(peer,index,latt_size);
+	    cout<<"Expect "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
+	  }
+	}}}}
+    }
+  }
+
+  Grid_finalize();
+}

tests/Grid_nersc_io.cc

@@ -0,0 +1,62 @@
+#include <Grid.h>
+#include <parallelIO/GridNerscIO.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> simd_layout({1,1,2,2});
+  std::vector<int> mpi_layout ({2,1,1,2});
+  std::vector<int> latt_size  ({16,16,16,32});
+    
+  GridCartesian     Fine(latt_size,simd_layout,mpi_layout);
+  GridRNG           FineRNG(&Fine);
+  LatticeGaugeField Umu(&Fine);
+
+  std::vector<LatticeColourMatrix> U(4,&Fine);
+  
+  NerscField header;
+  
+  std::string file("./ckpoint_lat.4000");
+  readNerscConfiguration(Umu,header,file);
+
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = peekIndex<3>(Umu,mu);
+  }
+
+  // Painful ; fix syntactical niceness
+  LatticeComplex LinkTrace(&Fine);
+  LinkTrace=zero;
+  for(int mu=0;mu<Nd;mu++){
+    LinkTrace = LinkTrace + trace(U[mu]);
+  }
+
+  // (1+2+3)=6 = N(N-1)/2 terms
+  LatticeComplex Plaq(&Fine);
+  Plaq = zero;
+  for(int mu=1;mu<Nd;mu++){
+    for(int nu=0;nu<mu;nu++){
+      Plaq = Plaq + trace(U[mu]*Cshift(U[nu],mu,1)*adj(Cshift(U[mu],nu,1))*adj(U[nu]));
+    }
+  }
+
+
+  double vol = Fine.gSites();
+
+  Complex PlaqScale(1.0/vol/6.0/3.0);
+  TComplex Tp = sum(Plaq);
+  Complex p  = TensorRemove(Tp);
+  std::cout << "calculated plaquettes " <<p*PlaqScale<<std::endl;
+
+  Complex LinkTraceScale(1.0/vol/4.0/3.0);
+  TComplex Tl = sum(LinkTrace);
+  Complex l  = TensorRemove(Tl);
+  std::cout << "calculated link trace " <<l*LinkTraceScale<<std::endl;
+
+  Grid_finalize();
+}

tests/Makefile.am

@@ -5,10 +5,16 @@ AM_LDFLAGS = -L$(top_srcdir)/lib
 #
 # Test code
 #
-bin_PROGRAMS = Grid_main test_Grid_stencil
+bin_PROGRAMS = Grid_main test_Grid_stencil Grid_nersc_io Grid_cshift
 
 Grid_main_SOURCES = Grid_main.cc
 Grid_main_LDADD = -lGrid
 
+Grid_nersc_io_SOURCES = Grid_nersc_io.cc
+Grid_nersc_io_LDADD = -lGrid
+
+Grid_cshift_SOURCES = Grid_cshift.cc
+Grid_cshift_LDADD = -lGrid
+
 test_Grid_stencil_SOURCES = test_Grid_stencil.cc
 test_Grid_stencil_LDADD = -lGrid