Updated

Grid/stencil/IcosahedralStencil.h

@@ -69,6 +69,7 @@ public:
 
 protected:
   GridBase *                        _grid;
+  GridBase *                        _vertexgrid;
 
 public: 
   GridBase *Grid(void) const { return _grid; }
@@ -572,11 +573,13 @@ public:
     std::cout << GridLogMessage<< " Icosahedral Stencil Geometry Test Complete"<<std::endl;
     std::cout << GridLogMessage<< "*************************************"<<std::endl;
   }
-  IcosahedralStencil(GridBase *grid) // Must be +1 or -1
+  IcosahedralStencil(GridBase *grid,GridBase *vertexgrid) 
   {
-    this->_grid    = grid;
+    this->_grid          = grid;
+    this->_vertexgrid    = vertexgrid;
     // Loop over L^2 x T x npatch and the
     assert(grid->isIcosahedral());
+    assert(grid->isIcosahedral());
   }
   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
   // VertexInputs = true  implies the neighbour has vertex support
@@ -589,14 +592,11 @@ public:
   //     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)
+  void NearestNeighbourStencil(int vertexInputs,int vertexOutputs)
   {
-    GridBase * grid = this->_grid;
-
-    int vertexInputs = grid->isIcosahedralVertex(); 
-
-    int osites  = grid->oSites();
-
+    GridBase * grid = this->_grid; // the edge grid
+    GridBase * vertexgrid = this->_vertexgrid;
+    
     uint64_t cart_sites = grid->CartesianOsites();
     uint64_t Npole_sites = grid->NorthPoleOsites();
     uint64_t Spole_sites = grid->SouthPoleOsites();
@@ -641,6 +641,7 @@ public:
 
       int Patch = Coor[nd-1];
       int HemiPatch = Patch%HemiPatches;
+      int Hemisphere= Patch/HemiPatches;
       int north = Patch/HemiPatches;
       int south = 1-north;
       int isPoleY;
@@ -681,7 +682,11 @@ public:
       int YmHemiPatch  = YmCoor[nd-1]%HemiPatches;
       int YmHemisphere = YmCoor[nd-1]/HemiPatches;
 
-      if ( vertexInputs ) {
+      int DmPatch  = DmCoor[nd-1];
+      int DmHemiPatch  = DmCoor[nd-1]%HemiPatches;
+      int DmHemisphere = DmCoor[nd-1]/HemiPatches;
+
+      if ( vertexInputs ) {// Neighbour will live on poles and peer point
 	////////////////////////////////////////////////
 	// XpCoor stencil entry; consider isPole case
 	////////////////////////////////////////////////
@@ -689,7 +694,8 @@ public:
 	SE._missing_link = false;
 	SE._offset       = grid->oIndex(XpCoor);
 	if ( isPoleX ) {
-	  SE._offset     = grid->PoleSiteForOcoor(Coor);
+	  SE._offset     = vertexgrid->PoleSiteForOcoor(Coor);
+	  //	  std::cout << site<<" setting X-Pole site "<<SE._offset<<" for coor "<<Coor<<std::endl;
 	}
 	SE._polarisation = IcosahedronPatchY;
 	SE._adjoint      = false;
@@ -701,12 +707,13 @@ public:
 	SE._missing_link = false;
 	SE._offset       = grid->oIndex(YpCoor);
 	if ( isPoleY ) {
-	  SE._offset     = grid->PoleSiteForOcoor(Coor); 
+	  SE._offset     = vertexgrid->PoleSiteForOcoor(Coor);
+	  //	  std::cout << site<<" setting Y-Pole site "<<SE._offset<<" for coor "<<Coor<<std::endl;
 	}
 	SE._polarisation = IcosahedronPatchX;
 	SE._adjoint      = false;
 	acceleratorPut(this->_entries[lexYp],SE);
-      } else { 
+      } else { // Neighbour will be a forward edge and connection may be more complicated
 	////////////////////////////////////////////////
 	// XpCoor stencil entry
 	// Store in look up table
@@ -718,7 +725,7 @@ public:
 	SE._offset       = grid->oIndex(XpCoor);
 	SE._polarisation = IcosahedronPatchY;
 	SE._adjoint      = false;
-	if ( DpHemiPatch != HemiPatch && south ) {
+	if ( DpHemiPatch != HemiPatch && south ) { // These are the sneaky redirect for edge / faces
 	  SE._offset       = grid->oIndex(DpCoor);
 	  SE._polarisation = IcosahedronPatchX;
 	  SE._adjoint      = true;
@@ -732,7 +739,7 @@ public:
 	SE._offset       = grid->oIndex(YpCoor);
 	SE._polarisation = IcosahedronPatchX;
 	SE._adjoint      = false;
-	if ( YpHemiPatch != HemiPatch && north ) {
+	if ( YpHemiPatch != HemiPatch && north ) {  // These are the sneaky redirect for edge / faces
 	  SE._offset       = grid->oIndex(DpCoor);
 	  SE._polarisation = IcosahedronPatchY;
 	  SE._adjoint      = true;
@@ -774,23 +781,34 @@ public:
     
       /////////////////////////////////////////////////////////////////////
       // for DmCoor ; never needed for staples, only for vertex diff ops
-      // no polarisation rotation
+      // No polarisation rotation.
+      // But polarisation rotation is needed for double storing.
       /////////////////////////////////////////////////////////////////////
       SE._offset       = grid->oIndex(DmCoor);
-      SE._polarisation = IcosahedronPatchDiagonal; // should ignore
+      SE._polarisation = IcosahedronPatchDiagonal; // default
+      if ( (DmHemiPatch != HemiPatch) && (DmHemisphere==Hemisphere)  && south ) {
+        SE._polarisation = IcosahedronPatchX;   // Basis rotates
+      }
+      if ( DmHemiPatch != HemiPatch && (DmHemisphere==Hemisphere) && north ) {
+        SE._polarisation = IcosahedronPatchY;   // Basis rotates
+      }
       SE._missing_link = missingLink;
       acceleratorPut(this->_entries[lexDm],SE);
     }
+
     if ( vertexOutputs ) {
       int ndm1 = grid->Nd()-1;
-      if ( grid->ownsSouthPole() ) {
+      if ( vertexgrid->ownsSouthPole() ) {
 	IcosahedralStencilEntry SE;
-	for(uint64_t site=0;site<cart_sites; site ++) {
+	for(uint64_t site=0;site<cart_sites; site ++) { // loops over volume
 	  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];
+	  vertexgrid->oCoorFromOindex(Coor,site);
+	  int north = Coor[ndm1]/HemiPatches;
+	  int south = 1-north;
+	  if( (Coor[0]==(L-1))&&(Coor[1]==0) &&south ) {
+	    int64_t pole_site = vertexgrid->PoleSiteForOcoor(Coor);
+	    int64_t lex       = pole_site*np+(Coor[ndm1]%HemiPatches);
+	    //	    std::cout << "Coor "<<Coor<<" connects to south pole_site "<<pole_site<<std::endl;
 	    SE._offset       = site;
 	    SE._is_local     = true;
 	    SE._polarisation = IcosahedronPatchX;  // ignored
@@ -804,17 +822,19 @@ public:
 	  }
 	}
       }
-      if ( grid->ownsNorthPole() ) {
+      if ( vertexgrid->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];
+	  vertexgrid->oCoorFromOindex(Coor,site);
+	  int north = Coor[ndm1]/HemiPatches;
+	  if( (Coor[0]==0)&&(Coor[1]==(L-1))&&north ) {
+	    int64_t pole_site = vertexgrid->PoleSiteForOcoor(Coor);
+	    int64_t lex       = pole_site*np+(Coor[ndm1]%HemiPatches);
+	    //	    std::cout << "Coor "<<Coor<<" connects to north pole_site "<<pole_site<<std::endl;
 	    SE._offset       = site;
 	    SE._is_local     = true;
-	    SE._polarisation = IcosahedronPatchX;  // ignored
+	    SE._polarisation = IcosahedronPatchY;  // ignored
 	    SE._adjoint      = false;              // ignored
 	    SE._missing_link = false;
 	    acceleratorPut(this->_entries[lex],SE);
@@ -986,3 +1006,4 @@ public:
    */
 };
 NAMESPACE_END(Grid);
+  

tests/debug/Test_icosahedron.cc

@@ -1,4 +1,4 @@
-    /*************************************************************************************
+/*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -33,30 +33,13 @@ using namespace Grid;
 
 const int MyNd=3;
 
-template<typename vtype> using iIcosahedralLorentzComplex      = iVector<iScalar<iScalar<vtype> >, MyNd+1 > ;
-template<typename vtype> using iIcosahedralLorentzColourMatrix = iVector<iScalar<iMatrix<vtype,Nc> >, MyNd+1 > ;
-template<typename vtype> using iIcosahedralColourMatrix        = iScalar<iScalar<iMatrix<vtype,Nc> > > ;
-
-typedef iIcosahedralLorentzComplex<Complex >  IcosahedralLorentzComplex;
-typedef iIcosahedralLorentzComplex<vComplex> vIcosahedralLorentzComplex;
-typedef Lattice<vIcosahedralLorentzComplex> LatticeIcosahedralLorentzComplex;
-
-
-typedef iIcosahedralLorentzColourMatrix<Complex >  IcosahedralLorentzColourMatrix;
-typedef iIcosahedralLorentzColourMatrix<vComplex> vIcosahedralLorentzColourMatrix;
-typedef Lattice<vIcosahedralLorentzColourMatrix>  LatticeIcosahedralLorentzColourMatrix;
-
-typedef iIcosahedralColourMatrix<Complex >  IcosahedralColourMatrix;
-typedef iIcosahedralColourMatrix<vComplex> vIcosahedralColourMatrix;
-typedef Lattice<vIcosahedralColourMatrix>  LatticeIcosahedralColourMatrix;
-
-
 class IcosahedralGimpl
 {
 public:
-  typedef LatticeIcosahedralLorentzColourMatrix  GaugeField;
-  typedef LatticeIcosahedralColourMatrix         GaugeLinkField;
-  typedef LatticeComplex                          ComplexField;
+  typedef LatticeLorentzColourMatrix  GaugeField;
+  typedef LatticeColourMatrix         GaugeLinkField;
+  typedef LatticeDoubledGaugeField    DoubledGaugeField;
+  typedef LatticeComplex              ComplexField;
 };
   
 template< class Gimpl>
@@ -66,6 +49,7 @@ public:
   typedef typename Gimpl::GaugeField GaugeField;
   typedef typename Gimpl::GaugeLinkField GaugeLinkField;
   typedef typename Gimpl::ComplexField ComplexField;
+  typedef typename Gimpl::DoubledGaugeField DoubledGaugeField;
   //
   GridBase *VertexGrid;
   GridBase *EdgeGrid;
@@ -74,18 +58,26 @@ public:
   IcosahedralStencil NNee; // edge neighbours with edge domain
   IcosahedralStencil NNev; // vertex neighbours but in edge domain
   IcosahedralStencil NNvv; // vertex neighbours with vertex domain
+  IcosahedralStencil NNve; // edge neighbours with vertex domain
   
   IcosahedralSupport(GridBase *_VertexGrid,GridBase *_EdgeGrid)
-    : FaceStencil (_EdgeGrid), NNee(_EdgeGrid), NNev(_VertexGrid), NNvv(_VertexGrid), VertexGrid(_VertexGrid), EdgeGrid(_EdgeGrid)
+    : FaceStencil (_EdgeGrid,_VertexGrid),
+      NNee(_EdgeGrid,_VertexGrid),
+      NNev(_EdgeGrid,_VertexGrid),
+      NNve(_EdgeGrid,_VertexGrid),
+      NNvv(_EdgeGrid,_VertexGrid),
+      VertexGrid(_VertexGrid), EdgeGrid(_EdgeGrid)
   {
     FaceStencil.FaceStencil();
-    NNee.NearestNeighbourStencil(false);// Edge nearest neighbour
-    NNev.NearestNeighbourStencil(false);// Edge result, vertex neighbour
-    NNvv.NearestNeighbourStencil(true); // vertex result and neighbour
+    std::cout << "NNee"<<std::endl;
+    // true/false is "vertexInput, VertexOutput"
+    NNee.NearestNeighbourStencil(false,false);// edge input + output      ; used by face stencil
+    NNev.NearestNeighbourStencil(true,false); // vertex input, edge ouput ; used by gauge transform
+    NNvv.NearestNeighbourStencil(true,true);  // vertex input + output    ; used by Laplacian
+    NNve.NearestNeighbourStencil(false,true); // edge input, vertex output; used by double store
   }
 
   ////////////////////////////////////////////////////////////////////////////////////
-  // Fixme: will need a version of "Gimpl" and a wrapper class following "WilsonLoops" style.
   // Gauge Link field GT is the gauge transform and lives on the VERTEX field
   ////////////////////////////////////////////////////////////////////////////////////
   void ForwardTriangles(GaugeField &Umu,LatticeComplex &plaq1,LatticeComplex &plaq2)
@@ -253,9 +245,6 @@ public:
       });
   }
   
-  /*
-   * Should be able to use the Vertex based stencil to do the GT, picking forward hops
-   */
   template<class MatterField>
   void Laplacian(MatterField &in,MatterField &out)
   {
@@ -312,8 +301,163 @@ public:
       coalescedWrite(out_v[ss](),o);
       });
   }
+
+
+  template<class MatterField>
+  void CovariantLaplacian(MatterField &in,MatterField &out,DoubledGaugeField &Uds)
+  {
+    autoView(out_v,out,AcceleratorWrite);
+    autoView(in_v,in,AcceleratorRead);
+    autoView(U_v,Uds,AcceleratorRead);
+    autoView(stencil_v,NNvv,AcceleratorRead);
+
+    const int np = NNvv._npoints;
+
+    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;
     
+    accelerator_for(ss,VertexGrid->oSites(),vComplex::Nsimd(),{
+	  
+      auto SE = stencil_v.GetEntry(ent_Xp,ss);
+      uint64_t xp_idx = SE->_offset;
+
+      SE = stencil_v.GetEntry(ent_Yp,ss);
+      uint64_t yp_idx = SE->_offset;
+
+      SE = stencil_v.GetEntry(ent_Dp,ss);
+      uint64_t dp_idx = SE->_offset;
+
+      SE = stencil_v.GetEntry(ent_Xm,ss);
+      uint64_t xm_idx = SE->_offset;
+
+      SE = stencil_v.GetEntry(ent_Ym,ss);
+      uint64_t ym_idx = SE->_offset;
+
+      SE = stencil_v.GetEntry(ent_Dm,ss);
+      uint64_t dm_idx = SE->_offset;
+      int missingLink = SE->_missing_link;
+      
+      auto i    = in_v(ss)();
+      auto inxp = in_v(xp_idx)();
+      auto inyp = in_v(yp_idx)();
+      auto indp = in_v(dp_idx)();
+      auto inxm = in_v(xm_idx)();
+      auto inym = in_v(ym_idx)();
+      auto indm = in_v(dm_idx)();
+
+      inxp = U_v(ss)(0)*inxp;
+      inyp = U_v(ss)(1)*inyp;
+      indp = U_v(ss)(2)*indp;
+      inxm = U_v(ss)(3)*inxm;
+      inym = U_v(ss)(4)*inym;
+      
+      auto o    = i;
+
+      if ( missingLink ) {
+	o = (1.0/5.0)*(inxp+inyp+indp+inxm+inym)-i;
+      } else {
+	indm = U_v(ss)(5)*indm;
+	o = (1.0/6.0)*(inxp+inyp+indp+inxm+inym+indm)-i;
+      }
+      
+      coalescedWrite(out_v[ss](),o);
+      });
+  }
   
+  void DoubleStore(GaugeField &U,DoubledGaugeField &Uds)
+  {
+    assert(U.Grid()==EdgeGrid);
+    assert(Uds.Grid()==VertexGrid);
+    autoView(Uds_v,Uds,AcceleratorWrite);
+    autoView(U_v,U,AcceleratorRead);
+    autoView(stencil_v,NNvv,AcceleratorRead);
+
+    // Vertex result
+    // Edge valued input
+    // Might need an extra case (?)
+    const int np = NNvv._npoints;
+
+    const int ent_Xm = 3;
+    const int ent_Ym = 4;
+    const int ent_Dm = 5;
+    
+    accelerator_for(ss,VertexGrid->CartesianOsites(),vComplex::Nsimd(),{
+	
+      // Three local links
+      {
+	auto Lx = U_v(ss)(IcosahedronPatchX);
+	auto Ly = U_v(ss)(IcosahedronPatchY);
+	auto Ld = U_v(ss)(IcosahedronPatchDiagonal);
+	coalescedWrite(Uds_v[ss](IcosahedronPatchX),Lx);
+	coalescedWrite(Uds_v[ss](IcosahedronPatchY),Ly);
+	coalescedWrite(Uds_v[ss](IcosahedronPatchDiagonal),Ld);
+      }	
+      // Three backwards links
+      {
+	auto SE = stencil_v.GetEntry(ent_Xm,ss);
+	auto pol = SE->_polarisation;
+	auto s   = SE->_offset;
+	int pol1 = IcosahedronPatchX;
+	//
+	// Xm link is given diagonal unless HemiPatch changes in Northern hemisphere
+	// But here for double storing we need the reverse link so
+	// return either the Xplus or Diag plus direction
+	//
+	if ( pol != IcosahedronPatchDiagonal ) {
+	  pol1 = IcosahedronPatchDiagonal;
+	}
+	auto Lx_at_xm = U_v(s)(pol1);
+	Lx_at_xm = adj(Lx_at_xm);
+	coalescedWrite(Uds_v[ss](3),Lx_at_xm );
+      }
+      {
+	auto SE   = stencil_v.GetEntry(ent_Ym,ss);
+	auto pol  = SE->_polarisation;
+	auto s    = SE->_offset;
+	//
+	// Ym link is given diagonal unless HemiPatch changes in the Southern hemisphere
+	// But here for double storing we need the reverse link so
+	// return either the Yplus or Diag plus direction
+	// 
+	int pol1 = IcosahedronPatchY;
+	if ( pol != IcosahedronPatchDiagonal ) {
+	  pol1 = IcosahedronPatchDiagonal;
+	}
+	auto Ly_at_ym = U_v(s)(pol1);
+	Ly_at_ym = adj(Ly_at_ym);
+	coalescedWrite(Uds_v[ss](4),Ly_at_ym );
+      }
+      int missingLink;
+      {      
+	auto SE = stencil_v.GetEntry(ent_Dm,ss);
+	auto s  = SE->_offset;
+	// Dm link given the correct value for this use case
+	auto pol= SE->_polarisation;
+	missingLink = SE->_missing_link;
+	if ( ! missingLink ) {
+	  auto Ld_at_dm  =  U_v(s)(pol);
+	  Ld_at_dm = adj(Ld_at_dm);
+	  coalescedWrite(Uds_v[ss](5),Ld_at_dm );
+	}
+      }
+    });
+    auto pole_sites = VertexGrid->oSites() - VertexGrid->CartesianOsites();
+    auto pole_offset= VertexGrid->CartesianOsites();
+    
+    accelerator_for(ss,pole_sites,vComplex::Nsimd(),{
+	for(int p=0;p<HemiPatches;p++){ // Neighbours of each pole
+	auto SE = stencil_v.GetEntry(p,pole_offset+ss);
+	auto s  = SE->_offset;
+	auto pol= SE->_polarisation;
+	auto Link =  adj(U_v(s)(pol));
+	coalescedWrite(Uds_v[pole_offset+ss](p),Link);
+      }
+    });
+  }
   void  GaugeTransform(GaugeLinkField &gt, GaugeField &Umu)
   {
     autoView(Umu_v,Umu,AcceleratorWrite);
@@ -321,8 +465,6 @@ public:
     autoView(stencil_v,NNev,AcceleratorRead);
 
     const int np = NNev._npoints;
-
-    std::cout << GridLogMessage<< "GaugeTransform via STENCIL "<<std::endl;
     
     const int ent_Xp = 0;
     const int ent_Yp = 1;
@@ -356,128 +498,14 @@ public:
       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  GaugeTransformCPU(GaugeLinkField &gt, GaugeField &Umu)
-  {
-    assert(gt.Grid()==VertexGrid);
-    assert(Umu.Grid()==EdgeGrid);
-    assert(VertexGrid->isIcosahedralVertex());
-    assert(EdgeGrid->isIcosahedralEdge());
-    
-    GridBase * vgrid = VertexGrid;
-    GridBase * grid  = EdgeGrid;
-
-    int osites  = grid->oSites();
-
-    uint64_t cart_sites = grid->CartesianOsites();
-    uint64_t Npole_sites = grid->NorthPoleOsites();
-    uint64_t Spole_sites = grid->SouthPoleOsites();
-    Coordinate pcoor = grid->ThisProcessorCoor();
-    Coordinate pgrid = grid->ProcessorGrid();
-
-    autoView(g_v,gt,CpuRead);
-    autoView(Umu_v,Umu,CpuWrite);
-    for(uint64_t site=0;site<cart_sites; site ++) {
-
-      Coordinate Coor;
-      Coordinate NbrCoor;
-
-      int nd = grid->Nd();
-      int L  = grid->LocalDimensions()[0];
-	
-      ////////////////////////////////////////////////
-      // Outer index of neighbour Offset calculation
-      ////////////////////////////////////////////////
-      grid->oCoorFromOindex(Coor,site);
-      NbrCoor = Coor;
-      assert( grid->LocalDimensions()[1]==grid->LocalDimensions()[0]);
-      assert( grid->_simd_layout[0]==1); // Cannot vectorise in these dims
-      assert( grid->_simd_layout[1]==1);
-      assert( grid->_processors[0]==1); // Cannot mpi distribute in these dims
-      assert( grid->_processors[1]==1);
-
-      int Patch = Coor[nd-1];
-      int HemiPatch = Patch%HemiPatches;
-      int north = Patch/HemiPatches;
-      int south = 1-north;
-      int isPoleY;
-      int isPoleX;
-      
-      assert(Patch<IcosahedralPatches);
-      assert((north==1)||(south==1));
-
-      Coordinate XpCoor;
-      Coordinate YpCoor;
-      Coordinate DpCoor;
-
-      FaceStencil.GetNbrForPlusDiagonal(grid,Coor,DpCoor);
-      FaceStencil.GetNbrForPlusX(grid,Coor,XpCoor,isPoleX);
-      FaceStencil.GetNbrForPlusY(grid,Coor,YpCoor,isPoleY);
-
-      int XpHemiPatch  = XpCoor[nd-1]%HemiPatches;
-      int XpHemisphere = XpCoor[nd-1]/HemiPatches;
-
-      int DpPatch  = DpCoor[nd-1];
-      int DpHemiPatch  = DpCoor[nd-1]%HemiPatches;
-      int DpHemisphere = DpCoor[nd-1]/HemiPatches;
-
-      // Work out the pole_osite
-      Coordinate rdims;
-      Coordinate ocoor;
-      int64_t pole_osite;
-      int Ndm1 = grid->Nd()-1;
-      for(int d=2;d<Ndm1;d++){
-	int dd=d-2;
-	rdims.push_back(grid->_rdimensions[d]);
-	ocoor.push_back(Coor[d]%grid->_rdimensions[d]);
-      }
-      Lexicographic::IndexFromCoor(ocoor,pole_osite,rdims);
-      
-      uint64_t xp_idx;
-      uint64_t yp_idx;
-      uint64_t dp_idx;
-      if ( isPoleX ) {
-	assert(vgrid->ownsSouthPole());
-	xp_idx = pole_osite + vgrid->SouthPoleOsite();
-      } else {
-	xp_idx = grid->oIndex(XpCoor);
-      }
-      if ( isPoleY ) {
-	assert(vgrid->ownsNorthPole());
-	yp_idx = pole_osite + vgrid->NorthPoleOsite();
-      } else {
-	yp_idx = grid->oIndex(YpCoor);
-      }
-      dp_idx = grid->oIndex(DpCoor);
-
-      auto g  = g_v(site)();
-      auto gx = g_v(xp_idx)();
-      auto gy = g_v(yp_idx)();
-      auto gd = g_v(dp_idx)();
-
-      auto lx = Umu_v(site)(IcosahedronPatchX);
-      auto ly = Umu_v(site)(IcosahedronPatchY);
-      auto ld = Umu_v(site)(IcosahedronPatchDiagonal);
-       
-      lx = g*lx*adj(gx);
-      ly = g*ly*adj(gy);
-      ld = g*ld*adj(gd);
-
-      coalescedWrite(Umu_v[site](IcosahedronPatchX),lx);
-      coalescedWrite(Umu_v[site](IcosahedronPatchY),ly);
-      coalescedWrite(Umu_v[site](IcosahedronPatchDiagonal),ld);
-    };
-  }
 
 };
-   */
 
 int main (int argc, char ** argv)
 {
@@ -499,8 +527,8 @@ int main (int argc, char ** argv)
   GridCartesianCrossIcosahedron VertexGrid(latt_size,simd_layout,mpi_layout,IcosahedralVertices);
 
   std::cout << GridLogMessage << " Created vertex grid "<<std::endl;
-  LatticeIcosahedralLorentzColourMatrix Umu(&EdgeGrid);
-  LatticeIcosahedralLorentzColourMatrix Umuck(&EdgeGrid);
+  LatticeLorentzColourMatrix Umu(&EdgeGrid);
+  LatticeLorentzColourMatrix Umuck(&EdgeGrid);
   LatticeComplex Phi(&VertexGrid);
   std::cout << GridLogMessage << " Created two fields "<<std::endl;
 
@@ -557,32 +585,28 @@ int main (int argc, char ** argv)
   std::vector<int> seeds({1,2,3,4});
   GridParallelRNG  vRNG(&EdgeGrid);   vRNG.SeedFixedIntegers(seeds);
   
-  //  SU<Nc>::LieRandomize(vRNG,g);
-  LatticeIcosahedralColourMatrix g(&VertexGrid);
+  LatticeColourMatrix g(&VertexGrid);
   LatticeReal gr(&VertexGrid);
   LatticeComplex gc(&VertexGrid);
-  gr = 1.0;
+  //  gr = Zero();
   gaussian(vRNG,gr);
   Complex ci(0.0,1.0);
   gc = toComplex(gr);
   g=one;
   g = g * exp(ci*gc);
-  
+ 
   std::cout << GridLogMessage << "****************************************"<<std::endl;
   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);
   std::cout << GridLogMessage << " applied gauge transform "<<std::endl;
-  //  std::cout << "Umu\n"<< Umu << std::endl;
+
   std::cout << GridLogMessage << " recalculating plaquette "<<std::endl;
   Support.ForwardTriangles(Umu,plaq1,plaq2);
   std::cout << GridLogMessage << " plaq1 "<< norm2(plaq1)<<std::endl;
   std::cout << GridLogMessage << " plaq2 "<< norm2(plaq2)<<std::endl;
 
-  //  std::cout << " plaq1 "<< plaq1<<std::endl;
-  //  std::cout << " plaq2 "<< plaq2<<std::endl;
-
   std::cout << GridLogMessage << " plaq1 err "<< norm2(plaq1-plaq_ref)<<std::endl;
   std::cout << GridLogMessage << " plaq2 err "<< norm2(plaq2-plaq_ref)<<std::endl;
 
@@ -629,17 +653,41 @@ int main (int argc, char ** argv)
   std::cout << GridLogMessage << " trace Y*StapleDX "<<norm2(trace(linkY * stapleDX))<<std::endl;
   std::cout << GridLogMessage << " err " << norm2(trace(linkY * stapleDX)-plaq_ref)<<std::endl;
 
-  //  std::cout << " D " << linkD<<std::endl;
-  //  std::cout << " X " << linkX<<std::endl;
-  //  std::cout << " Y " << linkY<<std::endl;
-  //  std::cout << " DXY\n " << closure(linkD * stapleYX) <<std::endl;
-  //  std::cout << " YXD\n " << closure(linkY * stapleXD) <<std::endl;
-  
   std::cout << GridLogMessage<< "Calling Laplacian" <<std::endl;
-  LatticeComplex in(&VertexGrid);
-  LatticeComplex out(&VertexGrid);
+  LatticeColourVector in(&VertexGrid);
+  LatticeColourVector out(&VertexGrid);
+  LatticeColourVector gout(&VertexGrid);
+  LatticeColourVector gin(&VertexGrid);
   gaussian(vRNG,in);
   Support.Laplacian(in,out);
-  
+
+  std::cout << GridLogMessage<< "Calling double storing gauge field" <<std::endl;
+  LatticeDoubledGaugeField Uds(&VertexGrid);
+  Support.DoubleStore(Umu,Uds);
+
+  Support.CovariantLaplacian(in,out,Uds);
+
+  auto ip = innerProduct(out, in);
+  std::cout << GridLogMessage<< "Applied covariant laplacian !" <<std::endl;
+  /*
+   * CovariantLaplacian testing -- check the laplacian is gauge invariant
+   *
+   * D[U_gt](gF) = g D[U] g^dag g F = g D[U] F
+   */
+  gout = g*out;
+  gin  = g*in;
+
+  Support.GaugeTransform(g,Umu);
+  Support.DoubleStore(Umu,Uds);
+  Support.CovariantLaplacian(gin,out,Uds);
+  std::cout << GridLogMessage<< "Applied gauge transformed covariant laplacian to transformed vector !" <<std::endl;
+  auto ipgt = innerProduct(out, gin);
+
+  std::cout << "Testing D[U_gt](gF) = g D[U] F : defect is "<<norm2(out-gout)<<std::endl;
+
+  ip = ip - ipgt;
+  std::cout << "Testing F D[U](F) = (gF) D[U_gt] gF : defect is "<<ip<<std::endl;
+
   Grid_finalize();
 }
+