From 7b110751021514e4ec2ef1aa79a36ed970475864 Mon Sep 17 00:00:00 2001
From: Christopher Kelly <ckelly@bnl.gov>
Date: Tue, 20 Jun 2023 17:09:56 -0400
Subject: [PATCH] The user can now specify the implementation of Cshift used by
 the PaddedCell class through a virtual base class API. Implementations for
 default (regular Cshift) and for gauge links (which respects the gauge BCs)
 Fixed const-correctness for PaddedCell and ConjugateGimpl::setDirections
 Modified test code for padded-cell implementation of staple, rect-staple to
 use cconj BCs

---
 Grid/lattice/PaddedCell.h                    | 28 +++++++++++++++-----
 Grid/qcd/action/gauge/GaugeImplementations.h |  2 +-
 tests/debug/Test_padded_cell_staple.cc       | 23 +++++++++-------
 3 files changed, 36 insertions(+), 17 deletions(-)

diff --git a/Grid/lattice/PaddedCell.h b/Grid/lattice/PaddedCell.h
index 863d361a..2efd2288 100644
--- a/Grid/lattice/PaddedCell.h
+++ b/Grid/lattice/PaddedCell.h
@@ -28,12 +28,28 @@ Author: Peter Boyle pboyle@bnl.gov
 
 NAMESPACE_BEGIN(Grid);
 
+//Allow the user to specify how the C-shift is performed, e.g. to respect the appropriate boundary conditions
+template<typename vobj>
+struct CshiftImplBase{
+  virtual Lattice<vobj> Cshift(const Lattice<vobj> &in, int dir, int shift) const = 0;
+  virtual ~CshiftImplBase(){}
+};
+template<typename vobj>
+struct CshiftImplDefault: public CshiftImplBase<vobj>{
+  Lattice<vobj> Cshift(const Lattice<vobj> &in, int dir, int shift) const override{ return Grid::Cshift(in,dir,shift); }
+};
+template<typename Gimpl>
+struct CshiftImplGauge: public CshiftImplBase<typename Gimpl::GaugeLinkField::vector_object>{
+  typename Gimpl::GaugeLinkField Cshift(const typename Gimpl::GaugeLinkField &in, int dir, int shift) const override{ return Gimpl::CshiftLink(in,dir,shift); }
+};  
+
 class PaddedCell {
 public:
   GridCartesian * unpadded_grid;
   int dims;
   int depth;
   std::vector<GridCartesian *> grids;
+
   ~PaddedCell()
   {
     DeleteGrids();
@@ -77,7 +93,7 @@ public:
     }
   };
   template<class vobj>
-  inline Lattice<vobj> Extract(Lattice<vobj> &in)
+  inline Lattice<vobj> Extract(const Lattice<vobj> &in)
   {
     Lattice<vobj> out(unpadded_grid);
 
@@ -88,19 +104,19 @@ public:
     return out;
   }
   template<class vobj>
-  inline Lattice<vobj> Exchange(Lattice<vobj> &in)
+  inline Lattice<vobj> Exchange(const Lattice<vobj> &in, const CshiftImplBase<vobj> &cshift = CshiftImplDefault<vobj>())
   {
     GridBase *old_grid = in.Grid();
     int dims = old_grid->Nd();
     Lattice<vobj> tmp = in;
     for(int d=0;d<dims;d++){
-      tmp = Expand(d,tmp); // rvalue && assignment
+      tmp = Expand(d,tmp,cshift); // rvalue && assignment
     }
     return tmp;
   }
   // expand up one dim at a time
   template<class vobj>
-  inline Lattice<vobj> Expand(int dim,Lattice<vobj> &in)
+  inline Lattice<vobj> Expand(int dim, const Lattice<vobj> &in, const CshiftImplBase<vobj> &cshift = CshiftImplDefault<vobj>())
   {
     GridBase *old_grid = in.Grid();
     GridCartesian *new_grid = grids[dim];//These are new grids
@@ -117,12 +133,12 @@ public:
       InsertSliceLocal(in,padded,x,depth+x,dim);
     }
     // High bit
-    shifted = Cshift(in,dim,depth);
+    shifted = cshift.Cshift(in,dim,depth);
     for(int x=0;x<depth;x++){
       InsertSliceLocal(shifted,padded,local[dim]-depth+x,depth+local[dim]+x,dim);
     }
     // Low bit
