From 4c6613d72cb6f576e8dd01958af372df661fe369 Mon Sep 17 00:00:00 2001 From: Christopher Kelly Date: Mon, 26 Jun 2023 10:20:23 -0400 Subject: [PATCH] Modified RectStapleDouble and RectStapleOptimised to use Gauge-BC respecting CshiftLink Added test code tests/debug/Test_optimized_staple_gaugebc demonstrating equivalence of above to RectStapleUnoptimised for cconj gauge BCs Removed optimized staple only being used for periodic gauge BCs; it is now always used --- Grid/qcd/utils/WilsonLoops.h | 34 +++---- tests/debug/Test_optimized_staple_gaugebc.cc | 94 ++++++++++++++++++++ 2 files changed, 107 insertions(+), 21 deletions(-) create mode 100644 tests/debug/Test_optimized_staple_gaugebc.cc diff --git a/Grid/qcd/utils/WilsonLoops.h b/Grid/qcd/utils/WilsonLoops.h index d6efbd5d..516158d0 100644 --- a/Grid/qcd/utils/WilsonLoops.h +++ b/Grid/qcd/utils/WilsonLoops.h @@ -707,15 +707,11 @@ public: // the sum over all staples on each site ////////////////////////////////////////////////// static void RectStapleDouble(GaugeMat &U2, const GaugeMat &U, int mu) { - U2 = U * Cshift(U, mu, 1); + U2 = U * CshiftLink(U, mu, 1); } //////////////////////////////////////////////////////////////////////////// - // Hop by two optimisation strategy does not work nicely with Gparity. (could - // do, - // but need to track two deep where cross boundary and apply a conjugation). - // Must differentiate this in Gimpl, and use Gimpl::isPeriodicGaugeField to do - // so . + // Hop by two optimisation strategy. Use RectStapleDouble to obtain 'U2' //////////////////////////////////////////////////////////////////////////// static void RectStapleOptimised(GaugeMat &Stap, std::vector &U2, std::vector &U, int mu) { @@ -732,9 +728,9 @@ public: // Up staple ___ ___ // | | - tmp = Cshift(adj(U[nu]), nu, -1); + tmp = CshiftLink(adj(U[nu]), nu, -1); tmp = adj(U2[mu]) * tmp; - tmp = Cshift(tmp, mu, -2); + tmp = CshiftLink(tmp, mu, -2); Staple2x1 = Gimpl::CovShiftForward(U[nu], nu, tmp); @@ -742,14 +738,14 @@ public: // |___ ___| // tmp = adj(U2[mu]) * U[nu]; - Staple2x1 += Gimpl::CovShiftBackward(U[nu], nu, Cshift(tmp, mu, -2)); + Staple2x1 += Gimpl::CovShiftBackward(U[nu], nu, CshiftLink(tmp, mu, -2)); // ___ ___ // | ___| // |___ ___| // - Stap += Cshift(Gimpl::CovShiftForward(U[mu], mu, Staple2x1), mu, 1); + Stap += CshiftLink(Gimpl::CovShiftForward(U[mu], mu, Staple2x1), mu, 1); // ___ ___ // |___ | @@ -758,7 +754,7 @@ public: // tmp= Staple2x1* Cshift(U[mu],mu,-2); // Stap+= Cshift(tmp,mu,1) ; - Stap += Cshift(Staple2x1, mu, 1) * Cshift(U[mu], mu, -1); + Stap += CshiftLink(Staple2x1, mu, 1) * CshiftLink(U[mu], mu, -1); ; // -- @@ -766,10 +762,10 @@ public: // // | | - tmp = Cshift(adj(U2[nu]), nu, -2); + tmp = CshiftLink(adj(U2[nu]), nu, -2); tmp = Gimpl::CovShiftBackward(U[mu], mu, tmp); - tmp = U2[nu] * Cshift(tmp, nu, 2); - Stap += Cshift(tmp, mu, 1); + tmp = U2[nu] * CshiftLink(tmp, nu, 2); + Stap += CshiftLink(tmp, mu, 1); // | | // @@ -778,8 +774,8 @@ public: tmp = Gimpl::CovShiftBackward(U[mu], mu, U2[nu]); tmp = adj(U2[nu]) * tmp; - tmp = Cshift(tmp, nu, -2); - Stap += Cshift(tmp, mu, 1); + tmp = CshiftLink(tmp, nu, -2); + Stap += CshiftLink(tmp, mu, 1); } } } @@ -790,11 +786,7 @@ public: static void RectStaple(const GaugeLorentz &Umu, GaugeMat &Stap, std::vector &U2, std::vector &U, int mu) { - if (Gimpl::isPeriodicGaugeField()) { - RectStapleOptimised(Stap, U2, U, mu); - } else { - RectStapleUnoptimised(Stap, Umu, mu); - } + RectStapleOptimised(Stap, U2, U, mu); } static void RectStapleUnoptimised(GaugeMat &Stap, const GaugeLorentz &Umu, diff --git a/tests/debug/Test_optimized_staple_gaugebc.cc b/tests/debug/Test_optimized_staple_gaugebc.cc new file mode 100644 index 00000000..63822378 --- /dev/null +++ b/tests/debug/Test_optimized_staple_gaugebc.cc @@ -0,0 +1,94 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_padded_cell_staple.cc + + Copyright (C) 2015 + +Author: Christopher Kelly +Author: Peter Boyle + + 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 + *************************************************************************************/ + /* END LEGAL */ +#include +#include +#include + +using namespace std; +using namespace Grid; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layout= GridDefaultSimd(Nd,vComplexD::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + std::cout << " mpi "<({45,12,81,9})); + LatticeGaugeField U(&GRID); + + SU::HotConfiguration(pRNG,U); + + //#define PRD +#ifdef PRD + typedef PeriodicGimplD Gimpl; +#else + typedef ConjugateGimplD Gimpl; + std::vector conj_dirs(Nd,0); conj_dirs[0]=1; conj_dirs[3]=1; + Gimpl::setDirections(conj_dirs); +#endif + + typedef typename WilsonLoops::GaugeMat GaugeMat; + typedef typename WilsonLoops::GaugeLorentz GaugeLorentz; + + int count = 0; + double torig=0, topt=0; + + std::vector Umu(Nd,&GRID), U2(Nd,&GRID); + for(int mu=0;mu(U,mu); + WilsonLoops::RectStapleDouble(U2[mu], Umu[mu], mu); + } + + std::cout << GridLogMessage << "Checking optimized vs unoptimized RectStaple" << std::endl; + for(int mu=0;mu::RectStapleUnoptimised(staple_orig,U,mu); + double t1 = usecond(); + WilsonLoops::RectStapleOptimised(staple_opt, U2, Umu, mu); + double t2 = usecond(); + torig += t1-t0; topt += t2-t1; + ++count; + + GaugeMat diff = staple_orig - staple_opt; + double n = norm2(diff); + std::cout << GridLogMessage << mu << " " << n << std::endl; + assert(n<1e-10); + } + std::cout << GridLogMessage << "RectStaple timings orig: " << torig/1000/count << "ms, optimized: " << topt/1000/count << "ms" << std::endl; + + Grid_finalize(); +}