From dbe210dd53405303cf57374c9b658902cbf8072a Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 25 Apr 2021 10:25:59 -0400 Subject: [PATCH 01/31] Open the ens_id --- Grid/parallelIO/NerscIO.h | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/Grid/parallelIO/NerscIO.h b/Grid/parallelIO/NerscIO.h index 3ebdf0cc..1ffda074 100644 --- a/Grid/parallelIO/NerscIO.h +++ b/Grid/parallelIO/NerscIO.h @@ -205,11 +205,22 @@ public: std::cout< + static inline void writeConfiguration(Lattice &Umu, + std::string file, + std::string ens_id = std::string("UKQCD"), + std::string ens_label = std::string("DWF")) + { + writeConfiguration(Umu,file,0,1,ens_id,ens_label); + } template static inline void writeConfiguration(Lattice &Umu, std::string file, int two_row, - int bits32) + int bits32, + std::string ens_id = std::string("UKQCD"), + std::string ens_label = std::string("DWF")) { typedef vLorentzColourMatrixD vobj; typedef typename vobj::scalar_object sobj; @@ -219,8 +230,8 @@ public: // Following should become arguments /////////////////////////////////////////// header.sequence_number = 1; - header.ensemble_id = "UKQCD"; - header.ensemble_label = "DWF"; + header.ensemble_id = ens_id; + header.ensemble_label = ens_label; typedef LorentzColourMatrixD fobj3D; typedef LorentzColour2x3D fobj2D; @@ -232,7 +243,7 @@ public: GaugeStats Stats; Stats(Umu,header); MachineCharacteristics(header); - uint64_t offset; + uint64_t offset; // Sod it -- always write 3x3 double header.floating_point = std::string("IEEE64BIG"); From 955a8113ded53fc55c204fc107832933ad120c2c Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 25 Apr 2021 10:36:38 -0400 Subject: [PATCH 02/31] Expose label only to reduce number of parameters --- Grid/parallelIO/NerscIO.h | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/Grid/parallelIO/NerscIO.h b/Grid/parallelIO/NerscIO.h index 1ffda074..04aae5d8 100644 --- a/Grid/parallelIO/NerscIO.h +++ b/Grid/parallelIO/NerscIO.h @@ -209,7 +209,6 @@ public: template static inline void writeConfiguration(Lattice &Umu, std::string file, - std::string ens_id = std::string("UKQCD"), std::string ens_label = std::string("DWF")) { writeConfiguration(Umu,file,0,1,ens_id,ens_label); @@ -219,7 +218,6 @@ public: std::string file, int two_row, int bits32, - std::string ens_id = std::string("UKQCD"), std::string ens_label = std::string("DWF")) { typedef vLorentzColourMatrixD vobj; @@ -230,7 +228,7 @@ public: // Following should become arguments /////////////////////////////////////////// header.sequence_number = 1; - header.ensemble_id = ens_id; + header.ensemble_id = std::string("UKQCD"); header.ensemble_label = ens_label; typedef LorentzColourMatrixD fobj3D; From d45c868656b4e28452946f90a6e71cdf12c21cf2 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 25 Apr 2021 10:53:34 -0400 Subject: [PATCH 03/31] Change interface --- Grid/parallelIO/NerscIO.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Grid/parallelIO/NerscIO.h b/Grid/parallelIO/NerscIO.h index 04aae5d8..99011e25 100644 --- a/Grid/parallelIO/NerscIO.h +++ b/Grid/parallelIO/NerscIO.h @@ -211,7 +211,7 @@ public: std::string file, std::string ens_label = std::string("DWF")) { - writeConfiguration(Umu,file,0,1,ens_id,ens_label); + writeConfiguration(Umu,file,0,1,ens_label); } template static inline void writeConfiguration(Lattice &Umu, From 8cd4263974060f9af3b002604d9036e2552cc307 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 25 Apr 2021 22:20:37 -0400 Subject: [PATCH 04/31] Tests compile --- tests/debug/Test_heatbath_dwf_eofa.cc | 6 ++++-- tests/debug/Test_heatbath_dwf_eofa_gparity.cc | 7 +++++-- tests/debug/Test_heatbath_mobius_eofa.cc | 6 ++++-- tests/debug/Test_heatbath_mobius_eofa_gparity.cc | 6 ++++-- tests/forces/Test_momentum_filter.cc | 4 +++- 5 files changed, 20 insertions(+), 9 deletions(-) diff --git a/tests/debug/Test_heatbath_dwf_eofa.cc b/tests/debug/Test_heatbath_dwf_eofa.cc index 9d453a96..e1c18021 100644 --- a/tests/debug/Test_heatbath_dwf_eofa.cc +++ b/tests/debug/Test_heatbath_dwf_eofa.cc @@ -66,7 +66,9 @@ int main(int argc, char** argv) // Set up RNGs std::vector seeds4({1, 2, 3, 4}); std::vector seeds5({5, 6, 7, 8}); + GridSerialRNG sRNG; GridParallelRNG RNG5(FGrid); + sRNG.SeedFixedIntegers(seeds5); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); @@ -84,7 +86,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu,sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -94,7 +96,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu,sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } diff --git a/tests/debug/Test_heatbath_dwf_eofa_gparity.cc b/tests/debug/Test_heatbath_dwf_eofa_gparity.cc index 22cc1e90..7eabfc65 100644 --- a/tests/debug/Test_heatbath_dwf_eofa_gparity.cc +++ b/tests/debug/Test_heatbath_dwf_eofa_gparity.cc @@ -74,6 +74,9 @@ int main(int argc, char** argv) RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridSerialRNG sRNG; + RNG4.SeedFixedIntegers(seeds4); + sRNG.SeedFixedIntegers(seeds5); // Random gauge field LatticeGaugeField Umu(UGrid); @@ -90,7 +93,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu,sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -100,7 +103,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu,sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } diff --git a/tests/debug/Test_heatbath_mobius_eofa.cc b/tests/debug/Test_heatbath_mobius_eofa.cc index 4cf4bf53..48806642 100644 --- a/tests/debug/Test_heatbath_mobius_eofa.cc +++ b/tests/debug/Test_heatbath_mobius_eofa.cc @@ -68,8 +68,10 @@ int main(int argc, char** argv) // Set up RNGs std::vector seeds4({1, 2, 3, 4}); std::vector seeds5({5, 6, 7, 8}); + GridSerialRNG sRNG; GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + sRNG.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); @@ -86,7 +88,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, sRNG,RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -96,7 +98,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, sRNG,RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } diff --git a/tests/debug/Test_heatbath_mobius_eofa_gparity.cc b/tests/debug/Test_heatbath_mobius_eofa_gparity.cc index 2fcb4b9f..52447e5e 100644 --- a/tests/debug/Test_heatbath_mobius_eofa_gparity.cc +++ b/tests/debug/Test_heatbath_mobius_eofa_gparity.cc @@ -73,7 +73,9 @@ int main(int argc, char** argv) std::vector seeds4({1, 2, 3, 4}); std::vector seeds5({5, 6, 7, 8}); GridParallelRNG RNG5(FGrid); + GridSerialRNG sRNG; RNG5.SeedFixedIntegers(seeds5); + sRNG.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); @@ -91,7 +93,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -101,7 +103,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, sRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } diff --git a/tests/forces/Test_momentum_filter.cc b/tests/forces/Test_momentum_filter.cc index 856ea0f2..794b5fa0 100644 --- a/tests/forces/Test_momentum_filter.cc +++ b/tests/forces/Test_momentum_filter.cc @@ -61,7 +61,9 @@ int main (int argc, char ** argv) std::vector seeds({1,2,3,4}); GridParallelRNG pRNG(&Grid); + GridSerialRNG sRNG; pRNG.SeedFixedIntegers(seeds); + sRNG.SeedFixedIntegers(seeds); typedef PeriodicGimplR Gimpl; typedef WilsonGaugeAction GaugeAction; @@ -115,7 +117,7 @@ int main (int argc, char ** argv) integrator.setMomentumFilter(filter); - integrator.refresh(U, pRNG); //doesn't actually change the gauge field + integrator.refresh(U, sRNG, pRNG); //doesn't actually change the gauge field //Check the momentum is zero on the boundary const auto &P = integrator.getMomentum(); From 009ccd581ede8faf0ba748fa49a1757419106e23 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Mon, 26 Apr 2021 10:36:33 +0100 Subject: [PATCH 05/31] bugfix 3D stout smearing --- Grid/qcd/smearing/StoutSmearing.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Grid/qcd/smearing/StoutSmearing.h b/Grid/qcd/smearing/StoutSmearing.h index ed2ccdb6..629f81e2 100644 --- a/Grid/qcd/smearing/StoutSmearing.h +++ b/Grid/qcd/smearing/StoutSmearing.h @@ -89,11 +89,12 @@ public: SmearBase->smear(C, U); for (int mu = 0; mu < Nd; mu++) { - if( mu == OrthogDim ) + Umu = peekLorentz(U, mu); + if( mu == OrthogDim ){ tmp = 1.0; // Don't smear in the orthogonal direction + } else { tmp = peekLorentz(C, mu); - Umu = peekLorentz(U, mu); iq_mu = Ta( tmp * adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper From cf2923d5ddb9c190ebb78efadab281e5a06ba247 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 27 Apr 2021 16:53:37 +0100 Subject: [PATCH 06/31] Jamie's fix --- Grid/qcd/smearing/StoutSmearing.h | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/Grid/qcd/smearing/StoutSmearing.h b/Grid/qcd/smearing/StoutSmearing.h index 629f81e2..6ee78e8c 100644 --- a/Grid/qcd/smearing/StoutSmearing.h +++ b/Grid/qcd/smearing/StoutSmearing.h @@ -85,22 +85,18 @@ public: std::cout << GridLogDebug << "Stout smearing started\n"; - // Smear the configurations + // C contains the staples multiplied by some rho + u_smr = U ; // set the smeared field to the current gauge field SmearBase->smear(C, U); for (int mu = 0; mu < Nd; mu++) { + if( mu == OrthogDim ) continue ; + // u_smr = exp(iQ_mu)*U_mu apart from Orthogdim Umu = peekLorentz(U, mu); - if( mu == OrthogDim ){ - tmp = 1.0; // Don't smear in the orthogonal direction - } - else { - tmp = peekLorentz(C, mu); - iq_mu = Ta( - tmp * - adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper - exponentiate_iQ(tmp, iq_mu); - } - pokeLorentz(u_smr, tmp * Umu, mu); // u_smr = exp(iQ_mu)*U_mu + tmp = peekLorentz(C, mu); + iq_mu = Ta( tmp * adj(Umu)); + exponentiate_iQ(tmp, iq_mu); + pokeLorentz(u_smr, tmp * Umu, mu); } std::cout << GridLogDebug << "Stout smearing completed\n"; }; From 834f536b5f426aa0b3a334a89b37d8da39fb4238 Mon Sep 17 00:00:00 2001 From: u61464 Date: Tue, 4 May 2021 08:40:18 -0700 Subject: [PATCH 07/31] Fastest option on SyCL is now std::complex --- Grid/tensors/Tensor_SIMT.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Grid/tensors/Tensor_SIMT.h b/Grid/tensors/Tensor_SIMT.h index 672f385f..0a7d3382 100644 --- a/Grid/tensors/Tensor_SIMT.h +++ b/Grid/tensors/Tensor_SIMT.h @@ -65,7 +65,8 @@ void coalescedWriteNonTemporal(vobj & __restrict__ vec,const vobj & __restrict__ #else -#ifndef GRID_SYCL +//#ifndef GRID_SYCL +#if 1 // Use the scalar as our own complex on GPU ... thrust::complex or std::complex template = 0> accelerator_inline typename vsimd::scalar_type From 8cfc7342cde4b93d1de8f41a6f909ddd6d5d351f Mon Sep 17 00:00:00 2001 From: u61464 Date: Wed, 5 May 2021 14:17:18 -0700 Subject: [PATCH 08/31] staggered hand unroll read coalesce --- Grid/qcd/action/fermion/Fermion.h | 6 - Grid/qcd/action/fermion/FermionOperatorImpl.h | 5 +- .../implementation/StaggeredKernelsAsm.h | 5 +- .../implementation/StaggeredKernelsHand.h | 203 +++++++++--------- 4 files changed, 106 insertions(+), 113 deletions(-) diff --git a/Grid/qcd/action/fermion/Fermion.h b/Grid/qcd/action/fermion/Fermion.h index 16252340..09777204 100644 --- a/Grid/qcd/action/fermion/Fermion.h +++ b/Grid/qcd/action/fermion/Fermion.h @@ -291,12 +291,6 @@ typedef ImprovedStaggeredFermion5D ImprovedStaggeredFermion5DR; typedef ImprovedStaggeredFermion5D ImprovedStaggeredFermion5DF; typedef ImprovedStaggeredFermion5D ImprovedStaggeredFermion5DD; -#ifndef GRID_CUDA -typedef ImprovedStaggeredFermion5D ImprovedStaggeredFermionVec5dR; -typedef ImprovedStaggeredFermion5D ImprovedStaggeredFermionVec5dF; -typedef ImprovedStaggeredFermion5D ImprovedStaggeredFermionVec5dD; -#endif - NAMESPACE_END(Grid); //////////////////// diff --git a/Grid/qcd/action/fermion/FermionOperatorImpl.h b/Grid/qcd/action/fermion/FermionOperatorImpl.h index 9345c0e6..56aaca12 100644 --- a/Grid/qcd/action/fermion/FermionOperatorImpl.h +++ b/Grid/qcd/action/fermion/FermionOperatorImpl.h @@ -183,7 +183,8 @@ NAMESPACE_CHECK(ImplStaggered); ///////////////////////////////////////////////////////////////////////////// // Single flavour one component spinors with colour index. 5d vec ///////////////////////////////////////////////////////////////////////////// -#include -NAMESPACE_CHECK(ImplStaggered5dVec); +// Deprecate Vec5d +//#include +//NAMESPACE_CHECK(ImplStaggered5dVec); diff --git a/Grid/qcd/action/fermion/implementation/StaggeredKernelsAsm.h b/Grid/qcd/action/fermion/implementation/StaggeredKernelsAsm.h index 63fd2a2f..e9cacbcf 100644 --- a/Grid/qcd/action/fermion/implementation/StaggeredKernelsAsm.h +++ b/Grid/qcd/action/fermion/implementation/StaggeredKernelsAsm.h @@ -680,7 +680,8 @@ void StaggeredKernels::DhopSiteAsm(StencilView &st, gauge2 =(uint64_t)&UU[sU]( Z ); \ gauge3 =(uint64_t)&UU[sU]( T ); - +#undef STAG_VEC5D +#ifdef STAG_VEC5D // This is the single precision 5th direction vectorised kernel #include template <> void StaggeredKernels::DhopSiteAsm(StencilView &st, @@ -790,7 +791,7 @@ template <> void StaggeredKernels::DhopSiteAsm(StencilView #endif } - +#endif #define PERMUTE_DIR3 __asm__ ( \ diff --git a/Grid/qcd/action/fermion/implementation/StaggeredKernelsHand.h b/Grid/qcd/action/fermion/implementation/StaggeredKernelsHand.h index 6bcb22b4..2b6087bc 100644 --- a/Grid/qcd/action/fermion/implementation/StaggeredKernelsHand.h +++ b/Grid/qcd/action/fermion/implementation/StaggeredKernelsHand.h @@ -32,25 +32,50 @@ Author: paboyle NAMESPACE_BEGIN(Grid); -#define LOAD_CHI(b) \ +#ifdef GRID_SIMT + +#define LOAD_CHI(ptype,b) \ + const SiteSpinor & ref (b[offset]); \ + Chi_0=coalescedReadPermute(ref()()(0),perm,lane); \ + Chi_1=coalescedReadPermute(ref()()(1),perm,lane); \ + Chi_2=coalescedReadPermute(ref()()(2),perm,lane); + +#define LOAD_CHI_COMMS(b) \ const SiteSpinor & ref (b[offset]); \ - Chi_0=ref()()(0);\ - Chi_1=ref()()(1);\ - Chi_2=ref()()(2); + Chi_0=coalescedRead(ref()()(0),lane); \ + Chi_1=coalescedRead(ref()()(1),lane); \ + Chi_2=coalescedRead(ref()()(2),lane); + +#define PERMUTE_DIR(dir) ; +#else +#define LOAD_CHI(ptype,b) LOAD_CHI_COMMS(b) + +#define LOAD_CHI_COMMS(b) \ + const SiteSpinor & ref (b[offset]); \ + Chi_0=ref()()(0); \ + Chi_1=ref()()(1); \ + Chi_2=ref()()(2); + +#define PERMUTE_DIR(dir) \ + permute##dir(Chi_0,Chi_0); \ + permute##dir(Chi_1,Chi_1); \ + permute##dir(Chi_2,Chi_2); + +#endif // To splat or not to splat depends on the implementation #define MULT(A,UChi) \ auto & ref(U[sU](A)); \ - Impl::loadLinkElement(U_00,ref()(0,0)); \ - Impl::loadLinkElement(U_10,ref()(1,0)); \ - Impl::loadLinkElement(U_20,ref()(2,0)); \ - Impl::loadLinkElement(U_01,ref()(0,1)); \ - Impl::loadLinkElement(U_11,ref()(1,1)); \ - Impl::loadLinkElement(U_21,ref()(2,1)); \ - Impl::loadLinkElement(U_02,ref()(0,2)); \ - Impl::loadLinkElement(U_12,ref()(1,2)); \ - Impl::loadLinkElement(U_22,ref()(2,2)); \ + U_00=coalescedRead(ref()(0,0),lane); \ + U_10=coalescedRead(ref()(1,0),lane); \ + U_20=coalescedRead(ref()(2,0),lane); \ + U_01=coalescedRead(ref()(0,1),lane); \ + U_11=coalescedRead(ref()(1,1),lane); \ + U_21=coalescedRead(ref()(2,1),lane); \ + U_02=coalescedRead(ref()(0,2),lane); \ + U_12=coalescedRead(ref()(1,2),lane); \ + U_22=coalescedRead(ref()(2,2),lane); \ UChi ## _0 = U_00*Chi_0; \ UChi ## _1 = U_10*Chi_0;\ UChi ## _2 = U_20*Chi_0;\ @@ -63,15 +88,15 @@ NAMESPACE_BEGIN(Grid); #define MULT_ADD(U,A,UChi) \ auto & ref(U[sU](A)); \ - Impl::loadLinkElement(U_00,ref()(0,0)); \ - Impl::loadLinkElement(U_10,ref()(1,0)); \ - Impl::loadLinkElement(U_20,ref()(2,0)); \ - Impl::loadLinkElement(U_01,ref()(0,1)); \ - Impl::loadLinkElement(U_11,ref()(1,1)); \ - Impl::loadLinkElement(U_21,ref()(2,1)); \ - Impl::loadLinkElement(U_02,ref()(0,2)); \ - Impl::loadLinkElement(U_12,ref()(1,2)); \ - Impl::loadLinkElement(U_22,ref()(2,2)); \ + U_00=coalescedRead(ref()(0,0),lane); \ + U_10=coalescedRead(ref()(1,0),lane); \ + U_20=coalescedRead(ref()(2,0),lane); \ + U_01=coalescedRead(ref()(0,1),lane); \ + U_11=coalescedRead(ref()(1,1),lane); \ + U_21=coalescedRead(ref()(2,1),lane); \ + U_02=coalescedRead(ref()(0,2),lane); \ + U_12=coalescedRead(ref()(1,2),lane); \ + U_22=coalescedRead(ref()(2,2),lane); \ UChi ## _0 += U_00*Chi_0; \ UChi ## _1 += U_10*Chi_0;\ UChi ## _2 += U_20*Chi_0;\ @@ -83,24 +108,18 @@ NAMESPACE_BEGIN(Grid); UChi ## _2 += U_22*Chi_2; -#define PERMUTE_DIR(dir) \ - permute##dir(Chi_0,Chi_0); \ - permute##dir(Chi_1,Chi_1); \ - permute##dir(Chi_2,Chi_2); - - #define HAND_STENCIL_LEG_BASE(Dir,Perm,skew) \ SE=st.GetEntry(ptype,Dir+skew,sF); \ offset = SE->_offset; \ local = SE->_is_local; \ perm = SE->_permute; \ if ( local ) { \ - LOAD_CHI(in); \ + LOAD_CHI(Perm,in); \ if ( perm) { \ PERMUTE_DIR(Perm); \ } \ } else { \ - LOAD_CHI(buf); \ + LOAD_CHI_COMMS(buf); \ } #define HAND_STENCIL_LEG_BEGIN(Dir,Perm,skew,even) \ @@ -116,19 +135,18 @@ NAMESPACE_BEGIN(Grid); } - #define HAND_STENCIL_LEG_INT(U,Dir,Perm,skew,even) \ SE=st.GetEntry(ptype,Dir+skew,sF); \ offset = SE->_offset; \ local = SE->_is_local; \ perm = SE->_permute; \ if ( local ) { \ - LOAD_CHI(in); \ + LOAD_CHI(Perm,in); \ if ( perm) { \ PERMUTE_DIR(Perm); \ } \ } else if ( st.same_node[Dir] ) { \ - LOAD_CHI(buf); \ + LOAD_CHI_COMMS(buf); \ } \ if (local || st.same_node[Dir] ) { \ MULT_ADD(U,Dir,even); \ @@ -140,10 +158,32 @@ NAMESPACE_BEGIN(Grid); local = SE->_is_local; \ if ((!local) && (!st.same_node[Dir]) ) { \ nmu++; \ - { LOAD_CHI(buf); } \ + { LOAD_CHI_COMMS(buf); } \ { MULT_ADD(U,Dir,even); } \ } +#define HAND_DECLARATIONS(Simd) \ + Simd even_0; \ + Simd even_1; \ + Simd even_2; \ + Simd odd_0; \ + Simd odd_1; \ + Simd odd_2; \ + \ + Simd Chi_0; \ + Simd Chi_1; \ + Simd Chi_2; \ + \ + Simd U_00; \ + Simd U_10; \ + Simd U_20; \ + Simd U_01; \ + Simd U_11; \ + Simd U_21; \ + Simd U_02; \ + Simd U_12; \ + Simd U_22; + template template accelerator_inline @@ -155,28 +195,14 @@ void StaggeredKernels::DhopSiteHand(StencilView &st, typedef typename Simd::scalar_type S; typedef typename Simd::vector_type V; - Simd even_0; // 12 regs on knc - Simd even_1; - Simd even_2; - Simd odd_0; // 12 regs on knc - Simd odd_1; - Simd odd_2; - Simd Chi_0; // two spinor; 6 regs - Simd Chi_1; - Simd Chi_2; - - Simd U_00; // two rows of U matrix - Simd U_10; - Simd U_20; - Simd U_01; - Simd U_11; - Simd U_21; // 2 reg left. - Simd U_02; - Simd U_12; - Simd U_22; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); + typedef decltype( coalescedRead( in[0]()()(0) )) Simt; + HAND_DECLARATIONS(Simt); - SiteSpinor result; + typedef decltype( coalescedRead( in[0] )) calcSiteSpinor; + calcSiteSpinor result; int offset,local,perm, ptype; StencilEntry *SE; @@ -215,7 +241,7 @@ void StaggeredKernels::DhopSiteHand(StencilView &st, result()()(1) = even_1 + odd_1; result()()(2) = even_2 + odd_2; } - vstream(out[sF],result); + coalescedWrite(out[sF],result); } } @@ -230,28 +256,13 @@ void StaggeredKernels::DhopSiteHandInt(StencilView &st, typedef typename Simd::scalar_type S; typedef typename Simd::vector_type V; - Simd even_0; // 12 regs on knc - Simd even_1; - Simd even_2; - Simd odd_0; // 12 regs on knc - Simd odd_1; - Simd odd_2; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); + typedef decltype( coalescedRead( in[0]()()(0) )) Simt; + HAND_DECLARATIONS(Simt); - Simd Chi_0; // two spinor; 6 regs - Simd Chi_1; - Simd Chi_2; - - Simd U_00; // two rows of U matrix - Simd U_10; - Simd U_20; - Simd U_01; - Simd U_11; - Simd U_21; // 2 reg left. - Simd U_02; - Simd U_12; - Simd U_22; - - SiteSpinor result; + typedef decltype( coalescedRead( in[0] )) calcSiteSpinor; + calcSiteSpinor result; int offset, ptype, local, perm; StencilEntry *SE; @@ -261,8 +272,8 @@ void StaggeredKernels::DhopSiteHandInt(StencilView &st, // int sF=s+LLs*sU; { - even_0 = Zero(); even_1 = Zero(); even_2 = Zero(); - odd_0 = Zero(); odd_1 = Zero(); odd_2 = Zero(); + zeroit(even_0); zeroit(even_1); zeroit(even_2); + zeroit(odd_0); zeroit(odd_1); zeroit(odd_2); skew = 0; HAND_STENCIL_LEG_INT(U,Xp,3,skew,even); @@ -294,7 +305,7 @@ void StaggeredKernels::DhopSiteHandInt(StencilView &st, result()()(1) = even_1 + odd_1; result()()(2) = even_2 + odd_2; } - vstream(out[sF],result); + coalescedWrite(out[sF],result); } } @@ -309,28 +320,13 @@ void StaggeredKernels::DhopSiteHandExt(StencilView &st, typedef typename Simd::scalar_type S; typedef typename Simd::vector_type V; - Simd even_0; // 12 regs on knc - Simd even_1; - Simd even_2; - Simd odd_0; // 12 regs on knc - Simd odd_1; - Simd odd_2; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); + typedef decltype( coalescedRead( in[0]()()(0) )) Simt; + HAND_DECLARATIONS(Simt); - Simd Chi_0; // two spinor; 6 regs - Simd Chi_1; - Simd Chi_2; - - Simd U_00; // two rows of U matrix - Simd U_10; - Simd U_20; - Simd U_01; - Simd U_11; - Simd U_21; // 2 reg left. - Simd U_02; - Simd U_12; - Simd U_22; - - SiteSpinor result; + typedef decltype( coalescedRead( in[0] )) calcSiteSpinor; + calcSiteSpinor result; int offset, ptype, local; StencilEntry *SE; @@ -340,8 +336,8 @@ void StaggeredKernels::DhopSiteHandExt(StencilView &st, // int sF=s+LLs*sU; { - even_0 = Zero(); even_1 = Zero(); even_2 = Zero(); - odd_0 = Zero(); odd_1 = Zero(); odd_2 = Zero(); + zeroit(even_0); zeroit(even_1); zeroit(even_2); + zeroit(odd_0); zeroit(odd_1); zeroit(odd_2); int nmu=0; skew = 0; HAND_STENCIL_LEG_EXT(U,Xp,3,skew,even); @@ -374,7 +370,7 @@ void StaggeredKernels::DhopSiteHandExt(StencilView &st, result()()(1) = even_1 + odd_1; result()()(2) = even_2 + odd_2; } - out[sF] = out[sF] + result; + coalescedWrite(out[sF] , out(sF)+ result); } } } @@ -397,6 +393,7 @@ void StaggeredKernels::DhopSiteHandExt(StencilView &st, const FermionFieldView &in, FermionFieldView &out, int dag); \ */ #undef LOAD_CHI +#undef HAND_DECLARATIONS NAMESPACE_END(Grid); From 0e27e3847d6252c7a950e59517a8af0bf1e15549 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 3 Jun 2021 04:24:19 +0000 Subject: [PATCH 09/31] Remove synch --- Grid/threads/Accelerator.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Grid/threads/Accelerator.h b/Grid/threads/Accelerator.h index 56b85c72..b76d6d1c 100644 --- a/Grid/threads/Accelerator.h +++ b/Grid/threads/Accelerator.h @@ -457,7 +457,7 @@ accelerator_inline void acceleratorSynchronise(void) __syncwarp(); #endif #ifdef GRID_SYCL - cl::sycl::detail::workGroupBarrier(); + //cl::sycl::detail::workGroupBarrier(); #endif #ifdef GRID_HIP __syncthreads(); From ca10bfa1c7f6615bcc13322ee5e050c6e4b37ad8 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 4 Jun 2021 11:12:22 +0100 Subject: [PATCH 10/31] removing Travis CI constantly failing due to overtime (no way we can compile Grid on free time anymore) --- .travis.yml | 56 ----------------------------------------------------- 1 file changed, 56 deletions(-) delete mode 100644 .travis.yml diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index 3a0e1e35..00000000 --- a/.travis.yml +++ /dev/null @@ -1,56 +0,0 @@ -language: cpp - -cache: - directories: - - clang - -matrix: - include: - - os: osx - osx_image: xcode8.3 - compiler: clang - -before_install: - - export GRIDDIR=`pwd` - - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]] && [ ! -e clang/bin ]; then wget $CLANG_LINK; tar -xf `basename $CLANG_LINK`; mkdir clang; mv clang+*/* clang/; fi - - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]]; then export PATH="${GRIDDIR}/clang/bin:${PATH}"; fi - - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]]; then export LD_LIBRARY_PATH="${GRIDDIR}/clang/lib:${LD_LIBRARY_PATH}"; fi - - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew update; fi - - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install libmpc openssl; fi - -install: - - export CWD=`pwd` - - echo $CWD - - export CC=$CC$VERSION - - export CXX=$CXX$VERSION - - echo $PATH - - which autoconf - - autoconf --version - - which automake - - automake --version - - which $CC - - $CC --version - - which $CXX - - $CXX --version - - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then export LDFLAGS='-L/usr/local/lib'; fi - - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then export EXTRACONF='--with-openssl=/usr/local/opt/openssl'; fi - -script: - - ./bootstrap.sh - - mkdir build - - cd build - - mkdir lime - - cd lime - - mkdir build - - cd build - - wget http://usqcd-software.github.io/downloads/c-lime/lime-1.3.2.tar.gz - - tar xf lime-1.3.2.tar.gz - - cd lime-1.3.2 - - ./configure --prefix=$CWD/build/lime/install - - make -j4 - - make install - - cd $CWD/build - - ../configure --enable-simd=SSE4 --enable-comms=none --with-lime=$CWD/build/lime/install ${EXTRACONF} - - make -j4 - - ./benchmarks/Benchmark_dwf --threads 1 --debug-signals - - make check From 92def28bd3331153da2b8a2414f471e4f7831a4c Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 6 Jun 2021 04:52:05 -0400 Subject: [PATCH 11/31] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index fff68dc6..88b922a5 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# Grid [![Teamcity status](http://ci.cliath.ph.ed.ac.uk/app/rest/builds/aggregated/strob:(buildType:(affectedProject(id:GridBasedSoftware_Grid)),branch:name:develop)/statusIcon.svg)](http://ci.cliath.ph.ed.ac.uk/project.html?projectId=GridBasedSoftware_Grid&tab=projectOverview) [![Travis status](https://travis-ci.org/paboyle/Grid.svg?branch=develop)](https://travis-ci.org/paboyle/Grid) +# Grid [![Teamcity status](http://ci.cliath.ph.ed.ac.uk/app/rest/builds/aggregated/strob:(buildType:(affectedProject(id:GridBasedSoftware_Grid)),branch:name:develop)/statusIcon.svg)](http://ci.cliath.ph.ed.ac.uk/project.html?projectId=GridBasedSoftware_Grid&tab=projectOverview) **Data parallel C++ mathematical object library.** From 4c5440fb0678b3a936ebe95f2c891d90b62feaaf Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 15 Jun 2021 21:45:07 +0000 Subject: [PATCH 12/31] const happy for sycl --- Grid/tensors/Tensor_extract_merge.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Grid/tensors/Tensor_extract_merge.h b/Grid/tensors/Tensor_extract_merge.h index f1ded209..ea619d0f 100644 --- a/Grid/tensors/Tensor_extract_merge.h +++ b/Grid/tensors/Tensor_extract_merge.h @@ -1,5 +1,5 @@ /************************************************************************************* - +n Grid physics library, www.github.com/paboyle/Grid Source file: ./lib/tensors/Tensor_extract_merge.h @@ -153,7 +153,7 @@ void insertLane(int lane, vobj & __restrict__ vec,const typename vobj::scalar_ob // Extract to a bunch of scalar object pointers of different scalar type, with offset. Useful for precision change //////////////////////////////////////////////////////////////////////// template accelerator -void extract(const vobj &vec,ExtractPointerArray &extracted, int offset) +void extract(const vobj &vec,const ExtractPointerArray &extracted, int offset) { typedef typename GridTypeMapper::scalar_type sobj_scalar_type; typedef typename GridTypeMapper::scalar_type scalar_type; @@ -181,7 +181,7 @@ void extract(const vobj &vec,ExtractPointerArray &extracted, int offset) // Merge bunch of scalar object pointers of different scalar type, with offset. Useful for precision change //////////////////////////////////////////////////////////////////////// template accelerator -void merge(vobj &vec,ExtractPointerArray &extracted, int offset) +void merge(vobj &vec,const ExtractPointerArray &extracted, int offset) { typedef typename GridTypeMapper::scalar_type sobj_scalar_type; typedef typename GridTypeMapper::scalar_type scalar_type; From 6cd9224dd78aca959d3997479287cf943832d79c Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 16 Jun 2021 17:10:55 +0000 Subject: [PATCH 13/31] SYCL comms buffer allocate --- Grid/communicator/SharedMemoryMPI.cc | 54 +++++++++++++++++++++++++--- 1 file changed, 49 insertions(+), 5 deletions(-) diff --git a/Grid/communicator/SharedMemoryMPI.cc b/Grid/communicator/SharedMemoryMPI.cc index 466f6a1e..786122fa 100644 --- a/Grid/communicator/SharedMemoryMPI.cc +++ b/Grid/communicator/SharedMemoryMPI.cc @@ -35,6 +35,9 @@ Author: Christoph Lehner #endif #ifdef GRID_HIP #include +#endif +#ifdef GRID_SYCl + #endif NAMESPACE_BEGIN(Grid); @@ -446,7 +449,46 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) //////////////////////////////////////////////////////////////////////////////////////////// // Hugetlbfs mapping intended //////////////////////////////////////////////////////////////////////////////////////////// -#if defined(GRID_CUDA) ||defined(GRID_HIP) +#if defined(GRID_CUDA) ||defined(GRID_HIP) || defined(GRID_SYCL) + +#if defined(GRID_SYCL) +void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) +{ + void * ShmCommBuf ; + assert(_ShmSetup==1); + assert(_ShmAlloc==0); + + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + // allocate the pointer array for shared windows for our group + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + MPI_Barrier(WorldShmComm); + WorldShmCommBufs.resize(WorldShmSize); + + /////////////////////////////////////////////////////////////////////////////////////////////////////////// + // Each MPI rank should allocate our own buffer + /////////////////////////////////////////////////////////////////////////////////////////////////////////// + ShmCommBuf = acceleratorAllocDevice(bytes); + + if (ShmCommBuf == (void *)NULL ) { + std::cerr << " SharedMemoryMPI.cc acceleratorAllocDevice failed NULL pointer for " << bytes<<" bytes " << std::endl; + exit(EXIT_FAILURE); + } + // if ( WorldRank == 0 ){ + if ( 1 ){ + std::cout << WorldRank << header " SharedMemoryMPI.cc acceleratorAllocDevice "<< bytes + << "bytes at "<< std::hex<< ShmCommBuf < Date: Tue, 22 Jun 2021 17:56:10 +0000 Subject: [PATCH 14/31] Force reqd subgroup size fo SYCL --- Grid/threads/Accelerator.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Grid/threads/Accelerator.h b/Grid/threads/Accelerator.h index b76d6d1c..c0af1019 100644 --- a/Grid/threads/Accelerator.h +++ b/Grid/threads/Accelerator.h @@ -257,11 +257,14 @@ accelerator_inline int acceleratorSIMTlane(int Nsimd) { unsigned long nt=acceleratorThreads(); \ unsigned long unum1 = num1; \ unsigned long unum2 = num2; \ + if(nt < 8)nt=8; \ cl::sycl::range<3> local {nt,1,nsimd}; \ cl::sycl::range<3> global{unum1,unum2,nsimd}; \ cgh.parallel_for( \ cl::sycl::nd_range<3>(global,local), \ - [=] (cl::sycl::nd_item<3> item) /*mutable*/ { \ + [=] (cl::sycl::nd_item<3> item) /*mutable*/ \ + [[intel::reqd_sub_group_size(8)]] \ + { \ auto iter1 = item.get_global_id(0); \ auto iter2 = item.get_global_id(1); \ auto lane = item.get_global_id(2); \ From 29a22ae603a3cf18d2ebeba2eb5aabcf27fe3e5d Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 22 Jun 2021 17:57:20 +0000 Subject: [PATCH 15/31] Simpler SYCL setup --- Grid/communicator/SharedMemoryMPI.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Grid/communicator/SharedMemoryMPI.cc b/Grid/communicator/SharedMemoryMPI.cc index 786122fa..caa03a60 100644 --- a/Grid/communicator/SharedMemoryMPI.cc +++ b/Grid/communicator/SharedMemoryMPI.cc @@ -473,13 +473,13 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) std::cerr << " SharedMemoryMPI.cc acceleratorAllocDevice failed NULL pointer for " << bytes<<" bytes " << std::endl; exit(EXIT_FAILURE); } - // if ( WorldRank == 0 ){ - if ( 1 ){ - std::cout << WorldRank << header " SharedMemoryMPI.cc acceleratorAllocDevice "<< bytes - << "bytes at "<< std::hex<< ShmCommBuf < Date: Thu, 5 Aug 2021 18:33:20 -0400 Subject: [PATCH 16/31] Check is wrong (HtoD / DtoH) --- Grid/communicator/SharedMemoryMPI.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Grid/communicator/SharedMemoryMPI.cc b/Grid/communicator/SharedMemoryMPI.cc index caa03a60..d7b41cee 100644 --- a/Grid/communicator/SharedMemoryMPI.cc +++ b/Grid/communicator/SharedMemoryMPI.cc @@ -844,7 +844,7 @@ void SharedMemory::SetCommunicator(Grid_MPI_Comm comm) } #endif - SharedMemoryTest(); + //SharedMemoryTest(); } ////////////////////////////////////////////////////////////////// // On node barrier From fe5aaf7677c11e917af03bd0cbbc14a29b67bd2d Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Mon, 9 Aug 2021 04:06:30 -0700 Subject: [PATCH 17/31] Make comms benchmark same as Benchmark_comms_host_device --- benchmarks/Benchmark_ITT.cc | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/benchmarks/Benchmark_ITT.cc b/benchmarks/Benchmark_ITT.cc index 032535b3..7311dfc4 100644 --- a/benchmarks/Benchmark_ITT.cc +++ b/benchmarks/Benchmark_ITT.cc @@ -133,15 +133,15 @@ public: std::vector xbuf(8); std::vector rbuf(8); - Grid.ShmBufferFreeAll(); + //Grid.ShmBufferFreeAll(); + uint64_t bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD); for(int d=0;d<8;d++){ - xbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD)); - rbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD)); + xbuf[d] = (HalfSpinColourVectorD *)acceleratorAllocDevice(bytes); + rbuf[d] = (HalfSpinColourVectorD *)acceleratorAllocDevice(bytes); // bzero((void *)xbuf[d],lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD)); // bzero((void *)rbuf[d],lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD)); } - int bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD); int ncomm; double dbytes; std::vector times(Nloop); @@ -152,7 +152,7 @@ public: dbytes=0; ncomm=0; - thread_for(dir,8,{ + for(int dir=0;dir<8;dir++) { double tbytes; int mu =dir % 4; @@ -168,15 +168,16 @@ public: int comm_proc = mpi_layout[mu]-1; Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank); } - tbytes= Grid.StencilSendToRecvFrom((void *)&xbuf[dir][0], xmit_to_rank, - (void *)&rbuf[dir][0], recv_from_rank, - bytes,dir); + Grid.SendToRecvFrom((void *)&xbuf[dir][0], xmit_to_rank, + (void *)&rbuf[dir][0], recv_from_rank, + bytes); + tbytes = bytes; thread_critical { ncomm++; dbytes+=tbytes; } } - }); + }; Grid.Barrier(); double stop=usecond(); t_time[i] = stop-start; // microseconds @@ -196,8 +197,12 @@ public: << "\t\t"<< bidibytes/timestat.mean<< " " << bidibytes*timestat.err/(timestat.mean*timestat.mean) << " " << bidibytes/timestat.max << " " << bidibytes/timestat.min << std::endl; - } - } + for(int d=0;d<8;d++){ + acceleratorFreeDevice(xbuf[d]); + acceleratorFreeDevice(rbuf[d]); + } + } + } return; } @@ -281,7 +286,6 @@ public: uint64_t lmax=32; -#define NLOOP (1000*lmax*lmax*lmax*lmax/lat/lat/lat/lat) GridSerialRNG sRNG; sRNG.SeedFixedIntegers(std::vector({45,12,81,9})); for(int lat=8;lat<=lmax;lat+=8){ From 75030637cca1f60173ca760b09dd37852a88299d Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 10 Aug 2021 05:16:30 -0700 Subject: [PATCH 18/31] Improved comms benchmark, same as benchmark_comms_host_device --- benchmarks/Benchmark_ITT.cc | 73 ++++++++++++++++--------------------- 1 file changed, 31 insertions(+), 42 deletions(-) diff --git a/benchmarks/Benchmark_ITT.cc b/benchmarks/Benchmark_ITT.cc index 7311dfc4..81d1acd4 100644 --- a/benchmarks/Benchmark_ITT.cc +++ b/benchmarks/Benchmark_ITT.cc @@ -144,23 +144,19 @@ public: int ncomm; double dbytes; - std::vector times(Nloop); - for(int i=0;i1 ) { - dbytes=0; - ncomm=0; + std::vector times(Nloop); + for(int i=0;i1 ) { - + dbytes=0; + double start=usecond(); int xmit_to_rank; int recv_from_rank; + if ( dir == mu ) { int comm_proc=1; Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank); @@ -171,42 +167,35 @@ public: Grid.SendToRecvFrom((void *)&xbuf[dir][0], xmit_to_rank, (void *)&rbuf[dir][0], recv_from_rank, bytes); - tbytes = bytes; - thread_critical { - ncomm++; - dbytes+=tbytes; - } + dbytes+=bytes; + + double stop=usecond(); + t_time[i] = stop-start; // microseconds + } - }; - Grid.Barrier(); - double stop=usecond(); - t_time[i] = stop-start; // microseconds + timestat.statistics(t_time); + + dbytes=dbytes*ppn; + double xbytes = dbytes*0.5; + double bidibytes = dbytes; + + std::cout<1) ) { + if ( do_comms ) { std::cout< Date: Tue, 10 Aug 2021 05:35:15 -0700 Subject: [PATCH 19/31] Level 0 IPC set up --- Grid/communicator/SharedMemoryMPI.cc | 56 +++++++++++++++++++++++++--- Grid/threads/Accelerator.cc | 5 ++- Grid/threads/Accelerator.h | 7 ++++ 3 files changed, 62 insertions(+), 6 deletions(-) diff --git a/Grid/communicator/SharedMemoryMPI.cc b/Grid/communicator/SharedMemoryMPI.cc index caa03a60..554f338b 100644 --- a/Grid/communicator/SharedMemoryMPI.cc +++ b/Grid/communicator/SharedMemoryMPI.cc @@ -73,6 +73,7 @@ void GlobalSharedMemory::Init(Grid_MPI_Comm comm) WorldNodes = WorldSize/WorldShmSize; assert( (WorldNodes * WorldShmSize) == WorldSize ); + // FIXME: Check all WorldShmSize are the same ? ///////////////////////////////////////////////////////////////////// @@ -451,7 +452,8 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) //////////////////////////////////////////////////////////////////////////////////////////// #if defined(GRID_CUDA) ||defined(GRID_HIP) || defined(GRID_SYCL) -#if defined(GRID_SYCL) +//if defined(GRID_SYCL) +#if 0 void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) { void * ShmCommBuf ; @@ -488,7 +490,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) } #endif -#if defined(GRID_CUDA) ||defined(GRID_HIP) +#if defined(GRID_CUDA) ||defined(GRID_HIP) ||defined(GRID_SYCL) void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) { void * ShmCommBuf ; @@ -511,8 +513,16 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) /////////////////////////////////////////////////////////////////////////////////////////////////////////// // Each MPI rank should allocate our own buffer /////////////////////////////////////////////////////////////////////////////////////////////////////////// + auto zeDevice = cl::sycl::get_native(theGridAccelerator->get_device()); + auto zeContext= cl::sycl::get_native(theGridAccelerator->get_context()); +#ifdef GRID_SYCL_LEVEL_ZERO_IPC + ze_device_mem_alloc_desc_t zeDesc = {}; + zeMemAllocDevice(zeContext,&zeDesc,bytes,2*1024*1024,zeDevice,&ShmCommBuf); + std::cout << WorldRank << header " SharedMemoryMPI.cc zeMemAllocDevice "<< bytes + << "bytes at "<< std::hex<< ShmCommBuf < #include + +#define GRID_SYCL_LEVEL_ZERO_IPC + +#ifdef GRID_SYCL_LEVEL_ZERO_IPC +#include +#include +#endif NAMESPACE_BEGIN(Grid); extern cl::sycl::queue *theGridAccelerator; From 5d29e175d82edc92af07ca9b708ec7f187fc006d Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 10 Aug 2021 18:25:43 +0100 Subject: [PATCH 20/31] Typo fix --- Grid/communicator/SharedMemoryMPI.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Grid/communicator/SharedMemoryMPI.cc b/Grid/communicator/SharedMemoryMPI.cc index d230b2a5..11788744 100644 --- a/Grid/communicator/SharedMemoryMPI.cc +++ b/Grid/communicator/SharedMemoryMPI.cc @@ -513,9 +513,9 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags) /////////////////////////////////////////////////////////////////////////////////////////////////////////// // Each MPI rank should allocate our own buffer /////////////////////////////////////////////////////////////////////////////////////////////////////////// +#ifdef GRID_SYCL_LEVEL_ZERO_IPC auto zeDevice = cl::sycl::get_native(theGridAccelerator->get_device()); auto zeContext= cl::sycl::get_native(theGridAccelerator->get_context()); -#ifdef GRID_SYCL_LEVEL_ZERO_IPC ze_device_mem_alloc_desc_t zeDesc = {}; zeMemAllocDevice(zeContext,&zeDesc,bytes,2*1024*1024,zeDevice,&ShmCommBuf); std::cout << WorldRank << header " SharedMemoryMPI.cc zeMemAllocDevice "<< bytes From ffbdd91e0e2770f0e12f9296e5e9991e52799f4d Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 20 Aug 2021 01:15:00 +0100 Subject: [PATCH 21/31] Apple happiness --- Grid/util/Init.cc | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Grid/util/Init.cc b/Grid/util/Init.cc index bfbc464d..3b661524 100644 --- a/Grid/util/Init.cc +++ b/Grid/util/Init.cc @@ -52,10 +52,12 @@ Author: paboyle #include -#ifdef __APPLE__ static int +#ifdef __APPLE__ feenableexcept (unsigned int excepts) { +#if 0 + // Fails on Apple M1 static fenv_t fenv; unsigned int new_excepts = excepts & FE_ALL_EXCEPT; unsigned int old_excepts; // previous masks @@ -70,6 +72,7 @@ feenableexcept (unsigned int excepts) iold_excepts = (int) old_excepts; return ( fesetenv (&fenv) ? -1 : iold_excepts ); +#endif } #endif From 7163b31a26bba90241480fe444a3584cd66f5893 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 20 Aug 2021 01:15:23 +0100 Subject: [PATCH 22/31] Examples --- scripts/filelist | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/scripts/filelist b/scripts/filelist index 8e29ba88..c3f0cf3f 100755 --- a/scripts/filelist +++ b/scripts/filelist @@ -82,3 +82,18 @@ for f in $TESTS; do echo >> Make.inc done cd .. + +# examples Make.inc +cd $home/examples/ +echo> Make.inc +TESTS=`ls *.cc` +TESTLIST=`echo ${TESTS} | sed s/.cc//g ` +echo bin_PROGRAMS = ${TESTLIST} > Make.inc +echo >> Make.inc +for f in $TESTS; do + BNAME=`basename $f .cc` + echo ${BNAME}_SOURCES=$f >> Make.inc + echo ${BNAME}_LDADD='$(top_builddir)/Grid/libGrid.a'>> Make.inc + echo >> Make.inc +done +cd .. From 40098424c75cb693cf92ebafa758efac8b25053a Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 22 Aug 2021 14:17:12 +0100 Subject: [PATCH 23/31] Examples --- Makefile.am | 2 +- configure.ac | 1 + examples/Example_Laplacian.cc | 395 ++++++++++++++++++++++++++++++++++ examples/Makefile.am | 6 + 4 files changed, 403 insertions(+), 1 deletion(-) create mode 100644 examples/Example_Laplacian.cc create mode 100644 examples/Makefile.am diff --git a/Makefile.am b/Makefile.am index 33b25026..d2a1a326 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,5 +1,5 @@ # additional include paths necessary to compile the C++ library -SUBDIRS = Grid HMC benchmarks tests +SUBDIRS = Grid HMC benchmarks tests examples include $(top_srcdir)/doxygen.inc diff --git a/configure.ac b/configure.ac index 4e5e33c8..721d890e 100644 --- a/configure.ac +++ b/configure.ac @@ -815,6 +815,7 @@ AC_CONFIG_FILES(tests/smearing/Makefile) AC_CONFIG_FILES(tests/qdpxx/Makefile) AC_CONFIG_FILES(tests/testu01/Makefile) AC_CONFIG_FILES(benchmarks/Makefile) +AC_CONFIG_FILES(examples/Makefile) AC_OUTPUT echo "" diff --git a/examples/Example_Laplacian.cc b/examples/Example_Laplacian.cc new file mode 100644 index 00000000..f6bbf7cf --- /dev/null +++ b/examples/Example_Laplacian.cc @@ -0,0 +1,395 @@ +#include +using namespace Grid; + +/* +///////////////////////////////////////////////////////////////////////////////////////////// +// Grid/algorithms/SparseMatrix.h: Interface defining what I expect of a general sparse matrix, such as a Fermion action +///////////////////////////////////////////////////////////////////////////////////////////// +template class SparseMatrixBase { +public: + virtual GridBase *Grid(void) =0; + + virtual void M (const Field &in, Field &out)=0; + virtual void Mdag (const Field &in, Field &out)=0; + virtual void MdagM(const Field &in, Field &out) { + Field tmp (in.Grid()); + M(in,tmp); + Mdag(tmp,out); + } + virtual void Mdiag (const Field &in, Field &out)=0; + virtual void Mdir (const Field &in, Field &out,int dir, int disp)=0; + virtual void MdirAll (const Field &in, std::vector &out)=0; +}; +*/ + +const std::vector directions ({Xdir,Ydir,Zdir,Xdir,Ydir,Zdir}); +const std::vector displacements({1,1,1,-1,-1,-1}); + +template class FreeLaplacianCshift : public SparseMatrixBase +{ +public: + GridBase *grid; + FreeLaplacianCshift(GridBase *_grid) + { + grid=_grid; + }; + virtual GridBase *Grid(void) { return grid; }; + + virtual void M (const Field &in, Field &out) + { + out = Zero(); + for(int mu=0;mu &out) {assert(0);}; // Unimplemented need only for multigrid +}; + +template class CovariantLaplacianCshift : public SparseMatrixBase +{ +public: + INHERIT_GIMPL_TYPES(Gimpl); + + GridBase *grid; + GaugeField U; + + CovariantLaplacianCshift(GaugeField &_U) : + grid(_U.Grid()), + U(_U) { }; + + virtual GridBase *Grid(void) { return grid; }; + + virtual void M (const Field &in, Field &out) + { + out=Zero(); + for(int mu=0;mu(U, mu); // NB: Inefficent + out = out + Gimpl::CovShiftForward(Umu,mu,in); + out = out + Gimpl::CovShiftBackward(Umu,mu,in); + out = out - 2.0*in; + } + }; + virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian + virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid + virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid + virtual void MdirAll (const Field &in, std::vector &out) {assert(0);}; // Unimplemented need only for multigrid +}; + + +#define LEG_LOAD(Dir) \ + SE = st.GetEntry(ptype, Dir, ss); \ + if (SE->_is_local ) { \ + int perm= SE->_permute; \ + chi = coalescedReadPermute(in[SE->_offset],ptype,perm,lane); \ + } else { \ + chi = coalescedRead(buf[SE->_offset],lane); \ + } \ + acceleratorSynchronise(); + +template class FreeLaplacianStencil : public SparseMatrixBase +{ +public: + typedef typename Field::vector_object siteObject; + typedef CartesianStencil StencilImpl; + + GridBase *grid; + StencilImpl Stencil; + SimpleCompressor Compressor; + + FreeLaplacianStencil(GridBase *_grid) + : Stencil (_grid,6,Even,directions,displacements,0), grid(_grid) + { }; + + virtual GridBase *Grid(void) { return grid; }; + + virtual void M (const Field &_in, Field &_out) + { + + /////////////////////////////////////////////// + // Halo exchange for this geometry of stencil + /////////////////////////////////////////////// + Stencil.HaloExchange(_in, Compressor); + + /////////////////////////////////// + // Arithmetic expressions + /////////////////////////////////// + StencilEntry *SE; + + // Views; device friendly/accessible pointers + auto st = Stencil.View(AcceleratorRead); + auto buf = st.CommBuf(); + autoView( in , _in , AcceleratorRead); + autoView( out , _out , AcceleratorWrite); + + typedef typename Field::vector_object vobj; + typedef decltype(coalescedRead(in[0])) calcObj; + + const int Nsimd = vobj::Nsimd(); + const uint64_t NN = grid->oSites(); + const int lane=acceleratorSIMTlane(Nsimd); + + accelerator_for( ss, NN, Nsimd, { + + const int lane=acceleratorSIMTlane(Nsimd); + + calcObj chi; + calcObj res; + int ptype; + + res = coalescedRead(in[ss])*(-6.0); + LEG_LOAD(0); res = res + chi; + LEG_LOAD(1); res = res + chi; + LEG_LOAD(2); res = res + chi; + LEG_LOAD(3); res = res + chi; + LEG_LOAD(4); res = res + chi; + LEG_LOAD(5); res = res + chi; + + coalescedWrite(out[ss], res,lane); + + }); + }; + virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian + virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid + virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid + virtual void MdirAll (const Field &in, std::vector &out) {assert(0);}; // Unimplemented need only for multigrid +}; + +template class CovariantLaplacianStencil : public SparseMatrixBase +{ +public: + INHERIT_GIMPL_TYPES(Gimpl); + + typedef typename Field::vector_object siteObject; + + template using iImplDoubledGaugeField = iVector >, Nds>; + typedef iImplDoubledGaugeField SiteDoubledGaugeField; + typedef Lattice DoubledGaugeField; + + typedef CartesianStencil StencilImpl; + + GridBase *grid; + StencilImpl Stencil; + SimpleCompressor Compressor; + DoubledGaugeField Uds; + CovariantLaplacianStencil(GaugeField &Umu) + : + grid(Umu.Grid()), + Stencil (grid,6,Even,directions,displacements,0), + Uds(grid) + { + for (int mu = 0; mu < Nd; mu++) { + auto U = PeekIndex(Umu, mu); + PokeIndex(Uds, U, mu ); + U = adj(Cshift(U, mu, -1)); + PokeIndex(Uds, U, mu + 4); + } + }; + + virtual GridBase *Grid(void) { return grid; }; + + virtual void M (const Field &_in, Field &_out) + { + /////////////////////////////////////////////// + // Halo exchange for this geometry of stencil + /////////////////////////////////////////////// + Stencil.HaloExchange(_in, Compressor); + + /////////////////////////////////// + // Arithmetic expressions + /////////////////////////////////// + auto st = Stencil.View(AcceleratorRead); + auto buf = st.CommBuf(); + StencilEntry *SE; + + autoView( in , _in , AcceleratorRead); + autoView( out , _out , AcceleratorWrite); + autoView( U , Uds , AcceleratorRead); + + typedef typename Field::vector_object vobj; + typedef decltype(coalescedRead(in[0])) calcObj; + typedef decltype(coalescedRead(U[0](0))) calcLink; + + const int Nsimd = vobj::Nsimd(); + const uint64_t NN = grid->oSites(); + const int lane=acceleratorSIMTlane(Nsimd); + accelerator_for( ss, NN, Nsimd, { + + const int lane=acceleratorSIMTlane(Nsimd); + + calcObj chi; + calcObj res; + calcObj Uchi; + calcLink UU; + int ptype; + + res = coalescedRead(in[ss])*(-6.0); + +#define LEG_LOAD_MULT(leg,polarisation) \ + UU = coalescedRead(U[ss](polarisation)); \ + LEG_LOAD(leg); \ + mult(&Uchi(), &UU, &chi()); \ + res = res + Uchi; + + LEG_LOAD_MULT(0,Xp); + LEG_LOAD_MULT(1,Yp); + LEG_LOAD_MULT(2,Zp); + LEG_LOAD_MULT(3,Xm); + LEG_LOAD_MULT(4,Ym); + LEG_LOAD_MULT(5,Zm); + + coalescedWrite(out[ss], res,lane); + }); + }; + virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian + virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid + virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid + virtual void MdirAll (const Field &in, std::vector &out) {assert(0);}; // Unimplemented need only for multigrid +}; + +#undef LEG_LOAD_MULT +#undef LEG_LOAD + +int main(int argc, char ** argv) +{ + Grid_init(&argc, &argv); + + typedef LatticeColourVector Field; + + auto latt_size = GridDefaultLatt(); + auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + auto mpi_layout = GridDefaultMpi(); + + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridParallelRNG RNG(&Grid); RNG.SeedFixedIntegers(std::vector({45,12,81,9})); + + FreeLaplacianCshift FLcs(&Grid); + FreeLaplacianStencil FLst(&Grid); + + LatticeGaugeField U(&Grid); + + SU::ColdConfiguration(RNG,U); + + std::cout << " Gauge field has norm " < CLcs(U); + CovariantLaplacianStencil CLst(U); + + Field in(&Grid); gaussian(RNG,in); + Field out_FLcs(&Grid); + Field out_FLst(&Grid); + Field out_CLcs(&Grid); + Field out_CLst(&Grid); + Field diff(&Grid); + + //////////////////////////////////////////////////////// + // First test: in free field these should all agree + //////////////////////////////////////////////////////// + FLcs.M(in,out_FLcs); + FLst.M(in,out_FLst); + CLcs.M(in,out_CLcs); + CLst.M(in,out_CLst); + + std:: cout << "******************************************************************" <::RandomGaugeTransform(RNG,U_GT,g); // Unit gauge + + Field in_GT(&Grid); + Field out_GT(&Grid); + + Field out_CLcs_GT(&Grid); + Field out_CLst_GT(&Grid); + + CovariantLaplacianCshift CLcs_GT(U_GT); + CovariantLaplacianStencil CLst_GT(U_GT); + + in_GT = g*in; + out_GT = g*out_FLcs; + + // Check M^GT_xy in_GT = g(x) M_xy g^dag(y) g(y) in = g(x) out(x) + CLcs_GT.M(in_GT,out_CLcs_GT); + CLst_GT.M(in_GT,out_CLst_GT); + + diff = out_CLcs_GT - out_GT; + std:: cout << " Difference between Gauge xformed result and covariant Cshift Laplacian in xformed gauge = " < dim_mask({1,1,1,0}); // 3d FFT + FFT theFFT(&Grid); + Field out(&Grid); + Field F_out(&Grid); + Field F_in(&Grid); + + // FFT the random input vector + theFFT.FFT_dim_mask(F_in,in,dim_mask,FFT::forward); + + // Convolution theorem: multiply by Fourier representation of (discrete) Laplacian to apply diff op + LatticeComplexD lap(&Grid); lap = Zero(); + LatticeComplexD kmu(&Grid); + ComplexD ci(0.0,1.0); + for(int mu=0;mu<3;mu++) { + + RealD TwoPiL = M_PI * 2.0/ latt_size[mu]; + + LatticeCoordinate(kmu,mu); + kmu = TwoPiL * kmu; + + // (e^ik_mu + e^-ik_mu - 2) = 2( cos kmu - 1) ~ 2 (1 - k_mu^2/2 -1 ) = - k_mu^2 + O(k^4) + lap = lap + 2.0*cos(kmu) - 2.0; + + } + F_out = lap * F_in; + + // Inverse FFT the result + theFFT.FFT_dim_mask(out,F_out,dim_mask,FFT::backward); + + std::cout<<"Fourier xformed (in) "< Date: Sun, 22 Aug 2021 18:25:07 +0100 Subject: [PATCH 24/31] Extra examples / solutions --- examples/Example_Laplacian_inverter.cc | 122 ++++++++++++++++++++++++ examples/Example_Laplacian_smearing.cc | 127 +++++++++++++++++++++++++ examples/Example_Laplacian_solver.cc | 127 +++++++++++++++++++++++++ 3 files changed, 376 insertions(+) create mode 100644 examples/Example_Laplacian_inverter.cc create mode 100644 examples/Example_Laplacian_smearing.cc create mode 100644 examples/Example_Laplacian_solver.cc diff --git a/examples/Example_Laplacian_inverter.cc b/examples/Example_Laplacian_inverter.cc new file mode 100644 index 00000000..1235d2b8 --- /dev/null +++ b/examples/Example_Laplacian_inverter.cc @@ -0,0 +1,122 @@ +#include +using namespace Grid; + +// Function used for Chebyshev smearing +// +Real MomentumSmearing(Real p2) +{ + return (1 - 4.0*p2) * exp(-p2/4); +} +Real DistillationSmearing(Real p2) +{ + if ( p2 > 0.5 ) return 0.0; + else return 1.0; +} + +// Flip sign to make prop to p^2, not -p^2 relative to last example +template class CovariantLaplacianCshift : public SparseMatrixBase +{ +public: + INHERIT_GIMPL_TYPES(Gimpl); + + GridBase *grid; + GaugeField U; + + CovariantLaplacianCshift(GaugeField &_U) : + grid(_U.Grid()), + U(_U) { }; + + virtual GridBase *Grid(void) { return grid; }; + + virtual void M (const Field &in, Field &out) + { + out=Zero(); + for(int mu=0;mu(U, mu); // NB: Inefficent + out = out - Gimpl::CovShiftForward(Umu,mu,in); + out = out - Gimpl::CovShiftBackward(Umu,mu,in); + out = out + 2.0*in; + } + }; + virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian + virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid + virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid + virtual void MdirAll (const Field &in, std::vector &out) {assert(0);}; // Unimplemented need only for multigrid +}; + + + +int main(int argc, char ** argv) +{ + Grid_init(&argc, &argv); + + typedef LatticeColourVector Field; + + auto latt_size = GridDefaultLatt(); + auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + auto mpi_layout = GridDefaultMpi(); + + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridParallelRNG RNG(&Grid); RNG.SeedFixedIntegers(std::vector({45,12,81,9})); + + + LatticeGaugeField U(&Grid); + + SU::ColdConfiguration(RNG,U); + + typedef CovariantLaplacianCshift Laplacian_t; + Laplacian_t Laplacian(U); + + ColourVector ColourKronecker; + ColourKronecker = Zero(); + ColourKronecker()()(0) = 1.0; + + Coordinate site({latt_size[0]/2, + latt_size[1]/2, + latt_size[2]/2, + 0}); + + Field kronecker(&Grid); + kronecker = Zero(); + pokeSite(ColourKronecker,kronecker,site); + + + Field psi(&Grid), chi(&Grid); + + ////////////////////////////////////// + // Classic Wuppertal smearing + ////////////////////////////////////// + + Integer Iterations = 80; + Real width = 2.0; + Real coeff = (width*width) / Real(4*Iterations); + + chi=kronecker; + // chi = (1-p^2/2N)^N kronecker + for(int n = 0; n < Iterations; ++n) { + Laplacian.M(chi,psi); + chi = chi - coeff*psi; + } + + std::cout << " Wuppertal smeared operator is chi = \n" << chi < HermOp(Laplacian); + + Chebyshev ChebySmear(lo,hi,20,DistillationSmearing); + // Chebyshev ChebySmear(lo,hi,20,MomentumSmearing); + { + std::ofstream of("chebysmear"); + ChebySmear.csv(of); + } + + ChebySmear(HermOp,kronecker,chi); + + std::cout << " Chebyshev smeared operator is chi = \n" << chi < +using namespace Grid; + +// Function used for Chebyshev smearing +// +Real MomentumSmearing(Real p2) +{ + return (1 - 4.0*p2) * exp(-p2/4); +} +Real DistillationSmearing(Real p2) +{ + if ( p2 > 0.5 ) return 0.0; + else return 1.0; +} + +// Flip sign to make prop to p^2, not -p^2 relative to last example +template class CovariantLaplacianCshift : public SparseMatrixBase +{ +public: + INHERIT_GIMPL_TYPES(Gimpl); + + GridBase *grid; + GaugeField U; + + CovariantLaplacianCshift(GaugeField &_U) : + grid(_U.Grid()), + U(_U) { }; + + virtual GridBase *Grid(void) { return grid; }; + + virtual void M (const Field &in, Field &out) + { + out=Zero(); + for(int mu=0;mu(U, mu); // NB: Inefficent + out = out - Gimpl::CovShiftForward(Umu,mu,in); + out = out - Gimpl::CovShiftBackward(Umu,mu,in); + out = out + 2.0*in; + } + }; + virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian + virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid + virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid + virtual void MdirAll (const Field &in, std::vector &out) {assert(0);}; // Unimplemented need only for multigrid +}; + + + +int main(int argc, char ** argv) +{ + Grid_init(&argc, &argv); + + typedef LatticeColourVector Field; + + auto latt_size = GridDefaultLatt(); + auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + auto mpi_layout = GridDefaultMpi(); + + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridParallelRNG RNG(&Grid); RNG.SeedFixedIntegers(std::vector({45,12,81,9})); + + + LatticeGaugeField U(&Grid); + + SU::ColdConfiguration(RNG,U); + + typedef CovariantLaplacianCshift Laplacian_t; + Laplacian_t Laplacian(U); + + + ColourVector ColourKronecker; + ColourKronecker = Zero(); + ColourKronecker()()(0) = 1.0; + + Coordinate site({latt_size[0]/2, + latt_size[1]/2, + latt_size[2]/2, + 0}); + + Field kronecker(&Grid); + kronecker = Zero(); + pokeSite(ColourKronecker,kronecker,site); + + + Field psi(&Grid), chi(&Grid); + + ////////////////////////////////////// + // Classic Wuppertal smearing + ////////////////////////////////////// + + Integer Iterations = 80; + Real width = 2.0; + Real coeff = (width*width) / Real(4*Iterations); + + chi=kronecker; + // chi = (1-p^2/2N)^N kronecker + for(int n = 0; n < Iterations; ++n) { + Laplacian.M(chi,psi); + chi = chi - coeff*psi; + } + + std::cout << " Wuppertal smeared operator is chi = \n" << chi < HermOp(Laplacian); + + std::cout << " Checking spectral range of our POSITIVE definite operator \n"; + PowerMethod PM; + PM(HermOp,kronecker); + + Chebyshev ChebySmear(lo,hi,20,DistillationSmearing); + // Chebyshev ChebySmear(lo,hi,20,MomentumSmearing); + { + std::ofstream of("chebysmear"); + ChebySmear.csv(of); + } + + ChebySmear(HermOp,kronecker,chi); + + std::cout << " Chebyshev smeared operator is chi = \n" << chi < +using namespace Grid; + +template +void SimpleConjugateGradient(LinearOperatorBase &HPDop,const Field &b, Field &x) +{ + RealD cp, c, alpha, d, beta, ssq, qq; + RealD Tolerance=1.0e-10; + int MaxIterations=10000; + + Field p(b), mmp(b), r(b); + + HPDop.HermOpAndNorm(x, mmp, d, beta); + + r = b - mmp; + p = r; + + alpha = norm2(p); + cp = alpha; + ssq = norm2(b); + + RealD rsq = Tolerance * Tolerance * ssq; + + for (int k = 1; k <= MaxIterations; k++) { + c = cp; + + HPDop.HermOp(p, mmp); + + ComplexD dc = innerProduct(p,mmp); + d = dc.real(); + alpha = c / d; + + cp = axpy_norm(r, -alpha, mmp, r); + beta = cp / c; + + x = x + alpha* p ; + p = r + beta* p ; + + std::cout << "iteration "< class CovariantLaplacianCshift : public SparseMatrixBase +{ +public: + INHERIT_GIMPL_TYPES(Gimpl); + + GridBase *grid; + GaugeField U; + RealD m2=1.0e-2; + CovariantLaplacianCshift(GaugeField &_U) : + grid(_U.Grid()), + U(_U) { }; + + virtual GridBase *Grid(void) { return grid; }; + + virtual void M (const Field &in, Field &out) + { + out=Zero(); + for(int mu=0;mu(U, mu); // NB: Inefficent + out = out - Gimpl::CovShiftForward(Umu,mu,in); + out = out - Gimpl::CovShiftBackward(Umu,mu,in); + out = out + 2.0*in + m2*in; + } + }; + virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian + virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid + virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid + virtual void MdirAll (const Field &in, std::vector &out) {assert(0);}; // Unimplemented need only for multigrid +}; + + + +int main(int argc, char ** argv) +{ + Grid_init(&argc, &argv); + + typedef LatticeColourVector Field; + + auto latt_size = GridDefaultLatt(); + auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + auto mpi_layout = GridDefaultMpi(); + + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridParallelRNG RNG(&Grid); RNG.SeedFixedIntegers(std::vector({45,12,81,9})); + + + LatticeGaugeField U(&Grid); + + SU::ColdConfiguration(RNG,U); + + typedef CovariantLaplacianCshift Laplacian_t; + Laplacian_t Laplacian(U); + + + ColourVector ColourKronecker; + ColourKronecker = Zero(); + ColourKronecker()()(0) = 1.0; + + Coordinate site({0,0,0,0}); // Point source at origin + + Field kronecker(&Grid); + kronecker = Zero(); + pokeSite(ColourKronecker,kronecker,site); + + Field psi(&Grid); psi=Zero(); + + HermitianLinearOperator HermOp(Laplacian); + SimpleConjugateGradient(HermOp, kronecker,psi); + + Field r(&Grid); + Laplacian.M(psi,r); + r=kronecker-r; + + std::cout << "True residual "<< norm2(r) < Date: Sun, 22 Aug 2021 18:28:39 +0100 Subject: [PATCH 25/31] Remove the file --- examples/Example_Laplacian_inverter.cc | 122 ------------------------- 1 file changed, 122 deletions(-) delete mode 100644 examples/Example_Laplacian_inverter.cc diff --git a/examples/Example_Laplacian_inverter.cc b/examples/Example_Laplacian_inverter.cc deleted file mode 100644 index 1235d2b8..00000000 --- a/examples/Example_Laplacian_inverter.cc +++ /dev/null @@ -1,122 +0,0 @@ -#include -using namespace Grid; - -// Function used for Chebyshev smearing -// -Real MomentumSmearing(Real p2) -{ - return (1 - 4.0*p2) * exp(-p2/4); -} -Real DistillationSmearing(Real p2) -{ - if ( p2 > 0.5 ) return 0.0; - else return 1.0; -} - -// Flip sign to make prop to p^2, not -p^2 relative to last example -template class CovariantLaplacianCshift : public SparseMatrixBase -{ -public: - INHERIT_GIMPL_TYPES(Gimpl); - - GridBase *grid; - GaugeField U; - - CovariantLaplacianCshift(GaugeField &_U) : - grid(_U.Grid()), - U(_U) { }; - - virtual GridBase *Grid(void) { return grid; }; - - virtual void M (const Field &in, Field &out) - { - out=Zero(); - for(int mu=0;mu(U, mu); // NB: Inefficent - out = out - Gimpl::CovShiftForward(Umu,mu,in); - out = out - Gimpl::CovShiftBackward(Umu,mu,in); - out = out + 2.0*in; - } - }; - virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian - virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid - virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid - virtual void MdirAll (const Field &in, std::vector &out) {assert(0);}; // Unimplemented need only for multigrid -}; - - - -int main(int argc, char ** argv) -{ - Grid_init(&argc, &argv); - - typedef LatticeColourVector Field; - - auto latt_size = GridDefaultLatt(); - auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); - auto mpi_layout = GridDefaultMpi(); - - GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridParallelRNG RNG(&Grid); RNG.SeedFixedIntegers(std::vector({45,12,81,9})); - - - LatticeGaugeField U(&Grid); - - SU::ColdConfiguration(RNG,U); - - typedef CovariantLaplacianCshift Laplacian_t; - Laplacian_t Laplacian(U); - - ColourVector ColourKronecker; - ColourKronecker = Zero(); - ColourKronecker()()(0) = 1.0; - - Coordinate site({latt_size[0]/2, - latt_size[1]/2, - latt_size[2]/2, - 0}); - - Field kronecker(&Grid); - kronecker = Zero(); - pokeSite(ColourKronecker,kronecker,site); - - - Field psi(&Grid), chi(&Grid); - - ////////////////////////////////////// - // Classic Wuppertal smearing - ////////////////////////////////////// - - Integer Iterations = 80; - Real width = 2.0; - Real coeff = (width*width) / Real(4*Iterations); - - chi=kronecker; - // chi = (1-p^2/2N)^N kronecker - for(int n = 0; n < Iterations; ++n) { - Laplacian.M(chi,psi); - chi = chi - coeff*psi; - } - - std::cout << " Wuppertal smeared operator is chi = \n" << chi < HermOp(Laplacian); - - Chebyshev ChebySmear(lo,hi,20,DistillationSmearing); - // Chebyshev ChebySmear(lo,hi,20,MomentumSmearing); - { - std::ofstream of("chebysmear"); - ChebySmear.csv(of); - } - - ChebySmear(HermOp,kronecker,chi); - - std::cout << " Chebyshev smeared operator is chi = \n" << chi < Date: Sun, 22 Aug 2021 18:40:55 +0100 Subject: [PATCH 26/31] Fail on non-apple --- Grid/util/Init.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Grid/util/Init.cc b/Grid/util/Init.cc index 3b661524..c0bfc906 100644 --- a/Grid/util/Init.cc +++ b/Grid/util/Init.cc @@ -52,8 +52,8 @@ Author: paboyle #include -static int #ifdef __APPLE__ +static int feenableexcept (unsigned int excepts) { #if 0 From 5b3c530aa778b23ae9383a5b7a63fd20f5d9c2cb Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Mon, 23 Aug 2021 15:30:45 +0100 Subject: [PATCH 27/31] Return value --- Grid/util/Init.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/Grid/util/Init.cc b/Grid/util/Init.cc index c0bfc906..ab2d2399 100644 --- a/Grid/util/Init.cc +++ b/Grid/util/Init.cc @@ -73,6 +73,7 @@ feenableexcept (unsigned int excepts) iold_excepts = (int) old_excepts; return ( fesetenv (&fenv) ? -1 : iold_excepts ); #endif + return 0; } #endif From 0d588b95f4e1847024c4eea0094f7fc0dc36ad28 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Mon, 23 Aug 2021 23:14:26 +0100 Subject: [PATCH 28/31] Bug fix to Example_Laplacian test --- Grid/threads/Accelerator.h | 6 ++++++ examples/Example_Laplacian.cc | 9 +++++---- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/Grid/threads/Accelerator.h b/Grid/threads/Accelerator.h index 2c9d15ba..a8c91aa8 100644 --- a/Grid/threads/Accelerator.h +++ b/Grid/threads/Accelerator.h @@ -39,6 +39,10 @@ Author: paboyle #ifdef HAVE_MM_MALLOC_H #include #endif +#ifdef __APPLE__ +// no memalign +inline void *memalign(size_t align, size_t bytes) { return malloc(bytes); } +#endif NAMESPACE_BEGIN(Grid); @@ -419,6 +423,8 @@ inline void acceleratorMemSet(void *base,int value,size_t bytes) { hipMemset(bas #undef GRID_SIMT + + #define accelerator #define accelerator_inline strong_inline #define accelerator_for(iterator,num,nsimd, ... ) thread_for(iterator, num, { __VA_ARGS__ }); diff --git a/examples/Example_Laplacian.cc b/examples/Example_Laplacian.cc index f6bbf7cf..fa8466cf 100644 --- a/examples/Example_Laplacian.cc +++ b/examples/Example_Laplacian.cc @@ -116,7 +116,6 @@ public: /////////////////////////////////// // Arithmetic expressions /////////////////////////////////// - StencilEntry *SE; // Views; device friendly/accessible pointers auto st = Stencil.View(AcceleratorRead); @@ -129,10 +128,11 @@ public: const int Nsimd = vobj::Nsimd(); const uint64_t NN = grid->oSites(); - const int lane=acceleratorSIMTlane(Nsimd); accelerator_for( ss, NN, Nsimd, { + StencilEntry *SE; + const int lane=acceleratorSIMTlane(Nsimd); calcObj chi; @@ -202,7 +202,6 @@ public: /////////////////////////////////// auto st = Stencil.View(AcceleratorRead); auto buf = st.CommBuf(); - StencilEntry *SE; autoView( in , _in , AcceleratorRead); autoView( out , _out , AcceleratorWrite); @@ -214,9 +213,11 @@ public: const int Nsimd = vobj::Nsimd(); const uint64_t NN = grid->oSites(); - const int lane=acceleratorSIMTlane(Nsimd); + accelerator_for( ss, NN, Nsimd, { + StencilEntry *SE; + const int lane=acceleratorSIMTlane(Nsimd); calcObj chi; From 114920b8de4048b250206b4bd3ac469cd7fffd95 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 25 Aug 2021 12:24:17 +0100 Subject: [PATCH 29/31] Some example clean up --- examples/Example_Laplacian_smearing.cc | 4 ++-- examples/Example_Laplacian_solver.cc | 14 ++++++++------ 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/examples/Example_Laplacian_smearing.cc b/examples/Example_Laplacian_smearing.cc index 327c5ca9..9780e8a0 100644 --- a/examples/Example_Laplacian_smearing.cc +++ b/examples/Example_Laplacian_smearing.cc @@ -112,8 +112,8 @@ int main(int argc, char ** argv) PowerMethod PM; PM(HermOp,kronecker); - Chebyshev ChebySmear(lo,hi,20,DistillationSmearing); - // Chebyshev ChebySmear(lo,hi,20,MomentumSmearing); + // Chebyshev ChebySmear(lo,hi,20,DistillationSmearing); + Chebyshev ChebySmear(lo,hi,20,MomentumSmearing); { std::ofstream of("chebysmear"); ChebySmear.csv(of); diff --git a/examples/Example_Laplacian_solver.cc b/examples/Example_Laplacian_solver.cc index 88275df8..4dc00280 100644 --- a/examples/Example_Laplacian_solver.cc +++ b/examples/Example_Laplacian_solver.cc @@ -15,8 +15,7 @@ void SimpleConjugateGradient(LinearOperatorBase &HPDop,const Field &b, Fi r = b - mmp; p = r; - alpha = norm2(p); - cp = alpha; + cp = alpha = norm2(p); ssq = norm2(b); RealD rsq = Tolerance * Tolerance * ssq; @@ -26,11 +25,12 @@ void SimpleConjugateGradient(LinearOperatorBase &HPDop,const Field &b, Fi HPDop.HermOp(p, mmp); - ComplexD dc = innerProduct(p,mmp); - d = dc.real(); + d = real(innerProduct(p,mmp)); + alpha = c / d; - cp = axpy_norm(r, -alpha, mmp, r); + r = r - alpha *mmp; + cp = norm2(r); beta = cp / c; x = x + alpha* p ; @@ -121,7 +121,9 @@ int main(int argc, char ** argv) r=kronecker-r; std::cout << "True residual "<< norm2(r) < Date: Mon, 30 Aug 2021 20:32:11 -0400 Subject: [PATCH 30/31] Some sample code --- examples/Example_Mobius_spectrum.cc | 334 ++++++++++++++++++++++++++++ 1 file changed, 334 insertions(+) create mode 100644 examples/Example_Mobius_spectrum.cc diff --git a/examples/Example_Mobius_spectrum.cc b/examples/Example_Mobius_spectrum.cc new file mode 100644 index 00000000..dd84a336 --- /dev/null +++ b/examples/Example_Mobius_spectrum.cc @@ -0,0 +1,334 @@ +/************************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Copyright (C) 2021 + +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 + +using namespace std; +using namespace Grid; + +// Flip sign to make prop to p^2, not -p^2 relative to last example +template class CovariantLaplacianCshift : public SparseMatrixBase +{ +public: + INHERIT_GIMPL_TYPES(Gimpl); + + GridBase *grid; + GaugeField U; + + CovariantLaplacianCshift(GaugeField &_U) : + grid(_U.Grid()), + U(_U) { }; + + virtual GridBase *Grid(void) { return grid; }; + + virtual void M (const Field &in, Field &out) + { + out=Zero(); + for(int mu=0;mu(U, mu); // NB: Inefficent + out = out - Gimpl::CovShiftForward(Umu,mu,in); + out = out - Gimpl::CovShiftBackward(Umu,mu,in); + out = out + 2.0*in; + } + }; + virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian + virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid + virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid + virtual void MdirAll (const Field &in, std::vector &out) {assert(0);}; // Unimplemented need only for multigrid +}; + +void MakePhase(Coordinate mom,LatticeComplex &phase) +{ + GridBase *grid = phase.Grid(); + auto latt_size = grid->GlobalDimensions(); + ComplexD ci(0.0,1.0); + phase=Zero(); + + LatticeComplex coor(phase.Grid()); + for(int mu=0;mu +void GaussianSmear(LatticeGaugeField &U,Field &unsmeared,Field &smeared) +{ + typedef CovariantLaplacianCshift Laplacian_t; + Laplacian_t Laplacian(U); + + Integer Iterations = 40; + Real width = 2.0; + Real coeff = (width*width) / Real(4*Iterations); + + Field tmp(U.Grid()); + smeared=unsmeared; + // chi = (1-p^2/2N)^N kronecker + for(int n = 0; n < Iterations; ++n) { + Laplacian.M(smeared,tmp); + smeared = smeared - coeff*tmp; + std::cout << " smear iter " << n<<" " < +void Solve(Action &D,LatticePropagator &source,LatticePropagator &propagator) +{ + GridBase *UGrid = D.GaugeGrid(); + GridBase *FGrid = D.FermionGrid(); + + LatticeFermion src4 (UGrid); + LatticeFermion src5 (FGrid); + LatticeFermion result5(FGrid); + LatticeFermion result4(UGrid); + + ConjugateGradient CG(1.0e-8,100000); + SchurRedBlackDiagMooeeSolve schur(CG); + ZeroGuesser ZG; // Could be a DeflatedGuesser if have eigenvectors + for(int s=0;s(src4,source,s,c); + + D.ImportPhysicalFermionSource(src4,src5); + + result5=Zero(); + schur(D,src5,result5,ZG); + std::cout<(propagator,result4,s,c); + } + } +} + +class MesonFile: Serializable { +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(MesonFile, std::vector >, data); +}; + +void MesonTrace(std::string file,LatticePropagator &q1,LatticePropagator &q2,LatticeComplex &phase) +{ + const int nchannel=4; + Gamma::Algebra Gammas[nchannel][2] = { + {Gamma::Algebra::Gamma5 ,Gamma::Algebra::Gamma5}, + {Gamma::Algebra::GammaTGamma5,Gamma::Algebra::GammaTGamma5}, + {Gamma::Algebra::GammaTGamma5,Gamma::Algebra::Gamma5}, + {Gamma::Algebra::Gamma5 ,Gamma::Algebra::GammaTGamma5} + }; + + Gamma G5(Gamma::Algebra::Gamma5); + + LatticeComplex meson_CF(q1.Grid()); + MesonFile MF; + + for(int ch=0;ch meson_T; + sliceSum(meson_CF,meson_T, Tdir); + + int nt=meson_T.size(); + + std::vector corr(nt); + for(int t=0;t seeds4({1,2,3,4}); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); + std::string config; + if( argc > 1 && argv[1][0] != '-' ) + { + std::cout<::ColdConfiguration(Umu); + // SU::HotConfiguration(RNG4,Umu); + config="HotConfig"; + } + + std::vector masses({ 0.03,0.04,0.45} ); // u/d, s, c ?? + + int nmass = masses.size(); + + std::vector FermActs; + + std::cout< PointProps(nmass,UGrid); + std::vector GaussProps(nmass,UGrid); + std::vector Z2Props (nmass,UGrid); + + for(int m=0;m Date: Mon, 30 Aug 2021 21:15:39 -0400 Subject: [PATCH 31/31] Comment update --- examples/Example_Mobius_spectrum.cc | 28 ++++------------------------ 1 file changed, 4 insertions(+), 24 deletions(-) diff --git a/examples/Example_Mobius_spectrum.cc b/examples/Example_Mobius_spectrum.cc index dd84a336..f4cd3335 100644 --- a/examples/Example_Mobius_spectrum.cc +++ b/examples/Example_Mobius_spectrum.cc @@ -1,33 +1,13 @@ -/************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid +/* + * Warning: This code illustrative only: not well tested, and not meant for production use + * without regression / tests being applied + */ - Copyright (C) 2021 - -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 using namespace std; using namespace Grid; -// Flip sign to make prop to p^2, not -p^2 relative to last example template class CovariantLaplacianCshift : public SparseMatrixBase { public: