Implemented gauge transform via stencil.

Now have ability to do Vertex AND Edge grids
Should now have no barriers to
a) Double Storing links for fermion operators / laplacian
b) Laplace or Wilson operators

Grid/cartesian/CartesianCrossIcosahedron.h

@@ -124,6 +124,34 @@ public:
   virtual int ownsNorthPole(void)   const override { return hasNorthPole; };
   virtual int ownsSouthPole(void)   const override { return hasSouthPole; };
   virtual int CartesianOsites(void) const override { return cartesianOsites; };
+  virtual int64_t PoleIdxForOcoor(Coordinate &Coor) override
+  {
+    // Work out the pole_osite. Pick the higher dims
+    Coordinate rdims;
+    Coordinate ocoor;
+    int64_t pole_idx;
+    int Ndm1 = this->Nd()-1;
+    for(int d=2;d<Ndm1;d++){
+      int dd=d-2;
+      rdims.push_back(this->_rdimensions[d]);
+      ocoor.push_back(Coor[d]%this->_rdimensions[d]);
+    }
+    Lexicographic::IndexFromCoor(ocoor,pole_idx,rdims);
+    return pole_idx;
+  }
+  virtual int64_t PoleSiteForOcoor(Coordinate &Coor) override
+  {
+    int Ndm1 = this->Nd()-1;
+    int64_t pole_idx = this->PoleIdxForOcoor(Coor);
+    int64_t pole_osite;
+    if ( Coor[Ndm1] >= HemiPatches ) {
+      pole_osite = pole_idx + this->NorthPoleOsite();
+    } else {
+      pole_osite = pole_idx + this->SouthPoleOsite();
+    }
+    return pole_osite;
+  }
+
 
   void InitPoles(void)
   {

Grid/cartesian/Cartesian_base.h

@@ -97,6 +97,8 @@ public:
   virtual int NorthPoleOsites(void) const { std::cout << "base osites" <<std::endl;return 0; };
   virtual int SouthPoleOsites(void) const { std::cout << "base osites" <<std::endl;return 0; };
   virtual int CartesianOsites(void) const { return this->oSites(); };
+  virtual int64_t PoleIdxForOcoor(Coordinate &Coor) { return 0;};
+  virtual int64_t PoleSiteForOcoor(Coordinate &Coor){ return 0;}
 
   ////////////////////////////////////////////////////////////////
   // Checkerboarding interface is virtual and overridden by 

Grid/stencil/IcosahedralStencil.h

@@ -578,11 +578,22 @@ public:
     // Loop over L^2 x T x npatch and the
     assert(grid->isIcosahedral());
   }
-  void NearestNeighbourStencil(void)
+  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+  // VertexInputs = true  implies the neighbour has vertex support
+  // VertexInputs = false implies the neighbout has edge support
+  // 
+  // isVertex implies must generate stencil entries to evaluate result on north/south pole
+  //
+  // These are independent:
+  //     can apply a vertex support gauge transform to edge supported gauge field
+  //     can apply a vertex supported link double store to edge supported gauge field
+  //     can apply a vertex supported laplace or dirac operator vertex supported matter field
+  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+  void NearestNeighbourStencil(int vertexOutputs)
   {
     GridBase * grid = this->_grid;
 
-    int isVertex = grid->isIcosahedralVertex();
+    int vertexInputs = grid->isIcosahedralVertex(); 
 
     int osites  = grid->oSites();
 
@@ -599,13 +610,14 @@ public:
     this->_entries.resize(this->_npoints * cart_sites);
     this->_entries_p = &_entries[0];
 
+    int nd = grid->Nd();
+    int L  = grid->LocalDimensions()[0];
+
     for(uint64_t site=0;site<cart_sites; site ++) {
 
       Coordinate Coor;
       Coordinate NbrCoor;
 
-      int nd = grid->Nd();
-      int L  = grid->LocalDimensions()[0];
 
       Integer lexXp = site*np  ;
       Integer lexYp = site*np+1;
@@ -669,42 +681,64 @@ public:
       int YmHemiPatch  = YmCoor[nd-1]%HemiPatches;
       int YmHemisphere = YmCoor[nd-1]/HemiPatches;
 
-      if ( isVertex ) assert(0);
+      if ( vertexInputs ) {
-
+	////////////////////////////////////////////////
-      ////////////////////////////////////////////////
+	// XpCoor stencil entry; consider isPole case
-      // XpCoor stencil entry
+	////////////////////////////////////////////////
-      // Store in look up table
+	SE._is_local     = true;
-      ////////////////////////////////////////////////
+	SE._missing_link = false;
-      // Basis rotates dictates BOTH adjoint and polarisation
+	SE._offset       = grid->oIndex(XpCoor);
-      // Could reduce the amount of information stored here
+	if ( isPoleX ) {
-      SE._adjoint      = false;
+	  SE._offset     = grid->PoleSiteForOcoor(Coor);
-      SE._is_local     = true;
+	}
-      SE._missing_link = false;
+	SE._polarisation = IcosahedronPatchY;
-      if ( DpHemiPatch != HemiPatch && south ) {
+	SE._adjoint      = false;
-	SE._offset       = grid->oIndex(DpCoor);
+	acceleratorPut(this->_entries[lexXp],SE);
+	////////////////////////////////////////////////
+	// for YpCoor
+	////////////////////////////////////////////////
+	SE._is_local     = true;
+	SE._missing_link = false;
+	SE._offset       = grid->oIndex(YpCoor);
+	if ( isPoleY ) {
+	  SE._offset     = grid->PoleSiteForOcoor(Coor); 
+	}
 	SE._polarisation = IcosahedronPatchX;
-	SE._adjoint      = true;
+	SE._adjoint      = false;
-      } else {
+	acceleratorPut(this->_entries[lexYp],SE);
+      } else { 
+	////////////////////////////////////////////////
+	// XpCoor stencil entry
+	// Store in look up table
+	////////////////////////////////////////////////
+	// Basis rotates dictates BOTH adjoint and polarisation
+	// Could reduce the amount of information stored here
+	SE._is_local     = true;
+	SE._missing_link = false;
 	SE._offset       = grid->oIndex(XpCoor);
 	SE._polarisation = IcosahedronPatchY;
-      }
+	SE._adjoint      = false;
-      acceleratorPut(this->_entries[lexXp],SE);
+	if ( DpHemiPatch != HemiPatch && south ) {
-      
+	  SE._offset       = grid->oIndex(DpCoor);
-      ////////////////////////////////////////////////
+	  SE._polarisation = IcosahedronPatchX;
-      // for YpCoor
+	  SE._adjoint      = true;
-      ////////////////////////////////////////////////
+	}
-      SE._adjoint      = false;
+	acceleratorPut(this->_entries[lexXp],SE);
-      SE._is_local     = true;
+	////////////////////////////////////////////////
-      SE._missing_link = false;
+	// for YpCoor
-      if ( YpHemiPatch != HemiPatch && north ) {
+	////////////////////////////////////////////////
-	SE._offset       = grid->oIndex(DpCoor);
+	SE._is_local     = true;
-	SE._polarisation = IcosahedronPatchY;
+	SE._missing_link = false;
-	SE._adjoint      = true;
-      } else {
 	SE._offset       = grid->oIndex(YpCoor);
 	SE._polarisation = IcosahedronPatchX;
+	SE._adjoint      = false;
+	if ( YpHemiPatch != HemiPatch && north ) {
+	  SE._offset       = grid->oIndex(DpCoor);
+	  SE._polarisation = IcosahedronPatchY;
+	  SE._adjoint      = true;
+	}
+	acceleratorPut(this->_entries[lexYp],SE);
       }
-      acceleratorPut(this->_entries[lexYp],SE);
 
       SE._adjoint      = false;
       SE._is_local     = true;
@@ -713,24 +747,20 @@ public:
       // XmCoor stencil entry
       // Store in look up table
       ////////////////////////////////////////////////
+      SE._offset       = grid->oIndex(XmCoor);
+      SE._polarisation = IcosahedronPatchDiagonal;
       if ( XmHemiPatch != HemiPatch && north ) {
-	SE._offset       = grid->oIndex(XmCoor);
 	SE._polarisation = IcosahedronPatchY; // nbrs Y instead of diagonal in North hemisphere exceptional case
-      } else {
-	SE._offset       = grid->oIndex(XmCoor);
-	SE._polarisation = IcosahedronPatchDiagonal;
       }
       acceleratorPut(this->_entries[lexXm],SE);
       
       ////////////////////////////////////////////////
       // for YmCoor
       ////////////////////////////////////////////////
+      SE._offset       = grid->oIndex(YmCoor);
+      SE._polarisation = IcosahedronPatchDiagonal;
       if ( YmHemiPatch != HemiPatch && south ) {
-	SE._offset       = grid->oIndex(YmCoor);
 	SE._polarisation = IcosahedronPatchX;   // Basis rotates
-      } else {
-	SE._offset       = grid->oIndex(YmCoor);
-	SE._polarisation = IcosahedronPatchDiagonal;
       }
       acceleratorPut(this->_entries[lexYm],SE);
       
@@ -751,6 +781,51 @@ public:
       SE._missing_link = missingLink;
       acceleratorPut(this->_entries[lexDm],SE);
     }
+    if ( vertexOutputs ) {
+      int ndm1 = grid->Nd()-1;
+      if ( grid->ownsSouthPole() ) {
+	IcosahedralStencilEntry SE;
+	for(uint64_t site=0;site<cart_sites; site ++) {
+	  Coordinate Coor;
+	  grid->oCoorFromOindex(Coor,site);
+	  if( (Coor[0]==L)&&(Coor[1]==0) ) {
+	    int64_t pole_site = grid->PoleSiteForOcoor(Coor);
+	    int64_t lex       = pole_site*np+Coor[ndm1];
+	    SE._offset       = site;
+	    SE._is_local     = true;
+	    SE._polarisation = IcosahedronPatchX;  // ignored
+	    SE._adjoint      = false;              // ignored
+	    SE._missing_link = false;
+	    acceleratorPut(this->_entries[lex],SE);
+	    
+	    int64_t lex5      = pole_site*np+5; // We miss the backwards link
+	    SE._missing_link = true;
+	    acceleratorPut(this->_entries[lex5],SE);
+	  }
+	}
+      }
+      if ( grid->ownsNorthPole() ) {
+	IcosahedralStencilEntry SE;
+	for(uint64_t site=0;site<cart_sites; site ++) {
+	  Coordinate Coor;
+	  grid->oCoorFromOindex(Coor,site);
+	  if( (Coor[0]==0)&&(Coor[1]==L) ) {
+	    int64_t pole_site = grid->PoleSiteForOcoor(Coor);
+	    int64_t lex       = pole_site*np+Coor[ndm1];
+	    SE._offset       = site;
+	    SE._is_local     = true;
+	    SE._polarisation = IcosahedronPatchX;  // ignored
+	    SE._adjoint      = false;              // ignored
+	    SE._missing_link = false;
+	    acceleratorPut(this->_entries[lex],SE);
+	    
+	    int64_t lex5     = pole_site*np+5; // We miss the backwards link
+	    SE._missing_link = true;
+	    acceleratorPut(this->_entries[lex5],SE);
+	  }
+	}
+      }
+    }
   }
     /*************************************************************
      * For gauge action implementation

Grid/util/Init.cc

@@ -189,8 +189,8 @@ void GridParseLayout(char **argv,int argc,
 {
   auto mpi =std::vector<int>(Nd,1);
   auto latt=std::vector<int>(Nd,8);
-  std::cout << "Default mpi "<<mpi<<std::endl;
+
-  std::cout << "Default latt"<<latt<<std::endl;
+
   GridThread::SetMaxThreads();
 
   std::string arg;

tests/debug/Test_icosahedron.cc

@@ -68,14 +68,19 @@ public:
   //
   GridBase *VertexGrid;
   GridBase *EdgeGrid;
+
   IcosahedralStencil FaceStencil;
-  IcosahedralStencil NNStencil;
+  IcosahedralStencil NNee; // edge neighbours with edge domain
+  IcosahedralStencil NNev; // vertex neighbours but in edge domain
+  IcosahedralStencil NNvv; // vertex neighbours with vertex domain
   
   IcosahedralEdgeSupport(GridBase *_VertexGrid,GridBase *_EdgeGrid)
-    : FaceStencil (EdgeGrid), NNStencil(EdgeGrid), VertexGrid(_VertexGrid), EdgeGrid(_EdgeGrid)
+    : FaceStencil (_EdgeGrid), NNee(_EdgeGrid), NNev(_VertexGrid), NNvv(_VertexGrid), VertexGrid(_VertexGrid), EdgeGrid(_EdgeGrid)
   {
     FaceStencil.FaceStencil();
-    NNStencil.NearestNeighbourStencil(); // Vertex nearest neighbour
+    NNee.NearestNeighbourStencil(false);// Edge nearest neighbour
+    NNev.NearestNeighbourStencil(false);// Edge result, vertex neighbour
+    NNvv.NearestNeighbourStencil(true); // vertex result and neighbour
   }
 
   ////////////////////////////////////////////////////////////////////////////////////
@@ -147,9 +152,9 @@ public:
       autoView(stapleDX_v,stapleDX,AcceleratorWrite);
       autoView(stapleYD_v,stapleYD,AcceleratorWrite);
       autoView(stapleDY_v,stapleDY,AcceleratorWrite);
-      autoView(stencil_v,NNStencil,AcceleratorRead);
+      autoView(stencil_v,NNee,AcceleratorRead);
 
-      const int np = NNStencil._npoints;
+      const int np = NNee._npoints;
       
       accelerator_for(ss,EdgeGrid->oSites(),vComplex::Nsimd(),{
 
@@ -265,10 +270,83 @@ public:
       });
   }
   
+  /*
+   * Should be able to use the Vertex based stencil to do the GT, picking forward hops
+   */
+  template<class MatterField>
+  void  CovariantLaplacian(MatterField &in,MatterField &out, GaugeField &Umu)
+  {
+  }
+  void  GaugeTransform(GaugeLinkField &gt, GaugeField &Umu)
+  {
+    autoView(Umu_v,Umu,AcceleratorWrite);
+    autoView(g_v,gt,AcceleratorRead);
+    autoView(stencil_v,NNev,AcceleratorRead);
+
+    const int np = NNev._npoints;
+
+    std::cout << GridLogMessage<< "GaugeTransform via STENCIL "<<std::endl;
+
+    accelerator_for(ss,EdgeGrid->oSites(),vComplex::Nsimd(),{
+
+      const int ent_Xp = 0;
+      const int ent_Yp = 1;
+      const int ent_Dp = 2;
+      const int ent_Xm = 3;
+      const int ent_Ym = 4;
+      const int ent_Dm = 5;
+      
+      Integer lexXp = ss*np+ent_Xp;
+      Integer lexYp = ss*np+ent_Yp;
+      Integer lexDp = ss*np+ent_Dp;
+      Integer lexXm = ss*np+ent_Xm;
+      Integer lexYm = ss*np+ent_Ym;
+      Integer lexDm = ss*np+ent_Dm;
+	  
+      const int x = IcosahedronPatchX;
+      const int y = IcosahedronPatchY;
+      const int d = IcosahedronPatchDiagonal;
+
+      // Three forward links from this site
+      auto Lx = Umu_v(ss)(x);
+      auto Ly = Umu_v(ss)(y);
+      auto Ld = Umu_v(ss)(d);
+
+      uint64_t xp_idx;
+      uint64_t yp_idx;
+      uint64_t dp_idx;
+
+      auto SE = stencil_v.GetEntry(ent_Xp,ss);
+      xp_idx = SE->_offset;
+
+      SE = stencil_v.GetEntry(ent_Yp,ss);
+      yp_idx = SE->_offset;
+
+      SE = stencil_v.GetEntry(ent_Dp,ss);
+      dp_idx = SE->_offset;
+
+      auto g  = g_v(ss)();
+      auto gx = g_v(xp_idx)();
+      auto gy = g_v(yp_idx)();
+      auto gd = g_v(dp_idx)();
+
+      auto lx = Umu_v(ss)(IcosahedronPatchX);
+      auto ly = Umu_v(ss)(IcosahedronPatchY);
+      auto ld = Umu_v(ss)(IcosahedronPatchDiagonal);
+
+      lx = g*lx*adj(gx);
+      ly = g*ly*adj(gy);
+      ld = g*ld*adj(gd);
+      
+      coalescedWrite(Umu_v[ss](IcosahedronPatchX),lx);
+      coalescedWrite(Umu_v[ss](IcosahedronPatchY),ly);
+      coalescedWrite(Umu_v[ss](IcosahedronPatchDiagonal),ld);
+      });
+  }
   /*
    * This routine is slow and single threaded on CPU
    */
-  void  GaugeTransform(GaugeLinkField &gt, GaugeField &Umu)
+  void  GaugeTransformCPU(GaugeLinkField &gt, GaugeField &Umu)
   {
     assert(gt.Grid()==VertexGrid);
     assert(Umu.Grid()==EdgeGrid);
@@ -372,6 +450,20 @@ public:
       auto ly = Umu_v(site)(IcosahedronPatchY);
       auto ld = Umu_v(site)(IcosahedronPatchDiagonal);
       
+      /*
+	std::cout << "site "<<site<<std::endl;
+	std::cout << " xp_idx "<<xp_idx<<std::endl;
+	std::cout << " yp_idx "<<yp_idx<<std::endl;
+	std::cout << " dp_idx "<<dp_idx<<std::endl;
+	std::cout << " g  "<<g<<std::endl;
+	std::cout << " gx "<<gx<<std::endl;
+	std::cout << " gy "<<gy<<std::endl;
+	std::cout << " gd "<<gd<<std::endl;
+	std::cout << " lx "<<lx<<std::endl;
+	std::cout << " ly "<<ly<<std::endl;
+	std::cout << " ld "<<ld<<std::endl;
+      */
+      
       lx = g*lx*adj(gx);
       ly = g*ly*adj(gy);
       ld = g*ld*adj(gd);
@@ -405,6 +497,7 @@ int main (int argc, char ** argv)
 
   std::cout << GridLogMessage << " Created vertex grid "<<std::endl;
   LatticeIcosahedralLorentzColourMatrix Umu(&EdgeGrid);
+  LatticeIcosahedralLorentzColourMatrix Umuck(&EdgeGrid);
   LatticeComplex Phi(&VertexGrid);
   std::cout << GridLogMessage << " Created two fields "<<std::endl;
 
@@ -476,7 +569,11 @@ int main (int argc, char ** argv)
   std::cout << GridLogMessage << " Check plaquette is gauge invariant "<<std::endl;
   std::cout << GridLogMessage << "****************************************"<<std::endl;
   std::cout << GridLogMessage << " applying gauge transform"<<std::endl;
-  Support.GaugeTransform(g,Umu);
+  Umuck = Umu;
+  Support.GaugeTransform    (g,Umuck);
+  Support.GaugeTransformCPU(g,Umu);
+  Umuck = Umuck - Umu;
+  std::cout << GridLogMessage <<"Diff between reference GT and stencil GT: "<<norm2(Umuck) <<std::endl;
   std::cout << GridLogMessage << " applied gauge transform "<<std::endl;
   //  std::cout << "Umu\n"<< Umu << std::endl;
   std::cout << GridLogMessage << " recalculating plaquette "<<std::endl;