From 026e736dfabdb256056a5fc6dc31d443556a5bfc Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Mon, 3 Apr 2023 16:31:19 +0100 Subject: [PATCH 01/20] Projection on algebra can now be templated. Fix #12 --- Grid/lattice/Lattice_ET.h | 2 + Grid/qcd/utils/Sp2n.impl | 24 +++- Grid/tensors/Tensor_Ta.h | 80 ++++++++++++ tests/sp2n/Test_project_on_Sp.cc | 211 +++++++++++++++++++++++-------- 4 files changed, 261 insertions(+), 56 deletions(-) diff --git a/Grid/lattice/Lattice_ET.h b/Grid/lattice/Lattice_ET.h index 9042eb8f..661d7441 100644 --- a/Grid/lattice/Lattice_ET.h +++ b/Grid/lattice/Lattice_ET.h @@ -345,6 +345,7 @@ GridUnopClass(UnaryNot, Not(a)); GridUnopClass(UnaryTrace, trace(a)); GridUnopClass(UnaryTranspose, transpose(a)); GridUnopClass(UnaryTa, Ta(a)); +GridUnopClass(UnaryProjectSp2nAlgebra, ProjectSp2nAlgebra(a)); GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a)); GridUnopClass(UnaryProjectOnSpGroup, ProjectOnSpGroup(a)); GridUnopClass(UnaryTimesI, timesI(a)); @@ -457,6 +458,7 @@ GRID_DEF_UNOP(operator!, UnaryNot); GRID_DEF_UNOP(trace, UnaryTrace); GRID_DEF_UNOP(transpose, UnaryTranspose); GRID_DEF_UNOP(Ta, UnaryTa); +GRID_DEF_UNOP(ProjectSp2nAlgebra, UnaryProjectSp2nAlgebra); GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup); GRID_DEF_UNOP(ProjectOnSpGroup, UnaryProjectOnSpGroup); GRID_DEF_UNOP(timesI, UnaryTimesI); diff --git a/Grid/qcd/utils/Sp2n.impl b/Grid/qcd/utils/Sp2n.impl index 77da9ad0..347eb3ce 100644 --- a/Grid/qcd/utils/Sp2n.impl +++ b/Grid/qcd/utils/Sp2n.impl @@ -274,8 +274,6 @@ static void OmegaInvariance(ColourMatrix &in) { Omega = Zero(); const int nsp=ncolour/2; - std::cout << GridLogMessage << "I am a ColourMatrix" << std::endl; - // for (int i = 0; i < ncolour; i++) wrong?! //{ // Omega()()(i, 2*ncolour-1-i) = 1.; @@ -314,8 +312,6 @@ static void OmegaInvariance(GaugeField &in) { Omega = Zero(); identity = 1.; - std::cout << GridLogMessage << "I am a GaugeField " << std::endl; - U = PeekIndex(in, 1); OmegaInvariance(U); @@ -333,8 +329,6 @@ static void OmegaInvariance(LatticeColourMatrixD &in) { Omega = Zero(); identity = 1.; - std::cout << GridLogMessage << "I am a LatticeColourMatrix " << std::endl; - for (int i = 0; i < nsp; i++) { Omega()()(i, nsp + i) = 1.; Omega()()(nsp + i, i) = -1; @@ -356,3 +350,21 @@ static void OmegaInvariance(LatticeColourMatrixD &in) { } } +template +static void Omega(LatticeColourMatrixD &in) { + const int nsp=ncolour/2; + LatticeColourMatrixD OmegaLatt(in.Grid()); + LatticeColourMatrixD identity(in.Grid()); + ColourMatrix Omega; + + OmegaLatt = Zero(); + Omega = Zero(); + identity = 1.; + + for (int i = 0; i < nsp; i++) { + Omega()()(i, nsp + i) = 1.; + Omega()()(nsp + i, i) = -1; + } + OmegaLatt = OmegaLatt + (identity * Omega); + in = OmegaLatt; +} diff --git a/Grid/tensors/Tensor_Ta.h b/Grid/tensors/Tensor_Ta.h index 2ee51fb4..87169ffc 100644 --- a/Grid/tensors/Tensor_Ta.h +++ b/Grid/tensors/Tensor_Ta.h @@ -66,6 +66,86 @@ template accelerator_inline iMatrix Ta(const iMatrix return ret; } +// for sp2n can't do something as simple as Ta. We do a Gram-Schmidt + +template accelerator_inline iScalar ProjectSp2nAlgebra(const iScalar&r) +{ + iScalar ret; + ret._internal = ProjectSp2nAlgebra(r._internal); + return ret; +} +template accelerator_inline iVector ProjectSp2nAlgebra(const iVector&r) +{ + iVector ret; + for(int i=0;i accelerator_inline iMatrix ProjectSp2nAlgebra(const iMatrix &arg) +{ + iMatrix ret; + vtype nrm; + vtype inner; + vtype tmp; + + for(int c1=0;c1 + prn += conjugate(ret._internal[c1][c])*ret._internal[b+N/2][c]; // + } + + + for(int c=0; c U_b + U_{b+N} ) + } + } + + zeroit(inner); + for(int c2=0;c2gSites(); std::vector pseeds({1,2,3,4,5}); - std::vector sseeds({6,7,8,9,10}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(pseeds); - GridSerialRNG sRNG; sRNG.SeedFixedIntegers(sseeds); SU::HotConfiguration(pRNG,Umu); - U = PeekIndex(Umu,2); + U = PeekIndex(Umu,0); aux = U*adj(U) - identity; - std::cout <::OmegaInvariance(U); + std::cout << GridLogMessage << "Checking Omega invariance after ProjectOnSpGroup" << std::endl; + Sp::OmegaInvariance(U); // no assertion here, but the next check will kill us if we are not simplectic // checks on elements - std::cout <::ProjectOnGaugeGroup(U); + aux = U*adj(U) - identity; + std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl; + assert( norm2(aux) < 1e-8); + Sp::OmegaInvariance(U); + + std::cout << GridLogMessage << "Checking the structure is " << std::endl; + std::cout << GridLogMessage << "U = ( W X ) " << std::endl; + std::cout << GridLogMessage << " ( -X^* W^* ) " << std::endl; + std::cout <(U,c1,c2); + auto Wstar = PeekIndex(U,c1+nsp,c2+nsp); + auto Ww = conjugate( Wstar ); + auto amizero = sum(W - Ww); + auto amizeroo = TensorRemove(amizero); + assert( amizeroo.real() < 10e-6 ); + amizeroo *= i; + assert( amizeroo.real() < 10e-6 ); + } + + } - Complex a(25041994., 12.); - Complex b(39., 0.22); - Complex d(10000., -2222.3333); + for (int c1 = 0; c1 < nsp ; c1++) + { + for (int c2 = 0; c2 < nsp; c2++) + { + auto X = PeekIndex(U,c1,c2+nsp); + auto minusXstar = PeekIndex(U,c1+nsp,c2); + auto minusXx = conjugate(minusXstar); + auto amizero = sum (X + minusXx); + auto amizeroo = TensorRemove(amizero); + assert( amizeroo.real() < 10e-6 ); + amizeroo *= i; + assert( amizeroo.real() < 10e-6 ); + } + } + + std::cout << GridLogMessage << "Testing ProjectGn" << std::endl; + U = U + 2932.111*identity; + std::cout << GridLogMessage << "Apply ProjectGn to deformed matrix" << std::endl; + Sp::ProjectGn(U); + aux = U*adj(U) - identity; + std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl; + assert( norm2(aux) < 1e-8); + std::cout << GridLogMessage << "Det after ProjectGn = " << norm2( Determinant(U) ) / vol << std::endl; + assert( norm2(aux) - 1 < 1e-8); + Sp::OmegaInvariance(U); + + std::cout << GridLogMessage << "Checking the structure is " << std::endl; + std::cout << GridLogMessage << "U = ( W X ) " << std::endl; + std::cout << GridLogMessage << " ( -X^* W^* ) " << std::endl; + std::cout <(U,c1,c2); + auto Wstar = PeekIndex(U,c1+nsp,c2+nsp); + auto Ww = conjugate( Wstar ); + auto amizero = sum(W - Ww); + auto amizeroo = TensorRemove(amizero); + assert( amizeroo.real() < 10e-6 ); + amizeroo *= i; + assert( amizeroo.real() < 10e-6 ); + } + + } + + for (int c1 = 0; c1 < nsp ; c1++) + { + for (int c2 = 0; c2 < nsp; c2++) + { + auto X = PeekIndex(U,c1,c2+nsp); + auto minusXstar = PeekIndex(U,c1+nsp,c2); + auto minusXx = conjugate(minusXstar); + auto amizero = sum (X + minusXx); + auto amizeroo = TensorRemove(amizero); + assert( amizeroo.real() < 10e-6 ); + amizeroo *= i; + assert( amizeroo.real() < 10e-6 ); + } + } + + std::cout << GridLogMessage << "Testing ProjectSp2nAlgebra" << std::endl; + + U = PeekIndex(Umu,1); + U = U + 666.*identity; + aux = U*adj(U) - identity; + std::cout << GridLogMessage << "Matrix deformed " << std::endl; + std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl; + std::cout << GridLogMessage << "Apply ProjectSp2nAlgebra to deformed matrix" << std::endl; + U = ProjectSp2nAlgebra(U); + aux = U*adj(U) - identity; + std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl; + assert( norm2(aux) < 1e-8); + std::cout << GridLogMessage << "Check that Omega U Omega = conj(U)" << std::endl; + + LatticeColourMatrixD Omega(&Grid); + + Sp::Omega(Omega); + aux = Omega*U*Omega - conjugate(U); + std::cout << GridLogMessage << "Omega U Omega - conj(U) = " << norm2(aux) << std::endl; + assert( norm2(aux) < 1e-8); + + + std::cout << GridLogMessage << "Checking the structure is " << std::endl; + std::cout << GridLogMessage << "U = ( W X ) " << std::endl; + std::cout << GridLogMessage << " ( X^* -W^* ) " << std::endl; + std::cout <(U,c1,c2); + auto Wstar = PeekIndex(U,c1+nsp,c2+nsp); + auto Ww = conjugate( Wstar ); + auto amizero = sum(W + Ww); + auto amizeroo = TensorRemove(amizero); + assert( amizeroo.real() < 10e-6 ); + amizeroo *= i; + assert( amizeroo.real() < 10e-6 ); + } + + } + + for (int c1 = 0; c1 < nsp ; c1++) + { + for (int c2 = 0; c2 < nsp; c2++) + { + auto X = PeekIndex(U,c1,c2+nsp); + auto minusXstar = PeekIndex(U,c1+nsp,c2); + auto minusXx = conjugate(minusXstar); + auto amizero = sum (X - minusXx); + auto amizeroo = TensorRemove(amizero); + assert( amizeroo.real() < 10e-6 ); + amizeroo *= i; + assert( amizeroo.real() < 10e-6 ); + } + } + - A()()(0,0) = a; - A()()(0,1) = b; - A()()(1,0) = i; - A()()(1,1) = d; - std::cout <(Umu,1); U = U + 666.*identity; - aux = U*adj(U) - identity; std::cout << GridLogMessage << "Matrix deformed " << std::endl; - std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl; - std::cout << GridLogMessage << "Apply ProjectSp2nAlgebra to deformed matrix" << std::endl; - U = ProjectSp2nAlgebra(U); - aux = U*adj(U) - identity; - std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl; - assert( norm2(aux) < 1e-8); + std::cout << GridLogMessage << "Apply SpTa to deformed matrix" << std::endl; + U = SpTa(U); + + aux = U - adj(U); + std::cout << GridLogMessage << "SpTa ::: T - Tda = " << norm2(aux) << std::endl; + aux = U + adj(U); + std::cout << GridLogMessage << "SpTa ::: T + Tda = " << norm2(aux) << std::endl; + std::cout << GridLogMessage << "Check that Omega U Omega = conj(U)" << std::endl; - + LatticeColourMatrixD Omega(&Grid); - Sp::Omega(Omega); aux = Omega*U*Omega - conjugate(U); std::cout << GridLogMessage << "Omega U Omega - conj(U) = " << norm2(aux) << std::endl; @@ -258,6 +268,71 @@ int main (int argc, char **argv) } } + //test Ta + /* + U = U + 666.*identity; + Up = Ta(U); + aux = Up - adj(Up); + std::cout << GridLogMessage << "TA !!! T - Tda = " << norm2(aux) << std::endl; + aux = Up + adj(Up); + std::cout << GridLogMessage << "TA !!! T + Tda = " << norm2(aux) << std::endl;*/ + + // test taProj + std::cout << GridLogMessage << "Testing taProj" << std::endl; + U = U + 666.*identity; + std::cout << GridLogMessage << "Matrix deformed " << std::endl; + std::cout << GridLogMessage << "Apply taProj to deformed matrix" << std::endl; + Sp::taProj(U, Up); + aux = Up - adj(Up); + std::cout << GridLogMessage << "taProj ::: T - Tda = " << norm2(aux) << std::endl; + aux = Up + adj(Up); + std::cout << GridLogMessage << "taProj ::: T + Tda = " << norm2(aux) << std::endl; + + std::cout << GridLogMessage << "Check that Omega U Omega = conj(U)" << std::endl; + Sp::Omega(Omega); + aux = Omega*Up*Omega - conjugate(Up); + std::cout << GridLogMessage << "Omega U Omega - conj(U) = " << norm2(aux) << std::endl; + assert( norm2(aux) < 1e-8); + // before it was + aux = Omega*U*Omega - conjugate(U); + std::cout << GridLogMessage << " before taProj Omega U Omega - conj(U) = " << norm2(aux) << std::endl; + + U = Up; + + std::cout << GridLogMessage << "Checking the structure is " << std::endl; + std::cout << GridLogMessage << "U = ( W X ) " << std::endl; + std::cout << GridLogMessage << " ( X^* -W^* ) " << std::endl; + std::cout <(U,c1,c2); + auto Wstar = PeekIndex(U,c1+nsp,c2+nsp); + auto Ww = conjugate( Wstar ); + auto amizero = sum(W + Ww); + auto amizeroo = TensorRemove(amizero); + assert( amizeroo.real() < 10e-6 ); + amizeroo *= i; + assert( amizeroo.real() < 10e-6 ); + } + + } + + for (int c1 = 0; c1 < nsp ; c1++) + { + for (int c2 = 0; c2 < nsp; c2++) + { + auto X = PeekIndex(U,c1,c2+nsp); + auto minusXstar = PeekIndex(U,c1+nsp,c2); + auto minusXx = conjugate(minusXstar); + auto amizero = sum (X - minusXx); + auto amizeroo = TensorRemove(amizero); + assert( amizeroo.real() < 10e-6 ); + amizeroo *= i; + assert( amizeroo.real() < 10e-6 ); + } + } Grid_finalize(); From 4ea29b8f0f1c619bebe1853b9b2dd39ced8f4762 Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Tue, 4 Apr 2023 17:49:28 +0100 Subject: [PATCH 03/20] Template group into GaugeImplTypes. Closing #2 --- Grid/qcd/action/gauge/GaugeImplTypes.h | 106 ++++++------------------- 1 file changed, 24 insertions(+), 82 deletions(-) diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h index 5a8887f2..17bee86a 100644 --- a/Grid/qcd/action/gauge/GaugeImplTypes.h +++ b/Grid/qcd/action/gauge/GaugeImplTypes.h @@ -62,7 +62,7 @@ NAMESPACE_BEGIN(Grid); // hardcodes the exponential approximation in the template //template class GaugeImplTypes { -template class GaugeImplTypes { +template > class GaugeImplTypes { public: typedef S Simd; typedef typename Simd::scalar_type scalar_type; @@ -80,7 +80,7 @@ public: typedef Lattice LinkField; typedef Lattice Field; - typedef SU Group; + //typedef SU Group; // Guido: we can probably separate the types from the HMC functions // this will create 2 kind of implementations @@ -124,53 +124,31 @@ public: for (int mu = 0; mu < Nd; mu++) { - if (isSp2n == true) - { - //const int nSp = Nrepresentation/2; - Sp::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); - } else - { - - - Group::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); - } - + Group::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR) ; Pmu = Pmu*scale; PokeIndex(P, Pmu, mu); } } - - static inline Field projectForce(Field &P) { - if (isSp2n == true) - { - P = Ta(P); - //const int nsp = Nc / 2; - - Sp::iGroupMatrix gen; - - - auto Psum = P; - - Psum = Zero(); - - for (int a = 0; a < Sp::AlgebraDimension; a++) + /* // this works + iScalar > > gen; + auto Psum = P; + Psum = Zero(); + for (int a = 0; a < Group::AlgebraDimension; a++) { - Sp::generator(a, gen); - + Group::generator(a, gen); auto coeff = 2. * trace(P * gen); Psum += coeff * gen; - } - return Psum; - - } else - { - return Ta(P); - } + + return Psum;}*/ + //this doesnt + Field ret(P.Grid()); + Group::taProj(P, ret); + return ret; } static inline void update_field(Field& P, Field& U, double ep){ @@ -181,14 +159,8 @@ public: autoView(P_v,P,AcceleratorRead); accelerator_for(ss, P.Grid()->oSites(),1,{ for (int mu = 0; mu < Nd; mu++) { - if (isSp2n == true) - { - U_v[ss](mu) = ProjectOnSpGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu)); - } else - { - U_v[ss](mu) = ProjectOnGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu)); - } - + U_v[ss](mu) = Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu); // + Group::ProjectOnGaugeGroup(U_v[ss](mu)); } }); //auto end = std::chrono::high_resolution_clock::now(); @@ -209,52 +181,22 @@ public: static inline void Project(Field &U) { - if (isSp2n == true) - { - ProjectSp2n(U); - } else - { - ProjectSUn(U); - } + Group::ProjectGn(U); } - static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) { - if (isSp2n == true) - { - //const int nSp = Nrepresentation/2; - Sp::HotConfiguration(pRNG, U); - } else - { - Group::HotConfiguration(pRNG, U); - } + Group::HotConfiguration(pRNG, U); } static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) { - if (isSp2n == true) - { - //const int nSp = Nrepresentation/2; - Sp::TepidConfiguration(pRNG, U); - } else - { - Group::TepidConfiguration(pRNG, U); - } + Group::TepidConfiguration(pRNG, U); } - static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) { - if (isSp2n == true) - { - //const int nSp = Nrepresentation/2; - Sp::ColdConfiguration(pRNG, U); - } else - { - Group::ColdConfiguration(pRNG, U); - } - + Group::ColdConfiguration(pRNG, U); } }; @@ -264,9 +206,9 @@ typedef GaugeImplTypes GimplTypesR; typedef GaugeImplTypes GimplTypesF; typedef GaugeImplTypes GimplTypesD; -typedef GaugeImplTypes SymplGimplTypesR; -typedef GaugeImplTypes SymplGimplTypesF; -typedef GaugeImplTypes SymplGimplTypesD; +typedef GaugeImplTypes> SymplGimplTypesR; +typedef GaugeImplTypes> SymplGimplTypesF; +typedef GaugeImplTypes> SymplGimplTypesD; typedef GaugeImplTypes::AdjointDimension> GimplAdjointTypesR; typedef GaugeImplTypes::AdjointDimension> GimplAdjointTypesF; From 178376f24bb99a99f268573a0239f88cd897e3b1 Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Thu, 6 Apr 2023 12:08:17 +0100 Subject: [PATCH 04/20] minor stylistic changes --- Grid/qcd/action/gauge/GaugeImplTypes.h | 21 ++++----------------- Grid/qcd/utils/GaugeGroup.h | 2 +- tests/core/Test_reunitarise.cc | 12 ++++-------- tests/sp2n/Test_project_on_Sp.cc | 14 ++++++-------- 4 files changed, 15 insertions(+), 34 deletions(-) diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h index 17bee86a..e2d98b85 100644 --- a/Grid/qcd/action/gauge/GaugeImplTypes.h +++ b/Grid/qcd/action/gauge/GaugeImplTypes.h @@ -133,19 +133,6 @@ public: static inline Field projectForce(Field &P) { - /* // this works - iScalar > > gen; - auto Psum = P; - Psum = Zero(); - for (int a = 0; a < Group::AlgebraDimension; a++) - { - Group::generator(a, gen); - auto coeff = 2. * trace(P * gen); - Psum += coeff * gen; - - - return Psum;}*/ - //this doesnt Field ret(P.Grid()); Group::taProj(P, ret); return ret; @@ -159,7 +146,7 @@ public: autoView(P_v,P,AcceleratorRead); accelerator_for(ss, P.Grid()->oSites(),1,{ for (int mu = 0; mu < Nd; mu++) { - U_v[ss](mu) = Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu); // + U_v[ss](mu) = Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu); Group::ProjectOnGaugeGroup(U_v[ss](mu)); } }); @@ -206,9 +193,9 @@ typedef GaugeImplTypes GimplTypesR; typedef GaugeImplTypes GimplTypesF; typedef GaugeImplTypes GimplTypesD; -typedef GaugeImplTypes> SymplGimplTypesR; -typedef GaugeImplTypes> SymplGimplTypesF; -typedef GaugeImplTypes> SymplGimplTypesD; +typedef GaugeImplTypes > SymplGimplTypesR; +typedef GaugeImplTypes > SymplGimplTypesF; +typedef GaugeImplTypes > SymplGimplTypesD; typedef GaugeImplTypes::AdjointDimension> GimplAdjointTypesR; typedef GaugeImplTypes::AdjointDimension> GimplAdjointTypesF; diff --git a/Grid/qcd/utils/GaugeGroup.h b/Grid/qcd/utils/GaugeGroup.h index 8e6c5ee4..0e8fbabe 100644 --- a/Grid/qcd/utils/GaugeGroup.h +++ b/Grid/qcd/utils/GaugeGroup.h @@ -170,7 +170,7 @@ class GaugeGroup { } } - template + template static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, double scale = 1.0) { GridBase *grid = out.Grid(); diff --git a/tests/core/Test_reunitarise.cc b/tests/core/Test_reunitarise.cc index e421a2a8..110ebabf 100644 --- a/tests/core/Test_reunitarise.cc +++ b/tests/core/Test_reunitarise.cc @@ -122,15 +122,16 @@ int main (int argc, char ** argv) detU=detU-1.0; std::cout << "Determinant defect before projection " <gSites(); @@ -118,7 +119,7 @@ int main (int argc, char **argv) std::cout << GridLogMessage << "Testing ProjectOnGaugeGroup" << std::endl; - U = U + 2932.111*identity; + U = U + Delta*identity; std::cout << GridLogMessage << "Apply ProjectOnGaugeGroup to deformed matrix" << std::endl; Sp::ProjectOnGaugeGroup(U); aux = U*adj(U) - identity; @@ -162,7 +163,7 @@ int main (int argc, char **argv) } std::cout << GridLogMessage << "Testing ProjectGn" << std::endl; - U = U + 2932.111*identity; + U = U + Delta*identity; std::cout << GridLogMessage << "Apply ProjectGn to deformed matrix" << std::endl; Sp::ProjectGn(U); aux = U*adj(U) - identity; @@ -214,7 +215,7 @@ int main (int argc, char **argv) std::cout << GridLogMessage << "Testing SpTa" << std::endl; U = PeekIndex(Umu,1); - U = U + 666.*identity; + U = U + Delta*identity; std::cout << GridLogMessage << "Matrix deformed " << std::endl; std::cout << GridLogMessage << "Apply SpTa to deformed matrix" << std::endl; U = SpTa(U); @@ -279,7 +280,7 @@ int main (int argc, char **argv) // test taProj std::cout << GridLogMessage << "Testing taProj" << std::endl; - U = U + 666.*identity; + U = U + Delta*identity; std::cout << GridLogMessage << "Matrix deformed " << std::endl; std::cout << GridLogMessage << "Apply taProj to deformed matrix" << std::endl; Sp::taProj(U, Up); @@ -333,10 +334,7 @@ int main (int argc, char **argv) assert( amizeroo.real() < 10e-6 ); } } - - + Grid_finalize(); - } - From be98d26610822bc39b847f18e52547978c54e79b Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Thu, 13 Apr 2023 17:48:43 +0100 Subject: [PATCH 05/20] small change I missed in previous commit --- Grid/qcd/action/gauge/GaugeImplTypes.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h index e2d98b85..036d93ef 100644 --- a/Grid/qcd/action/gauge/GaugeImplTypes.h +++ b/Grid/qcd/action/gauge/GaugeImplTypes.h @@ -80,8 +80,6 @@ public: typedef Lattice LinkField; typedef Lattice Field; - //typedef SU Group; - // Guido: we can probably separate the types from the HMC functions // this will create 2 kind of implementations // probably confusing the users From dace904c10c32d9513a0298313dfe17e1f064457 Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Fri, 14 Apr 2023 18:06:18 +0100 Subject: [PATCH 06/20] fix typo --- Grid/tensors/Tensor_Ta.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Grid/tensors/Tensor_Ta.h b/Grid/tensors/Tensor_Ta.h index a7103343..c9132854 100644 --- a/Grid/tensors/Tensor_Ta.h +++ b/Grid/tensors/Tensor_Ta.h @@ -72,7 +72,7 @@ template accelerator_inline iScalar SpTa(const iScalar ret; ret._internal = SpTa(r._internal); - return Ta(ret); + return ret; } template accelerator_inline iVector SpTa(const iVector&r) { @@ -80,11 +80,11 @@ template accelerator_inline iVector SpTa(const iVect for(int i=0;i accelerator_inline iMatrix SpTa(const iMatrix &arg) { - iMatrix ret; + iMatrix ret(arg); vtype nrm; vtype inner; vtype tmp; From 5aabe074fe36f24a18ef8933f10c14281f7e1318 Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Tue, 18 Apr 2023 11:50:20 +0100 Subject: [PATCH 07/20] Rename Sympl* to Sp* --- Grid/qcd/action/gauge/Gauge.h | 6 +++--- Grid/qcd/action/gauge/GaugeImplTypes.h | 6 +++--- Grid/qcd/action/gauge/GaugeImplementations.h | 6 +++--- Grid/qcd/hmc/GenericHMCrunner.h | 4 ++-- tests/sp2n/Test_hmc_Sp_WF_2_Fund_3_2AS.cc | 2 +- tests/sp2n/Test_hmc_Sp_Wilson2ASFermionGauge.cc | 2 +- tests/sp2n/Test_hmc_Sp_WilsonFundFermionGauge.cc | 2 +- tests/sp2n/Test_hmc_Sp_pureGaugeWilson.cc | 2 +- 8 files changed, 15 insertions(+), 15 deletions(-) diff --git a/Grid/qcd/action/gauge/Gauge.h b/Grid/qcd/action/gauge/Gauge.h index 535f0845..7a72d087 100644 --- a/Grid/qcd/action/gauge/Gauge.h +++ b/Grid/qcd/action/gauge/Gauge.h @@ -39,9 +39,9 @@ NAMESPACE_BEGIN(Grid); typedef WilsonGaugeAction WilsonGaugeActionR; typedef WilsonGaugeAction WilsonGaugeActionF; typedef WilsonGaugeAction WilsonGaugeActionD; -typedef WilsonGaugeAction SymplWilsonGaugeActionR; -typedef WilsonGaugeAction SymplWilsonGaugeActionF; -typedef WilsonGaugeAction SymplWilsonGaugeActionD; +typedef WilsonGaugeAction SpWilsonGaugeActionR; +typedef WilsonGaugeAction SpWilsonGaugeActionF; +typedef WilsonGaugeAction SpWilsonGaugeActionD; typedef PlaqPlusRectangleAction PlaqPlusRectangleActionR; typedef PlaqPlusRectangleAction PlaqPlusRectangleActionF; typedef PlaqPlusRectangleAction PlaqPlusRectangleActionD; diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h index 036d93ef..c3e42f80 100644 --- a/Grid/qcd/action/gauge/GaugeImplTypes.h +++ b/Grid/qcd/action/gauge/GaugeImplTypes.h @@ -191,9 +191,9 @@ typedef GaugeImplTypes GimplTypesR; typedef GaugeImplTypes GimplTypesF; typedef GaugeImplTypes GimplTypesD; -typedef GaugeImplTypes > SymplGimplTypesR; -typedef GaugeImplTypes > SymplGimplTypesF; -typedef GaugeImplTypes > SymplGimplTypesD; +typedef GaugeImplTypes > SpGimplTypesR; +typedef GaugeImplTypes > SpGimplTypesF; +typedef GaugeImplTypes > SpGimplTypesD; typedef GaugeImplTypes::AdjointDimension> GimplAdjointTypesR; typedef GaugeImplTypes::AdjointDimension> GimplAdjointTypesF; diff --git a/Grid/qcd/action/gauge/GaugeImplementations.h b/Grid/qcd/action/gauge/GaugeImplementations.h index bc056263..4a5a6651 100644 --- a/Grid/qcd/action/gauge/GaugeImplementations.h +++ b/Grid/qcd/action/gauge/GaugeImplementations.h @@ -155,9 +155,9 @@ typedef ConjugateGaugeImpl ConjugateGimplR; // Real.. whichever pre typedef ConjugateGaugeImpl ConjugateGimplF; // Float typedef ConjugateGaugeImpl ConjugateGimplD; // Double -typedef PeriodicGaugeImpl SymplPeriodicGimplR; // Real.. whichever prec -typedef PeriodicGaugeImpl SymplPeriodicGimplF; // Float -typedef PeriodicGaugeImpl SymplPeriodicGimplD; // Double +typedef PeriodicGaugeImpl SpPeriodicGimplR; // Real.. whichever prec +typedef PeriodicGaugeImpl SpPeriodicGimplF; // Float +typedef PeriodicGaugeImpl SpPeriodicGimplD; // Double NAMESPACE_END(Grid); diff --git a/Grid/qcd/hmc/GenericHMCrunner.h b/Grid/qcd/hmc/GenericHMCrunner.h index 16143455..a9b51a4d 100644 --- a/Grid/qcd/hmc/GenericHMCrunner.h +++ b/Grid/qcd/hmc/GenericHMCrunner.h @@ -207,12 +207,12 @@ using GenericHMCRunnerHirep = // sp2n template