From bf91778550218008230e56aff7c4c2f6b9574589 Mon Sep 17 00:00:00 2001 From: david clarke Date: Wed, 17 May 2023 15:15:54 -0600 Subject: [PATCH 01/40] verbose plaquette example; fat link test frame --- examples/Example_plaquette.cc | 180 ++++++++++++++++++++++++++++++++ tests/smearing/Test_fatLinks.cc | 40 +++++++ 2 files changed, 220 insertions(+) create mode 100644 examples/Example_plaquette.cc create mode 100644 tests/smearing/Test_fatLinks.cc diff --git a/examples/Example_plaquette.cc b/examples/Example_plaquette.cc new file mode 100644 index 00000000..2aec72ff --- /dev/null +++ b/examples/Example_plaquette.cc @@ -0,0 +1,180 @@ +/* + * 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 +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 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. + // WARNING: The first time you include this or take it out, the compile time will increase a lot. + 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 &U, const int mu, const int nu) { + // These CovShift calls seem to carry out the multiplication already. A positive shift moves the lattice + // site x_mu = 1 in the RHS to x_mu = 0 in the result. + plaq = Gimpl::CovShiftForward(U[mu],mu, + Gimpl::CovShiftForward(U[nu],nu, + Gimpl::CovShiftBackward(U[mu],mu, + Gimpl::CovShiftIdentityBackward(U[nu], nu)))); + } + + // tr U_mu_nu(x) + static void traceDirPlaquette(ComplexField &plaq, const std::vector &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 &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 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(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 + ConfParameters(Reader& Reader){ + // If we are reading an XML file, it should be structured like: + // + // + // l20t20b06498a_nersc.302500 + // + // + 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. + Coordinate simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + Coordinate latt_size = GridDefaultLatt(); + + // Instantiate the Grid on which everything will be built. + GridCartesian spacetime(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(&spacetime); + + // 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("params.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::avgPlaquette(U); + + // This is how you make log messages. + std::cout << GridLogMessage << std::setprecision(std::numeric_limits::digits10 + 1) << "Plaquette = " << plaq << std::endl; + + // To wrap things up. + Grid_finalize(); +} \ No newline at end of file diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc new file mode 100644 index 00000000..4f9d608d --- /dev/null +++ b/tests/smearing/Test_fatLinks.cc @@ -0,0 +1,40 @@ +/* + * Test_fatLinks.cc + * + * D. Clarke + * + * Test the various constructs used to make fat links. + * + */ + + +#include +using namespace Grid; + +template class : public Gimpl { +public: + + INHERIT_GIMPL_TYPES(Gimpl); + + typedef typename Gimpl::GaugeLinkField GaugeMat; + typedef typename Gimpl::GaugeField GaugeLorentz; + + static void staple(GaugeMat &plaq, const std::vector &U, const int mu, const int nu) { + } + +} + +int main (int argc, char **argv) +{ + Grid_init(&argc, &argv); + + Coordinate simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + Coordinate latt_size = GridDefaultLatt(); + + GridCartesian spacetime(latt_size,simd_layout,mpi_layout); + + PeriodicGimplD::Field U(&spacetime); + + Grid_finalize(); +} \ No newline at end of file From c7bdf2c0e475f5b6d54af8db7b01075efc94485a Mon Sep 17 00:00:00 2001 From: david clarke Date: Sun, 21 May 2023 04:33:20 -0600 Subject: [PATCH 02/40] 3-link test at least gives an answer --- examples/Example_plaquette.cc | 13 +-- tests/debug/Test_padded_cell.cc | 104 +++++++++++++--------- tests/smearing/Test_fatLinks.cc | 151 ++++++++++++++++++++++++++++---- 3 files changed, 202 insertions(+), 66 deletions(-) diff --git a/examples/Example_plaquette.cc b/examples/Example_plaquette.cc index 2aec72ff..17de4762 100644 --- a/examples/Example_plaquette.cc +++ b/examples/Example_plaquette.cc @@ -148,21 +148,22 @@ int main (int argc, char **argv) // --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. - Coordinate simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); + // 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 Grid on which everything will be built. - GridCartesian spacetime(latt_size,simd_layout,mpi_layout); + // 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(&spacetime); + 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("params.xml",false, "grid"); + XmlReader Reader("Example_plaquette.xml",false, "grid"); ConfParameters param(Reader); // Load a lattice from SIMULATeQCD into U. SIMULATeQCD finds plaquette = 0.6381995717 diff --git a/tests/debug/Test_padded_cell.cc b/tests/debug/Test_padded_cell.cc index 4fb461fe..f110df46 100644 --- a/tests/debug/Test_padded_cell.cc +++ b/tests/debug/Test_padded_cell.cc @@ -32,6 +32,7 @@ Author: Peter Boyle using namespace std; using namespace Grid; +// This is to optimize the SIMD template void gpermute(vobj & inout,int perm){ vobj tmp=inout; if (perm & 0x1 ) { permute(inout,tmp,0); tmp=inout;} @@ -39,7 +40,8 @@ template 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 "<({45,12,81,9})); LatticeGaugeField Umu(&GRID); - + pRNG.SeedFixedIntegers(std::vector({45,12,81,9})); SU::HotConfiguration(pRNG,Umu); Real plaq=WilsonLoops::avgPlaquette(Umu); LatticeComplex trplaq(&GRID); + // Store Umu in U. Peek/Poke mean respectively getElement/setElement. std::vector U(Nd, Umu.Grid()); for (int mu = 0; mu < Nd; mu++) { U[mu] = PeekIndex(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 shifts; for(int mu=0;mu_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_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; diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index 4f9d608d..85a60ee3 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -7,34 +7,151 @@ * */ - #include +#include +#include using namespace Grid; -template class : public Gimpl { -public: - - INHERIT_GIMPL_TYPES(Gimpl); - - typedef typename Gimpl::GaugeLinkField GaugeMat; - typedef typename Gimpl::GaugeField GaugeLorentz; - - static void staple(GaugeMat &plaq, const std::vector &U, const int mu, const int nu) { - } - +// This is to optimize the SIMD +template void gpermute(vobj & inout,int perm) { + vobj tmp=inout; + if (perm & 0x1) {permute(inout,tmp,0); tmp=inout;} + if (perm & 0x2) {permute(inout,tmp,1); tmp=inout;} + if (perm & 0x4) {permute(inout,tmp,2); tmp=inout;} + if (perm & 0x8) {permute(inout,tmp,3); tmp=inout;} } +// Make the logger work like Python print() +template +inline std::string sjoin(Args&&... args) noexcept { + std::ostringstream msg; + (msg << ... << args); + return msg.str(); +} +template +inline void Grid_log(Args&&... args) { + std::string msg = sjoin(std::forward(args)...); + std::cout << GridLogMessage << msg << std::endl; +} + +struct fatParams: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS( + fatParams, + std::string, conf_in, + std::string, conf_out); + + template + fatParams(Reader& Reader){ + read(Reader, "parameters", *this); + } +}; + + + int main (int argc, char **argv) { - Grid_init(&argc, &argv); + Grid_init(&argc,&argv); - Coordinate simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); - Coordinate mpi_layout = GridDefaultMpi(); Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); - GridCartesian spacetime(latt_size,simd_layout,mpi_layout); + Grid_log("mpi = ",mpi_layout); + Grid_log("simd = ",simd_layout); + Grid_log("latt = ",latt_size); - PeriodicGimplD::Field U(&spacetime); + GridCartesian GRID(latt_size,simd_layout,mpi_layout); + + XmlReader Reader("fatParams.xml",false, "grid"); + fatParams param(Reader); + + LatticeGaugeField Umu(&GRID); + FieldMetaData header; + NerscIO::readConfiguration(Umu, header, param.conf_in); + + // Create a padded cell of extra padding depth=1 + int depth = 1; + PaddedCell Ghost(depth,&GRID); + LatticeGaugeField Ughost = Ghost.Exchange(Umu); + + // Array for (x) + GridBase *GhostGrid = Ughost.Grid(); + LatticeComplex gplaq(GhostGrid); + + // This is where the 3-link constructs will be stored + LatticeGaugeField Ughost_3link(Ughost.Grid()); + + // Create 3-link stencil + std::vector shifts; + for(int mu=0;mu_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); + + auto W = U1*U2*U3; + + // We add together contributions coming from each orientation. + U_3link_v[ss](mu) = U_3link_v[ss](mu) + W; + + 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. + Umu = Ghost.Extract(Ughost_3link); + + NerscIO::writeConfiguration(Umu,param.conf_out,"HISQ"); Grid_finalize(); } \ No newline at end of file From ab56ad8d7a23a5662e1045dbedad4d06e57453a8 Mon Sep 17 00:00:00 2001 From: david clarke Date: Wed, 7 Jun 2023 21:14:58 -0600 Subject: [PATCH 03/40] fix 3-link stencil --- tests/smearing/Test_fatLinks.cc | 51 +++++++++++++++++++++++++-------- 1 file changed, 39 insertions(+), 12 deletions(-) diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index 85a60ee3..99081fea 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -12,7 +12,7 @@ #include using namespace Grid; -// This is to optimize the SIMD +// This is to optimize the SIMD (will also need to be in the class, at least for now) template void gpermute(vobj & inout,int perm) { vobj tmp=inout; if (perm & 0x1) {permute(inout,tmp,0); tmp=inout;} @@ -81,19 +81,27 @@ int main (int argc, char **argv) // This is where the 3-link constructs will be stored LatticeGaugeField Ughost_3link(Ughost.Grid()); - // Create 3-link stencil + // Create 3-link stencil (class will build its own stencils) + // writing your own stencil, you're hard-coding the periodic BCs, so you don't need + // the policy-based stuff, at least for now std::vector shifts; for(int mu=0;mu_offset; int o1 = SE1->_offset; int o2 = SE2->_offset; int o3 = SE3->_offset; + int o4 = SE4->_offset; + int o5 = SE5->_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)); + auto U0 = U_v[o0](nu); + auto U1 = adj(U_v[o1](mu)); + auto U2 = adj(U_v[o2](nu)); gpermute(U0,SE0->_permute); gpermute(U1,SE1->_permute); gpermute(U2,SE2->_permute); + + auto U3 = adj(U_v[o3](nu)); + auto U4 = adj(U_v[o4](mu)); + auto U5 = U_v[o5](nu); + gpermute(U3,SE3->_permute); + gpermute(U4,SE4->_permute); + gpermute(U5,SE5->_permute); - auto W = U1*U2*U3; - - // We add together contributions coming from each orientation. + // Forward contribution from this orientation + auto W = U0*U1*U2; U_3link_v[ss](mu) = U_3link_v[ss](mu) + W; - s=s+4; + // Backward contribution from this orientation + W = U3*U4*U5; + U_3link_v[ss](mu) = U_3link_v[ss](mu) + W; + + s=s+6; } } } From 4b994a1bc75901f124c886523c278c7ec771c899 Mon Sep 17 00:00:00 2001 From: david clarke Date: Thu, 8 Jun 2023 17:37:25 -0600 Subject: [PATCH 04/40] trouble with compilation --- Grid/qcd/smearing/HISQSmearing.h | 190 +++++++++++++++++++++++++++++++ tests/smearing/Test_fatLinks.cc | 127 ++------------------- 2 files changed, 200 insertions(+), 117 deletions(-) create mode 100644 Grid/qcd/smearing/HISQSmearing.h diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h new file mode 100644 index 00000000..53480f25 --- /dev/null +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -0,0 +1,190 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/smearing/StoutSmearing.h + +Copyright (C) 2019 + +Author: D. A. Clarke + +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 +*/ + +// things like @brief are seen by things like doxygen and javadocs + +#pragma once + +#include +#include +#include + +NAMESPACE_BEGIN(Grid); + + +// This is to optimize the SIMD (will also need to be in the class, at least for now) +template void gpermute(vobj & inout,int perm) { + vobj tmp=inout; + if (perm & 0x1) {permute(inout,tmp,0); tmp=inout;} + if (perm & 0x2) {permute(inout,tmp,1); tmp=inout;} + if (perm & 0x4) {permute(inout,tmp,2); tmp=inout;} + if (perm & 0x8) {permute(inout,tmp,3); tmp=inout;} +} + + +/*! @brief 3-link smearing of link variable. */ +//template +//class Smear_HISQ_3link : public Smear { +class Smear_HISQ_3link { +// TODO: I'm not using Gimpl so I don't know how to inherit + +private: +// std::vector _linkTreatment; + GridBase* const _grid; + +public: +// INHERIT_GIMPL_TYPES(Gimpl) + + // Eventually this will take, e.g., coefficients as argument + Smear_HISQ_3link(GridBase* grid) : _grid(grid) { + assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3"); + } + + ~Smear_HISQ_3link() {} + + void smear(LatticeGaugeField& u_smr, const LatticeGaugeField& U) const { + + // Create a padded cell of extra padding depth=1 + int depth = 1; + PaddedCell Ghost(depth,this->_grid); + LatticeGaugeField Ughost = Ghost.Exchange(u_smr); + + // Array for (x) + GridBase *GhostGrid = Ughost.Grid(); + LatticeComplex gplaq(GhostGrid); + + // This is where the 3-link constructs will be stored + LatticeGaugeField Ughost_3link(Ughost.Grid()); + + // Create 3-link stencil (class will build its own stencils) + // writing your own stencil, you're hard-coding the periodic BCs, so you don't need + // the policy-based stuff, at least for now + std::vector shifts; + for(int mu=0;mu_offset; + int o1 = SE1->_offset; + int o2 = SE2->_offset; + int o3 = SE3->_offset; + int o4 = SE4->_offset; + int o5 = SE5->_offset; + + auto U0 = U_v[o0](nu); + auto U1 = adj(U_v[o1](mu)); + auto U2 = adj(U_v[o2](nu)); + + gpermute(U0,SE0->_permute); + gpermute(U1,SE1->_permute); + gpermute(U2,SE2->_permute); + + auto U3 = adj(U_v[o3](nu)); + auto U4 = adj(U_v[o4](mu)); + auto U5 = U_v[o5](nu); + + gpermute(U3,SE3->_permute); + gpermute(U4,SE4->_permute); + gpermute(U5,SE5->_permute); + + // Forward contribution from this orientation + auto W = U0*U1*U2; + U_3link_v[ss](mu) = U_3link_v[ss](mu) + W; + + // Backward contribution from this orientation + W = U3*U4*U5; + U_3link_v[ss](mu) = U_3link_v[ss](mu) + W; + + s=s+6; + } + } + } + + // 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. + u_smr = Ghost.Extract(Ughost_3link); + }; + +// void derivative(const GaugeField& Gauge) const { +// }; +}; + +NAMESPACE_END(Grid); diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index 99081fea..124cfac3 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -10,16 +10,9 @@ #include #include #include +#include using namespace Grid; -// This is to optimize the SIMD (will also need to be in the class, at least for now) -template void gpermute(vobj & inout,int perm) { - vobj tmp=inout; - if (perm & 0x1) {permute(inout,tmp,0); tmp=inout;} - if (perm & 0x2) {permute(inout,tmp,1); tmp=inout;} - if (perm & 0x4) {permute(inout,tmp,2); tmp=inout;} - if (perm & 0x8) {permute(inout,tmp,3); tmp=inout;} -} // Make the logger work like Python print() template @@ -46,7 +39,11 @@ struct fatParams: Serializable { } }; - +// +// one method: input --> fat +// another : input --> long (naik) +// another : input --> unitarize +// int main (int argc, char **argv) { @@ -66,119 +63,15 @@ int main (int argc, char **argv) fatParams param(Reader); LatticeGaugeField Umu(&GRID); + LatticeGaugeField U_smr(&GRID); FieldMetaData header; NerscIO::readConfiguration(Umu, header, param.conf_in); - // Create a padded cell of extra padding depth=1 - int depth = 1; - PaddedCell Ghost(depth,&GRID); - LatticeGaugeField Ughost = Ghost.Exchange(Umu); + Smear_HISQ_3link hisq_3link(&GRID); - // Array for (x) - GridBase *GhostGrid = Ughost.Grid(); - LatticeComplex gplaq(GhostGrid); + hisq_3link.smear(U_smr,Umu); - // This is where the 3-link constructs will be stored - LatticeGaugeField Ughost_3link(Ughost.Grid()); - - // Create 3-link stencil (class will build its own stencils) - // writing your own stencil, you're hard-coding the periodic BCs, so you don't need - // the policy-based stuff, at least for now - std::vector shifts; - for(int mu=0;mu_offset; - int o1 = SE1->_offset; - int o2 = SE2->_offset; - int o3 = SE3->_offset; - int o4 = SE4->_offset; - int o5 = SE5->_offset; - - auto U0 = U_v[o0](nu); - auto U1 = adj(U_v[o1](mu)); - auto U2 = adj(U_v[o2](nu)); - - gpermute(U0,SE0->_permute); - gpermute(U1,SE1->_permute); - gpermute(U2,SE2->_permute); - - auto U3 = adj(U_v[o3](nu)); - auto U4 = adj(U_v[o4](mu)); - auto U5 = U_v[o5](nu); - - gpermute(U3,SE3->_permute); - gpermute(U4,SE4->_permute); - gpermute(U5,SE5->_permute); - - // Forward contribution from this orientation - auto W = U0*U1*U2; - U_3link_v[ss](mu) = U_3link_v[ss](mu) + W; - - // Backward contribution from this orientation - W = U3*U4*U5; - U_3link_v[ss](mu) = U_3link_v[ss](mu) + W; - - s=s+6; - } - } - } - - // 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. - Umu = Ghost.Extract(Ughost_3link); - - NerscIO::writeConfiguration(Umu,param.conf_out,"HISQ"); + NerscIO::writeConfiguration(U_smr,param.conf_out,"HISQ"); Grid_finalize(); } \ No newline at end of file From 1cf9ec1cce35748763df382532a611bc5386e1f5 Mon Sep 17 00:00:00 2001 From: david clarke Date: Fri, 9 Jun 2023 16:27:45 -0600 Subject: [PATCH 05/40] now compiles --- Grid/qcd/smearing/HISQSmearing.h | 164 +++++++++++++++++++++++++++---- tests/smearing/Test_fatLinks.cc | 4 +- 2 files changed, 149 insertions(+), 19 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 53480f25..95187d48 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -51,25 +51,21 @@ template void gpermute(vobj & inout,int perm) { } -/*! @brief 3-link smearing of link variable. */ -//template -//class Smear_HISQ_3link : public Smear { -class Smear_HISQ_3link { -// TODO: I'm not using Gimpl so I don't know how to inherit + +/*! @brief create fat links from link variables. */ +class Smear_HISQ_fat { private: -// std::vector _linkTreatment; - GridBase* const _grid; + GridCartesian* const _grid; public: -// INHERIT_GIMPL_TYPES(Gimpl) // Eventually this will take, e.g., coefficients as argument - Smear_HISQ_3link(GridBase* grid) : _grid(grid) { + Smear_HISQ_fat(GridCartesian* grid) : _grid(grid) { assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3"); } - ~Smear_HISQ_3link() {} + ~Smear_HISQ_fat() {} void smear(LatticeGaugeField& u_smr, const LatticeGaugeField& U) const { @@ -83,7 +79,7 @@ public: LatticeComplex gplaq(GhostGrid); // This is where the 3-link constructs will be stored - LatticeGaugeField Ughost_3link(Ughost.Grid()); + LatticeGaugeField Ughost_fat(Ughost.Grid()); // Create 3-link stencil (class will build its own stencils) // writing your own stencil, you're hard-coding the periodic BCs, so you don't need @@ -110,11 +106,11 @@ public: } GeneralLocalStencil gStencil(GhostGrid,shifts); - Ughost_3link=Zero(); + Ughost_fat=Zero(); - // Create the accessors, here U_v and U_3link_v + // Create the accessors, here U_v and U_fat_v autoView(U_v , Ughost , CpuRead); - autoView(U_3link_v, Ughost_3link, CpuWrite); + autoView(U_fat_v, Ughost_fat, CpuWrite); // This is a loop over local sites. for(int ss=0;ss_grid); + LatticeGaugeField Ughost = Ghost.Exchange(u_smr); + + GridBase *GhostGrid = Ughost.Grid(); + LatticeComplex gplaq(GhostGrid); + + LatticeGaugeField Ughost_naik(Ughost.Grid()); + + std::vector shifts; + for(int mu=0;mu_offset; + int o1 = SE1->_offset; + int o2 = SE2->_offset; + int o3 = SE3->_offset; + int o4 = SE4->_offset; + int o5 = SE5->_offset; + + auto U0 = U_v[o0](nu); + auto U1 = adj(U_v[o1](mu)); + auto U2 = adj(U_v[o2](nu)); + + gpermute(U0,SE0->_permute); + gpermute(U1,SE1->_permute); + gpermute(U2,SE2->_permute); + + auto U3 = adj(U_v[o3](nu)); + auto U4 = adj(U_v[o4](mu)); + auto U5 = U_v[o5](nu); + + gpermute(U3,SE3->_permute); + gpermute(U4,SE4->_permute); + gpermute(U5,SE5->_permute); + + // Forward contribution from this orientation + auto W = U0*U1*U2; + U_naik_v[ss](mu) = U_naik_v[ss](mu) + W; + + // Backward contribution from this orientation + W = U3*U4*U5; + U_naik_v[ss](mu) = U_naik_v[ss](mu) + W; + + s=s+6; + } + } + } + + // 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. + u_smr = Ghost.Extract(Ughost_naik); }; // void derivative(const GaugeField& Gauge) const { // }; }; + NAMESPACE_END(Grid); diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index 124cfac3..e87c48e0 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -67,9 +67,9 @@ int main (int argc, char **argv) FieldMetaData header; NerscIO::readConfiguration(Umu, header, param.conf_in); - Smear_HISQ_3link hisq_3link(&GRID); + Smear_HISQ_fat hisq_fat(&GRID); - hisq_3link.smear(U_smr,Umu); + hisq_fat.smear(U_smr,Umu); NerscIO::writeConfiguration(U_smr,param.conf_out,"HISQ"); From 14d352ea4f5a5d513fd21165a31f599e59680012 Mon Sep 17 00:00:00 2001 From: david clarke Date: Mon, 12 Jun 2023 16:55:44 -0600 Subject: [PATCH 06/40] added smearParams struct --- Grid/qcd/smearing/HISQSmearing.h | 33 +++++++++++++++++++++++++++++--- tests/smearing/Test_fatLinks.cc | 21 +++++++++++++------- 2 files changed, 44 insertions(+), 10 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 95187d48..9fae32fd 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -51,17 +51,44 @@ template void gpermute(vobj & inout,int perm) { } +/*! @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. */ + +/*! @brief create fat links from link variables */ class Smear_HISQ_fat { private: GridCartesian* const _grid; + SmearingParameters _LVL1; public: - // Eventually this will take, e.g., coefficients as argument - Smear_HISQ_fat(GridCartesian* grid) : _grid(grid) { + Smear_HISQ_fat(GridCartesian* grid, Real c1=1/8., Real cnaik=0., Real c3=1/16., Real c5=1/64., Real c7=1/384., Real clp=0.) + : _grid(grid), + _LVL1(c1,cnaik,c3,c5,c7,clp) { + assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3"); + } + + // Allow to pass a pointer to a C-style, double array for MILC convenience + Smear_HISQ_fat(GridCartesian* grid, double* coeff) + : _grid(grid), + _LVL1(coeff[0],coeff[1],coeff[2],coeff[3],coeff[4],coeff[5]) { assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3"); } diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index e87c48e0..1854ea71 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -47,31 +47,38 @@ struct fatParams: Serializable { int main (int argc, char **argv) { + // Initialize the Grid Grid_init(&argc,&argv); - Coordinate latt_size = GridDefaultLatt(); 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); - GridCartesian GRID(latt_size,simd_layout,mpi_layout); - XmlReader Reader("fatParams.xml",false, "grid"); - fatParams param(Reader); - + // Instantiate the LatticeGaugeField objects holding thin (Umu) and fat (U_smr) links LatticeGaugeField Umu(&GRID); LatticeGaugeField U_smr(&GRID); + + // Read in the parameter file + XmlReader Reader("fatParams.xml",false, "grid"); + fatParams param(Reader); FieldMetaData header; + + // Read the configuration into Umu NerscIO::readConfiguration(Umu, header, param.conf_in); + // Smear Umu and store result in U_smr Smear_HISQ_fat hisq_fat(&GRID); - hisq_fat.smear(U_smr,Umu); NerscIO::writeConfiguration(U_smr,param.conf_out,"HISQ"); + // In the following, we test instantiation of Smear_HISQ_fat in different ways: + Smear_HISQ_fat hisq_fat_typical(&GRID,1,2,3,4,5,6); + double path_coeff[6] = {1, 2, 3, 4, 5, 6}; + Smear_HISQ_fat hisq_fat_Cstyle(&GRID,path_coeff); + Grid_finalize(); } \ No newline at end of file From 26b2caf5706b482c0461d52fb03a60e0c20b1127 Mon Sep 17 00:00:00 2001 From: david clarke Date: Tue, 20 Jun 2023 15:37:54 -0600 Subject: [PATCH 07/40] add template parameter to Smear_HISQ_fat for MILC interfacing --- Grid/qcd/smearing/HISQSmearing.h | 38 ++++++++++++++------------------ tests/smearing/Test_fatLinks.cc | 7 +++--- 2 files changed, 20 insertions(+), 25 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 9fae32fd..b8fc5aef 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -71,6 +71,7 @@ struct SmearingParameters{ /*! @brief create fat links from link variables */ +template // TODO: change to Gimpl? class Smear_HISQ_fat { private: @@ -79,7 +80,8 @@ private: public: - Smear_HISQ_fat(GridCartesian* grid, Real c1=1/8., Real cnaik=0., Real c3=1/16., Real c5=1/64., Real c7=1/384., Real clp=0.) + // Don't allow default values here. + Smear_HISQ_fat(GridCartesian* grid, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) : _grid(grid), _LVL1(c1,cnaik,c3,c5,c7,clp) { assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3"); @@ -94,19 +96,21 @@ public: ~Smear_HISQ_fat() {} - void smear(LatticeGaugeField& u_smr, const LatticeGaugeField& U) const { - + void smear(LGF& u_smr, LGF& u_thin) const { + + SmearingParameters lvl1 = this->_LVL1; + // Create a padded cell of extra padding depth=1 int depth = 1; PaddedCell Ghost(depth,this->_grid); - LatticeGaugeField Ughost = Ghost.Exchange(u_smr); + LGF Ughost = Ghost.Exchange(u_thin); // Array for (x) GridBase *GhostGrid = Ughost.Grid(); LatticeComplex gplaq(GhostGrid); // This is where the 3-link constructs will be stored - LatticeGaugeField Ughost_fat(Ughost.Grid()); + LGF Ughost_fat(Ughost.Grid()); // Create 3-link stencil (class will build its own stencils) // writing your own stencil, you're hard-coding the periodic BCs, so you don't need @@ -136,7 +140,7 @@ public: Ughost_fat=Zero(); // Create the accessors, here U_v and U_fat_v - autoView(U_v , Ughost , CpuRead); + autoView(U_v , Ughost , CpuRead); autoView(U_fat_v, Ughost_fat, CpuWrite); // This is a loop over local sites. @@ -186,12 +190,8 @@ public: gpermute(U4,SE4->_permute); gpermute(U5,SE5->_permute); - // Forward contribution from this orientation - auto W = U0*U1*U2; - U_fat_v[ss](mu) = U_fat_v[ss](mu) + W; - - // Backward contribution from this orientation - W = U3*U4*U5; + // forward backward + auto W = U0*U1*U2 + U3*U4*U5; U_fat_v[ss](mu) = U_fat_v[ss](mu) + W; s=s+6; @@ -199,14 +199,9 @@ public: } } - // 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. - u_smr = Ghost.Extract(Ughost_fat); + u_smr = lvl1.c_3*Ghost.Extract(Ughost_fat) + lvl1.c_1*u_thin; }; - // I guess the way this will go is: // 1. 3-link smear // 2. exchange @@ -219,6 +214,7 @@ public: /*! @brief create long links from link variables. */ +template class Smear_HISQ_Naik { private: @@ -233,16 +229,16 @@ public: ~Smear_HISQ_Naik() {} - void smear(LatticeGaugeField& u_smr, const LatticeGaugeField& U) const { + void smear(LGF& u_smr, const LGF& U) const { int depth = 1; PaddedCell Ghost(depth,this->_grid); - LatticeGaugeField Ughost = Ghost.Exchange(u_smr); + LGF Ughost = Ghost.Exchange(u_smr); GridBase *GhostGrid = Ughost.Grid(); LatticeComplex gplaq(GhostGrid); - LatticeGaugeField Ughost_naik(Ughost.Grid()); + LGF Ughost_naik(Ughost.Grid()); std::vector shifts; for(int mu=0;mu hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,0.); hisq_fat.smear(U_smr,Umu); NerscIO::writeConfiguration(U_smr,param.conf_out,"HISQ"); - // In the following, we test instantiation of Smear_HISQ_fat in different ways: - Smear_HISQ_fat hisq_fat_typical(&GRID,1,2,3,4,5,6); + // Test a C-style instantiation double path_coeff[6] = {1, 2, 3, 4, 5, 6}; - Smear_HISQ_fat hisq_fat_Cstyle(&GRID,path_coeff); + Smear_HISQ_fat hisq_fat_Cstyle(&GRID,path_coeff); Grid_finalize(); } \ No newline at end of file From f44f005dad99325f52e2390033f38fb6cbc30181 Mon Sep 17 00:00:00 2001 From: david clarke Date: Tue, 20 Jun 2023 15:48:27 -0600 Subject: [PATCH 08/40] rename _lvl1 --> _linkTreatment --- Grid/qcd/smearing/HISQSmearing.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index b8fc5aef..1e7e93ab 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -76,21 +76,21 @@ class Smear_HISQ_fat { private: GridCartesian* const _grid; - SmearingParameters _LVL1; + SmearingParameters _linkTreatment; public: // Don't allow default values here. Smear_HISQ_fat(GridCartesian* grid, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) : _grid(grid), - _LVL1(c1,cnaik,c3,c5,c7,clp) { + _linkTreatment(c1,cnaik,c3,c5,c7,clp) { assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3"); } // Allow to pass a pointer to a C-style, double array for MILC convenience Smear_HISQ_fat(GridCartesian* grid, double* coeff) : _grid(grid), - _LVL1(coeff[0],coeff[1],coeff[2],coeff[3],coeff[4],coeff[5]) { + _linkTreatment(coeff[0],coeff[1],coeff[2],coeff[3],coeff[4],coeff[5]) { assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3"); } @@ -98,7 +98,7 @@ public: void smear(LGF& u_smr, LGF& u_thin) const { - SmearingParameters lvl1 = this->_LVL1; + SmearingParameters lt = this->_linkTreatment; // Create a padded cell of extra padding depth=1 int depth = 1; @@ -199,7 +199,7 @@ public: } } - u_smr = lvl1.c_3*Ghost.Extract(Ughost_fat) + lvl1.c_1*u_thin; + u_smr = lt.c_3*Ghost.Extract(Ughost_fat) + lt.c_1*u_thin; }; // I guess the way this will go is: From d536c67b9d407617bf8d2f17c5b9efa95bad8d00 Mon Sep 17 00:00:00 2001 From: david clarke Date: Tue, 20 Jun 2023 16:04:48 -0600 Subject: [PATCH 09/40] add HISQSmearing to Smearing.h --- Grid/qcd/smearing/Smearing.h | 1 + 1 file changed, 1 insertion(+) diff --git a/Grid/qcd/smearing/Smearing.h b/Grid/qcd/smearing/Smearing.h index da5ede72..41a305ae 100644 --- a/Grid/qcd/smearing/Smearing.h +++ b/Grid/qcd/smearing/Smearing.h @@ -5,4 +5,5 @@ #include #include #include +#include From df99f227c16a152b648ec3b86122c02443bb04e7 Mon Sep 17 00:00:00 2001 From: david clarke Date: Thu, 22 Jun 2023 14:57:10 -0600 Subject: [PATCH 10/40] include missing staple orientations; invert path direction, which was backwards --- Grid/qcd/smearing/HISQSmearing.h | 185 ++++++++++++++++--------------- 1 file changed, 95 insertions(+), 90 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 1e7e93ab..94b47e52 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -104,101 +104,106 @@ public: int depth = 1; PaddedCell Ghost(depth,this->_grid); LGF Ughost = Ghost.Exchange(u_thin); - + // Array for (x) GridBase *GhostGrid = Ughost.Grid(); LatticeComplex gplaq(GhostGrid); - + // This is where the 3-link constructs will be stored LGF Ughost_fat(Ughost.Grid()); - // Create 3-link stencil (class will build its own stencils) - // writing your own stencil, you're hard-coding the periodic BCs, so you don't need - // the policy-based stuff, at least for now + // Create 3-link stencil. Writing your own stencil, you're hard-coding the + // periodic BCs, so you don't need the policy-based stuff, at least for now. + // Loop over all orientations, i.e. demand mu != nu. std::vector shifts; - for(int mu=0;mu_offset; - int o1 = SE1->_offset; - int o2 = SE2->_offset; - int o3 = SE3->_offset; - int o4 = SE4->_offset; - int o5 = SE5->_offset; - - auto U0 = U_v[o0](nu); - auto U1 = adj(U_v[o1](mu)); - auto U2 = adj(U_v[o2](nu)); - - gpermute(U0,SE0->_permute); - gpermute(U1,SE1->_permute); - gpermute(U2,SE2->_permute); - - auto U3 = adj(U_v[o3](nu)); - auto U4 = adj(U_v[o4](mu)); - auto U5 = U_v[o5](nu); - - gpermute(U3,SE3->_permute); - gpermute(U4,SE4->_permute); - gpermute(U5,SE5->_permute); - - // forward backward - auto W = U0*U1*U2 + U3*U4*U5; - U_fat_v[ss](mu) = U_fat_v[ss](mu) + W; - - s=s+6; - } + + for(int mu=0;mu_offset; + int o1 = SE1->_offset; + int o2 = SE2->_offset; + int o3 = SE3->_offset; + int o4 = SE4->_offset; + int o5 = 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. + auto U0 = adj(U_v[o0](nu)); + auto U1 = U_v[o1](mu); + auto U2 = U_v[o2](nu); + + gpermute(U0,SE0->_permute); + gpermute(U1,SE1->_permute); + gpermute(U2,SE2->_permute); + + auto U3 = U_v[o3](nu); + auto U4 = U_v[o4](mu); + auto U5 = adj(U_v[o5](nu)); + + gpermute(U3,SE3->_permute); + gpermute(U4,SE4->_permute); + gpermute(U5,SE5->_permute); + + // "left" "right" + auto W = U2*U1*U0 + U5*U4*U3; + U_fat_v[ss](mu) = U_fat_v[ss](mu) + W; + + s=s+6; } } - + u_smr = lt.c_3*Ghost.Extract(Ughost_fat) + lt.c_1*u_thin; }; @@ -261,23 +266,23 @@ public: } } GeneralLocalStencil gStencil(GhostGrid,shifts); - + Ughost_naik=Zero(); - + // Create the accessors, here U_v and U_fat_v autoView(U_v , Ughost , CpuRead); autoView(U_naik_v, Ughost_naik, CpuWrite); - + // This is a loop over local sites. for(int ss=0;ss_offset; int o1 = SE1->_offset; @@ -298,36 +303,36 @@ public: int o3 = SE3->_offset; int o4 = SE4->_offset; int o5 = SE5->_offset; - + auto U0 = U_v[o0](nu); auto U1 = adj(U_v[o1](mu)); auto U2 = adj(U_v[o2](nu)); - + gpermute(U0,SE0->_permute); gpermute(U1,SE1->_permute); gpermute(U2,SE2->_permute); - + auto U3 = adj(U_v[o3](nu)); auto U4 = adj(U_v[o4](mu)); auto U5 = U_v[o5](nu); - + gpermute(U3,SE3->_permute); gpermute(U4,SE4->_permute); gpermute(U5,SE5->_permute); - + // Forward contribution from this orientation auto W = U0*U1*U2; U_naik_v[ss](mu) = U_naik_v[ss](mu) + W; - + // Backward contribution from this orientation W = U3*U4*U5; U_naik_v[ss](mu) = U_naik_v[ss](mu) + W; - + s=s+6; } } } - + // 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 From eeb4703b84b723722ad654e834ad0e7af8a637c3 Mon Sep 17 00:00:00 2001 From: david clarke Date: Mon, 26 Jun 2023 17:45:35 -0600 Subject: [PATCH 11/40] develop wrappers to make the stencils easier to construct --- Grid/qcd/smearing/HISQSmearing.h | 49 ++++++++++++++++---------------- tests/smearing/Test_fatLinks.cc | 7 +++++ 2 files changed, 31 insertions(+), 25 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 94b47e52..d36fd85f 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -51,8 +51,17 @@ template void gpermute(vobj & inout,int perm) { } +void appendShift(std::vector& shifts, int mu, int steps=1) { + Coordinate shift(Nd,0); + shift[mu]=steps; + // push_back creates an element at the end of shifts and + // assigns the data in the argument to it. + shifts.push_back(shift); +} + + /*! @brief structure holding the link treatment */ -struct SmearingParameters{ +struct SmearingParameters { SmearingParameters(){} Real c_1; // 1 link Real c_naik; // Naik term @@ -100,8 +109,10 @@ public: SmearingParameters lt = this->_linkTreatment; - // Create a padded cell of extra padding depth=1 - int depth = 1; + // We create a cell with extra padding 2. This allows us to capture the LePage + // term without needing to save intermediate gauge fields or extra halo exchanges. + // The tradeoff is that we compute extra constructs in the padding. + int depth = 2; PaddedCell Ghost(depth,this->_grid); LGF Ughost = Ghost.Exchange(u_thin); @@ -112,28 +123,20 @@ public: // This is where the 3-link constructs will be stored LGF Ughost_fat(Ughost.Grid()); - // Create 3-link stencil. Writing your own stencil, you're hard-coding the + // Next we make the stencils. Writing your own stencil, you're hard-coding the // periodic BCs, so you don't need the policy-based stuff, at least for now. // Loop over all orientations, i.e. demand mu != nu. std::vector shifts; for(int mu=0;mu hisq_fat_Cstyle(&GRID,path_coeff); + // Make sure result doesn't change w.r.t. a trusted lattice + NerscIO::readConfiguration(Umu, header, "nersc.l8t4b3360.3link.control"); + LatticeGaugeField diff(&GRID); + diff = Umu-U_smr; + auto absDiff = norm2(diff)/norm2(Umu); + Grid_log(" |Umu-U|/|Umu| =",absDiff); + Grid_finalize(); } \ No newline at end of file From a7eabaad566cd2b740ad98e370eb7800f9394039 Mon Sep 17 00:00:00 2001 From: david clarke Date: Mon, 26 Jun 2023 23:59:28 -0600 Subject: [PATCH 12/40] rudimentary appendShift convenience method, which allows the user to append an arbitrary shift in one line --- Grid/qcd/smearing/HISQSmearing.h | 76 +++++++++++++++++++++++--------- tests/smearing/Test_fatLinks.cc | 4 +- 2 files changed, 56 insertions(+), 24 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index d36fd85f..e8587c9b 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -30,7 +30,6 @@ directory @brief Declares classes related to HISQ smearing */ -// things like @brief are seen by things like doxygen and javadocs #pragma once @@ -38,6 +37,9 @@ directory #include #include +#define BACKWARD_CONST 16 +#define NO_SHIFT -1 + NAMESPACE_BEGIN(Grid); @@ -51,9 +53,47 @@ template void gpermute(vobj & inout,int perm) { } -void appendShift(std::vector& shifts, int mu, int steps=1) { - Coordinate shift(Nd,0); - shift[mu]=steps; +/*! @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. Should work as long as BACKWARD_CONST > Nd. + return dir + BACKWARD_CONST; +} + + +/*! @brief shift one unit in direction dir */ +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 +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...); +} + + +/*! @brief append arbitrary shift path to shifts */ +template +void appendShift(std::vector& 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); @@ -123,22 +163,17 @@ public: // This is where the 3-link constructs will be stored LGF Ughost_fat(Ughost.Grid()); - // Next we make the stencils. Writing your own stencil, you're hard-coding the - // periodic BCs, so you don't need the policy-based stuff, at least for now. - // Loop over all orientations, i.e. demand mu != nu. + // Create a stencil, which is a collection of sites neighboring some initial site. std::vector shifts; for(int mu=0;mu_offset; @@ -177,7 +210,6 @@ public: int o2 = SE2->_offset; int o3 = SE3->_offset; int o4 = SE4->_offset; - int o5 = 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 @@ -193,17 +225,17 @@ public: auto U3 = U_v[o3](nu); auto U4 = U_v[o4](mu); - auto U5 = adj(U_v[o5](nu)); + auto U5 = adj(U_v[o4](nu)); gpermute(U3,SE3->_permute); gpermute(U4,SE4->_permute); - gpermute(U5,SE5->_permute); + gpermute(U4,SE4->_permute); // "left" "right" auto W = U2*U1*U0 + U5*U4*U3; U_fat_v[ss](mu) = U_fat_v[ss](mu) + W; - s=s+6; + s=s+5; } } diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index aa8e8f92..de1796b7 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -15,13 +15,13 @@ using namespace Grid; // Make the logger work like Python print() -template +template inline std::string sjoin(Args&&... args) noexcept { std::ostringstream msg; (msg << ... << args); return msg.str(); } -template +template inline void Grid_log(Args&&... args) { std::string msg = sjoin(std::forward(args)...); std::cout << GridLogMessage << msg << std::endl; From 9015c229dc194fdf2b6552e6137b7ffbe23c6e78 Mon Sep 17 00:00:00 2001 From: david clarke Date: Tue, 27 Jun 2023 21:28:26 -0600 Subject: [PATCH 13/40] add benchmark to see whether matrix multiplication is slower than read from object --- Grid/qcd/smearing/HISQSmearing.h | 149 +++------------------- benchmarks/Benchmark_su3mult_vs_lookup.cc | 125 ++++++++++++++++++ tests/smearing/Test_fatLinks.cc | 78 ++++++----- 3 files changed, 188 insertions(+), 164 deletions(-) create mode 100644 benchmarks/Benchmark_su3mult_vs_lookup.cc diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index e8587c9b..0c11bde4 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -2,9 +2,9 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: ./lib/qcd/smearing/StoutSmearing.h +Source file: ./lib/qcd/smearing/HISQSmearing.h -Copyright (C) 2019 +Copyright (C) 2023 Author: D. A. Clarke @@ -56,7 +56,7 @@ template void gpermute(vobj & inout,int perm) { /*! @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. Should work as long as BACKWARD_CONST > Nd. + // backward. Should work as long as BACKWARD_CONST > Nd. Trick inspired by SIMULATeQCD. return dir + BACKWARD_CONST; } @@ -191,45 +191,33 @@ public: for(int mu=0;mu_offset; - int o1 = SE1->_offset; - int o2 = SE2->_offset; - int o3 = SE3->_offset; - int o4 = SE4->_offset; + auto SE0 = gStencil.GetEntry(s+0,ss); int x_p_mu = SE0->_offset; + auto SE1 = gStencil.GetEntry(s+1,ss); int x_p_nu = SE1->_offset; + auto SE2 = gStencil.GetEntry(s+2,ss); int x = SE2->_offset; + auto SE3 = gStencil.GetEntry(s+3,ss); int x_p_mu_m_nu = SE3->_offset; + auto SE4 = gStencil.GetEntry(s+4,ss); int x_m_nu = SE4->_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. - auto U0 = adj(U_v[o0](nu)); - auto U1 = U_v[o1](mu); - auto U2 = U_v[o2](nu); + auto U0 = adj(U_v[x_p_mu](nu)); + auto U1 = U_v[x_p_nu](mu) ; + auto U2 = U_v[x ](nu) ; gpermute(U0,SE0->_permute); gpermute(U1,SE1->_permute); gpermute(U2,SE2->_permute); - auto U3 = U_v[o3](nu); - auto U4 = U_v[o4](mu); - auto U5 = adj(U_v[o4](nu)); + auto U3 = U_v[x_p_mu_m_nu](nu) ; + auto U4 = U_v[x_m_nu ](mu) ; + auto U5 = adj(U_v[x_m_nu ](nu)); gpermute(U3,SE3->_permute); gpermute(U4,SE4->_permute); - gpermute(U4,SE4->_permute); // "left" "right" auto W = U2*U1*U0 + U5*U4*U3; @@ -265,111 +253,8 @@ public: ~Smear_HISQ_Naik() {} - void smear(LGF& u_smr, const LGF& U) const { - - int depth = 1; - PaddedCell Ghost(depth,this->_grid); - LGF Ughost = Ghost.Exchange(u_smr); - - GridBase *GhostGrid = Ughost.Grid(); - LatticeComplex gplaq(GhostGrid); - - LGF Ughost_naik(Ughost.Grid()); - - std::vector shifts; - for(int mu=0;mu_offset; - int o1 = SE1->_offset; - int o2 = SE2->_offset; - int o3 = SE3->_offset; - int o4 = SE4->_offset; - int o5 = SE5->_offset; - - auto U0 = U_v[o0](nu); - auto U1 = adj(U_v[o1](mu)); - auto U2 = adj(U_v[o2](nu)); - - gpermute(U0,SE0->_permute); - gpermute(U1,SE1->_permute); - gpermute(U2,SE2->_permute); - - auto U3 = adj(U_v[o3](nu)); - auto U4 = adj(U_v[o4](mu)); - auto U5 = U_v[o5](nu); - - gpermute(U3,SE3->_permute); - gpermute(U4,SE4->_permute); - gpermute(U5,SE5->_permute); - - // Forward contribution from this orientation - auto W = U0*U1*U2; - U_naik_v[ss](mu) = U_naik_v[ss](mu) + W; - - // Backward contribution from this orientation - W = U3*U4*U5; - U_naik_v[ss](mu) = U_naik_v[ss](mu) + W; - - s=s+6; - } - } - } - - // 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. - u_smr = Ghost.Extract(Ughost_naik); - }; +// void smear(LGF& u_smr, const LGF& U) const { +// }; // void derivative(const GaugeField& Gauge) const { // }; diff --git a/benchmarks/Benchmark_su3mult_vs_lookup.cc b/benchmarks/Benchmark_su3mult_vs_lookup.cc new file mode 100644 index 00000000..5f9c98ba --- /dev/null +++ b/benchmarks/Benchmark_su3mult_vs_lookup.cc @@ -0,0 +1,125 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./benchmarks/Benchmark_su3mult_vs_lookup.cc + + Copyright (C) 2023 + + Author: D. A. Clarke + + 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 +using namespace Grid; + +/*! @brief make the logger work like python print */ +template +inline std::string sjoin(Args&&... args) noexcept { + std::ostringstream msg; + (msg << ... << args); + return msg.str(); +} +template +inline void Grid_log(Args&&... args) { + std::string msg = sjoin(std::forward(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 + ConfParameters(Reader& 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 + +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 #include @@ -14,7 +38,7 @@ using namespace Grid; -// Make the logger work like Python print() +/*! @brief make the logger work like python print */ template inline std::string sjoin(Args&&... args) noexcept { std::ostringstream msg; @@ -27,32 +51,26 @@ inline void Grid_log(Args&&... args) { std::cout << GridLogMessage << msg << std::endl; } -struct fatParams: Serializable { - GRID_SERIALIZABLE_CLASS_MEMBERS( - fatParams, - std::string, conf_in, - std::string, conf_out); - - template - fatParams(Reader& Reader){ - read(Reader, "parameters", *this); - } -}; - // // one method: input --> fat // another : input --> long (naik) // another : input --> unitarize // -int main (int argc, char **argv) -{ +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"; + std::string conf_out = "nersc.l8t4b3360.3link"; + // Initialize the Grid Grid_init(&argc,&argv); - Coordinate latt_size = GridDefaultLatt(); 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("latt = ",latt_size); GridCartesian GRID(latt_size,simd_layout,mpi_layout); @@ -61,19 +79,15 @@ int main (int argc, char **argv) LatticeGaugeField Umu(&GRID); LatticeGaugeField U_smr(&GRID); - // Read in the parameter file - XmlReader Reader("fatParams.xml",false, "grid"); - fatParams param(Reader); - FieldMetaData header; - // Read the configuration into Umu - NerscIO::readConfiguration(Umu, header, param.conf_in); + FieldMetaData header; + NerscIO::readConfiguration(Umu, header, conf_in); // Smear Umu and store result in U_smr Smear_HISQ_fat hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,0.); hisq_fat.smear(U_smr,Umu); - NerscIO::writeConfiguration(U_smr,param.conf_out,"HISQ"); + NerscIO::writeConfiguration(U_smr,conf_out,"HISQ"); // Test a C-style instantiation double path_coeff[6] = {1, 2, 3, 4, 5, 6}; From 9d263d9a7d7268fd8e8e909e76cd19ead9a6fa0c Mon Sep 17 00:00:00 2001 From: david clarke Date: Wed, 28 Jun 2023 10:05:34 -0600 Subject: [PATCH 14/40] fix bug in HISQSmearing; move benchmark b/c i don't understand how makefiles work --- Grid/qcd/smearing/HISQSmearing.h | 42 ++++++++----------- .../smearing}/Benchmark_su3mult_vs_lookup.cc | 0 tests/smearing/Test_fatLinks.cc | 2 +- 3 files changed, 18 insertions(+), 26 deletions(-) rename {benchmarks => tests/smearing}/Benchmark_su3mult_vs_lookup.cc (100%) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 0c11bde4..cba761f4 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -62,6 +62,7 @@ inline int Back(const int dir) { /*! @brief shift one unit in direction dir */ +template void generalShift(Coordinate& shift, int dir) { if (dir >= BACKWARD_CONST) { dir -= BACKWARD_CONST; @@ -101,7 +102,7 @@ void appendShift(std::vector& shifts, int dir, Args... args) { /*! @brief structure holding the link treatment */ -struct SmearingParameters { +struct SmearingParameters{ SmearingParameters(){} Real c_1; // 1 link Real c_naik; // Naik term @@ -149,10 +150,8 @@ public: SmearingParameters lt = this->_linkTreatment; - // We create a cell with extra padding 2. This allows us to capture the LePage - // term without needing to save intermediate gauge fields or extra halo exchanges. - // The tradeoff is that we compute extra constructs in the padding. - int depth = 2; + // Create a padded cell of extra padding depth=1 + int depth = 1; PaddedCell Ghost(depth,this->_grid); LGF Ughost = Ghost.Exchange(u_thin); @@ -163,7 +162,9 @@ public: // This is where the 3-link constructs will be stored LGF Ughost_fat(Ughost.Grid()); - // Create a stencil, which is a collection of sites neighboring some initial site. + // Create 3-link stencil. Writing your own stencil, you're hard-coding the + // periodic BCs, so you don't need the policy-based stuff, at least for now. + // Loop over all orientations, i.e. demand mu != nu. std::vector shifts; for(int mu=0;mu_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); - gpermute(U0,SE0->_permute); - gpermute(U1,SE1->_permute); - gpermute(U2,SE2->_permute); - - auto U3 = U_v[x_p_mu_m_nu](nu) ; - auto U4 = U_v[x_m_nu ](mu) ; - auto U5 = adj(U_v[x_m_nu ](nu)); - - gpermute(U3,SE3->_permute); - gpermute(U4,SE4->_permute); - - // "left" "right" - auto W = U2*U1*U0 + U5*U4*U3; + // "left" "right" + auto W = U2*U1*adj(U0) + adj(U5)*U4*U3; U_fat_v[ss](mu) = U_fat_v[ss](mu) + W; s=s+5; @@ -235,8 +229,6 @@ public: // }; }; - - /*! @brief create long links from link variables. */ template class Smear_HISQ_Naik { @@ -261,4 +253,4 @@ public: }; -NAMESPACE_END(Grid); +NAMESPACE_END(Grid); \ No newline at end of file diff --git a/benchmarks/Benchmark_su3mult_vs_lookup.cc b/tests/smearing/Benchmark_su3mult_vs_lookup.cc similarity index 100% rename from benchmarks/Benchmark_su3mult_vs_lookup.cc rename to tests/smearing/Benchmark_su3mult_vs_lookup.cc diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index 5bfdd891..5a19e835 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -98,7 +98,7 @@ int main (int argc, char** argv) { LatticeGaugeField diff(&GRID); diff = Umu-U_smr; auto absDiff = norm2(diff)/norm2(Umu); - Grid_log(" |Umu-U|/|Umu| =",absDiff); + Grid_log(" |Umu-U|/|Umu| = ",absDiff); Grid_finalize(); } \ No newline at end of file From 99d879ea7ff4dda0309311529b58e185cba3ac7c Mon Sep 17 00:00:00 2001 From: david clarke Date: Fri, 11 Aug 2023 22:56:30 -0600 Subject: [PATCH 15/40] 5-link first attempt --- Grid/qcd/smearing/HISQSmearing.h | 125 ++++++++++++------ tests/smearing/Benchmark_su3mult_vs_lookup.cc | 125 ------------------ tests/smearing/Test_fatLinks.cc | 65 ++++++++- 3 files changed, 145 insertions(+), 170 deletions(-) delete mode 100644 tests/smearing/Benchmark_su3mult_vs_lookup.cc diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index cba761f4..e6bab0f1 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/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 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 (x) - GridBase *GhostGrid = Ughost.Grid(); - LatticeComplex gplaq(GhostGrid); - - // This is where the 3-link constructs will be stored + // This is where auxiliary N-link fields and the final smear 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 - // periodic BCs, so you don't need the policy-based stuff, at least for now. - // Loop over all orientations, i.e. demand mu != nu. + // Create 3-link stencil. We allow mu==nu just to make the indexing easier. + // Shifts with mu==nu will not be used. std::vector shifts; for(int mu=0;mu_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; - auto SE1 = gStencil.GetEntry(s+1,ss); int x_p_nu = SE1->_offset; - auto SE2 = gStencil.GetEntry(s+2,ss); int x = SE2->_offset; - auto SE3 = gStencil.GetEntry(s+3,ss); int x_p_mu_m_nu = SE3->_offset; - auto SE4 = gStencil.GetEntry(s+4,ss); int x_m_nu = SE4->_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. + 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 - // 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. - 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" + auto W = U2*U1*adj(U0) + adj(U5)*U4*U3; - // "left" "right" - auto W = U2*U1*adj(U0) + adj(U5)*U4*U3; - U_fat_v[ss](mu) = U_fat_v[ss](mu) + W; + U_3link_v[site](nu) = 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_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 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() {} diff --git a/tests/smearing/Benchmark_su3mult_vs_lookup.cc b/tests/smearing/Benchmark_su3mult_vs_lookup.cc deleted file mode 100644 index 5f9c98ba..00000000 --- a/tests/smearing/Benchmark_su3mult_vs_lookup.cc +++ /dev/null @@ -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 - - 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 -using namespace Grid; - -/*! @brief make the logger work like python print */ -template -inline std::string sjoin(Args&&... args) noexcept { - std::ostringstream msg; - (msg << ... << args); - return msg.str(); -} -template -inline void Grid_log(Args&&... args) { - std::string msg = sjoin(std::forward(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 - ConfParameters(Reader& 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 + ConfParameters(Reader& 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("simd = ",simd_layout); - Grid_log("latt = ",latt_size); + 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); + // 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 hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,0.); +// Smear_HISQ_fat 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 Date: Sat, 16 Sep 2023 23:18:16 -0600 Subject: [PATCH 16/40] try 7-link --- Grid/qcd/smearing/HISQSmearing.h | 70 +++++++++++++++++++++++++++++--- tests/smearing/Test_fatLinks.cc | 5 +-- 2 files changed, 67 insertions(+), 8 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index e6bab0f1..9a69b660 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -129,7 +129,7 @@ struct SmearingParameters{ /*! @brief create fat links from link variables */ -template // TODO: change to Gimpl? +template class Smear_HISQ_fat { private: @@ -168,6 +168,8 @@ public: // This is where auxiliary N-link fields and the final smear will be stored. LGF Ughost_fat(Ughost.Grid()); LGF Ughost_3link(Ughost.Grid()); + LGF Ughost_5linkA(Ughost.Grid()); + LGF Ughost_5linkB(Ughost.Grid()); // Create 3-link stencil. We allow mu==nu just to make the indexing easier. // Shifts with mu==nu will not be used. @@ -188,13 +190,18 @@ public: Ughost_fat=Zero(); // Create the accessors - autoView(U_v , Ughost , CpuRead); - autoView(U_fat_v , Ughost_fat , CpuWrite); - autoView(U_3link_v, Ughost_3link, CpuWrite); + autoView(U_v , Ughost , CpuRead); + autoView(U_fat_v , Ughost_fat , CpuWrite); + autoView(U_3link_v , Ughost_3link , CpuWrite); + autoView(U_5linkA_v, Ughost_5linkA, CpuWrite); + autoView(U_5linkB_v, Ughost_5linkB, CpuWrite); for(int mu=0;mu_offset; @@ -253,7 +263,57 @@ public: auto W = U2*U1*adj(U0) + adj(U5)*U4*U3; + if(sigmaIndex<3) { + U_5linkA_v[site](rho) = W; + } else { + U_5linkB_v[site](rho) = W; + } + U_fat_v[site](mu) = U_fat_v[site](mu) + lt.c_5*W; + + sigmaIndex++; + } + } + } + + // 7-link + for(int site=0;site_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 = U0; + if(sigmaIndex<3) { + U1 = U_5linkB_v[x_p_nu](rho); gpermute(U1,SE1->_permute); + } else { + U1 = U_5linkA_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 = U0; + if(sigmaIndex<3) { + U4 = U_5linkB_v[x_m_nu](rho); gpermute(U4,SE4->_permute); + } else { + U4 = U_5linkA_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_7*W; + + sigmaIndex++; } } } diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index 97e02400..f7c422c1 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -78,7 +78,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.3link"; + std::string conf_out = "nersc.l8t4b3360.35link"; int threads = GridThread::GetThreads(); // Initialize the Grid @@ -105,7 +105,6 @@ int main (int argc, char** argv) { // Smear Umu and store result in U_smr Smear_HISQ_fat hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,0.); -// Smear_HISQ_fat 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"); @@ -115,7 +114,7 @@ int main (int argc, char** argv) { Smear_HISQ_fat hisq_fat_Cstyle(&GRID,path_coeff); // Make sure result doesn't change w.r.t. a trusted lattice - NerscIO::readConfiguration(Umu, header, "nersc.l8t4b3360.3link.control"); + NerscIO::readConfiguration(Umu, header, "nersc.l8t4b3360.35link.control"); LatticeGaugeField diff(&GRID); diff = Umu-U_smr; auto absDiff = norm2(diff)/norm2(Umu); From 0cfd13d18b41978afbdd8ba4b5ccbbce30e9dc09 Mon Sep 17 00:00:00 2001 From: david clarke Date: Tue, 10 Oct 2023 22:41:52 -0600 Subject: [PATCH 17/40] 7-link working --- Grid/qcd/smearing/HISQSmearing.h | 29 ++++++++++++++++------------- tests/smearing/Test_fatLinks.cc | 4 ++-- 2 files changed, 18 insertions(+), 15 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 9a69b660..8c60d874 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -82,7 +82,7 @@ void generalShift(Coordinate& shift, int dir) { } } - +// Move into general stencil header, beneath definition of general stencil /*! @brief follow a path of directions, shifting one unit in each direction */ template void generalShift(Coordinate& shift, int dir, Args... args) { @@ -118,13 +118,13 @@ struct SmearingParameters{ 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){} + 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){} }; @@ -198,11 +198,12 @@ public: for(int mu=0;mu_offset; auto SE1 = gStencil.GetEntry(s+1,site); int x_p_nu = SE1->_offset; @@ -254,6 +254,8 @@ public: 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; + // gpermutes will be replaced with single line of code, combines load and permute + // into one step. still in pull request stage 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); @@ -283,8 +285,7 @@ public: if(nu==mu) continue; int s = stencilIndex(mu,nu); for(int rho=0;rho_offset; auto SE1 = gStencil.GetEntry(s+1,site); int x_p_nu = SE1->_offset; @@ -293,6 +294,7 @@ public: 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); + // decltype, or auto U1 = { ? ... } auto U1 = U0; if(sigmaIndex<3) { U1 = U_5linkB_v[x_p_nu](rho); gpermute(U1,SE1->_permute); @@ -311,6 +313,7 @@ public: auto W = U2*U1*adj(U0) + adj(U5)*U4*U3; + // std::vector(3) ? U_fat_v[site](mu) = U_fat_v[site](mu) + lt.c_7*W; sigmaIndex++; diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index f7c422c1..200e2af6 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -78,7 +78,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.35link"; + std::string conf_out = "nersc.l8t4b3360.357link"; int threads = GridThread::GetThreads(); // Initialize the Grid @@ -114,7 +114,7 @@ int main (int argc, char** argv) { Smear_HISQ_fat hisq_fat_Cstyle(&GRID,path_coeff); // Make sure result doesn't change w.r.t. a trusted lattice - NerscIO::readConfiguration(Umu, header, "nersc.l8t4b3360.35link.control"); + NerscIO::readConfiguration(Umu, header, "nersc.l8t4b3360.357link.control"); LatticeGaugeField diff(&GRID); diff = Umu-U_smr; auto absDiff = norm2(diff)/norm2(Umu); From 36600899e21094e130ec21b7aad74deab6550a62 Mon Sep 17 00:00:00 2001 From: david clarke Date: Thu, 12 Oct 2023 11:11:39 -0600 Subject: [PATCH 18/40] working 7-link; Grid_log; generalShift --- Grid/log/Log.h | 21 ++++++++ Grid/qcd/smearing/HISQSmearing.h | 81 ++++++++---------------------- Grid/stencil/GeneralLocalStencil.h | 44 ++++++++++++++++ tests/smearing/Test_fatLinks.cc | 14 ------ 4 files changed, 85 insertions(+), 75 deletions(-) diff --git a/Grid/log/Log.h b/Grid/log/Log.h index 2d663a3c..b88bf61f 100644 --- a/Grid/log/Log.h +++ b/Grid/log/Log.h @@ -191,6 +191,27 @@ extern Colours GridLogColours; std::string demangle(const char* name) ; +template +inline std::string sjoin(Args&&... args) noexcept { + std::ostringstream msg; + (msg << ... << args); + return msg.str(); +} + +/*! @brief make log messages work like python print */ +template +inline void Grid_log(Args&&... args) { + std::string msg = sjoin(std::forward(args)...); + std::cout << GridLogMessage << msg << std::endl; +} + +/*! @brief make warning messages work like python print */ +template +inline void Grid_warn(Args&&... args) { + std::string msg = sjoin(std::forward(args)...); + std::cout << GridLogWarning << msg << std::endl; +} + #define _NBACKTRACE (256) extern void * Grid_backtrace_buffer[_NBACKTRACE]; diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 8c60d874..1ea1b7b9 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -32,17 +32,25 @@ directory #pragma once - #include #include #include -#define BACKWARD_CONST 16 -#define NO_SHIFT -1 NAMESPACE_BEGIN(Grid); +/*! @brief append arbitrary shift path to shifts */ +template +void appendShift(std::vector& 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); +} + + // This is to optimize the SIMD (will also need to be in the class, at least for now) template void gpermute(vobj & inout,int perm) { vobj tmp=inout; @@ -53,14 +61,6 @@ template void gpermute(vobj & inout,int perm) { } -/*! @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. Should work as long as BACKWARD_CONST > Nd. Trick inspired by SIMULATeQCD. - return dir + BACKWARD_CONST; -} - - /*! @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 @@ -69,46 +69,6 @@ inline int stencilIndex(int mu, int nu) { } -/*! @brief shift one unit in direction dir */ -template -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; - } -} - -// Move into general stencil header, beneath definition of general stencil -/*! @brief follow a path of directions, shifting one unit in each direction */ -template -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...); -} - - -/*! @brief append arbitrary shift path to shifts */ -template -void appendShift(std::vector& 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 structure holding the link treatment */ struct SmearingParameters{ SmearingParameters(){} @@ -189,17 +149,16 @@ public: // This is where contributions from the smearing get added together Ughost_fat=Zero(); - // Create the accessors - autoView(U_v , Ughost , CpuRead); - autoView(U_fat_v , Ughost_fat , CpuWrite); - autoView(U_3link_v , Ughost_3link , CpuWrite); - autoView(U_5linkA_v, Ughost_5linkA, CpuWrite); - autoView(U_5linkB_v, Ughost_5linkB, CpuWrite); - for(int mu=0;mu 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 +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 +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); diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index 200e2af6..f5c7b5ca 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -38,20 +38,6 @@ directory using namespace Grid; -/*! @brief make the logger work like python print */ -template -inline std::string sjoin(Args&&... args) noexcept { - std::ostringstream msg; - (msg << ... << args); - return msg.str(); -} -template -inline void Grid_log(Args&&... args) { - std::string msg = sjoin(std::forward(args)...); - std::cout << GridLogMessage << msg << std::endl; -} - - /*! @brief parameter file to easily adjust Nloop */ struct ConfParameters: Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS( From bf4369f72de9bd241a0bde9e5382c3dc57e38f22 Mon Sep 17 00:00:00 2001 From: david clarke Date: Thu, 12 Oct 2023 12:41:06 -0600 Subject: [PATCH 19/40] clean up HISQSmear with decltypes --- Grid/qcd/smearing/HISQSmearing.h | 101 ++++++++++++++----------------- 1 file changed, 47 insertions(+), 54 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 1ea1b7b9..432184e0 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -40,6 +40,8 @@ directory 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 void appendShift(std::vector& shifts, int dir, Args... args) { @@ -51,16 +53,6 @@ void appendShift(std::vector& shifts, int dir, Args... args) { } -// This is to optimize the SIMD (will also need to be in the class, at least for now) -template void gpermute(vobj & inout,int perm) { - vobj tmp=inout; - if (perm & 0x1) {permute(inout,tmp,0); tmp=inout;} - if (perm & 0x2) {permute(inout,tmp,1); tmp=inout;} - if (perm & 0x4) {permute(inout,tmp,2); tmp=inout;} - if (perm & 0x8) {permute(inout,tmp,3); tmp=inout;} -} - - /*! @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 @@ -163,6 +155,12 @@ public: Ughost_5linkA=Zero(); Ughost_5linkB=Zero(); + // 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; + U3matrix U0, U1, U2, U3, U4, U5, W; + // 3-link for(int site=0;site_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; + SE0 = gStencil.GetEntry(s+0,site); int x_p_mu = SE0->_offset; + SE1 = gStencil.GetEntry(s+1,site); int x_p_nu = SE1->_offset; + 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; // 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. - 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); + 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" - auto W = U2*U1*adj(U0) + adj(U5)*U4*U3; + // "left" "right" + W = U2*U1*adj(U0) + adj(U5)*U4*U3; U_3link_v[site](nu) = W; @@ -197,7 +195,6 @@ public: } } - // 5-link for(int site=0;site_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; + SE0 = gStencil.GetEntry(s+0,site); int x_p_mu = SE0->_offset; + SE1 = gStencil.GetEntry(s+1,site); int x_p_nu = SE1->_offset; + 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; // gpermutes will be replaced with single line of code, combines load and permute // into one step. still in pull request stage - 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); + 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); - auto W = U2*U1*adj(U0) + adj(U5)*U4*U3; + W = U2*U1*adj(U0) + adj(U5)*U4*U3; if(sigmaIndex<3) { U_5linkA_v[site](rho) = W; @@ -246,33 +243,29 @@ public: for(int rho=0;rho_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; + SE0 = gStencil.GetEntry(s+0,site); int x_p_mu = SE0->_offset; + SE1 = gStencil.GetEntry(s+1,site); int x_p_nu = SE1->_offset; + 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; - auto U0 = U_v[x_p_mu ](nu) ; gpermute(U0,SE0->_permute); - // decltype, or auto U1 = { ? ... } - auto U1 = U0; + U0 = coalescedReadGeneralPermute(U_v[x_p_mu](nu),SE0->_permute,Nd); if(sigmaIndex<3) { - U1 = U_5linkB_v[x_p_nu](rho); gpermute(U1,SE1->_permute); + U1 = coalescedReadGeneralPermute(U_5linkB_v[x_p_nu](rho),SE1->_permute,Nd); } else { - U1 = U_5linkA_v[x_p_nu](rho); gpermute(U1,SE1->_permute); + U1 = coalescedReadGeneralPermute(U_5linkA_v[x_p_nu](rho),SE1->_permute,Nd); } - 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 = U0; + 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 = U_5linkB_v[x_m_nu](rho); gpermute(U4,SE4->_permute); + U4 = coalescedReadGeneralPermute(U_5linkB_v[x_m_nu](rho),SE4->_permute,Nd); } else { - U4 = U_5linkA_v[x_m_nu](rho); gpermute(U4,SE4->_permute); + U4 = coalescedReadGeneralPermute(U_5linkA_v[x_m_nu](rho),SE4->_permute,Nd); } - auto U5 = U_v[x_m_nu ](nu) ; gpermute(U5,SE4->_permute); + U5 = coalescedReadGeneralPermute(U_v[x_m_nu](nu),SE4->_permute,Nd); - auto W = U2*U1*adj(U0) + adj(U5)*U4*U3; + W = U2*U1*adj(U0) + adj(U5)*U4*U3; - // std::vector(3) ? U_fat_v[site](mu) = U_fat_v[site](mu) + lt.c_7*W; sigmaIndex++; From 391fd9cc6a4c4d826eb71c2d8d761afaea3b7a69 Mon Sep 17 00:00:00 2001 From: david clarke Date: Tue, 17 Oct 2023 14:57:15 -0600 Subject: [PATCH 20/40] try lepage term --- Grid/qcd/smearing/HISQSmearing.h | 52 +++++++++++++++++++++++++------- examples/Example_plaquette.cc | 14 +++++---- tests/smearing/Test_fatLinks.cc | 4 +-- 3 files changed, 51 insertions(+), 19 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 432184e0..44e14e85 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -81,8 +81,8 @@ struct SmearingParameters{ /*! @brief create fat links from link variables */ -template -class Smear_HISQ_fat { +template +class Smear_HISQ_fat : public Gimpl { private: GridCartesian* const _grid; @@ -90,6 +90,8 @@ private: public: + INHERIT_GIMPL_TYPES(Gimpl); + // Don't allow default values here. Smear_HISQ_fat(GridCartesian* grid, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) : _grid(grid), @@ -141,6 +143,7 @@ public: // 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_offset; SE4 = gStencil.GetEntry(s+4,site); int x_m_nu = SE4->_offset; - // gpermutes will be replaced with single line of code, combines load and permute - // into one step. still in pull request stage 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); @@ -234,8 +233,7 @@ public: } } - // 7-link - for(int site=0;site U(Nd, u_thin.Grid()); + std::vector V(Nd, u_smr.Grid()); + for (int mu = 0; mu < Nd; mu++) { + U[mu] = PeekIndex(u_thin, mu); + V[mu] = PeekIndex(u_smr, mu); + } + + // Compute LePage term from U_thin: + for(int mu=0;mu(u_smr, V[mu], mu); + } + }; // void derivative(const GaugeField& Gauge) const { // }; diff --git a/examples/Example_plaquette.cc b/examples/Example_plaquette.cc index 17de4762..faf17d82 100644 --- a/examples/Example_plaquette.cc +++ b/examples/Example_plaquette.cc @@ -29,7 +29,6 @@ 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. - // WARNING: The first time you include this or take it out, the compile time will increase a lot. INHERIT_GIMPL_TYPES(Gimpl); // Some example Gimpls can be found in GaugeImplementations.h, at the bottom. These are in turn built @@ -53,12 +52,15 @@ public: // U_mu_nu(x) static void dirPlaquette(GaugeMat &plaq, const std::vector &U, const int mu, const int nu) { - // These CovShift calls seem to carry out the multiplication already. A positive shift moves the lattice - // site x_mu = 1 in the RHS to x_mu = 0 in the result. + // 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, - Gimpl::CovShiftForward(U[nu],nu, - Gimpl::CovShiftBackward(U[mu],mu, - Gimpl::CovShiftIdentityBackward(U[nu], nu)))); + // 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) diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index f5c7b5ca..1b24c7ca 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -90,14 +90,14 @@ int main (int argc, char** argv) { NerscIO::readConfiguration(Umu, header, conf_in); // Smear Umu and store result in U_smr - Smear_HISQ_fat hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,0.); + Smear_HISQ_fat hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,-1/8.); hisq_fat.smear(U_smr,Umu); NerscIO::writeConfiguration(U_smr,conf_out,"HISQ"); // Test a C-style instantiation double path_coeff[6] = {1, 2, 3, 4, 5, 6}; - Smear_HISQ_fat hisq_fat_Cstyle(&GRID,path_coeff); + Smear_HISQ_fat 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"); From 7bb8ab7000962dd9da9e56254c4550be4cbe0c44 Mon Sep 17 00:00:00 2001 From: david clarke Date: Fri, 20 Oct 2023 08:41:02 -0600 Subject: [PATCH 21/40] improve smearing templating --- Grid/qcd/smearing/HISQSmearing.h | 4 +++- tests/smearing/Test_fatLinks.cc | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 44e14e85..82eedc1d 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -81,7 +81,8 @@ struct SmearingParameters{ /*! @brief create fat links from link variables */ -template +//template +template class Smear_HISQ_fat : public Gimpl { private: @@ -91,6 +92,7 @@ private: public: INHERIT_GIMPL_TYPES(Gimpl); + typedef typename Gimpl::GaugeField LGF; // Don't allow default values here. Smear_HISQ_fat(GridCartesian* grid, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index 1b24c7ca..2ff3c116 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -90,14 +90,14 @@ int main (int argc, char** argv) { NerscIO::readConfiguration(Umu, header, conf_in); // Smear Umu and store result in U_smr - Smear_HISQ_fat hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,-1/8.); + Smear_HISQ_fat hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,-1/8.); hisq_fat.smear(U_smr,Umu); NerscIO::writeConfiguration(U_smr,conf_out,"HISQ"); // Test a C-style instantiation double path_coeff[6] = {1, 2, 3, 4, 5, 6}; - Smear_HISQ_fat hisq_fat_Cstyle(&GRID,path_coeff); + Smear_HISQ_fat 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"); From 21ed6ac0f4c442564da142098952bf6b226c7d52 Mon Sep 17 00:00:00 2001 From: david clarke Date: Fri, 20 Oct 2023 13:54:26 -0600 Subject: [PATCH 22/40] added floating-point support --- Grid/qcd/smearing/HISQSmearing.h | 25 +++++++++++++------------ tests/smearing/Test_fatLinks.cc | 10 ++++++---- 2 files changed, 19 insertions(+), 16 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 82eedc1d..9ac6eddd 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -81,7 +81,7 @@ struct SmearingParameters{ /*! @brief create fat links from link variables */ -//template +//template template class Smear_HISQ_fat : public Gimpl { @@ -92,7 +92,8 @@ private: public: INHERIT_GIMPL_TYPES(Gimpl); - typedef typename Gimpl::GaugeField LGF; + typedef typename Gimpl::GaugeField GF; + typedef typename Gimpl::GaugeLinkField LF; // Don't allow default values here. Smear_HISQ_fat(GridCartesian* grid, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) @@ -112,20 +113,20 @@ public: ~Smear_HISQ_fat() {} - void smear(LGF& u_smr, LGF& u_thin) const { + void smear(GF& u_smr, GF& u_thin) const { SmearingParameters lt = this->_linkTreatment; // 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); + GF Ughost = Ghost.Exchange(u_thin); // This is where auxiliary N-link fields and the final smear will be stored. - LGF Ughost_fat(Ughost.Grid()); - LGF Ughost_3link(Ughost.Grid()); - LGF Ughost_5linkA(Ughost.Grid()); - LGF Ughost_5linkB(Ughost.Grid()); + GF Ughost_fat(Ughost.Grid()); + GF Ughost_3link(Ughost.Grid()); + GF Ughost_5linkA(Ughost.Grid()); + GF Ughost_5linkB(Ughost.Grid()); // Create 3-link stencil. We allow mu==nu just to make the indexing easier. // Shifts with mu==nu will not be used. @@ -279,8 +280,8 @@ public: 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 U(Nd, u_thin.Grid()); - std::vector V(Nd, u_smr.Grid()); + std::vector U(Nd, u_thin.Grid()); + std::vector V(Nd, u_smr.Grid()); for (int mu = 0; mu < Nd; mu++) { U[mu] = PeekIndex(u_thin, mu); V[mu] = PeekIndex(u_smr, mu); @@ -317,7 +318,7 @@ public: /*! @brief create long links from link variables. */ -template +template class Smear_HISQ_Naik { private: @@ -333,7 +334,7 @@ public: ~Smear_HISQ_Naik() {} -// void smear(LGF& u_smr, const LGF& U) const { +// void smear(GF& u_smr, const GF& U) const { // }; // void derivative(const GaugeField& Gauge) const { diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index 2ff3c116..acfb626c 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -67,6 +67,8 @@ int main (int argc, char** argv) { std::string conf_out = "nersc.l8t4b3360.357link"; int threads = GridThread::GetThreads(); + typedef LatticeGaugeFieldD LGF; + // Initialize the Grid Grid_init(&argc,&argv); Coordinate simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd()); @@ -81,9 +83,9 @@ int main (int argc, char** argv) { 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); + // Instantiate the LGF objects holding thin (Umu) and fat (U_smr) links + LGF Umu(&GRID); + LGF U_smr(&GRID); // Read the configuration into Umu FieldMetaData header; @@ -101,7 +103,7 @@ int main (int argc, char** argv) { // Make sure result doesn't change w.r.t. a trusted lattice NerscIO::readConfiguration(Umu, header, "nersc.l8t4b3360.357link.control"); - LatticeGaugeField diff(&GRID); + LGF diff(&GRID); diff = Umu-U_smr; auto absDiff = norm2(diff)/norm2(Umu); Grid_log(" |Umu-U|/|Umu| = ",absDiff); From 3d3376d1a321b007daf52dcbfc8746da59611b96 Mon Sep 17 00:00:00 2001 From: david clarke Date: Fri, 27 Oct 2023 16:26:31 -0600 Subject: [PATCH 23/40] LePage works, trying Naik --- Grid/qcd/smearing/HISQSmearing.h | 46 +++++++++++--------------------- tests/smearing/Test_fatLinks.cc | 6 ++--- 2 files changed, 19 insertions(+), 33 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 9ac6eddd..4ec61141 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -56,7 +56,7 @@ void appendShift(std::vector& 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. - // Shifts with mu==nu will not be used. + // mu-nu plane stencil. We allow mu==nu to make indexing the stencil easier, + // but these entries will not be used. std::vector shifts; for(int mu=0;mu_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_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; - 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 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); \ No newline at end of file diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index acfb626c..61668b66 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/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 hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,-1/8.); + Smear_HISQ_fat 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 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); From df9b958c40e920ff8e3ae5c565be9be57b377735 Mon Sep 17 00:00:00 2001 From: david clarke Date: Mon, 30 Oct 2023 17:40:53 -0600 Subject: [PATCH 24/40] naik now returns separately --- Grid/log/Log.h | 22 +++++++++++++--- Grid/qcd/smearing/HISQSmearing.h | 27 ++++++++++++-------- tests/smearing/Test_fatLinks.cc | 44 ++++++++++++++++++++------------ 3 files changed, 63 insertions(+), 30 deletions(-) diff --git a/Grid/log/Log.h b/Grid/log/Log.h index b88bf61f..370b0428 100644 --- a/Grid/log/Log.h +++ b/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; @@ -209,7 +209,21 @@ inline void Grid_log(Args&&... args) { template inline void Grid_warn(Args&&... args) { std::string msg = sjoin(std::forward(args)...); - std::cout << GridLogWarning << msg << std::endl; + std::cout << "\033[33m" << GridLogWarning << msg << "\033[0m" << std::endl; +} + +/*! @brief make error messages work like python print */ +template +inline void Grid_error(Args&&... args) { + std::string msg = sjoin(std::forward(args)...); + std::cout << "\033[31m" << GridLogError << msg << "\033[0m" << std::endl; +} + +/*! @brief make pass messages work like python print */ +template +inline void Grid_pass(Args&&... args) { + std::string msg = sjoin(std::forward(args)...); + std::cout << "\033[32m" << GridLogMessage << msg << "\033[0m" << std::endl; } #define _NBACKTRACE (256) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 4ec61141..a0b60dcd 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -81,9 +81,9 @@ struct SmearingParameters{ /*! @brief create fat links from link variables */ -//template template -class Smear_HISQ_fat : public Gimpl { +class Smear_HISQ : public Gimpl { +// TODO: this needs to be renamed, becaues the Naik guy is not part of the fat smear private: GridCartesian* const _grid; @@ -96,7 +96,7 @@ public: typedef typename Gimpl::GaugeLinkField LF; // Don't allow default values here. - Smear_HISQ_fat(GridCartesian* grid, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) + 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"); @@ -104,16 +104,18 @@ public: } // Allow to pass a pointer to a C-style, double array for MILC convenience - Smear_HISQ_fat(GridCartesian* grid, double* coeff) + 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_fat() {} + ~Smear_HISQ() {} - void smear(GF& u_smr, GF& u_thin) const { + // 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; @@ -127,6 +129,7 @@ public: GF Ughost_3link(Ughost.Grid()); GF Ughost_5linkA(Ughost.Grid()); GF Ughost_5linkB(Ughost.Grid()); + GF Ughost_naik(Ughost.Grid()); // mu-nu plane stencil. We allow mu==nu to make indexing the stencil easier, // but these entries will not be used. @@ -146,6 +149,7 @@ public: // This is where contributions from the smearing get added together Ughost_fat=Zero(); + Ughost_naik=Zero(); // This loop handles 3-, 5-, and 7-link constructs, minus Lepage and Naik. for(int mu=0;mu_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; + U_naik_v[site](mu) = U_fat_v[site](mu) + lt.c_naik*W; } } @@ -290,6 +294,9 @@ public: // c1, c3, c5, c7 construct contributions u_smr = Ghost.Extract(Ughost_fat) + lt.c_1*u_thin; + // Naik contribution + u_naik = Ghost.Extract(Ughost_naik); + // Load up U and V std::vectors to access thin and smeared links. std::vector U(Nd, u_thin.Grid()); std::vector V(Nd, u_smr.Grid()); diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index 61668b66..05e80929 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -57,6 +57,22 @@ struct ConfParameters: Serializable { // another : input --> unitarize // + +void 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 hisq_fat(&GRID,c1,cnaik,c3,c5,c7,clp); + hisq_fat.smear(Usmr, Unaik, Umu); + LatticeGaugeFieldD diff(&GRID); + diff = Ucontrol-Usmr; + auto absDiff = norm2(diff)/norm2(Ucontrol); + if (absDiff < 1e-30) { + Grid_pass(" |Umu-U|/|Umu| = ",absDiff); + } else { + Grid_error(" |Umu-U|/|Umu| = ",absDiff); + } +} + + int main (int argc, char** argv) { // Params for the test. @@ -83,30 +99,26 @@ int main (int argc, char** argv) { ConfParameters param(Reader); if(param.benchmark) Grid_log(" Nloop = ",param.Nloop); - // Instantiate the LGF objects holding thin (Umu) and fat (U_smr) links - LGF Umu(&GRID); - LGF U_smr(&GRID); + LGF Umu(&GRID), Usmr(&GRID), Unaik(&GRID), Ucontrol(&GRID); // Read the configuration into Umu FieldMetaData header; NerscIO::readConfiguration(Umu, header, conf_in); - // Smear Umu and store result in U_smr - Smear_HISQ_fat 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"); + // Carry out various tests + NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.357lplink.control"); + testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,1/384.,-1/8.); + NerscIO::writeConfiguration(Usmr,conf_out,"HISQ"); + NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.357link.control"); + testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,1/384.,0.); + NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.35link.control"); + testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,0.,0.); + NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.3link.control"); + testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,0.,0.,0.); // Test a C-style instantiation double path_coeff[6] = {1, 2, 3, 4, 5, 6}; - Smear_HISQ_fat hisq_fat_Cstyle(&GRID,path_coeff); - - // Make sure result doesn't change w.r.t. a trusted lattice - NerscIO::readConfiguration(Umu, header, "nersc.l8t4b3360.357lplink.control"); - LGF diff(&GRID); - diff = Umu-U_smr; - auto absDiff = norm2(diff)/norm2(Umu); - Grid_log(" |Umu-U|/|Umu| = ",absDiff); + Smear_HISQ hisq_fat_Cstyle(&GRID,path_coeff); if (param.benchmark) { From 69c869d345dc0b842138cab082cf01d84b5e8dac Mon Sep 17 00:00:00 2001 From: david clarke Date: Mon, 30 Oct 2023 17:41:52 -0600 Subject: [PATCH 25/40] fixed stupid typo --- Grid/qcd/smearing/HISQSmearing.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index a0b60dcd..46cc527d 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -211,7 +211,7 @@ public: U1 = coalescedReadGeneralPermute(U_v[x ](mu),SE2->_permute,Nd); U2 = coalescedReadGeneralPermute(U_v[x_p_mu](mu),SE0->_permute,Nd); W = U0*U1*U2; - U_naik_v[site](mu) = U_fat_v[site](mu) + lt.c_naik*W; + U_naik_v[site](mu) = lt.c_naik*W; } } From 2ae2a81e859976996894d289f198b53dfe19c637 Mon Sep 17 00:00:00 2001 From: david clarke Date: Tue, 31 Oct 2023 13:54:55 -0600 Subject: [PATCH 26/40] attempt to fix Naik --- Grid/qcd/smearing/HISQSmearing.h | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 46cc527d..43fe0f77 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -203,15 +203,16 @@ public: 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; + U_3link_v[x](nu) = W; + U_fat_v[x](mu) = U_fat_v[x](mu) + lt.c_3*W; - // TODO: May need to be shifted by 1? + // Naik term starts at x-mu, save at x-mu. The idea will be to keep track + // of this shift, and then take into account when we use Naik later. 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; - U_naik_v[site](mu) = lt.c_naik*W; + U_naik_v[x_m_mu](mu) = lt.c_naik*W; } } @@ -239,12 +240,12 @@ public: W = U2*U1*adj(U0) + adj(U5)*U4*U3; if(sigmaIndex<3) { - U_5linkA_v[site](rho) = W; + U_5linkA_v[x](rho) = W; } else { - U_5linkB_v[site](rho) = W; + U_5linkB_v[x](rho) = W; } - U_fat_v[site](mu) = U_fat_v[site](mu) + lt.c_5*W; + U_fat_v[x](mu) = U_fat_v[x](mu) + lt.c_5*W; sigmaIndex++; } @@ -282,7 +283,7 @@ public: W = U2*U1*adj(U0) + adj(U5)*U4*U3; - U_fat_v[site](mu) = U_fat_v[site](mu) + lt.c_7*W; + U_fat_v[x](mu) = U_fat_v[x](mu) + lt.c_7*W; sigmaIndex++; } From c8b17c9526679541443141559604636fe0657e66 Mon Sep 17 00:00:00 2001 From: david clarke Date: Thu, 2 Nov 2023 12:43:22 -0600 Subject: [PATCH 27/40] Naik to CShift --- Grid/qcd/smearing/HISQSmearing.h | 28 +++++++++++----------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 43fe0f77..36f66ee5 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -129,7 +129,6 @@ public: GF Ughost_3link(Ughost.Grid()); GF Ughost_5linkA(Ughost.Grid()); GF Ughost_5linkB(Ughost.Grid()); - GF Ughost_naik(Ughost.Grid()); // mu-nu plane stencil. We allow mu==nu to make indexing the stencil easier, // but these entries will not be used. @@ -149,7 +148,6 @@ public: // This is where contributions from the smearing get added together Ughost_fat=Zero(); - Ughost_naik=Zero(); // This loop handles 3-, 5-, and 7-link constructs, minus Lepage and Naik. for(int mu=0;mu_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; - U_naik_v[x_m_mu](mu) = lt.c_naik*W; } } @@ -295,19 +284,24 @@ public: // c1, c3, c5, c7 construct contributions u_smr = Ghost.Extract(Ughost_fat) + lt.c_1*u_thin; - // Naik contribution - u_naik = Ghost.Extract(Ughost_naik); - // Load up U and V std::vectors to access thin and smeared links. std::vector U(Nd, u_thin.Grid()); std::vector V(Nd, u_smr.Grid()); + std::vector Vnaik(Nd, u_naik.Grid()); for (int mu = 0; mu < Nd; mu++) { - U[mu] = PeekIndex(u_thin, mu); - V[mu] = PeekIndex(u_smr, mu); + U[mu] = PeekIndex(u_thin, mu); + V[mu] = PeekIndex(u_smr, mu); + Vnaik[mu] = PeekIndex(u_naik, mu); } - // Compute LePage term from U_thin: for(int mu=0;mu Date: Fri, 3 Nov 2023 14:11:38 -0600 Subject: [PATCH 28/40] fix naik bug --- Grid/qcd/smearing/HISQSmearing.h | 10 +++++----- tests/smearing/Test_fatLinks.cc | 23 ++++++++++++++++------- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 36f66ee5..959a6cf0 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -291,15 +291,14 @@ public: for (int mu = 0; mu < Nd; mu++) { U[mu] = PeekIndex(u_thin, mu); V[mu] = PeekIndex(u_smr, mu); - Vnaik[mu] = PeekIndex(u_naik, mu); } for(int mu=0;mu(u_smr, V[mu], mu); + PokeIndex(u_smr , V[mu] , mu); + PokeIndex(u_naik, Vnaik[mu], mu); } }; diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index 05e80929..e3aa0e33 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -62,13 +62,20 @@ void testSmear(GridCartesian& GRID, LatticeGaugeFieldD Umu, LatticeGaugeFieldD U LatticeGaugeFieldD Ucontrol, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) { Smear_HISQ hisq_fat(&GRID,c1,cnaik,c3,c5,c7,clp); hisq_fat.smear(Usmr, Unaik, Umu); - LatticeGaugeFieldD diff(&GRID); - diff = Ucontrol-Usmr; - auto absDiff = norm2(diff)/norm2(Ucontrol); - if (absDiff < 1e-30) { - Grid_pass(" |Umu-U|/|Umu| = ",absDiff); + LatticeGaugeFieldD diff1(&GRID), diff2(&GRID); + diff1 = Ucontrol-Usmr; + diff2 = Ucontrol-Unaik; + auto absDiff1 = norm2(diff1)/norm2(Ucontrol); + auto absDiff2 = norm2(diff2)/norm2(Ucontrol); + if (absDiff1 < 1e-30) { + Grid_pass(" |Umu-Usmr|/|Umu| = ",absDiff1); } else { - Grid_error(" |Umu-U|/|Umu| = ",absDiff); + Grid_error(" |Umu-Usmr|/|Umu| = ",absDiff1); + } + if (absDiff2 < 1e-30) { + Grid_pass(" |Umu-Unaik|/|Umu| = ",absDiff2); + } else { + Grid_error(" |Umu-Unaik|/|Umu| = ",absDiff2); } } @@ -114,7 +121,9 @@ int main (int argc, char** argv) { NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.35link.control"); testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,0.,0.); NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.3link.control"); - testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,0.,0.,0.); + testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,1.,1/16.,0.,0.,0.); + NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.3link.control"); + testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,2.,1/16.,0.,0.,0.); // Test a C-style instantiation double path_coeff[6] = {1, 2, 3, 4, 5, 6}; From 981c93d67ac2ca2fbe66bc238c4c538492db8ad6 Mon Sep 17 00:00:00 2001 From: david clarke Date: Sun, 21 Jan 2024 21:09:19 -0700 Subject: [PATCH 29/40] update Test_fatLinks to accept Naik --- tests/smearing/Test_fatLinks.cc | 42 +++++++++++++++++---------------- 1 file changed, 22 insertions(+), 20 deletions(-) diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index e3aa0e33..d9d35e69 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -58,24 +58,28 @@ struct ConfParameters: Serializable { // -void testSmear(GridCartesian& GRID, LatticeGaugeFieldD Umu, LatticeGaugeFieldD Usmr, LatticeGaugeFieldD Unaik, +void 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 hisq_fat(&GRID,c1,cnaik,c3,c5,c7,clp); + LatticeGaugeFieldD diff(&GRID); hisq_fat.smear(Usmr, Unaik, Umu); - LatticeGaugeFieldD diff1(&GRID), diff2(&GRID); - diff1 = Ucontrol-Usmr; - diff2 = Ucontrol-Unaik; - auto absDiff1 = norm2(diff1)/norm2(Ucontrol); - auto absDiff2 = norm2(diff2)/norm2(Ucontrol); - if (absDiff1 < 1e-30) { - Grid_pass(" |Umu-Usmr|/|Umu| = ",absDiff1); - } else { - Grid_error(" |Umu-Usmr|/|Umu| = ",absDiff1); - } - if (absDiff2 < 1e-30) { - Grid_pass(" |Umu-Unaik|/|Umu| = ",absDiff2); - } else { - Grid_error(" |Umu-Unaik|/|Umu| = ",absDiff2); + 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); + } else { + Grid_error(" |Umu-Usmr|/|Umu| = ",absDiff); + } + } else { // Testing Naik specifically + diff = Ucontrol-Unaik; + auto absDiff = norm2(diff)/norm2(Ucontrol); + if (absDiff < 1e-30) { + Grid_pass(" |Umu-Unaik|/|Umu| = ",absDiff); + } else { + Grid_error(" |Umu-Unaik|/|Umu| = ",absDiff); + } +// NerscIO::writeConfiguration(Unaik,"nersc.l8t4b3360.naik"); } } @@ -87,7 +91,6 @@ 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.357lplink"; int threads = GridThread::GetThreads(); typedef LatticeGaugeFieldD LGF; @@ -115,15 +118,14 @@ int main (int argc, char** argv) { // Carry out various tests NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.357lplink.control"); testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,1/384.,-1/8.); - NerscIO::writeConfiguration(Usmr,conf_out,"HISQ"); NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.357link.control"); testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,1/384.,0.); NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.35link.control"); testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,0.,0.); NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.3link.control"); - testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,1.,1/16.,0.,0.,0.); - NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.3link.control"); - testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,2.,1/16.,0.,0.,0.); + testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,0.,0.,0.); + NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.naik.control"); + testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,0.,0.8675309,0.,0.,0.,0.); // Test a C-style instantiation double path_coeff[6] = {1, 2, 3, 4, 5, 6}; From f5b3d582b044ea420256f2470b4efe934c7012bb Mon Sep 17 00:00:00 2001 From: david clarke Date: Mon, 22 Jan 2024 02:49:40 -0700 Subject: [PATCH 30/40] first attempt at U3 projection --- Grid/qcd/smearing/HISQSmearing.h | 46 +++++++++++++++++++++++++++----- tests/smearing/Test_fatLinks.cc | 3 ++- 2 files changed, 42 insertions(+), 7 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 959a6cf0..43b06534 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -83,7 +83,6 @@ struct SmearingParameters{ /*! @brief create fat links from link variables */ template class Smear_HISQ : public Gimpl { -// TODO: this needs to be renamed, becaues the Naik guy is not part of the fat smear private: GridCartesian* const _grid; @@ -289,16 +288,16 @@ public: std::vector V(Nd, u_smr.Grid()); std::vector Vnaik(Nd, u_naik.Grid()); for (int mu = 0; mu < Nd; mu++) { - U[mu] = PeekIndex(u_thin, mu); - V[mu] = PeekIndex(u_smr, mu); + U[mu] = PeekIndex(u_thin, mu); + V[mu] = PeekIndex(u_smr, mu); } for(int mu=0;mu(u_mu, mu); + Q = adj(V[mu])*V[mu]; + c1 = trace(Q*Q)/2.; // SU(N) matrices are traceless, so c0=0. + c2 = trace(Q*Q*Q)/3.; + S = c1/3.; + R = c2/2.; + theta = std::acos(R/std::pow(S,1.5)); + g0 = 2.*std::sqrt(S)*std::cos(theta/3.-2*M_PI/3.); + g1 = 2.*std::sqrt(S)*std::cos(theta/3. ); + g2 = 2.*std::sqrt(S)*std::cos(theta/3.+2*M_PI/3.); +// if (fabs(Q.determinant()/(g0*g1*g2)-1.0) > 1e-5) {} + u = std::sqrt(g0) + std::sqrt(g1) + std::sqrt(g2); + v = std::sqrt(g0*g1) + std::sqrt(g0*g2) + std::sqrt(g1*g2); + w = std::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; + + sqrtQinv = f0 + f1*Q + f2*Q*Q; + PokeIndex(u_proj, V*sqrtQinv, mu); + } + }; + + // void derivative(const GaugeField& Gauge) const { // }; }; diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index d9d35e69..e2dc5d6d 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -61,7 +61,7 @@ struct ConfParameters: Serializable { void 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 hisq_fat(&GRID,c1,cnaik,c3,c5,c7,clp); - LatticeGaugeFieldD diff(&GRID); + LatticeGaugeFieldD diff(&GRID), Uproj(&GRID); hisq_fat.smear(Usmr, Unaik, Umu); if (cnaik < 1e-30) { // Testing anything but Naik term diff = Ucontrol-Usmr; @@ -79,6 +79,7 @@ void testSmear(GridCartesian& GRID, LatticeGaugeFieldD Umu, LatticeGaugeFieldD U } else { Grid_error(" |Umu-Unaik|/|Umu| = ",absDiff); } + hisq_fat.projectU3(Uproj,Usmr); // NerscIO::writeConfiguration(Unaik,"nersc.l8t4b3360.naik"); } } From 00f24f87653d92276f9ffeeaf403f78b5adef4cc Mon Sep 17 00:00:00 2001 From: david clarke Date: Mon, 22 Jan 2024 05:50:16 -0700 Subject: [PATCH 31/40] already found some bugs in projection, still needs testing --- Grid/qcd/smearing/HISQSmearing.h | 38 ++++++++++++++++++-------------- 1 file changed, 22 insertions(+), 16 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 43b06534..17959495 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -93,6 +93,7 @@ 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) @@ -329,31 +330,36 @@ public: // IN--u_mu void projectU3(GF& u_proj, GF& u_mu) const { - LF V, Q, sqrtQinv; - Real c1, c2, g0, g1, g2, S, R, theta, u, v, w, den, f0, f1, f2; + LF V(u_mu.Grid()), Q(u_mu.Grid()), sqrtQinv(u_mu.Grid()), id_3(u_mu.Grid()); + CF c0(u_mu.Grid()), c1(u_mu.Grid()), c2(u_mu.Grid()), g0(u_mu.Grid()), g1(u_mu.Grid()), + g2(u_mu.Grid()), S(u_mu.Grid()), R(u_mu.Grid()), theta(u_mu.Grid()), u(u_mu.Grid()), + v(u_mu.Grid()), w(u_mu.Grid()), den(u_mu.Grid()), f0(u_mu.Grid()), f1(u_mu.Grid()), + f2(u_mu.Grid()); // Follow MILC 10.1103/PhysRevD.82.074501, eqs (B2-B3) and (C1-C8) for (int mu = 0; mu < Nd; mu++) { V = PeekIndex(u_mu, mu); - Q = adj(V[mu])*V[mu]; - c1 = trace(Q*Q)/2.; // SU(N) matrices are traceless, so c0=0. - c2 = trace(Q*Q*Q)/3.; - S = c1/3.; - R = c2/2.; - theta = std::acos(R/std::pow(S,1.5)); - g0 = 2.*std::sqrt(S)*std::cos(theta/3.-2*M_PI/3.); - g1 = 2.*std::sqrt(S)*std::cos(theta/3. ); - g2 = 2.*std::sqrt(S)*std::cos(theta/3.+2*M_PI/3.); + 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; + 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) {} - u = std::sqrt(g0) + std::sqrt(g1) + std::sqrt(g2); - v = std::sqrt(g0*g1) + std::sqrt(g0*g2) + std::sqrt(g1*g2); - w = std::sqrt(g0*g1*g2); + 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; + f1 = (-w-u*u*u+2.*u*v)/den; f2 = u/den; + id_3 = 1.; - sqrtQinv = f0 + f1*Q + f2*Q*Q; + sqrtQinv = f0*id_3 + f1*Q + f2*Q*Q; PokeIndex(u_proj, V*sqrtQinv, mu); } }; From 4924b3209e9bb5c358699ebfd8c3d929eaec87c2 Mon Sep 17 00:00:00 2001 From: david clarke Date: Tue, 23 Jan 2024 14:43:58 -0700 Subject: [PATCH 32/40] projectU3 yields a unitary matrix --- .gitignore | 4 ++++ Grid/qcd/smearing/HISQSmearing.h | 31 ++++++++++++++++++------------- tests/smearing/Test_fatLinks.cc | 8 +------- 3 files changed, 23 insertions(+), 20 deletions(-) diff --git a/.gitignore b/.gitignore index 40156f9d..94e866e2 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,7 @@ +# Doxygen stuff +html/* +latex/* + # Compiled Object files # ######################### *.slo diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 17959495..c8255acc 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -330,7 +330,7 @@ public: // IN--u_mu void projectU3(GF& u_proj, GF& u_mu) const { - LF V(u_mu.Grid()), Q(u_mu.Grid()), sqrtQinv(u_mu.Grid()), id_3(u_mu.Grid()); + LF V(u_mu.Grid()), Q(u_mu.Grid()), sqrtQinv(u_mu.Grid()), id_3(u_mu.Grid()), diff(u_mu.Grid()); CF c0(u_mu.Grid()), c1(u_mu.Grid()), c2(u_mu.Grid()), g0(u_mu.Grid()), g1(u_mu.Grid()), g2(u_mu.Grid()), S(u_mu.Grid()), R(u_mu.Grid()), theta(u_mu.Grid()), u(u_mu.Grid()), v(u_mu.Grid()), w(u_mu.Grid()), den(u_mu.Grid()), f0(u_mu.Grid()), f1(u_mu.Grid()), @@ -338,18 +338,22 @@ public: // Follow MILC 10.1103/PhysRevD.82.074501, eqs (B2-B3) and (C1-C8) for (int mu = 0; mu < Nd; mu++) { - V = PeekIndex(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; - 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) {} + V = PeekIndex(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); @@ -360,6 +364,7 @@ public: id_3 = 1.; sqrtQinv = f0*id_3 + f1*Q + f2*Q*Q; + PokeIndex(u_proj, V*sqrtQinv, mu); } }; diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc index e2dc5d6d..742cb205 100644 --- a/tests/smearing/Test_fatLinks.cc +++ b/tests/smearing/Test_fatLinks.cc @@ -51,12 +51,6 @@ struct ConfParameters: Serializable { } }; -// -// one method: input --> fat -// another : input --> long (naik) -// another : input --> unitarize -// - void testSmear(GridCartesian& GRID, LatticeGaugeFieldD Umu, LatticeGaugeFieldD Usmr, LatticeGaugeFieldD Unaik, LatticeGaugeFieldD Ucontrol, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) { @@ -79,7 +73,7 @@ void testSmear(GridCartesian& GRID, LatticeGaugeFieldD Umu, LatticeGaugeFieldD U } else { Grid_error(" |Umu-Unaik|/|Umu| = ",absDiff); } - hisq_fat.projectU3(Uproj,Usmr); + hisq_fat.projectU3(Uproj,Ucontrol); // NerscIO::writeConfiguration(Unaik,"nersc.l8t4b3360.naik"); } } From 0a6e2f42c5b8382bf59e24a37aa1adc8da0c7577 Mon Sep 17 00:00:00 2001 From: david clarke Date: Tue, 6 Feb 2024 16:32:07 -0700 Subject: [PATCH 33/40] small amount of cleanup --- Grid/qcd/smearing/HISQSmearing.h | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index c8255acc..f053afcf 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -118,10 +118,11 @@ public: 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,this->_grid); + PaddedCell Ghost(depth,grid); GF Ughost = Ghost.Exchange(u_thin); // This is where auxiliary N-link fields and the final smear will be stored. @@ -285,9 +286,9 @@ public: 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 U(Nd, u_thin.Grid()); - std::vector V(Nd, u_smr.Grid()); - std::vector Vnaik(Nd, u_naik.Grid()); + std::vector U(Nd, grid); + std::vector V(Nd, grid); + std::vector Vnaik(Nd, grid); for (int mu = 0; mu < Nd; mu++) { U[mu] = PeekIndex(u_thin, mu); V[mu] = PeekIndex(u_smr, mu); @@ -330,11 +331,11 @@ public: // IN--u_mu void projectU3(GF& u_proj, GF& u_mu) const { - LF V(u_mu.Grid()), Q(u_mu.Grid()), sqrtQinv(u_mu.Grid()), id_3(u_mu.Grid()), diff(u_mu.Grid()); - CF c0(u_mu.Grid()), c1(u_mu.Grid()), c2(u_mu.Grid()), g0(u_mu.Grid()), g1(u_mu.Grid()), - g2(u_mu.Grid()), S(u_mu.Grid()), R(u_mu.Grid()), theta(u_mu.Grid()), u(u_mu.Grid()), - v(u_mu.Grid()), w(u_mu.Grid()), den(u_mu.Grid()), f0(u_mu.Grid()), f1(u_mu.Grid()), - f2(u_mu.Grid()); + 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++) { From a38fb0e04a01e751c87a719b79294086f8864435 Mon Sep 17 00:00:00 2001 From: david clarke Date: Tue, 6 Feb 2024 18:24:55 -0700 Subject: [PATCH 34/40] first effort toward accelerators --- Grid/qcd/smearing/HISQSmearing.h | 34 ++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index f053afcf..0deb080d 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -154,7 +154,7 @@ public: for(int mu=0;mu Date: Tue, 6 Feb 2024 18:40:13 -0700 Subject: [PATCH 35/40] acceleration compiles and doesn't break scalar mode --- Grid/qcd/smearing/HISQSmearing.h | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 0deb080d..2ae11bbb 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -153,18 +153,18 @@ public: // This loop handles 3-, 5-, and 7-link constructs, minus Lepage and Naik. for(int mu=0;mu_permute,Nd)) U3matrix; From db420525b3f872e8dab0b72ff1d5696be575590d Mon Sep 17 00:00:00 2001 From: david clarke Date: Mon, 12 Feb 2024 15:03:53 -0700 Subject: [PATCH 36/40] fix Simd::Nsimd typo --- Grid/qcd/smearing/HISQSmearing.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 2ae11bbb..7635ef06 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -173,7 +173,7 @@ public: int Nsites = U_v.size(); - accelerator_for(site,Nsites,Simd:Nsimd(),{ // ----------- 3-link constructs + accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 3-link constructs // for(int site=0;site Date: Wed, 14 Feb 2024 13:56:57 -0700 Subject: [PATCH 37/40] accelerator_inline bug --- Grid/qcd/smearing/HISQSmearing.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 7635ef06..d2091806 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -54,7 +54,7 @@ void appendShift(std::vector& shifts, int dir, Args... args) { /*! @brief figure out the stencil index from mu and nu */ -inline int stencilIndex(int mu, int 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; From 94581e3c7ad7d3b065f16959160b2e6e8669f4f0 Mon Sep 17 00:00:00 2001 From: david clarke Date: Fri, 23 Feb 2024 15:58:33 -0700 Subject: [PATCH 38/40] accelerator_for is broken --- Grid/qcd/smearing/HISQSmearing.h | 24 ++++++++++++------------ tests/smearing/Test_fatLinks.cc | 29 +++++++++++++++++++++-------- 2 files changed, 33 insertions(+), 20 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index d2091806..ac4cb8b6 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -173,8 +173,8 @@ public: int Nsites = U_v.size(); - accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 3-link constructs -// for(int site=0;site 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); + hisq_fat.projectU3(Uproj,Ucontrol); // NerscIO::writeConfiguration(Unaik,"nersc.l8t4b3360.naik"); } + return result; } @@ -110,23 +116,30 @@ int main (int argc, char** argv) { FieldMetaData header; NerscIO::readConfiguration(Umu, header, conf_in); + bool pass=true; + // Carry out various tests NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.357lplink.control"); - testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,1/384.,-1/8.); + 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"); - testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,1/384.,0.); + 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"); - testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,0.,0.); + pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,0.,0.); NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.3link.control"); - testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,0.,0.,0.); + pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,0.,0.,0.); NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.naik.control"); - testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,0.,0.8675309,0.,0.,0.,0.); + 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 hisq_fat_Cstyle(&GRID,path_coeff); - if (param.benchmark) { autoView(U_v, Umu, CpuRead); // Gauge accessor From b02d022993031f5fe58658a0ddb63a381cecb92f Mon Sep 17 00:00:00 2001 From: david clarke Date: Fri, 23 Feb 2024 17:14:28 -0700 Subject: [PATCH 39/40] fixed race condition (thx michael) --- Grid/qcd/smearing/HISQSmearing.h | 59 +++++++++++++++++--------------- 1 file changed, 32 insertions(+), 27 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index ac4cb8b6..6fc6993e 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -168,25 +168,26 @@ 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, SE5; - U3matrix U0, U1, U2, U3, U4, U5, W; int Nsites = U_v.size(); + auto gStencil_v = gStencil.View(); -// accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 3-link constructs - for(int site=0;site_offset; - SE1 = gStencil.GetEntry(s+1,site); int x_p_nu = SE1->_offset; - 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; + 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 @@ -212,10 +213,12 @@ public: // 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_offset; - SE1 = gStencil.GetEntry(s+1,site); int x_p_nu = SE1->_offset; - 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; + 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); @@ -248,10 +251,12 @@ public: sigmaIndex++; } } - }//) + }) -// accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 7-link - for(int site=0;site_offset; - SE1 = gStencil.GetEntry(s+1,site); int x_p_nu = SE1->_offset; - 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; + 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) { @@ -286,7 +291,7 @@ public: sigmaIndex++; } } - }//) + }) } // end mu loop From f70df6e1955107d83b62df323608f969a53c7c0e Mon Sep 17 00:00:00 2001 From: david clarke Date: Thu, 29 Feb 2024 12:29:30 -0700 Subject: [PATCH 40/40] changed NO_SHIFT and BACKWARD_CONST from define to enum --- Grid/qcd/smearing/HISQSmearing.h | 5 +---- Grid/stencil/GeneralLocalStencil.h | 24 +++++++++++++++--------- 2 files changed, 16 insertions(+), 13 deletions(-) diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h index 6fc6993e..529ea090 100644 --- a/Grid/qcd/smearing/HISQSmearing.h +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -138,7 +138,7 @@ public: for(int nu=0;nu Nd! @@ -150,16 +156,16 @@ public: 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; + return dir + shiftSignal::BACKWARD_CONST; } /*! @brief shift one unit in direction dir */ template void generalShift(Coordinate& shift, int dir) { - if (dir >= BACKWARD_CONST) { - dir -= BACKWARD_CONST; + if (dir >= shiftSignal::BACKWARD_CONST) { + dir -= shiftSignal::BACKWARD_CONST; shift[dir]+=-1; - } else if (dir == NO_SHIFT) { + } else if (dir == shiftSignal::NO_SHIFT) { ; // do nothing } else { shift[dir]+=1; @@ -169,10 +175,10 @@ void generalShift(Coordinate& shift, int dir) { /*! @brief follow a path of directions, shifting one unit in each direction */ template void generalShift(Coordinate& shift, int dir, Args... args) { - if (dir >= BACKWARD_CONST) { - dir -= BACKWARD_CONST; + if (dir >= shiftSignal::BACKWARD_CONST) { + dir -= shiftSignal::BACKWARD_CONST; shift[dir]+=-1; - } else if (dir == NO_SHIFT) { + } else if (dir == shiftSignal::NO_SHIFT) { ; // do nothing } else { shift[dir]+=1;