Clean up

tests/debug/Test_icosahedron.cc

@@ -41,6 +41,7 @@ typedef iIcosahedralLorentzComplex<Complex >  IcosahedralLorentzComplex;
 typedef iIcosahedralLorentzComplex<vComplex> vIcosahedralLorentzComplex;
 typedef Lattice<vIcosahedralLorentzComplex> LatticeIcosahedralLorentzComplex;
 
+
 typedef iIcosahedralLorentzColourMatrix<Complex >  IcosahedralLorentzColourMatrix;
 typedef iIcosahedralLorentzColourMatrix<vComplex> vIcosahedralLorentzColourMatrix;
 typedef Lattice<vIcosahedralLorentzColourMatrix>  LatticeIcosahedralLorentzColourMatrix;
@@ -96,14 +97,10 @@ public:
       autoView(stencil_v,FaceStencil,AcceleratorRead);
 
       accelerator_for(ss,EdgeGrid->oSites(),vComplex::Nsimd(),{
-
-	const int x = IcosahedronPatchX;
-	const int y = IcosahedronPatchY;
-	const int d = IcosahedronPatchDiagonal;
 	
-	auto Lx = Umu_v(ss)(x);
-	auto Ly = Umu_v(ss)(y);
-	auto Ld = Umu_v(ss)(d);
+	auto Lx = Umu_v(ss)(IcosahedronPatchX);
+	auto Ly = Umu_v(ss)(IcosahedronPatchY);
+	auto Ld = Umu_v(ss)(IcosahedronPatchDiagonal);
 
 	// for trace [ U_x(z) U_y(z+\hat x) adj(U_d(z)) ]
 	{
@@ -118,8 +115,6 @@ public:
 	  coalescedWrite(plaq1_v[ss](),trace(Lx*L1*adj(Ld) ) );
 	}
 
-	// This was wrong after GT
-	// Could be EITHER the GT or the the plaq / stencil
 	// for trace [  U_y(z) U_x(z+\hat y) adj(U_d(z))   ]
 	{
 	  auto SE2 = stencil_v.GetEntry(1,ss);
@@ -130,7 +125,6 @@ public:
 	  if(doAdj)
 	    L2 = adj(L2);
 	  coalescedWrite(plaq2_v[ss](),trace(Ly*L2*adj(Ld) ) );
-	  //	  std::cout << "site "<< ss<<" plaq "<< plaq2_v[ss] << " doAdj "<< (int) doAdj<<" pol "<<(int) pol <<std::endl;
 	}
       });
     }
@@ -156,28 +150,17 @@ public:
 
       const int np = NNee._npoints;
       
+      const int ent_Xp = 0;
+      const int ent_Yp = 1;
+      const int ent_Xm = 3;
+      const int ent_Ym = 4;
+
       accelerator_for(ss,EdgeGrid->oSites(),vComplex::Nsimd(),{
 
-	const int ent_Xp = 0;
-	const int ent_Yp = 1;
-	const int ent_Xm = 3;
-	const int ent_Ym = 4;
-
-	Integer lexXp = ss*np+ent_Xp;
-	Integer lexYp = ss*np+ent_Yp;
-	Integer lexXm = ss*np+ent_Xm;
-	Integer lexYm = ss*np+ent_Ym;
-	//	Integer lexDp = ss*np+2; // Not touched by staples.
-	//	Integer lexDm = ss*np+5;
-	  
-	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);
+	auto Lx = Umu_v(ss)(IcosahedronPatchX);
+	auto Ly = Umu_v(ss)(IcosahedronPatchY);
+	auto Ld = Umu_v(ss)(IcosahedronPatchDiagonal);
 
 	///////////////////////////////////////////////////////////////////
 	// Terms for the staple orthog to PlusDiagonal
@@ -274,9 +257,63 @@ 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 Laplacian(MatterField &in,MatterField &out)
   {
+    autoView(out_v,out,AcceleratorWrite);
+    autoView(in_v,in,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)();
+
+      auto o    = out_v(ss)();
+
+      if ( missingLink ) {
+	o = (1.0/5.0)*(inxp+inyp+indp+inxm+inym)-o;
+      } else {
+	o = (1.0/6.0)*(inxp+inyp+indp+inxm+inym+indm)-o;
+      }
+      
+      coalescedWrite(out_v[ss](),o);
+      });
   }
+    
+  
   void  GaugeTransform(GaugeLinkField &gt, GaugeField &Umu)
   {
     autoView(Umu_v,Umu,AcceleratorWrite);
@@ -286,44 +323,26 @@ public:
     const int np = NNev._npoints;
 
     std::cout << GridLogMessage<< "GaugeTransform via STENCIL "<<std::endl;
+    
+    const int ent_Xp = 0;
+    const int ent_Yp = 1;
+    const int ent_Dp = 2;
 
     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 Lx = Umu_v(ss)(IcosahedronPatchX);
+      auto Ly = Umu_v(ss)(IcosahedronPatchY);
+      auto Ld = Umu_v(ss)(IcosahedronPatchDiagonal);
 
       auto SE = stencil_v.GetEntry(ent_Xp,ss);
-      xp_idx = SE->_offset;
+      uint64_t xp_idx = SE->_offset;
 
       SE = stencil_v.GetEntry(ent_Yp,ss);
-      yp_idx = SE->_offset;
+      uint64_t yp_idx = SE->_offset;
 
       SE = stencil_v.GetEntry(ent_Dp,ss);
-      dp_idx = SE->_offset;
+      uint64_t dp_idx = SE->_offset;
 
       auto g  = g_v(ss)();
       auto gx = g_v(xp_idx)();
@@ -345,7 +364,6 @@ public:
   }
   /*
    * This routine is slow and single threaded on CPU
-   */
   void  GaugeTransformCPU(GaugeLinkField &gt, GaugeField &Umu)
   {
     assert(gt.Grid()==VertexGrid);
@@ -363,9 +381,7 @@ public:
     uint64_t Spole_sites = grid->SouthPoleOsites();
     Coordinate pcoor = grid->ThisProcessorCoor();
     Coordinate pgrid = grid->ProcessorGrid();
-    /*
-     * resize the stencil entries array and set npoints
-     */
+
     autoView(g_v,gt,CpuRead);
     autoView(Umu_v,Umu,CpuWrite);
     for(uint64_t site=0;site<cart_sites; site ++) {
@@ -449,21 +465,7 @@ public:
       auto lx = Umu_v(site)(IcosahedronPatchX);
       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);
@@ -475,6 +477,7 @@ public:
   }
 
 };
+   */
 
 int main (int argc, char ** argv)
 {
@@ -569,11 +572,7 @@ 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;
-  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;
+  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;
@@ -586,7 +585,6 @@ int main (int argc, char ** argv)
 
   std::cout << GridLogMessage << " plaq1 err "<< norm2(plaq1-plaq_ref)<<std::endl;
   std::cout << GridLogMessage << " plaq2 err "<< norm2(plaq2-plaq_ref)<<std::endl;
-  
 
   typedef IcosahedralGimpl::GaugeLinkField GaugeLinkField;
   typedef IcosahedralGimpl::GaugeField     GaugeField;
@@ -636,8 +634,12 @@ int main (int argc, char ** argv)
   //  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);
+  gaussian(vRNG,in);
+  Support.Laplacian(in,out);
   
   Grid_finalize();
 }