Merge b02d022993 into ee1b8bbdbd

.gitignore

@@ -1,3 +1,7 @@
+# Doxygen stuff
+html/*
+latex/*
+
 # Compiled Object files #
 #########################
 *.slo

Grid/log/Log.h

@@ -179,11 +179,11 @@ extern GridLogger GridLogSolver;
 extern GridLogger GridLogError;
 extern GridLogger GridLogWarning;
 extern GridLogger GridLogMessage;
-extern GridLogger GridLogDebug  ;
+extern GridLogger GridLogDebug;
 extern GridLogger GridLogPerformance;
 extern GridLogger GridLogDslash;
-extern GridLogger GridLogIterative  ;
-extern GridLogger GridLogIntegrator  ;
+extern GridLogger GridLogIterative;
+extern GridLogger GridLogIntegrator;
 extern GridLogger GridLogHMC;
 extern GridLogger GridLogMemory;
 extern GridLogger GridLogTracing;
@@ -191,6 +191,41 @@ extern Colours    GridLogColours;
 
 std::string demangle(const char* name) ;
 
+template<typename... Args>
+inline std::string sjoin(Args&&... args) noexcept {
+    std::ostringstream msg;
+    (msg << ... << args);
+    return msg.str();
+}
+
+/*!  @brief make log messages work like python print */
+template <typename... Args>
+inline void Grid_log(Args&&... args) {
+    std::string msg = sjoin(std::forward<Args>(args)...);
+    std::cout << GridLogMessage << msg << std::endl;
+}
+
+/*!  @brief make warning messages work like python print */
+template <typename... Args>
+inline void Grid_warn(Args&&... args) {
+    std::string msg = sjoin(std::forward<Args>(args)...);
+    std::cout << "\033[33m" << GridLogWarning << msg << "\033[0m" << std::endl;
+}
+
+/*!  @brief make error messages work like python print */
+template <typename... Args>
+inline void Grid_error(Args&&... args) {
+    std::string msg = sjoin(std::forward<Args>(args)...);
+    std::cout << "\033[31m" << GridLogError << msg << "\033[0m" << std::endl;
+}
+
+/*!  @brief make pass messages work like python print */
+template <typename... Args>
+inline void Grid_pass(Args&&... args) {
+    std::string msg = sjoin(std::forward<Args>(args)...);
+    std::cout << "\033[32m" << GridLogMessage << msg << "\033[0m" << std::endl;
+}
+
 #define _NBACKTRACE (256)
 extern void * Grid_backtrace_buffer[_NBACKTRACE];
 

Grid/qcd/smearing/HISQSmearing.h

