From 3358a77c7a111193b6310ed7b68b2cee48d34775 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Mon, 25 May 2015 13:45:08 +0100
Subject: [PATCH] Better checkerboard tracking.

---
 lib/lattice/Grid_lattice_ET.h    | 41 ++++++++++++++++++++++++++
 lib/lattice/Grid_lattice_arith.h | 28 ++++++++++++++++--
 lib/lattice/Grid_lattice_base.h  | 50 ++++++++++++++++++++++++++++++++
 3 files changed, 117 insertions(+), 2 deletions(-)

diff --git a/lib/lattice/Grid_lattice_ET.h b/lib/lattice/Grid_lattice_ET.h
index b0205f2f..2a7172ea 100644
--- a/lib/lattice/Grid_lattice_ET.h
+++ b/lib/lattice/Grid_lattice_ET.h
@@ -91,6 +91,47 @@ inline void GridFromExpression( GridBase * &grid,const LatticeTrinaryExpression<
   GridFromExpression(grid,std::get<2>(expr.second));
 }
 
+
+//////////////////////////////////////////////////////////////////////////
+// Obtain the CB from an expression, ensuring conformable. This must follow a tree recursion
+//////////////////////////////////////////////////////////////////////////
+template<class T1, typename std::enable_if<is_lattice<T1>::value, T1>::type * =nullptr >
+inline void CBFromExpression(int &cb,const T1& lat)   // Lattice leaf
+{
+  if ( (cb==Odd) || (cb==Even) ) {
+    assert(cb==lat.checkerboard);
+  } 
+  cb=lat.checkerboard;
+  //  std::cout<<"Lattice leaf cb "<<cb<<std::endl;
+}
+template<class T1,typename std::enable_if<!is_lattice<T1>::value, T1>::type * = nullptr >
+inline void CBFromExpression(int &cb,const T1& notlat)   // non-lattice leaf
+{
+  //  std::cout<<"Non lattice leaf cb"<<cb<<std::endl;
+}
+template <typename Op, typename T1>
+inline void CBFromExpression(int &cb,const LatticeUnaryExpression<Op,T1 > &expr)
+{
+  CBFromExpression(cb,std::get<0>(expr.second));// recurse 
+  //  std::cout<<"Unary node cb "<<cb<<std::endl;
+}
+
+template <typename Op, typename T1, typename T2>
+inline void CBFromExpression(int &cb,const LatticeBinaryExpression<Op,T1,T2> &expr) 
+{
+  CBFromExpression(cb,std::get<0>(expr.second));// recurse
+  CBFromExpression(cb,std::get<1>(expr.second));
+  //  std::cout<<"Binary node cb "<<cb<<std::endl;
+}
+template <typename Op, typename T1, typename T2, typename T3>
+inline void CBFromExpression( int &cb,const LatticeTrinaryExpression<Op,T1,T2,T3 > &expr) 
+{
+  CBFromExpression(cb,std::get<0>(expr.second));// recurse
+  CBFromExpression(cb,std::get<1>(expr.second));
+  CBFromExpression(cb,std::get<2>(expr.second));
+  //  std::cout<<"Trinary node cb "<<cb<<std::endl;
+}
+
 ////////////////////////////////////////////
 // Unary operators and funcs
 ////////////////////////////////////////////