-    shifted = Cshift(in,dim,-depth);
+    shifted = cshift.Cshift(in,dim,-depth);
     for(int x=0;x<depth;x++){
       InsertSliceLocal(shifted,padded,x,x,dim);
     }
diff --git a/Grid/qcd/action/gauge/GaugeImplementations.h b/Grid/qcd/action/gauge/GaugeImplementations.h
index f518b236..ae9a5c4b 100644
--- a/Grid/qcd/action/gauge/GaugeImplementations.h
+++ b/Grid/qcd/action/gauge/GaugeImplementations.h
@@ -176,7 +176,7 @@ public:
       return PeriodicBC::CshiftLink(Link,mu,shift);
   }
 
-  static inline void       setDirections(std::vector<int> &conjDirs) { _conjDirs=conjDirs; }
+  static inline void       setDirections(const std::vector<int> &conjDirs) { _conjDirs=conjDirs; }
   static inline std::vector<int> getDirections(void) { return _conjDirs; }
   static inline bool isPeriodicGaugeField(void) { return false; }
 };
diff --git a/tests/debug/Test_padded_cell_staple.cc b/tests/debug/Test_padded_cell_staple.cc
index 5a7f412a..bed00cad 100644
--- a/tests/debug/Test_padded_cell_staple.cc
+++ b/tests/debug/Test_padded_cell_staple.cc
@@ -106,9 +106,10 @@ public:
 
     GaugeMat U_mu = PeekIndex<LorentzIndex>(U, mu);
     GaugeMat U_nu = PeekIndex<LorentzIndex>(U, nu);
-  
-    GaugeMat Ug_mu = Ghost.Exchange(U_mu);
-    GaugeMat Ug_nu = Ghost.Exchange(U_nu);
+
+    CshiftImplGauge<Gimpl> cshift_impl;
+    GaugeMat Ug_mu = Ghost.Exchange(U_mu, cshift_impl);
+    GaugeMat Ug_nu = Ghost.Exchange(U_nu, cshift_impl);
     GridBase *ggrid = Ug_mu.Grid();
 
     GaugeMat gStaple(ggrid);
@@ -299,11 +300,11 @@ public:
   static void RectStaplePadded(GaugeMat &Stap, const GaugeLorentz &U,
 			       int mu) {
     PaddedCell Ghost(2,U.Grid());
-    GaugeLorentz Ug = Ghost.Exchange(const_cast<GaugeLorentz&>(U));
-    GridBase *ggrid = Ug.Grid();
-
+    GridBase *ggrid = Ghost.grids.back();
+    
+    CshiftImplGauge<Gimpl> cshift_impl;
     std::vector<GaugeMat> Ug_dirs(Nd,ggrid);
-    for(int i=0;i<Nd;i++) Ug_dirs[i] = PeekIndex<LorentzIndex>(Ug, i);
+    for(int i=0;i<Nd;i++) Ug_dirs[i] = Ghost.Exchange(PeekIndex<LorentzIndex>(U, i), cshift_impl);
 
     GaugeMat gStaple(ggrid);
 
@@ -495,11 +496,14 @@ int main (int argc, char ** argv)
 
   SU<Nc>::HotConfiguration(pRNG,U);
 
-  typedef PeriodicGimplD Gimpl;
+  //typedef PeriodicGimplD Gimpl;
+  typedef ConjugateGimplD Gimpl;
+  std::vector<int> conj_dirs(Nd,0); conj_dirs[0]=1; conj_dirs[3]=1;
+  Gimpl::setDirections(conj_dirs);
+
   typedef typename WilsonLoopsTest<Gimpl>::GaugeMat GaugeMat;
   typedef typename WilsonLoopsTest<Gimpl>::GaugeLorentz GaugeLorentz;
 
-#if 0
   std::cout << "Checking Staple" << std::endl;
   for(int mu=0;mu<Nd;mu++){
     for(int nu=0;nu<Nd;nu++){
@@ -514,7 +518,6 @@ int main (int argc, char ** argv)
       }
     }
   }
-#endif
 
   std::cout << "Checking RectStaple" << std::endl;
   for(int mu=0;mu<Nd;mu++){