@@ -0,0 +1,392 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/smearing/HISQSmearing.h
+
+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 HISQSmearing.h
+    @brief Declares classes related to HISQ smearing 
+*/
+
+
+#pragma once
+#include <Grid/Grid.h>
+#include <Grid/lattice/PaddedCell.h>
+#include <Grid/stencil/GeneralLocalStencil.h>
+
+
+NAMESPACE_BEGIN(Grid);
+
+
+// TODO: find a way to fold this into the stencil header. need to access grid to get
+// Nd, since you don't want to inherit from QCD.h
+/*!  @brief append arbitrary shift path to shifts */
+template<typename... Args>
+void appendShift(std::vector<Coordinate>& shifts, int dir, Args... args) {
+    Coordinate shift(Nd,0);
+    generalShift(shift, dir, args...); 
+    // push_back creates an element at the end of shifts and
+    // assigns the data in the argument to it.
+    shifts.push_back(shift);
+}
+
+
+/*!  @brief figure out the stencil index from mu and nu */
+accelerator_inline int stencilIndex(int mu, int nu) {
+    // Nshifts depends on how you built the stencil
+    int Nshifts = 6;
+    return Nshifts*nu + Nd*Nshifts*mu;
+}
+
+
+/*!  @brief structure holding the link treatment */
+struct SmearingParameters{
+    SmearingParameters(){}
+    Real c_1;               // 1 link
+    Real c_naik;            // Naik term
+    Real c_3;               // 3 link
+    Real c_5;               // 5 link
+    Real c_7;               // 7 link
+    Real c_lp;              // 5 link Lepage
+    SmearingParameters(Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) 
+        : c_1(c1),
+          c_naik(cnaik),
+          c_3(c3),
+          c_5(c5),
+          c_7(c7),
+          c_lp(clp){}
+};
+
+
+/*!  @brief create fat links from link variables */
+template<class Gimpl> 
+class Smear_HISQ : public Gimpl {
+
+private:
+    GridCartesian* const _grid;
+    SmearingParameters _linkTreatment;
+
+public:
+
+    INHERIT_GIMPL_TYPES(Gimpl);
+    typedef typename Gimpl::GaugeField     GF;
+    typedef typename Gimpl::GaugeLinkField LF;
+    typedef typename Gimpl::ComplexField   CF;
+
+    // Don't allow default values here.
+    Smear_HISQ(GridCartesian* grid, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) 
+        : _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
+    Smear_HISQ(GridCartesian* grid, double* coeff) 
+        : _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() {}
+
+    // Intent: OUT--u_smr, u_naik
+    //          IN--u_thin
+    void smear(GF& u_smr, GF& u_naik, GF& u_thin) const {
+
+        SmearingParameters lt = this->_linkTreatment;
+        auto grid = this->_grid;
+
+        // Create a padded cell of extra padding depth=1 and fill the padding.
+        int depth = 1;
+        PaddedCell Ghost(depth,grid);
+        GF Ughost = Ghost.Exchange(u_thin);
+
+        // This is where auxiliary N-link fields and the final smear will be stored. 
+        GF Ughost_fat(Ughost.Grid());
+        GF Ughost_3link(Ughost.Grid());
+        GF Ughost_5linkA(Ughost.Grid());
+        GF Ughost_5linkB(Ughost.Grid());
+
+        // mu-nu plane stencil. We allow mu==nu to make indexing the stencil easier,
+        // 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++) {
+            appendShift(shifts,mu);
+            appendShift(shifts,nu);
+            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 
+        GeneralLocalStencil gStencil(Ughost.Grid(),shifts);
+
+        // This is where contributions from the smearing get added together
+        Ughost_fat=Zero();
+
+        // This loop handles 3-, 5-, and 7-link constructs, minus Lepage and Naik.
+        for(int mu=0;mu<Nd;mu++) {
+
+            // TODO: This approach is slightly memory inefficient. It uses 25% extra memory 
+            Ughost_3link =Zero();
+            Ughost_5linkA=Zero();
+            Ughost_5linkB=Zero();
+
+            // Create the accessors
+            autoView(U_v       , Ughost       , AcceleratorRead);
+            autoView(U_fat_v   , Ughost_fat   , AcceleratorWrite);
+            autoView(U_3link_v , Ughost_3link , AcceleratorWrite);
+            autoView(U_5linkA_v, Ughost_5linkA, AcceleratorWrite);
+            autoView(U_5linkB_v, Ughost_5linkB, AcceleratorWrite);
+
+            // 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;
+
+            int Nsites = U_v.size();
+            auto gStencil_v = gStencil.View(); 
+
+            accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 3-link constructs
+//            for(int site=0;site<Nsites;site++){ // ----------- 3-link constructs
+                stencilElement SE0, SE1, SE2, SE3, SE4, SE5;
+                U3matrix U0, U1, U2, U3, U4, U5, W;
+                for(int nu=0;nu<Nd;nu++) {
+                    if(nu==mu) continue;
+                    int s = stencilIndex(mu,nu);
+
+                    // The stencil gives us support points in the mu-nu plane that we will use to
+                    // grab the links we need.
+                    SE0 = gStencil_v.GetEntry(s+0,site); int x_p_mu      = SE0->_offset;
+                    SE1 = gStencil_v.GetEntry(s+1,site); int x_p_nu      = SE1->_offset;
+                    SE2 = gStencil_v.GetEntry(s+2,site); int x           = SE2->_offset;
+                    SE3 = gStencil_v.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset;
+                    SE4 = gStencil_v.GetEntry(s+4,site); int x_m_nu      = SE4->_offset;
+                    SE5 = gStencil_v.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
+                    // 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.
+                    U0 = coalescedReadGeneralPermute(U_v[x_p_mu     ](nu),SE0->_permute,Nd);
+                    U1 = coalescedReadGeneralPermute(U_v[x_p_nu     ](mu),SE1->_permute,Nd);
+                    U2 = coalescedReadGeneralPermute(U_v[x          ](nu),SE2->_permute,Nd);
+                    U3 = coalescedReadGeneralPermute(U_v[x_p_mu_m_nu](nu),SE3->_permute,Nd);
+                    U4 = coalescedReadGeneralPermute(U_v[x_m_nu     ](mu),SE4->_permute,Nd);
+                    U5 = coalescedReadGeneralPermute(U_v[x_m_nu     ](nu),SE4->_permute,Nd);
+
+                    //  "left"          "right"
+                    W = U2*U1*adj(U0) + adj(U5)*U4*U3;
+
+                    // Save 3-link construct for later and add to smeared field.
+                    coalescedWrite(U_3link_v[x](nu), W);
+
+                    // The index operator (x) returns the coalesced read on GPU. The view [] index returns 
+                    // a reference to the vector object. The [x](mu) returns a reference to the densely 
+                    // packed (contiguous in memory) mu-th element of the vector object. On CPU, 
+                    // coalescedRead/Write is the identity mapping assigning vector object to vector object.
+                    // But on GPU it's non-trivial and maps scalar object to vector object and vice versa.
+                    coalescedWrite(U_fat_v[x](mu), U_fat_v(x)(mu) + lt.c_3*W);
+                }
+            })
+
+            accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 5-link 
+//            for(int site=0;site<Nsites;site++){ // ----------- 5-link
+                stencilElement SE0, SE1, SE2, SE3, SE4, SE5;
+                U3matrix U0, U1, U2, U3, U4, U5, W;
+                int sigmaIndex = 0;
+                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 == mu || rho == nu) continue;
+
+                        SE0 = gStencil_v.GetEntry(s+0,site); int x_p_mu      = SE0->_offset;
+                        SE1 = gStencil_v.GetEntry(s+1,site); int x_p_nu      = SE1->_offset;
+                        SE2 = gStencil_v.GetEntry(s+2,site); int x           = SE2->_offset;
+                        SE3 = gStencil_v.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset;
+                        SE4 = gStencil_v.GetEntry(s+4,site); int x_m_nu      = SE4->_offset;
+
+                        U0 = coalescedReadGeneralPermute(      U_v[x_p_mu     ](nu ),SE0->_permute,Nd);
+                        U1 = coalescedReadGeneralPermute(U_3link_v[x_p_nu     ](rho),SE1->_permute,Nd);
+                        U2 = coalescedReadGeneralPermute(      U_v[x          ](nu ),SE2->_permute,Nd);
+                        U3 = coalescedReadGeneralPermute(      U_v[x_p_mu_m_nu](nu ),SE3->_permute,Nd);
+                        U4 = coalescedReadGeneralPermute(U_3link_v[x_m_nu     ](rho),SE4->_permute,Nd);
+                        U5 = coalescedReadGeneralPermute(      U_v[x_m_nu     ](nu ),SE4->_permute,Nd);
+
+                        W  = U2*U1*adj(U0) + adj(U5)*U4*U3;
+
+                        if(sigmaIndex<3) {
+                            coalescedWrite(U_5linkA_v[x](rho), W);
+                        } else {
+                            coalescedWrite(U_5linkB_v[x](rho), W);
+                        }    
+
+                        coalescedWrite(U_fat_v[x](mu), U_fat_v(x)(mu) + lt.c_5*W);
+                        sigmaIndex++;
+                    }
+                }
+            })
+
+            accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 7-link
+//            for(int site=0;site<Nsites;site++){ // ----------- 7-link
+                stencilElement SE0, SE1, SE2, SE3, SE4, SE5;
+                U3matrix U0, U1, U2, U3, U4, U5, W;
+                int sigmaIndex = 0;
+                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 == mu || rho == nu) continue;
+
+                        SE0 = gStencil_v.GetEntry(s+0,site); int x_p_mu      = SE0->_offset;
+                        SE1 = gStencil_v.GetEntry(s+1,site); int x_p_nu      = SE1->_offset;
+                        SE2 = gStencil_v.GetEntry(s+2,site); int x           = SE2->_offset;
+                        SE3 = gStencil_v.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset;
+                        SE4 = gStencil_v.GetEntry(s+4,site); int x_m_nu      = SE4->_offset;
+
+                        U0 = coalescedReadGeneralPermute(U_v[x_p_mu](nu),SE0->_permute,Nd);
+                        if(sigmaIndex<3) {
+                            U1 = coalescedReadGeneralPermute(U_5linkB_v[x_p_nu](rho),SE1->_permute,Nd);
+                        } else {
+                            U1 = coalescedReadGeneralPermute(U_5linkA_v[x_p_nu](rho),SE1->_permute,Nd);
+                        }  
+                        U2 = coalescedReadGeneralPermute(U_v[x](nu),SE2->_permute,Nd);
+                        U3 = coalescedReadGeneralPermute(U_v[x_p_mu_m_nu](nu),SE3->_permute,Nd);
+                        if(sigmaIndex<3) {
+                            U4 = coalescedReadGeneralPermute(U_5linkB_v[x_m_nu](rho),SE4->_permute,Nd);
+                        } else {
+                            U4 = coalescedReadGeneralPermute(U_5linkA_v[x_m_nu](rho),SE4->_permute,Nd);
+                        }  
+                        U5 = coalescedReadGeneralPermute(U_v[x_m_nu](nu),SE4->_permute,Nd);
+
+                        W  = U2*U1*adj(U0) + adj(U5)*U4*U3;
+
+                        coalescedWrite(U_fat_v[x](mu), U_fat_v(x)(mu) + lt.c_7*W);
+                        sigmaIndex++;
+                    }
+                }
+            })
+
+        } // end mu loop
+
+        // c1, c3, c5, c7 construct contributions
+        u_smr = Ghost.Extract(Ughost_fat) + lt.c_1*u_thin;
+
+        // Load up U and V std::vectors to access thin and smeared links.
+        std::vector<LF> U(Nd, grid);
+        std::vector<LF> V(Nd, grid);
+        std::vector<LF> Vnaik(Nd, grid);
+        for (int mu = 0; mu < Nd; mu++) {
+            U[mu] = PeekIndex<LorentzIndex>(u_thin, mu);
+            V[mu] = PeekIndex<LorentzIndex>(u_smr, mu);
+        }
+
+        for(int mu=0;mu<Nd;mu++) {
+
+            // Naik
+            Vnaik[mu] = lt.c_naik*Gimpl::CovShiftForward(U[mu],mu,
+                                    Gimpl::CovShiftForward(U[mu],mu,
+                                      Gimpl::CovShiftIdentityForward(U[mu],mu)));
+
+            // LePage
+            for (int nu_h=1;nu_h<Nd;nu_h++) {
+                int nu=(mu+nu_h)%Nd;
+                                // nu, nu, mu, Back(nu), Back(nu)
+                V[mu] = V[mu] + lt.c_lp*Gimpl::CovShiftForward(U[nu],nu,
+                                          Gimpl::CovShiftForward(U[nu],nu,
+                                            Gimpl::CovShiftForward(U[mu],mu,
+                                              Gimpl::CovShiftBackward(U[nu],nu,
+                                                Gimpl::CovShiftIdentityBackward(U[nu],nu)))))
+                                // Back(nu), Back(nu), mu, nu, nu
+                              + lt.c_lp*Gimpl::CovShiftBackward(U[nu],nu,
+                                          Gimpl::CovShiftBackward(U[nu],nu,
+                                            Gimpl::CovShiftForward(U[mu],mu,
+                                              Gimpl::CovShiftForward(U[nu],nu,
+                                                Gimpl::CovShiftIdentityForward(U[nu],nu)))));
+            }
+        }
+
+        // Put V back into u_smr.
+        for (int mu = 0; mu < Nd; mu++) {
+            PokeIndex<LorentzIndex>(u_smr , V[mu]    , mu);
+            PokeIndex<LorentzIndex>(u_naik, Vnaik[mu], mu);
+        }
+    };
+
+
+    // Intent: OUT--u_proj
+    //          IN--u_mu
+    void projectU3(GF& u_proj, GF& u_mu) const {
+
+        auto grid = this->_grid;
+
+        LF V(grid), Q(grid), sqrtQinv(grid), id_3(grid), diff(grid);
+        CF c0(grid), c1(grid), c2(grid), g0(grid), g1(grid), g2(grid), S(grid), R(grid), theta(grid), 
+           u(grid), v(grid), w(grid), den(grid), f0(grid), f1(grid), f2(grid);
+
+        // Follow MILC 10.1103/PhysRevD.82.074501, eqs (B2-B3) and (C1-C8)
+        for (int mu = 0; mu < Nd; mu++) {
+            V  = PeekIndex<LorentzIndex>(u_mu, mu);
+            Q  = adj(V)*V;
+            c0 =        real(trace(Q));
+            c1 = (1/2.)*real(trace(Q*Q));
+            c2 = (1/3.)*real(trace(Q*Q*Q));
+            S  = (1/3.)*c1-(1/18.)*c0*c0;
+            if (norm2(S)<1e-28) {
+                g0 = (1/3.)*c0; g1 = g0; g2 = g1;
+            } else {
+                R     = (1/2.)*c2-(1/3. )*c0*c1+(1/27.)*c0*c0*c0;
+                theta = acos(R*pow(S,-1.5));
+                g0    = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta-2*M_PI/3.);
+                g1    = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta          );
+                g2    = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta+2*M_PI/3.);
+            }
+//            if (fabs(Q.determinant()/(g0*g1*g2)-1.0) > 1e-5) { SVD }
+            u     = sqrt(g0) + sqrt(g1) + sqrt(g2);
+            v     = sqrt(g0*g1) + sqrt(g0*g2) + sqrt(g1*g2);
+            w     = sqrt(g0*g1*g2);
+            den   = w*(u*v-w);
+            f0    = (-w*(u*u+v)+u*v*v)/den;
+            f1    = (-w-u*u*u+2.*u*v)/den;
+            f2    = u/den;
+            id_3  = 1.;
+
+            sqrtQinv = f0*id_3 + f1*Q + f2*Q*Q;
+
+            PokeIndex<LorentzIndex>(u_proj, V*sqrtQinv, mu);
+        }
+    };
+
+
+//    void derivative(const GaugeField& Gauge) const {
+//    };
+};
+
+
+NAMESPACE_END(Grid);