diff --git a/lib/lattice/Grid_lattice_arith.h b/lib/lattice/Grid_lattice_arith.h
index a6f59b56..dc5e57e8 100644
--- a/lib/lattice/Grid_lattice_arith.h
+++ b/lib/lattice/Grid_lattice_arith.h
@@ -9,6 +9,8 @@ namespace Grid {
   //////////////////////////////////////////////////////////////////////////////////////////////////////
   template<class obj1,class obj2,class obj3> strong_inline
     void mult(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const Lattice<obj3> &rhs){
+    ret.checkerboard = lhs.checkerboard;
+    conformable(ret,rhs);
     conformable(lhs,rhs);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
@@ -24,6 +26,8 @@ PARALLEL_FOR_LOOP
   
   template<class obj1,class obj2,class obj3> strong_inline
     void mac(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const Lattice<obj3> &rhs){
+    ret.checkerboard = lhs.checkerboard;
+    conformable(ret,rhs);
     conformable(lhs,rhs);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
@@ -39,6 +43,8 @@ PARALLEL_FOR_LOOP
   
   template<class obj1,class obj2,class obj3> strong_inline
     void sub(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const Lattice<obj3> &rhs){
+    ret.checkerboard = lhs.checkerboard;
+    conformable(ret,rhs);
     conformable(lhs,rhs);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
@@ -53,6 +59,8 @@ PARALLEL_FOR_LOOP
   }
   template<class obj1,class obj2,class obj3> strong_inline
     void add(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const Lattice<obj3> &rhs){
+    ret.checkerboard = lhs.checkerboard;
+    conformable(ret,rhs);
     conformable(lhs,rhs);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
@@ -71,6 +79,7 @@ PARALLEL_FOR_LOOP
   //////////////////////////////////////////////////////////////////////////////////////////////////////
   template<class obj1,class obj2,class obj3> strong_inline
     void mult(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const obj3 &rhs){
+    ret.checkerboard = lhs.checkerboard;
     conformable(lhs,ret);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
@@ -82,7 +91,8 @@ PARALLEL_FOR_LOOP
   
   template<class obj1,class obj2,class obj3> strong_inline
     void mac(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const obj3 &rhs){
-    conformable(lhs,ret);
+    ret.checkerboard = lhs.checkerboard;
+    conformable(ret,lhs);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
       obj1 tmp;
@@ -93,7 +103,8 @@ PARALLEL_FOR_LOOP
   
   template<class obj1,class obj2,class obj3> strong_inline
     void sub(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const obj3 &rhs){
-    conformable(lhs,ret);
+    ret.checkerboard = lhs.checkerboard;
+    conformable(ret,lhs);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
 #ifdef STREAMING_STORES
@@ -107,6 +118,7 @@ PARALLEL_FOR_LOOP
   }
   template<class obj1,class obj2,class obj3> strong_inline
     void add(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const obj3 &rhs){
+    ret.checkerboard = lhs.checkerboard;
     conformable(lhs,ret);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
@@ -125,6 +137,7 @@ PARALLEL_FOR_LOOP
   //////////////////////////////////////////////////////////////////////////////////////////////////////
     template<class obj1,class obj2,class obj3> strong_inline
     void mult(Lattice<obj1> &ret,const obj2 &lhs,const Lattice<obj3> &rhs){
+    ret.checkerboard = rhs.checkerboard;
     conformable(ret,rhs);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<rhs._grid->oSites();ss++){
@@ -140,6 +153,7 @@ PARALLEL_FOR_LOOP
   
   template<class obj1,class obj2,class obj3> strong_inline
     void mac(Lattice<obj1> &ret,const obj2 &lhs,const Lattice<obj3> &rhs){
+    ret.checkerboard = rhs.checkerboard;
     conformable(ret,rhs);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<rhs._grid->oSites();ss++){
@@ -155,6 +169,7 @@ PARALLEL_FOR_LOOP
   
   template<class obj1,class obj2,class obj3> strong_inline
     void sub(Lattice<obj1> &ret,const obj2 &lhs,const Lattice<obj3> &rhs){
+    ret.checkerboard = rhs.checkerboard;
     conformable(ret,rhs);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<rhs._grid->oSites();ss++){
@@ -169,6 +184,7 @@ PARALLEL_FOR_LOOP
   }
   template<class obj1,class obj2,class obj3> strong_inline
     void add(Lattice<obj1> &ret,const obj2 &lhs,const Lattice<obj3> &rhs){
+    ret.checkerboard = rhs.checkerboard;
     conformable(ret,rhs);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<rhs._grid->oSites();ss++){
@@ -184,6 +200,8 @@ PARALLEL_FOR_LOOP
   
   template<class sobj,class vobj> strong_inline
   void axpy(Lattice<vobj> &ret,sobj a,const Lattice<vobj> &x,const Lattice<vobj> &y){
+    ret.checkerboard = x.checkerboard;
+    conformable(ret,x);
     conformable(x,y);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<x._grid->oSites();ss++){
@@ -197,6 +215,8 @@ PARALLEL_FOR_LOOP
   }
   template<class sobj,class vobj> strong_inline
   void axpby(Lattice<vobj> &ret,sobj a,sobj b,const Lattice<vobj> &x,const Lattice<vobj> &y){
+    ret.checkerboard = x.checkerboard;
+    conformable(ret,x);
     conformable(x,y);
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<x._grid->oSites();ss++){
@@ -211,12 +231,16 @@ PARALLEL_FOR_LOOP
 
   template<class sobj,class vobj> strong_inline
   RealD axpy_norm(Lattice<vobj> &ret,sobj a,const Lattice<vobj> &x,const Lattice<vobj> &y){
+    ret.checkerboard = x.checkerboard;
+    conformable(ret,x);
     conformable(x,y);
     axpy(ret,a,x,y);
     return norm2(ret);
   }
   template<class sobj,class vobj> strong_inline
   RealD axpby_norm(Lattice<vobj> &ret,sobj a,sobj b,const Lattice<vobj> &x,const Lattice<vobj> &y){
+    ret.checkerboard = x.checkerboard;
+    conformable(ret,x);
     conformable(x,y);
     axpby(ret,a,b,x,y);
     return norm2(ret); // FIXME implement parallel norm in ss loop
diff --git a/lib/lattice/Grid_lattice_base.h b/lib/lattice/Grid_lattice_base.h
index 41f349b7..1d3b1efb 100644
--- a/lib/lattice/Grid_lattice_base.h
+++ b/lib/lattice/Grid_lattice_base.h
@@ -66,6 +66,16 @@ public:
   ////////////////////////////////////////////////////////////////////////////////
   template <typename Op, typename T1>                         strong_inline Lattice<vobj> & operator=(const LatticeUnaryExpression<Op,T1> &expr)
   {
+    GridBase *egrid(nullptr);
+    GridFromExpression(egrid,expr);
+    assert(egrid!=nullptr);
+    conformable(_grid,egrid);
+
+    int cb=-1;
+    CBFromExpression(cb,expr);
+    assert( (cb==Odd) || (cb==Even));
+    checkerboard=cb;
+
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
 #ifdef STREAMING_STORES
@@ -79,6 +89,16 @@ PARALLEL_FOR_LOOP
   }
   template <typename Op, typename T1,typename T2>             strong_inline Lattice<vobj> & operator=(const LatticeBinaryExpression<Op,T1,T2> &expr)
   {
+    GridBase *egrid(nullptr);
+    GridFromExpression(egrid,expr);
+    assert(egrid!=nullptr);
+    conformable(_grid,egrid);
+
+    int cb=-1;
+    CBFromExpression(cb,expr);
+    assert( (cb==Odd) || (cb==Even));
+    checkerboard=cb;
+
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
 #ifdef STREAMING_STORES
@@ -92,6 +112,16 @@ PARALLEL_FOR_LOOP
   }
   template <typename Op, typename T1,typename T2,typename T3> strong_inline Lattice<vobj> & operator=(const LatticeTrinaryExpression<Op,T1,T2,T3> &expr)
   {
+    GridBase *egrid(nullptr);
+    GridFromExpression(egrid,expr);
+    assert(egrid!=nullptr);
+    conformable(_grid,egrid);
+
+    int cb=-1;
+    CBFromExpression(cb,expr);
+    assert( (cb==Odd) || (cb==Even));
+    checkerboard=cb;
+
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
 #ifdef STREAMING_STORES
@@ -106,8 +136,15 @@ PARALLEL_FOR_LOOP
   //GridFromExpression is tricky to do
   template<class Op,class T1>
     Lattice(const LatticeUnaryExpression<Op,T1> & expr):    _grid(nullptr){
+
     GridFromExpression(_grid,expr);
     assert(_grid!=nullptr);
+
+    int cb=-1;
+    CBFromExpression(cb,expr);
+    assert( (cb==Odd) || (cb==Even));
+    checkerboard=cb;
+
     _odata.resize(_grid->oSites());
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
@@ -123,6 +160,12 @@ PARALLEL_FOR_LOOP
   Lattice(const LatticeBinaryExpression<Op,T1,T2> & expr):    _grid(nullptr){
     GridFromExpression(_grid,expr);
     assert(_grid!=nullptr);
+
+    int cb=-1;
+    CBFromExpression(cb,expr);
+    assert( (cb==Odd) || (cb==Even));
+    checkerboard=cb;
+
     _odata.resize(_grid->oSites());
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
@@ -138,6 +181,12 @@ PARALLEL_FOR_LOOP
   Lattice(const LatticeTrinaryExpression<Op,T1,T2,T3> & expr):    _grid(nullptr){
     GridFromExpression(_grid,expr);
     assert(_grid!=nullptr);
+
+    int cb=-1;
+    CBFromExpression(cb,expr);
+    assert( (cb==Odd) || (cb==Even));
+    checkerboard=cb;
+
     _odata.resize(_grid->oSites());
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
@@ -169,6 +218,7 @@ PARALLEL_FOR_LOOP
         return *this;
     }
     template<class robj> strong_inline Lattice<vobj> & operator = (const Lattice<robj> & r){
+      this->checkerboard = r.checkerboard;
       conformable(*this,r);
       std::cout<<"Lattice operator ="<<std::endl;
 PARALLEL_FOR_LOOP