5-link first attempt

Grid/qcd/smearing/HISQSmearing.h

@@ -61,6 +61,14 @@ inline int Back(const int dir) {
 }
 
 
+/*!  @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;
+    return Nshifts*nu + Nd*Nshifts*mu;
+}
+
+
 /*!  @brief shift one unit in direction dir */
 template<typename... Args>
 void generalShift(Coordinate& shift, int dir) {
@@ -135,6 +143,7 @@ public:
         : _grid(grid), 
           _linkTreatment(c1,cnaik,c3,c5,c7,clp) {
         assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3");
+        assert(Nd == 4 && "HISQ smearing only defined for Nd==4");
     }
 
     // Allow to pass a pointer to a C-style, double array for MILC convenience
@@ -142,6 +151,7 @@ public:
         : _grid(grid), 
           _linkTreatment(coeff[0],coeff[1],coeff[2],coeff[3],coeff[4],coeff[5]) {
         assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3");
+        assert(Nd == 4 && "HISQ smearing only defined for Nd==4");
     }
 
     ~Smear_HISQ_fat() {}
@@ -150,25 +160,20 @@ public:
 
         SmearingParameters lt = this->_linkTreatment;
 
-        // Create a padded cell of extra padding depth=1
+        // Create a padded cell of extra padding depth=1 and fill the padding.
         int depth = 1;
         PaddedCell Ghost(depth,this->_grid);
         LGF Ughost = Ghost.Exchange(u_thin);
 
-        // Array for <tr U_mu_nu>(x)
+        // This is where auxiliary N-link fields and the final smear will be stored. 
-        GridBase *GhostGrid = Ughost.Grid();
-        LatticeComplex gplaq(GhostGrid); 
-
-        // This is where the 3-link constructs will be stored
         LGF Ughost_fat(Ughost.Grid());
+        LGF Ughost_3link(Ughost.Grid());
 
-        // Create 3-link stencil. Writing your own stencil, you're hard-coding the 
+        // Create 3-link stencil. We allow mu==nu just to make the indexing easier.
-        // periodic BCs, so you don't need the policy-based stuff, at least for now.
+        // Shifts with mu==nu will not be used. 
-        // Loop over all orientations, i.e. demand mu != nu.
         std::vector<Coordinate> shifts;
         for(int mu=0;mu<Nd;mu++)
         for(int nu=0;nu<Nd;nu++) {
-            if(mu==nu) continue;
             appendShift(shifts,mu);
             appendShift(shifts,nu);
             appendShift(shifts,NO_SHIFT);
@@ -176,52 +181,86 @@ public:
             appendShift(shifts,Back(nu));
         }
 
-        GeneralLocalStencil gStencil(GhostGrid,shifts);
+        // A GeneralLocalStencil has two indices: a site and stencil index 
+        GeneralLocalStencil gStencil(Ughost.Grid(),shifts);
 
+        // This is where contributions from the smearing get added together
         Ughost_fat=Zero();
 
-        // Create the accessors, here U_v and U_fat_v 
+        // Create the accessors
-        autoView(U_v    , Ughost    , CpuRead);
+        autoView(U_v      , Ughost      , CpuRead);
-        autoView(U_fat_v, Ughost_fat, CpuWrite);
+        autoView(U_fat_v  , Ughost_fat  , CpuWrite);
+        autoView(U_3link_v, Ughost_3link, CpuWrite);
 
-        // This is a loop over local sites.
+        for(int mu=0;mu<Nd;mu++) {
-        for(int ss=0;ss<U_v.size();ss++){
 
-            // This is the stencil index. It increases as we make our way through the spacetime sites,
+            Ughost_3link=Zero();
-            // plaquette orientations, and as we travel around a plaquette.
-            int s=0;
 
-            for(int mu=0;mu<Nd;mu++)
+            // 3-link
-            for(int nu=0;nu<Nd;nu++) {
+            for(int site=0;site<U_v.size();site++){
+                for(int nu=0;nu<Nd;nu++) {
+                    if(nu==mu) continue;
+                    int s = stencilIndex(mu,nu);
 
-                if(mu==nu) continue;
+                    // The stencil gives us support points in the mu-nu plane that we will use to
+                    // grab the links we need.
+                    auto SE0 = gStencil.GetEntry(s+0,site); int x_p_mu      = SE0->_offset;
+                    auto SE1 = gStencil.GetEntry(s+1,site); int x_p_nu      = SE1->_offset;
+                    auto SE2 = gStencil.GetEntry(s+2,site); int x           = SE2->_offset;
+                    auto SE3 = gStencil.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset;
+                    auto SE4 = gStencil.GetEntry(s+4,site); int x_m_nu      = SE4->_offset;
 
-                auto SE0 = gStencil.GetEntry(s+0,ss); int x_p_mu      = SE0->_offset;
+                    // When you're deciding whether to take an adjoint, the question is: how is the
-                auto SE1 = gStencil.GetEntry(s+1,ss); int x_p_nu      = SE1->_offset;
+                    // stored link oriented compared to the one you want? If I imagine myself travelling
-                auto SE2 = gStencil.GetEntry(s+2,ss); int x           = SE2->_offset;
+                    // with the to-be-updated link, I have two possible, alternative 3-link paths I can
-                auto SE3 = gStencil.GetEntry(s+3,ss); int x_p_mu_m_nu = SE3->_offset;
+                    // take, one starting by going to the left, the other starting by going to the right.
-                auto SE4 = gStencil.GetEntry(s+4,ss); int x_m_nu      = SE4->_offset;
+                    auto U0 = U_v[x_p_mu     ](nu); gpermute(U0,SE0->_permute);
+                    auto U1 = U_v[x_p_nu     ](mu); gpermute(U1,SE1->_permute);
+                    auto U2 = U_v[x          ](nu); gpermute(U2,SE2->_permute);
+                    auto U3 = U_v[x_p_mu_m_nu](nu); gpermute(U3,SE3->_permute);
+                    auto U4 = U_v[x_m_nu     ](mu); gpermute(U4,SE4->_permute);
+                    auto U5 = U_v[x_m_nu     ](nu); gpermute(U5,SE4->_permute);
 
-                // When you're deciding whether to take an adjoint, the question is: how is the
+                    //       "left"          "right"
-                // stored link oriented compared to the one you want? If I imagine myself travelling
+                    auto W = U2*U1*adj(U0) + adj(U5)*U4*U3;
-                // with the to-be-updated link, I have two possible, alternative 3-link paths I can
-                // take, one starting by going to the left, the other starting by going to the right.
-                auto U0 = U_v[x_p_mu     ](nu); gpermute(U0,SE0->_permute);
-                auto U1 = U_v[x_p_nu     ](mu); gpermute(U1,SE1->_permute);
-                auto U2 = U_v[x          ](nu); gpermute(U2,SE2->_permute);
-                auto U3 = U_v[x_p_mu_m_nu](nu); gpermute(U3,SE3->_permute);
-                auto U4 = U_v[x_m_nu     ](mu); gpermute(U4,SE4->_permute);
-                auto U5 = U_v[x_m_nu     ](nu); gpermute(U5,SE4->_permute);
 
-                //       "left"          "right"
+                    U_3link_v[site](nu) = W;
-                auto W = U2*U1*adj(U0) + adj(U5)*U4*U3;
-                U_fat_v[ss](mu) = U_fat_v[ss](mu) + W;
 
-                s=s+5;
+                    U_fat_v[site](mu) = U_fat_v[site](mu) + lt.c_3*W;
+                }
             }
+
+            // 5-link
+            for(int site=0;site<U_v.size();site++){
+                for(int nu=0;nu<Nd;nu++) {
+                    if(nu==mu) continue;
+                    int s = stencilIndex(mu,nu);
+                    for(int rho=0;rho<Nd;rho++) {
+                        if(rho==nu) continue;
+
+                        auto SE0 = gStencil.GetEntry(s+0,site); int x_p_mu      = SE0->_offset;
+                        auto SE1 = gStencil.GetEntry(s+1,site); int x_p_nu      = SE1->_offset;
+                        auto SE2 = gStencil.GetEntry(s+2,site); int x           = SE2->_offset;
+                        auto SE3 = gStencil.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset;
+                        auto SE4 = gStencil.GetEntry(s+4,site); int x_m_nu      = SE4->_offset;
+
+                        auto U0 =       U_v[x_p_mu     ](nu) ; gpermute(U0,SE0->_permute);
+                        auto U1 = U_3link_v[x_p_nu     ](rho); gpermute(U1,SE1->_permute);
+                        auto U2 =       U_v[x          ](nu) ; gpermute(U2,SE2->_permute);
+                        auto U3 =       U_v[x_p_mu_m_nu](nu) ; gpermute(U3,SE3->_permute);
+                        auto U4 = U_3link_v[x_m_nu     ](rho); gpermute(U4,SE4->_permute);
+                        auto U5 =       U_v[x_m_nu     ](nu) ; gpermute(U5,SE4->_permute);
+
+                        auto W  = U2*U1*adj(U0) + adj(U5)*U4*U3;
+
+                        U_fat_v[site](mu) = U_fat_v[site](mu) + lt.c_5*W;
+                    }
+                }
+            }
+
         }
 
-        u_smr = lt.c_3*Ghost.Extract(Ughost_fat) + lt.c_1*u_thin;
+        u_smr = Ghost.Extract(Ughost_fat) + lt.c_1*u_thin;
     };
 
 
@@ -229,6 +268,7 @@ public:
 //    };
 };
 
+
 /*!  @brief create long links from link variables. */
 template<class LGF>
 class Smear_HISQ_Naik {
@@ -241,6 +281,7 @@ 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() {}

tests/smearing/Benchmark_su3mult_vs_lookup.cc

@@ -1,125 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./benchmarks/Benchmark_su3mult_vs_lookup.cc
-
-    Copyright (C) 2023
-
-    Author: D. A. Clarke <clarke.davida@gmail.com>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-/*
-    @file Benchmark_su3mult_vs_lookup.cc
-    @brief check to see whether su3 multiplication or lookup tables is faster 
-*/
-
-#include <Grid/Grid.h>
-using namespace Grid;
-
-/*!  @brief make the logger work like python print */
-template<typename... Args>
-inline std::string sjoin(Args&&... args) noexcept {
-    std::ostringstream msg;
-    (msg << ... << args);
-    return msg.str();
-}
-template <typename... Args>
-inline void Grid_log(Args&&... args) {
-    std::string msg = sjoin(std::forward<Args>(args)...);
-    std::cout << GridLogMessage << msg << std::endl;
-}
-
-/*!  @brief parameter file to easily adjust Nloop */
-struct ConfParameters: Serializable {
-    GRID_SERIALIZABLE_CLASS_MEMBERS(
-        ConfParameters,
-        int, Nloop);
-
-    template <class ReaderClass>
-    ConfParameters(Reader<ReaderClass>& Reader){
-        read(Reader, "parameters", *this);
-    }
-};
-
-int main (int argc, char** argv) {
-
-    // Params for the test.
-    int Ns               = 8;
-    int Nt               = 4;
-    int threads          = GridThread::GetThreads();
-    std::string conf_in  = "nersc.l8t4b3360";
-    Coordinate latt_size(Nd,0); latt_size[0]=Ns; latt_size[1]=Ns; latt_size[2]=Ns; latt_size[3]=Nt;
-
-    Grid_init(&argc,&argv);
-
-    Coordinate simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd());
-    Coordinate mpi_layout  = GridDefaultMpi();
-    GridCartesian GRID(latt_size,simd_layout,mpi_layout);
-
-    Grid_log("    mpi = ",mpi_layout);
-    Grid_log("   simd = ",simd_layout);
-    Grid_log("   latt = ",latt_size);
-    Grid_log("threads = ",threads);
-
-    XmlReader Reader("mult_vs_lookup.xml",false, "grid");
-    ConfParameters param(Reader); 
-    Grid_log("  Nloop = ",param.Nloop);
-
-    // Gauge field and accessor
-    LatticeGaugeField Umu(&GRID);
-    autoView(U_v, Umu, CpuRead);
-
-    // Read the configuration into Umu
-    FieldMetaData header;
-    NerscIO::readConfiguration(Umu, header, conf_in);
-
-    // Read in lattice sequentially, Nloop times 
-    double lookupTime = 0.; 
-    for(int i=0;i<param.Nloop;i++) {
-        double start = usecond();
-        for(int ss=0;ss<U_v.size();ss++)
-            for(int mu=0;mu<Nd;mu++) {
-                auto U1 = U_v[ss](mu);
-        }
-        double stop  = usecond();
-    	lookupTime += stop-start; // microseconds
-    }
-    Grid_log("Time to lookup: ",lookupTime,"[ms]");
-
-    // Raise a matrix to the power nmat, for each link. 
-    auto U1 = U_v[0](0);
-    for(int nmat=1;nmat<8;nmat++) {
-        double multTime = 0.; 
-        for(int i=0;i<param.Nloop;i++) {
-            double start=usecond();
-            for(int ss=0;ss<U_v.size();ss++)
-                for(int mu=0;mu<Nd;mu++) {
-                    auto U2 = U1;
-                    for(int j=1;j<nmat;j++) {
-                        U2 *= U1;
-                    }
-            }
-            double stop=usecond();
-            multTime += stop-start;
-        }
-        Grid_log("Time to multiply ",nmat," matrices: ",lookupTime," [ms]");
-    }
-
-    Grid_finalize();
-}

tests/smearing/Test_fatLinks.cc

@@ -51,6 +51,20 @@ inline void Grid_log(Args&&... args) {
     std::cout << GridLogMessage << msg << std::endl;
 }
 
+
+/*!  @brief parameter file to easily adjust Nloop */
+struct ConfParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(
+        ConfParameters,
+        int, benchmark, 
+        int, Nloop);
+
+    template <class ReaderClass>
+    ConfParameters(Reader<ReaderClass>& Reader){
+        read(Reader, "parameters", *this);
+    }
+};
+
 //
 // one method: input --> fat
 // another   : input --> long (naik)