Grid/qcd/smearing/Smearing.h

@@ -5,4 +5,5 @@
 #include <Grid/qcd/smearing/StoutSmearing.h>
 #include <Grid/qcd/smearing/GaugeConfiguration.h>
 #include <Grid/qcd/smearing/WilsonFlow.h>
+#include <Grid/qcd/smearing/HISQSmearing.h>
 

Grid/stencil/GeneralLocalStencil.h

@@ -137,5 +137,49 @@ public:
   
 };
 
+
+////////////////////////////////////////////////
+// Some machinery to streamline making a stencil 
+////////////////////////////////////////////////
+#define BACKWARD_CONST 16
+#define NO_SHIFT -1
+
+// TODO: put a check somewhere that BACKWARD_CONST > Nd!
+
+/*!  @brief signals that you want to go backwards in direction dir */
+inline int Back(const int dir) {
+    // generalShift will use BACKWARD_CONST to determine whether we step forward or 
+    // backward. Trick inspired by SIMULATeQCD. 
+    return dir + BACKWARD_CONST;
+}
+
+/*!  @brief shift one unit in direction dir */
+template<typename... Args>
+void generalShift(Coordinate& shift, int dir) {
+    if (dir >= BACKWARD_CONST) {
+        dir -= BACKWARD_CONST;
+        shift[dir]+=-1;
+    } else if (dir == NO_SHIFT) {
+        ; // do nothing
+    } else {
+        shift[dir]+=1;
+    }
+}
+
+/*!  @brief follow a path of directions, shifting one unit in each direction */
+template<typename... Args>
+void generalShift(Coordinate& shift, int dir, Args... args) {
+    if (dir >= BACKWARD_CONST) {
+        dir -= BACKWARD_CONST;
+        shift[dir]+=-1;
+    } else if (dir == NO_SHIFT) {
+        ; // do nothing
+    } else {
+        shift[dir]+=1;
+    }
+    generalShift(shift, args...);
+}
+
+
 NAMESPACE_END(Grid);
 

