LePage works, trying Naik

Grid/qcd/smearing/HISQSmearing.h

@@ -56,7 +56,7 @@ void appendShift(std::vector<Coordinate>& shifts, int dir, Args... args) {
 /*!  @brief figure out the stencil index from mu and nu */
 inline int stencilIndex(int mu, int nu) {
     // Nshifts depends on how you built the stencil
-    int Nshifts = 5;
+    int Nshifts = 6;
     return Nshifts*nu + Nd*Nshifts*mu;
 }
 
@@ -128,8 +128,8 @@ public:
         GF Ughost_5linkA(Ughost.Grid());
         GF Ughost_5linkB(Ughost.Grid());
 
-        // Create 3-link stencil. We allow mu==nu just to make the indexing easier.
+        // mu-nu plane stencil. We allow mu==nu to make indexing the stencil easier,
-        // Shifts with mu==nu will not be used. 
+        // but these entries will not be used. 
         std::vector<Coordinate> shifts;
         for(int mu=0;mu<Nd;mu++)
         for(int nu=0;nu<Nd;nu++) {
@@ -138,6 +138,7 @@ public:
             appendShift(shifts,NO_SHIFT);
             appendShift(shifts,mu,Back(nu));
             appendShift(shifts,Back(nu));
+            appendShift(shifts,Back(mu));
         }
 
         // A GeneralLocalStencil has two indices: a site and stencil index 
@@ -164,9 +165,10 @@ public:
             // We infer some types that will be needed in the calculation.
             typedef decltype(gStencil.GetEntry(0,0)) stencilElement;
             typedef decltype(coalescedReadGeneralPermute(U_v[0](0),gStencil.GetEntry(0,0)->_permute,Nd)) U3matrix;
-            stencilElement SE0, SE1, SE2, SE3, SE4;
+            stencilElement SE0, SE1, SE2, SE3, SE4, SE5;
             U3matrix U0, U1, U2, U3, U4, U5, W;
 
+
             for(int site=0;site<U_v.size();site++){ // ----------- 3-link
                 for(int nu=0;nu<Nd;nu++) {
                     if(nu==mu) continue;
@@ -179,6 +181,7 @@ public:
                     SE2 = gStencil.GetEntry(s+2,site); int x           = SE2->_offset;
                     SE3 = gStencil.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset;
                     SE4 = gStencil.GetEntry(s+4,site); int x_m_nu      = SE4->_offset;
+                    SE5 = gStencil.GetEntry(s+5,site); int x_m_mu      = SE5->_offset;
 
                     // When you're deciding whether to take an adjoint, the question is: how is the
                     // stored link oriented compared to the one you want? If I imagine myself travelling
@@ -194,9 +197,17 @@ public:
                     //  "left"          "right"
                     W = U2*U1*adj(U0) + adj(U5)*U4*U3;
 
+                    // Save 3-link construct for later and add to smeared field.
                     U_3link_v[site](nu) = W;
-
                     U_fat_v[site](mu)   = U_fat_v[site](mu) + lt.c_3*W;
+
+                    U0 = coalescedReadGeneralPermute(U_v[x_m_mu](mu),SE5->_permute,Nd);
+                    U1 = coalescedReadGeneralPermute(U_v[x     ](mu),SE2->_permute,Nd);
+                    U2 = coalescedReadGeneralPermute(U_v[x_p_mu](mu),SE0->_permute,Nd);
+                    W  = U0*U1*U2;
+
+                    // Add Naik term to smeared field.
+                    U_fat_v[site](mu) = U_fat_v[site](mu) + lt.c_naik*W;
                 }
             }
 
@@ -317,29 +328,4 @@ public:
 };
 
 
-/*!  @brief create long links from link variables. */
-template<class GF>
-class Smear_HISQ_Naik {
-
-private:
-    GridCartesian* const _grid;
-
-public:
-
-    // Eventually this will take, e.g., coefficients as argument 
-    Smear_HISQ_Naik(GridCartesian* grid) : _grid(grid) {
-        assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3");
-        assert(Nd == 4 && "HISQ smearing only defined for Nd==4");
-    }
-
-    ~Smear_HISQ_Naik() {}
-
-//    void smear(GF& u_smr, const GF& U) const {
-//    };
-
-//    void derivative(const GaugeField& Gauge) const {
-//    };
-};
-
-
 NAMESPACE_END(Grid);

tests/smearing/Test_fatLinks.cc

@@ -64,7 +64,7 @@ int main (int argc, char** argv) {
     int Nt = 4;
     Coordinate latt_size(Nd,0); latt_size[0]=Ns; latt_size[1]=Ns; latt_size[2]=Ns; latt_size[3]=Nt;
     std::string conf_in  = "nersc.l8t4b3360";
-    std::string conf_out = "nersc.l8t4b3360.357link";
+    std::string conf_out = "nersc.l8t4b3360.357lplink";
     int threads          = GridThread::GetThreads();
 
     typedef LatticeGaugeFieldD LGF;
@@ -92,7 +92,7 @@ int main (int argc, char** argv) {
     NerscIO::readConfiguration(Umu, header, conf_in);
 
     // Smear Umu and store result in U_smr
-    Smear_HISQ_fat<PeriodicGimplD> hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,-1/8.);
+    Smear_HISQ_fat<PeriodicGimplD> hisq_fat(&GRID,1/8.,-1/24.,1/16.,1/64.,1/384.,-1/8.);
     hisq_fat.smear(U_smr,Umu);
 
     NerscIO::writeConfiguration(U_smr,conf_out,"HISQ");
@@ -102,7 +102,7 @@ int main (int argc, char** argv) {
     Smear_HISQ_fat<PeriodicGimplD> hisq_fat_Cstyle(&GRID,path_coeff);
 
     // Make sure result doesn't change w.r.t. a trusted lattice
-    NerscIO::readConfiguration(Umu, header, "nersc.l8t4b3360.357link.control");
+    NerscIO::readConfiguration(Umu, header, "nersc.l8t4b3360.357lplink.control");
     LGF diff(&GRID);
     diff = Umu-U_smr;
     auto absDiff = norm2(diff)/norm2(Umu);