@@ -65,16 +79,22 @@ int main (int argc, char** argv) {
     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.3link";
+    int threads          = GridThread::GetThreads();
 
     // Initialize the Grid
     Grid_init(&argc,&argv);
     Coordinate simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd());
     Coordinate mpi_layout  = GridDefaultMpi();
-    Grid_log(" mpi = ",mpi_layout);
+    Grid_log("mpi     = ",mpi_layout);
-    Grid_log("simd = ",simd_layout);
+    Grid_log("simd    = ",simd_layout);
-    Grid_log("latt = ",latt_size);
+    Grid_log("latt    = ",latt_size);
+    Grid_log("threads = ",threads);
     GridCartesian GRID(latt_size,simd_layout,mpi_layout);
 
+    XmlReader Reader("fatParams.xml",false,"grid");
+    ConfParameters param(Reader);
+    if(param.benchmark) Grid_log("  Nloop = ",param.Nloop);
+
     // Instantiate the LatticeGaugeField objects holding thin (Umu) and fat (U_smr) links
     LatticeGaugeField Umu(&GRID);
     LatticeGaugeField U_smr(&GRID);
@@ -85,6 +105,7 @@ int main (int argc, char** argv) {
 
     // Smear Umu and store result in U_smr
     Smear_HISQ_fat<LatticeGaugeField> hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,0.);