examples/Example_plaquette.cc

@@ -0,0 +1,183 @@
+/* 
+ * Example_plaquette.cc                                                               
+ * 
+ * D. Clarke 
+ * 
+ * Here I just want to create an incredibly simple main to get started with GRID and get used
+ * to its syntax. If the reader is like me, they vaguely understand something about lattice coding,
+ * they don't know a ton of C++, don't know much of the fine details, and certainly know nothing about GRID.
+ *
+ * Once you've made a new executable, like this one, you can bootstrap.sh again. At this point,
+ * the code should be able to find your new executable. You can tell that bootstrap.sh worked by
+ * having a look at Make.inc. You should see your executable inside there.
+ *
+ * Warning: This code illustrative only, not well tested, and not meant for production use. The best
+ * way to read this code is to start at the main.
+ * 
+ */
+
+
+// All your mains should have this
+#include <Grid/Grid.h>
+using namespace Grid;
+
+
+// This copies what already exists in WilsonLoops.h. The point here is to be pedagogical and explain in
+// detail what everything does so we can see how GRID works.
+template <class Gimpl> class WLoops : public Gimpl {
+public:
+    // Gimpl seems to be an arbitrary class. Within this class, it is expected that certain types are
+    // already defined, things like Scalar and Field. This macro includes a bunch of #typedefs that
+    // implement this equivalence at compile time.
+    INHERIT_GIMPL_TYPES(Gimpl);
+
+    // Some example Gimpls can be found in GaugeImplementations.h, at the bottom. These are in turn built
+    // out of GaugeImplTypes, which can be found in GaugeImplTypes.h. The GaugeImplTypes contain the base
+    // field/vector/link/whatever types. These inherit from iScalar, iVector, and iMatrix objects, which
+    // are sort of the building blocks for gerenal math objects. The "i" at the beginning of these names
+    // indicates that they should be for internal use only. It seems like these base types have the
+    // acceleration, e.g. SIMD or GPU or what-have-you, abstracted away. How you accelerate these things
+    // appears to be controlled through a template parameter called vtype.
+
+    // The general math/physics objects, such as a color matrix, are built up by nesting these objects.
+    // For instance a general color matrix has two color indices, so it's built up like
+    //     iScalar<iScalar<iMatrix<vtype ...
+    // where the levels going from the inside out are color, spin, then Lorentz indices. Scalars have
+    // no indices, so it's what we use when such an index isn't needed. Lattice objects are made by one
+    // higher level of indexing using iVector.
+
+    // These types will be used for U and U_mu objects, respectively.
+    typedef typename Gimpl::GaugeLinkField GaugeMat;
+    typedef typename Gimpl::GaugeField GaugeLorentz;
+
+    // U_mu_nu(x)
+    static void dirPlaquette(GaugeMat &plaq, const std::vector<GaugeMat> &U, const int mu, const int nu) {
+        // Calls like CovShiftForward and CovShiftBackward have 3 arguments, and they multiply together
+        // the first and last argument. (Second arg gives the shift direction.) The CovShiftIdentityBackward
+        // has meanwhile only two arguments; it just returns the shifted (adjoint since backward) link. 
+        plaq = Gimpl::CovShiftForward(U[mu],mu,
+                   // Means Link*Cshift(field,mu,1), arguments are Link, mu, field in that order.
+                   Gimpl::CovShiftForward(U[nu],nu,
+                       Gimpl::CovShiftBackward(U[mu],mu,
+                           // This means Cshift(adj(Link), mu, -1)
+                           Gimpl::CovShiftIdentityBackward(U[nu], nu))));
+    }
+
+    // tr U_mu_nu(x)
+    static void traceDirPlaquette(ComplexField &plaq, const std::vector<GaugeMat> &U, const int mu, const int nu) {
+        // This .Grid() syntax seems to get the pointer to the GridBase. Apparently this is needed as argument
+        // to instantiate a Lattice object.
+        GaugeMat sp(U[0].Grid());
+        dirPlaquette(sp, U, mu, nu);
+        plaq = trace(sp);
+    }
+
+    // sum_mu_nu tr U_mu_nu(x)
+    static void sitePlaquette(ComplexField &Plaq, const std::vector<GaugeMat> &U) {
+        ComplexField sitePlaq(U[0].Grid());
+        Plaq = Zero();
+        // Nd=4 and Nc=3 are set as global constants in QCD.h
+        for (int mu = 1; mu < Nd; mu++) {
+            for (int nu = 0; nu < mu; nu++) {
+                traceDirPlaquette(sitePlaq, U, mu, nu);
+                Plaq = Plaq + sitePlaq;
+            }
+        }
+    }
+
+    // sum_mu_nu_x Re tr U_mu_nu(x)
+    static RealD sumPlaquette(const GaugeLorentz &Umu) {
+        std::vector<GaugeMat> U(Nd, Umu.Grid());
+        for (int mu = 0; mu < Nd; mu++) {
+            // Umu is a GaugeLorentz object, and as such has a non-trivial Lorentz index. We can
+            // access the element in the mu Lorentz index with this PeekIndex syntax.
+            U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
+        }
+        ComplexField Plaq(Umu.Grid());
+        sitePlaquette(Plaq, U);
+        // I guess this should be the line that sums over all space-time sites.
+        auto Tp = sum(Plaq);
+        // Until now, we have been working with objects inside the tensor nest. This TensorRemove gets
+        // rid of the tensor nest to return whatever is inside.
+        auto p  = TensorRemove(Tp);
+        return p.real();
+    }
+
+    // < Re tr U_mu_nu(x) >
+    static RealD avgPlaquette(const GaugeLorentz &Umu) {
+        // Real double type
+        RealD sumplaq = sumPlaquette(Umu);
+        // gSites() is the number of global sites. there is also lSites() for local sites.
+        double vol = Umu.Grid()->gSites();
+        // The number of orientations. 4*3/2=6 for Nd=4, as known.
+        double faces = (1.0 * Nd * (Nd - 1)) / 2.0;
+        return sumplaq / vol / faces / Nc;
+    }
+};
+
+
+// Next we show an example of how to construct an input parameter class. We first inherit
+// from Serializable. Then all class data members have to be defined using the
+// GRID_SERIALIZABLE_CLASS_MEMBERS macro. This variadic macro allows for arbitrarily many
+// class data members. In the below case, we make a parameter file holding the configuration
+// name. Here, it expects the name to be labeled with "conf_name" in the configuration file. 
+struct ConfParameters: Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(
+        ConfParameters,
+        std::string, conf_name);
+
+    template <class ReaderClass>
+    ConfParameters(Reader<ReaderClass>& Reader){
+        // If we are reading an XML file, it should be structured like:
+        // <grid>
+        //   <parameters>
+        //     <conf_name>l20t20b06498a_nersc.302500</conf_name>
+        //   </parameters>
+        // </grid>
+        read(Reader, "parameters", *this);
+    }
+};
+
+
+
+// This syntax lets you pass command line arguments to main. An asterisk means that what follows is
+// a pointer. Two asterisks means what follows is a pointer to an array. 
+int main (int argc, char **argv)
+{
+    // This initializes Grid. Some command line options include
+    //   --mpi n.n.n.n
+    //   --threads n
+    //   --grid n.n.n.n
+    Grid_init(&argc, &argv);
+
+    // This is where you would specify a custom lattice size, if not from the command line. Here
+    // Nd is a global quantity that is currently set to 4.
+    Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+    Coordinate mpi_layout  = GridDefaultMpi();
+    Coordinate latt_size   = GridDefaultLatt();
+
+    // Instantiate the spacetime Grid on which everything will be built.
+    GridCartesian GRID(latt_size,simd_layout,mpi_layout);
+
+    // The PeriodicGimplD type is what you want for gauge matrices. There is also a LatticeGaugeFieldD
+    // type that you can use, which will work perfectly with what follows. 
+    PeriodicGimplD::Field U(&GRID);
+
+    // Here we read in the parameter file params.json to get conf_name. The last argument is what the
+    // top organizational level is called in the param file. 
+    XmlReader Reader("Example_plaquette.xml",false, "grid");
+    ConfParameters param(Reader);  
+
+    // Load a lattice from SIMULATeQCD into U. SIMULATeQCD finds plaquette = 0.6381995717
+    FieldMetaData header;
+    NerscIO::readConfiguration(U, header, param.conf_name);
+
+    // Let's see what we find.
+    RealD plaq = WLoops<PeriodicGimplD>::avgPlaquette(U);
+
+    // This is how you make log messages.
+    std::cout << GridLogMessage << std::setprecision(std::numeric_limits<Real>::digits10 + 1) << "Plaquette = " << plaq << std::endl;
+
+    // To wrap things up.
+    Grid_finalize();
+}

tests/debug/Test_padded_cell.cc

@@ -32,6 +32,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 using namespace std;
 using namespace Grid;
 
+// This is to optimize the SIMD
 template<class vobj> void gpermute(vobj & inout,int perm){
   vobj tmp=inout;
   if (perm & 0x1 ) { permute(inout,tmp,0); tmp=inout;}
@@ -39,7 +40,8 @@ template<class vobj> void gpermute(vobj & inout,int perm){
   if (perm & 0x4 ) { permute(inout,tmp,2); tmp=inout;}
   if (perm & 0x8 ) { permute(inout,tmp,3); tmp=inout;}
 }
-  
+
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
@@ -47,20 +49,21 @@ int main (int argc, char ** argv)
   Coordinate latt_size  = GridDefaultLatt();
   Coordinate simd_layout= GridDefaultSimd(Nd,vComplexD::Nsimd());
   Coordinate mpi_layout = GridDefaultMpi();
-  std::cout << " mpi "<<mpi_layout<<std::endl;
-  std::cout << " simd "<<simd_layout<<std::endl;
-  std::cout << " latt "<<latt_size<<std::endl;
+  std::cout << GridLogMessage << " mpi "<<mpi_layout<<std::endl;
+  std::cout << GridLogMessage << " simd "<<simd_layout<<std::endl;
+  std::cout << GridLogMessage << " latt "<<latt_size<<std::endl;
   GridCartesian GRID(latt_size,simd_layout,mpi_layout);
 
+  // Initialize configuration as hot start.
   GridParallelRNG   pRNG(&GRID);
-  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
   LatticeGaugeField Umu(&GRID);
-
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
   SU<Nc>::HotConfiguration(pRNG,Umu);
 
   Real plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
   LatticeComplex trplaq(&GRID);
 
+  // Store Umu in U. Peek/Poke mean respectively getElement/setElement.
   std::vector<LatticeColourMatrix> U(Nd, Umu.Grid());
   for (int mu = 0; mu < Nd; mu++) {
     U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
@@ -70,9 +73,7 @@ int main (int argc, char ** argv)
 
   LatticeComplex cplaq(&GRID); cplaq=Zero();
 
-  /////////////////////////////////////////////////
   // Create a padded cell of extra padding depth=1
-  /////////////////////////////////////////////////
   int depth = 1;
   PaddedCell Ghost(depth,&GRID);
   LatticeGaugeField Ughost = Ghost.Exchange(Umu);
@@ -114,18 +115,25 @@ int main (int argc, char ** argv)
   }
 #endif
 
-  ///// Array for the site plaquette
+  // Array for the site plaquette
   GridBase *GhostGrid = Ughost.Grid();
   LatticeComplex gplaq(GhostGrid); 
-  
+
+  // Now we're going to put together the "stencil" that will be useful to us when
+  // calculating the plaquette. Our eventual goal is to make the product
+  //    Umu(x) Unu(x+mu) Umu^dag(x+nu) Unu^dag(x),
+  // which requires, in order, the sites x, x+mu, x+nu, and x. We arrive at these
+  // sites relative to x through "shifts", which is represented here by a 4-d
+  // vector of 0s (no movement) and 1s (shift one unit) at each site. The
+  // "stencil" is the set of all these shifts.
   std::vector<Coordinate> shifts;
   for(int mu=0;mu<Nd;mu++){
     for(int nu=mu+1;nu<Nd;nu++){
-  
-      //    Umu(x) Unu(x+mu) Umu^dag(x+nu) Unu^dag(x)
       Coordinate shift_0(Nd,0);
       Coordinate shift_mu(Nd,0); shift_mu[mu]=1;
       Coordinate shift_nu(Nd,0); shift_nu[nu]=1;
+      // push_back creates an element at the end of shifts and
+      // assigns the data in the argument to it.
       shifts.push_back(shift_0);
       shifts.push_back(shift_mu);
       shifts.push_back(shift_nu);
@@ -135,41 +143,51 @@ int main (int argc, char ** argv)
   GeneralLocalStencil gStencil(GhostGrid,shifts);
 
   gplaq=Zero();
-  {
-    autoView( gp_v , gplaq, CpuWrite);
-    autoView( t_v , trplaq, CpuRead);
-    autoView( U_v , Ughost, CpuRead);
-    for(int ss=0;ss<gp_v.size();ss++){
-      int s=0;
-      for(int mu=0;mu<Nd;mu++){
-	for(int nu=mu+1;nu<Nd;nu++){
 
-	  auto SE0 = gStencil.GetEntry(s+0,ss);
-	  auto SE1 = gStencil.GetEntry(s+1,ss);
-	  auto SE2 = gStencil.GetEntry(s+2,ss);
-	  auto SE3 = gStencil.GetEntry(s+3,ss);
-	
-	  int o0 = SE0->_offset;
-	  int o1 = SE1->_offset;
-	  int o2 = SE2->_offset;
-	  int o3 = SE3->_offset;
-	  
-	  auto U0 = U_v[o0](mu);
-	  auto U1 = U_v[o1](nu);
-	  auto U2 = adj(U_v[o2](mu));
-	  auto U3 = adj(U_v[o3](nu));
+  // Before doing accelerator stuff, there is an opening and closing of "Views". I guess the
+  // "Views" are stored in *_v variables listed below.
+  autoView( gp_v , gplaq, CpuWrite);
+  autoView( t_v , trplaq, CpuRead);
+  autoView( U_v , Ughost, CpuRead);
 
-	  gpermute(U0,SE0->_permute);
-	  gpermute(U1,SE1->_permute);
-	  gpermute(U2,SE2->_permute);
-	  gpermute(U3,SE3->_permute);
-	  
-	  gp_v[ss]() =gp_v[ss]() + trace( U0*U1*U2*U3 );
-	  s=s+4;
-	}
-      }
+  // This is now a loop over stencil shift elements. That is, s increases as we make our
+  // way through the spacetimes sites, but also as we make our way around the plaquette.
+  for(int ss=0;ss<gp_v.size();ss++){
+    int s=0;
+    for(int mu=0;mu<Nd;mu++){
+    	for(int nu=mu+1;nu<Nd;nu++){
+    
+    	  auto SE0 = gStencil.GetEntry(s+0,ss);
+    	  auto SE1 = gStencil.GetEntry(s+1,ss);
+    	  auto SE2 = gStencil.GetEntry(s+2,ss);
+    	  auto SE3 = gStencil.GetEntry(s+3,ss);
+
+        // Due to our strategy, each offset corresponds to a site.
+    	  int o0 = SE0->_offset;
+    	  int o1 = SE1->_offset;
+    	  int o2 = SE2->_offset;
+    	  int o3 = SE3->_offset;
+    	  
+    	  auto U0 = U_v[o0](mu);
+    	  auto U1 = U_v[o1](nu);
+    	  auto U2 = adj(U_v[o2](mu));
+    	  auto U3 = adj(U_v[o3](nu));
+    
+    	  gpermute(U0,SE0->_permute);
+    	  gpermute(U1,SE1->_permute);
+    	  gpermute(U2,SE2->_permute);
+    	  gpermute(U3,SE3->_permute);
+    	  
+    	  gp_v[ss]() =gp_v[ss]() + trace( U0*U1*U2*U3 );
+    	  s=s+4;
+    	}
     }
   }
+
+  // Here is my understanding of this part: The padded cell has its own periodic BCs, so
+  // if I take a step to the right at the right-most side of the cell, I end up on the
+  // left-most side. This means that the plaquettes in the padding are wrong. Luckily
+  // all we care about are the plaquettes in the cell, which we obtain from Extract.
   cplaq = Ghost.Extract(gplaq);
   RealD vol = cplaq.Grid()->gSites();
   RealD faces = (Nd * (Nd-1))/2;

tests/smearing/Test_fatLinks.cc

@@ -0,0 +1,181 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/smearing/Test_fatLinks.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 Test_fatLinks.cc
+    @brief test of the HISQ smearing 
+*/
+
+
+#include <Grid/Grid.h>
+#include <Grid/lattice/PaddedCell.h>
+#include <Grid/stencil/GeneralLocalStencil.h>
+#include <Grid/qcd/smearing/HISQSmearing.h>
+using namespace Grid;
+
+
+/*!  @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);
+    }
+};
+
+
+bool testSmear(GridCartesian& GRID, LatticeGaugeFieldD Umu, LatticeGaugeFieldD Usmr, LatticeGaugeFieldD Unaik, 
+               LatticeGaugeFieldD Ucontrol, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) {
+    Smear_HISQ<PeriodicGimplD> hisq_fat(&GRID,c1,cnaik,c3,c5,c7,clp);
+    LatticeGaugeFieldD diff(&GRID), Uproj(&GRID);
+    hisq_fat.smear(Usmr, Unaik, Umu);
+    bool result;
+    if (cnaik < 1e-30) { // Testing anything but Naik term
+        diff = Ucontrol-Usmr;
+        auto absDiff = norm2(diff)/norm2(Ucontrol);
+        if (absDiff < 1e-30) {
+            Grid_pass(" |Umu-Usmr|/|Umu| = ",absDiff);
+            result = true;
+        } else {
+            Grid_error(" |Umu-Usmr|/|Umu| = ",absDiff);
+            result = false;
+        }
+    } else { // Testing Naik specifically
+        diff = Ucontrol-Unaik;
+        auto absDiff = norm2(diff)/norm2(Ucontrol);
+        if (absDiff < 1e-30) {
+            Grid_pass(" |Umu-Unaik|/|Umu| = ",absDiff);
+            result = true;
+        } else {
+            Grid_error(" |Umu-Unaik|/|Umu| = ",absDiff);
+            result = false;
+        }
+        hisq_fat.projectU3(Uproj,Ucontrol);
+//        NerscIO::writeConfiguration(Unaik,"nersc.l8t4b3360.naik");
+    }
+    return result;
+}
+
+
+int main (int argc, char** argv) {
+
+    // Params for the test.
+    int Ns = 8;
+    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";
+    int threads          = GridThread::GetThreads();
+
+    typedef LatticeGaugeFieldD LGF;
+
+    // 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("simd    = ",simd_layout);
+    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);
+
+    LGF Umu(&GRID), Usmr(&GRID), Unaik(&GRID), Ucontrol(&GRID);
+
+    // Read the configuration into Umu
+    FieldMetaData header;
+    NerscIO::readConfiguration(Umu, header, conf_in);
+
+    bool pass=true;
+
+    // Carry out various tests    
+    NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.357lplink.control");
+    pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,1/384.,-1/8.);
+    NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.357link.control");
+    pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,1/384.,0.);
+    NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.35link.control");
+    pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,0.,0.);
+    NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.3link.control");
+    pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,0.,0.,0.);
+    NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.naik.control");
+    pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,0.,0.8675309,0.,0.,0.,0.);
+
+    if(pass){
+        Grid_pass("All tests passed.");
+    } else {
+        Grid_error("At least one test failed.");
+    }
+
+    // Test a C-style instantiation 
+    double path_coeff[6] = {1, 2, 3, 4, 5, 6};
+    Smear_HISQ<PeriodicGimplD> hisq_fat_Cstyle(&GRID,path_coeff);
+
+    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();
+}