+//    Smear_HISQ_fat<LatticeGaugeField> hisq_fat(&GRID,1/8.,0.,1/16.,0.,1/384.,0.);
     hisq_fat.smear(U_smr,Umu);
 
     NerscIO::writeConfiguration(U_smr,conf_out,"HISQ");
@@ -100,5 +121,43 @@ int main (int argc, char** argv) {
     auto absDiff = norm2(diff)/norm2(Umu);
     Grid_log(" |Umu-U|/|Umu| = ",absDiff);
 
+
+    if (param.benchmark) {
+
+        autoView(U_v, Umu, CpuRead); // Gauge accessor
+
+        // Read in lattice sequentially, Nloop times 
+        double lookupTime = 0.; 
+        for(int i=0;i<param.Nloop;i++) {
+            double start = usecond();
+            for(int ss=0;ss<U_v.size();ss++)
+                for(int mu=0;mu<Nd;mu++) {
+                    auto U1 = U_v[ss](mu);
+            }
+            double stop  = usecond();
+        	lookupTime += stop-start; // microseconds
+        }
+        Grid_log("Time to lookup: ",lookupTime,"[ms]");
+
+        // Raise a matrix to the power nmat, for each link. 
+        auto U1 = U_v[0](0);
+        for(int nmat=1;nmat<8;nmat++) {
+            double multTime = 0.; 
+            for(int i=0;i<param.Nloop;i++) {
+                double start=usecond();
+                for(int ss=0;ss<U_v.size();ss++)
+                    for(int mu=0;mu<Nd;mu++) {
+                        auto U2 = U1;
+                        for(int j=1;j<nmat;j++) {
+                            U2 *= U1;
+                        }
+                }
+                double stop=usecond();
+                multTime += stop-start;
+            }
+            Grid_log("Time to multiply ",nmat," matrices: ",multTime," [ms]");
+        }
+    }
+
     Grid_finalize();
 }