From f51222086c3cc795808a0adaa820db0245d6c48b Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Thu, 9 Mar 2023 16:22:20 +0000 Subject: [PATCH 1/9] Move functions from GaugeGroup to group specific implementations --- Grid/qcd/utils/GaugeGroup.h | 94 +++--------------------------- Grid/qcd/utils/SUn.impl | 113 ++++++++++++++++++++++++++++++++++++ Grid/qcd/utils/Sp2n.impl | 113 ++++++++++++++++++++++++++++++++++++ 3 files changed, 235 insertions(+), 85 deletions(-) diff --git a/Grid/qcd/utils/GaugeGroup.h b/Grid/qcd/utils/GaugeGroup.h index ecd92780..f0442e2f 100644 --- a/Grid/qcd/utils/GaugeGroup.h +++ b/Grid/qcd/utils/GaugeGroup.h @@ -170,118 +170,40 @@ class GaugeGroup { } } - // reunitarise?? template static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, double scale = 1.0) { - GridBase *grid = out.Grid(); - - typedef typename LatticeMatrixType::vector_type vector_type; - typedef typename LatticeMatrixType::scalar_type scalar_type; - - typedef iSinglet vTComplexType; - - typedef Lattice LatticeComplexType; - typedef typename GridTypeMapper< - typename LatticeMatrixType::vector_object>::scalar_object MatrixType; - - LatticeComplexType ca(grid); - LatticeMatrixType lie(grid); - LatticeMatrixType la(grid); - ComplexD ci(0.0, scale); - // ComplexD cone(1.0, 0.0); - MatrixType ta; - - lie = Zero(); - - for (int a = 0; a < AdjointDimension; a++) { - random(pRNG, ca); - - ca = (ca + conjugate(ca)) * 0.5; - ca = ca - 0.5; - - generator(a, ta); - - la = ci * ca * ta; - - lie = lie + la; // e^{i la ta} - } - taExp(lie, out); + LieRandomize(pRNG, out, group_name(), scale); } static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG, LatticeMatrix &out, Real scale = 1.0) { - GridBase *grid = out.Grid(); - LatticeReal ca(grid); - LatticeMatrix la(grid); - Complex ci(0.0, scale); - Matrix ta; - - out = Zero(); - for (int a = 0; a < AdjointDimension; a++) { - gaussian(pRNG, ca); - generator(a, ta); - la = toComplex(ca) * ta; - out += la; - } - out *= ci; + GaussianFundamentalLieAlgebraMatrix(pRNG, out, group_name(), scale); } static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h, LatticeMatrix &out, Real scale = 1.0) { - conformable(h, out); - GridBase *grid = out.Grid(); - LatticeMatrix la(grid); - Matrix ta; - - out = Zero(); - for (int a = 0; a < AdjointDimension; a++) { - generator(a, ta); - la = peekColour(h, a) * timesI(ta) * scale; - out += la; - } + FundamentalLieAlgebraMatrix(h, out, group_name(), scale); } // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 // ) inverse operation: FundamentalLieAlgebraMatrix static void projectOnAlgebra(LatticeAlgebraVector &h_out, const LatticeMatrix &in, Real scale = 1.0) { - conformable(h_out, in); - h_out = Zero(); - Matrix Ta; - - for (int a = 0; a < AdjointDimension; a++) { - generator(a, Ta); - pokeColour(h_out, -2.0 * (trace(timesI(Ta) * in)) * scale, a); - } + projectOnAlgebra(h_out, in, group_name(), scale); } template static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) { - typedef typename GaugeField::vector_type vector_type; - typedef iGroupMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType Umu(out.Grid()); - for (int mu = 0; mu < Nd; mu++) { - LieRandomize(pRNG, Umu, 1.0); - PokeIndex(out, Umu, mu); - } + HotConfiguration(pRNG, out, group_name()); } template static void TepidConfiguration(GridParallelRNG &pRNG, GaugeField &out) { - typedef typename GaugeField::vector_type vector_type; - typedef iGroupMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType Umu(out.Grid()); - for (int mu = 0; mu < Nd; mu++) { - LieRandomize(pRNG, Umu, 0.01); - PokeIndex(out, Umu, mu); - } + TepidConfiguration(pRNG, out, group_name()); } + template static void ColdConfiguration(GaugeField &out) { typedef typename GaugeField::vector_type vector_type; @@ -294,6 +216,7 @@ class GaugeGroup { PokeIndex(out, Umu, mu); } } + template static void ColdConfiguration(GridParallelRNG &pRNG, GaugeField &out) { ColdConfiguration(out); @@ -303,6 +226,7 @@ class GaugeGroup { static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) { out = Ta(in); } + template static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { typedef typename LatticeMatrixType::scalar_type ComplexType; diff --git a/Grid/qcd/utils/SUn.impl b/Grid/qcd/utils/SUn.impl index eaf82877..eef039c7 100644 --- a/Grid/qcd/utils/SUn.impl +++ b/Grid/qcd/utils/SUn.impl @@ -508,6 +508,119 @@ static void testGenerators(GroupName::SU) { std::cout << GridLogMessage << std::endl; } +template +static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, GroupName::SU, + double scale = 1.0) { + GridBase *grid = out.Grid(); + + typedef typename LatticeMatrixType::vector_type vector_type; + typedef typename LatticeMatrixType::scalar_type scalar_type; + + typedef iSinglet vTComplexType; + + typedef Lattice LatticeComplexType; + typedef typename GridTypeMapper< + typename LatticeMatrixType::vector_object>::scalar_object MatrixType; + + LatticeComplexType ca(grid); + LatticeMatrixType lie(grid); + LatticeMatrixType la(grid); + ComplexD ci(0.0, scale); + // ComplexD cone(1.0, 0.0); + MatrixType ta; + + lie = Zero(); + + for (int a = 0; a < AdjointDimension; a++) { + random(pRNG, ca); + + ca = (ca + conjugate(ca)) * 0.5; + ca = ca - 0.5; + + generator(a, ta); + + la = ci * ca * ta; + + lie = lie + la; // e^{i la ta} + } + taExp(lie, out); +} + +static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG, + LatticeMatrix &out, GroupName::SU, + Real scale = 1.0) { + GridBase *grid = out.Grid(); + LatticeReal ca(grid); + LatticeMatrix la(grid); + Complex ci(0.0, scale); + Matrix ta; + + out = Zero(); + for (int a = 0; a < AdjointDimension; a++) { + gaussian(pRNG, ca); + generator(a, ta); + la = toComplex(ca) * ta; + out += la; + } + out *= ci; +} + +static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h, + LatticeMatrix &out, GroupName::SU, + Real scale = 1.0) { + conformable(h, out); + GridBase *grid = out.Grid(); + LatticeMatrix la(grid); + Matrix ta; + + out = Zero(); + for (int a = 0; a < AdjointDimension; a++) { + generator(a, ta); + la = peekColour(h, a) * timesI(ta) * scale; + out += la; + } +} + +template +static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out, GroupName::SU) { + typedef typename GaugeField::vector_type vector_type; + typedef iGroupMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType Umu(out.Grid()); + for (int mu = 0; mu < Nd; mu++) { + LieRandomize(pRNG, Umu, 1.0); + PokeIndex(out, Umu, mu); + } +} + +template +static void TepidConfiguration(GridParallelRNG &pRNG, GaugeField &out, GroupName::SU) { + typedef typename GaugeField::vector_type vector_type; + typedef iGroupMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType Umu(out.Grid()); + for (int mu = 0; mu < Nd; mu++) { + LieRandomize(pRNG, Umu, 0.01); + PokeIndex(out, Umu, mu); + } +} + +// Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 +// ) inverse operation: FundamentalLieAlgebraMatrix +static void projectOnAlgebra(LatticeAlgebraVector &h_out, + const LatticeMatrix &in, GroupName::SU, Real scale = 1.0) { + conformable(h_out, in); + h_out = Zero(); + Matrix Ta; + + for (int a = 0; a < AdjointDimension; a++) { + generator(a, Ta); + pokeColour(h_out, -2.0 * (trace(timesI(Ta) * in)) * scale, a); + } +} + /* * Fundamental rep gauge xform */ diff --git a/Grid/qcd/utils/Sp2n.impl b/Grid/qcd/utils/Sp2n.impl index c8ed2ddc..411235ec 100644 --- a/Grid/qcd/utils/Sp2n.impl +++ b/Grid/qcd/utils/Sp2n.impl @@ -253,6 +253,119 @@ static void testGenerators(GroupName::Sp) { } } +template + static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, GroupName::Sp, + double scale = 1.0) { + GridBase *grid = out.Grid(); + + typedef typename LatticeMatrixType::vector_type vector_type; + typedef typename LatticeMatrixType::scalar_type scalar_type; + + typedef iSinglet vTComplexType; + + typedef Lattice LatticeComplexType; + typedef typename GridTypeMapper< + typename LatticeMatrixType::vector_object>::scalar_object MatrixType; + + LatticeComplexType ca(grid); + LatticeMatrixType lie(grid); + LatticeMatrixType la(grid); + ComplexD ci(0.0, scale); + // ComplexD cone(1.0, 0.0); + MatrixType ta; + + lie = Zero(); + + for (int a = 0; a < AdjointDimension; a++) { + random(pRNG, ca); + + ca = (ca + conjugate(ca)) * 0.5; + ca = ca - 0.5; + + generator(a, ta); + + la = ci * ca * ta; + + lie = lie + la; // e^{i la ta} + } + taExp(lie, out); + } + +static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG, + LatticeMatrix &out, GroupName::Sp, + Real scale = 1.0) { + GridBase *grid = out.Grid(); + LatticeReal ca(grid); + LatticeMatrix la(grid); + Complex ci(0.0, scale); + Matrix ta; + + out = Zero(); + for (int a = 0; a < AdjointDimension; a++) { + gaussian(pRNG, ca); + generator(a, ta); + la = toComplex(ca) * ta; + out += la; + } + out *= ci; +} + +static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h, + LatticeMatrix &out, GroupName::Sp, + Real scale = 1.0) { + conformable(h, out); + GridBase *grid = out.Grid(); + LatticeMatrix la(grid); + Matrix ta; + + out = Zero(); + for (int a = 0; a < AdjointDimension; a++) { + generator(a, ta); + la = peekColour(h, a) * timesI(ta) * scale; + out += la; + } +} + +// Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 +// ) inverse operation: FundamentalLieAlgebraMatrix +static void projectOnAlgebra(LatticeAlgebraVector &h_out, + const LatticeMatrix &in, GroupName::Sp, Real scale = 1.0) { + conformable(h_out, in); + h_out = Zero(); + Matrix Ta; + + for (int a = 0; a < AdjointDimension; a++) { + generator(a, Ta); + pokeColour(h_out, -2.0 * (trace(timesI(Ta) * in)) * scale, a); + } +} + +template +static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out, GroupName::Sp) { + typedef typename GaugeField::vector_type vector_type; + typedef iGroupMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType Umu(out.Grid()); + for (int mu = 0; mu < Nd; mu++) { + LieRandomize(pRNG, Umu, 1.0); + PokeIndex(out, Umu, mu); + } +} + +template +static void TepidConfiguration(GridParallelRNG &pRNG, GaugeField &out, GroupName::Sp) { + typedef typename GaugeField::vector_type vector_type; + typedef iGroupMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType Umu(out.Grid()); + for (int mu = 0; mu < Nd; mu++) { + LieRandomize(pRNG, Umu, 0.01); + PokeIndex(out, Umu, mu); + } +} + public: template static void OmegaInvariance(ColourMatrix &in) { From fd057c838f4788f809b49ce137b33acd6b3a4377 Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Fri, 10 Mar 2023 12:10:46 +0000 Subject: [PATCH 2/9] add ProjectOnGaugeGroup and ProjectGn to allow future templating in GaugeImplTypes --- Grid/qcd/utils/GaugeGroup.h | 27 ++++++++++++++++++++++++ Grid/qcd/utils/SUn.impl | 41 +++++++++++++++++++++++++++++++++++++ Grid/qcd/utils/Sp2n.impl | 41 +++++++++++++++++++++++++++++++++++++ 3 files changed, 109 insertions(+) diff --git a/Grid/qcd/utils/GaugeGroup.h b/Grid/qcd/utils/GaugeGroup.h index f0442e2f..dea8f985 100644 --- a/Grid/qcd/utils/GaugeGroup.h +++ b/Grid/qcd/utils/GaugeGroup.h @@ -244,6 +244,29 @@ class GaugeGroup { ex = ex + xn * nfac; // x2/2!, x3/3!.... } } + + template // reunitarise, resimplectify... + static void ProjectOnGaugeGroup(Lattice >, Nd> > &U) { + ProjectOnGaugeGroup(U, group_name()); + } + + template // reunitarise, resimplectify... + static void ProjectOnGaugeGroup(Lattice > > > &Umu) { + ProjectOnGaugeGroup(Umu, group_name()); + } + + + template // reunitarise, resimplectify... previously ProjectSUn + static void ProjectGn(Lattice > > > &Umu) { + ProjectGn(Umu, group_name()); + } + + template // reunitarise, resimplectify... previously ProjectSUn + static void ProjectGn(Lattice >, Nd> > &U) { + ProjectGn(U, group_name()); + + } + }; template @@ -271,6 +294,8 @@ LatticeComplexD Determinant( }); return ret; } + + template static void ProjectSUn( Lattice > > > &Umu) { @@ -285,6 +310,7 @@ static void ProjectSUn( PokeIndex(Umu, element, Nc - 1, i); } } + template static void ProjectSUn( Lattice >, Nd> > &U) { @@ -297,6 +323,7 @@ static void ProjectSUn( PokeIndex(U, Umu, mu); } } + // Explicit specialisation for SU(3). // Explicit specialisation for SU(3). static void ProjectSU3( diff --git a/Grid/qcd/utils/SUn.impl b/Grid/qcd/utils/SUn.impl index eef039c7..18be96b1 100644 --- a/Grid/qcd/utils/SUn.impl +++ b/Grid/qcd/utils/SUn.impl @@ -621,6 +621,47 @@ static void projectOnAlgebra(LatticeAlgebraVector &h_out, } } +template +static void ProjectOnGaugeGroup(Lattice > > > &Umu, GroupName::SU) { + Umu = ProjectOnGroup(Umu); + } + +template +static void ProjectOnGaugeGroup(Lattice >, Nd> > &U, GroupName::SU) { + GridBase *grid = U.Grid(); + // Reunitarise + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnGroup(Umu); + } +} + +template +static void ProjectGn(Lattice > > > &Umu, GroupName::SU) { + Umu = ProjectOnGroup(Umu); + auto det = Determinant(Umu); + + det = conjugate(det); + + for (int i = 0; i < N; i++) { + auto element = PeekIndex(Umu, N - 1, i); + element = element * det; + PokeIndex(Umu, element, Nc - 1, i); + } +} + +template +static void ProjectGn(Lattice >, Nd> > &U, GroupName::SU) { + GridBase *grid = U.Grid(); + // Reunitarise + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnGroup(Umu); + ProjectSUn(Umu); + PokeIndex(U, Umu, mu); + } +} + /* * Fundamental rep gauge xform */ diff --git a/Grid/qcd/utils/Sp2n.impl b/Grid/qcd/utils/Sp2n.impl index 411235ec..44771a14 100644 --- a/Grid/qcd/utils/Sp2n.impl +++ b/Grid/qcd/utils/Sp2n.impl @@ -366,6 +366,47 @@ static void TepidConfiguration(GridParallelRNG &pRNG, GaugeField &out, GroupName } } +template +static void ProjectOnGaugeGroup(Lattice > > > &Umu, GroupName::Sp) { + Umu = ProjectOnSpGroup(Umu); + } + +template +static void ProjectOnGaugeGroup(Lattice >, Nd> > &U, GroupName::Sp) { + GridBase *grid = U.Grid(); + // Reunitarise + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnSpGroup(Umu); + } +} + +template +static void ProjectGn(Lattice > > > &Umu, GroupName::Sp) { + Umu = ProjectOnSpGroup(Umu); + auto det = Determinant(Umu); + + det = conjugate(det); + + for (int i = 0; i < N; i++) { + auto element = PeekIndex(Umu, N - 1, i); + element = element * det; + PokeIndex(Umu, element, Nc - 1, i); + } +} + +template +static void ProjectGn(Lattice >, Nd> > &U, GroupName::Sp) { + GridBase *grid = U.Grid(); + // Reunitarise + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnSpGroup(Umu); + ProjectSUn(Umu); + PokeIndex(U, Umu, mu); + } +} + public: template static void OmegaInvariance(ColourMatrix &in) { From 29586f6b5ee09bc21b61f00230ec16d9b881e83a Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Mon, 13 Mar 2023 08:17:14 +0000 Subject: [PATCH 3/9] Deactivate some tests for Nc!=3 --- tests/core/Test_lie_generators.cc | 2 ++ tests/core/Test_reunitarise.cc | 2 ++ 2 files changed, 4 insertions(+) diff --git a/tests/core/Test_lie_generators.cc b/tests/core/Test_lie_generators.cc index 372e1408..a49dbca1 100644 --- a/tests/core/Test_lie_generators.cc +++ b/tests/core/Test_lie_generators.cc @@ -28,6 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#if Nc==3 #include #include @@ -530,3 +531,4 @@ int main(int argc, char** argv) { Grid_finalize(); } +#endif diff --git a/tests/core/Test_reunitarise.cc b/tests/core/Test_reunitarise.cc index 6644be1a..0ff8ec7e 100644 --- a/tests/core/Test_reunitarise.cc +++ b/tests/core/Test_reunitarise.cc @@ -26,6 +26,7 @@ Author: Peter Boyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#if Nc==3 #include using namespace std; @@ -143,3 +144,4 @@ int main (int argc, char ** argv) +#endif From d6ff644aab51e0e33a8efbfc3b7a862622fd2598 Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Tue, 14 Mar 2023 10:43:25 +0000 Subject: [PATCH 4/9] Towards the day all tests compile --- tests/core/Test_lie_generators.cc | 113 +++++++++++++-------------- tests/core/Test_reunitarise.cc | 10 +-- tests/solver/Test_coarse_even_odd.cc | 2 +- 3 files changed, 62 insertions(+), 63 deletions(-) diff --git a/tests/core/Test_lie_generators.cc b/tests/core/Test_lie_generators.cc index a49dbca1..5b044061 100644 --- a/tests/core/Test_lie_generators.cc +++ b/tests/core/Test_lie_generators.cc @@ -28,7 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#if Nc==3 + #include #include @@ -43,11 +43,11 @@ directory using namespace std; using namespace Grid; - ; + int main(int argc, char** argv) { Grid_init(&argc, &argv); - + std::vector latt({4, 4, 4, 8}); GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid( latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); @@ -67,7 +67,7 @@ int main(int argc, char** argv) { std::cout << GridLogMessage << "*********************************************" << std::endl; - std::cout << GridLogMessage << "* Generators for SU(Nc" << std::endl; + std::cout << GridLogMessage << "* Generators for SU(3)" << std::endl; std::cout << GridLogMessage << "*********************************************" << std::endl; SU3::printGenerators(); @@ -88,22 +88,22 @@ int main(int argc, char** argv) { // Projectors GridParallelRNG gridRNG(grid); gridRNG.SeedFixedIntegers(std::vector({45,12,81,9})); - SU3Adjoint::LatticeAdjMatrix Gauss(grid); - SU3::LatticeAlgebraVector ha(grid); - SU3::LatticeAlgebraVector hb(grid); + SU_Adjoint::LatticeAdjMatrix Gauss(grid); + SU::LatticeAlgebraVector ha(grid); + SU::LatticeAlgebraVector hb(grid); random(gridRNG,Gauss); std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl; - SU3Adjoint::projectOnAlgebra(ha, Gauss); + SU_Adjoint::projectOnAlgebra(ha, Gauss); std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "Start projector" << std::endl; - SU3Adjoint::projector(hb, Gauss); + SU_Adjoint::projector(hb, Gauss); std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "ReStart projector" << std::endl; - SU3Adjoint::projector(hb, Gauss); + SU_Adjoint::projector(hb, Gauss); std::cout << GridLogMessage << "end projector" << std::endl; - SU3::LatticeAlgebraVector diff = ha -hb; + SU::LatticeAlgebraVector diff = ha -hb; std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl; @@ -115,8 +115,8 @@ int main(int argc, char** argv) { LatticeGaugeField U(grid), V(grid); - SU3::HotConfiguration(gridRNG, U); - SU3::HotConfiguration(gridRNG, V); + SU::HotConfiguration(gridRNG, U); + SU::HotConfiguration(gridRNG, V); // Adjoint representation // Test group structure @@ -124,8 +124,8 @@ int main(int argc, char** argv) { LatticeGaugeField UV(grid); UV = Zero(); for (int mu = 0; mu < Nd; mu++) { - SU3::LatticeMatrix Umu = peekLorentz(U,mu); - SU3::LatticeMatrix Vmu = peekLorentz(V,mu); + SU::LatticeMatrix Umu = peekLorentz(U,mu); + SU::LatticeMatrix Vmu = peekLorentz(V,mu); pokeLorentz(UV,Umu*Vmu, mu); } @@ -152,16 +152,16 @@ int main(int argc, char** argv) { // Check correspondence of algebra and group transformations // Create a random vector - SU3::LatticeAlgebraVector h_adj(grid); + SU::LatticeAlgebraVector h_adj(grid); typename AdjointRep::LatticeMatrix Ar(grid); random(gridRNG,h_adj); h_adj = real(h_adj); SU_Adjoint::AdjointLieAlgebraMatrix(h_adj,Ar); // Re-extract h_adj - SU3::LatticeAlgebraVector h_adj2(grid); + SU::LatticeAlgebraVector h_adj2(grid); SU_Adjoint::projectOnAlgebra(h_adj2, Ar); - SU3::LatticeAlgebraVector h_diff = h_adj - h_adj2; + SU::LatticeAlgebraVector h_diff = h_adj - h_adj2; std::cout << GridLogMessage << "Projections structure check vector difference (Adjoint representation) : " << norm2(h_diff) << std::endl; // Exponentiate @@ -184,14 +184,14 @@ int main(int argc, char** argv) { // Construct the fundamental matrix in the group - SU3::LatticeMatrix Af(grid); - SU3::FundamentalLieAlgebraMatrix(h_adj,Af); - SU3::LatticeMatrix Ufund(grid); + SU::LatticeMatrix Af(grid); + SU::FundamentalLieAlgebraMatrix(h_adj,Af); + SU::LatticeMatrix Ufund(grid); Ufund = expMat(Af, 1.0, 16); // Check unitarity - SU3::LatticeMatrix uno_f(grid); + SU::LatticeMatrix uno_f(grid); uno_f = 1.0; - SU3::LatticeMatrix UnitCheck(grid); + SU::LatticeMatrix UnitCheck(grid); UnitCheck = Ufund * adj(Ufund) - uno_f; std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck) << std::endl; @@ -261,20 +261,20 @@ int main(int argc, char** argv) { std::cout << GridLogMessage << "Test for the Two Index Symmetric projectors" << std::endl; // Projectors - SU3TwoIndexSymm::LatticeTwoIndexMatrix Gauss2(grid); + SU_TwoIndex::LatticeTwoIndexMatrix Gauss2(grid); random(gridRNG,Gauss2); std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl; - SU3TwoIndexSymm::projectOnAlgebra(ha, Gauss2); + SU_TwoIndex::projectOnAlgebra(ha, Gauss2); std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "Start projector" << std::endl; - SU3TwoIndexSymm::projector(hb, Gauss2); + SU_TwoIndex::projector(hb, Gauss2); std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "ReStart projector" << std::endl; - SU3TwoIndexSymm::projector(hb, Gauss2); + SU_TwoIndex::projector(hb, Gauss2); std::cout << GridLogMessage << "end projector" << std::endl; - SU3::LatticeAlgebraVector diff2 = ha - hb; + SU::LatticeAlgebraVector diff2 = ha - hb; std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl; std::cout << GridLogMessage << "*********************************************" << std::endl; @@ -285,20 +285,20 @@ int main(int argc, char** argv) { std::cout << GridLogMessage << "Test for the Two index anti-Symmetric projectors" << std::endl; // Projectors - SU3TwoIndexAntiSymm::LatticeTwoIndexMatrix Gauss2a(grid); + SU_TwoIndex::LatticeTwoIndexMatrix Gauss2a(grid); random(gridRNG,Gauss2a); std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl; - SU3TwoIndexAntiSymm::projectOnAlgebra(ha, Gauss2a); + SU_TwoIndex::projectOnAlgebra(ha, Gauss2a); std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "Start projector" << std::endl; - SU3TwoIndexAntiSymm::projector(hb, Gauss2a); + SU_TwoIndex::projector(hb, Gauss2a); std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "ReStart projector" << std::endl; - SU3TwoIndexAntiSymm::projector(hb, Gauss2a); + SU_TwoIndex::projector(hb, Gauss2a); std::cout << GridLogMessage << "end projector" << std::endl; - SU3::LatticeAlgebraVector diff2a = ha - hb; + SU::LatticeAlgebraVector diff2a = ha - hb; std::cout << GridLogMessage << "Difference: " << norm2(diff2a) << std::endl; std::cout << GridLogMessage << "*********************************************" << std::endl; @@ -312,14 +312,14 @@ int main(int argc, char** argv) { // Test group structure // (U_f * V_f)_r = U_r * V_r LatticeGaugeField U2(grid), V2(grid); - SU3::HotConfiguration(gridRNG, U2); - SU3::HotConfiguration(gridRNG, V2); + SU::HotConfiguration(gridRNG, U2); + SU::HotConfiguration(gridRNG, V2); LatticeGaugeField UV2(grid); UV2 = Zero(); for (int mu = 0; mu < Nd; mu++) { - SU3::LatticeMatrix Umu2 = peekLorentz(U2,mu); - SU3::LatticeMatrix Vmu2 = peekLorentz(V2,mu); + SU::LatticeMatrix Umu2 = peekLorentz(U2,mu); + SU::LatticeMatrix Vmu2 = peekLorentz(V2,mu); pokeLorentz(UV2,Umu2*Vmu2, mu); } @@ -346,16 +346,16 @@ int main(int argc, char** argv) { // Check correspondence of algebra and group transformations // Create a random vector - SU3::LatticeAlgebraVector h_sym(grid); + SU::LatticeAlgebraVector h_sym(grid); typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ar_sym(grid); random(gridRNG,h_sym); h_sym = real(h_sym); SU_TwoIndex::TwoIndexLieAlgebraMatrix(h_sym,Ar_sym); // Re-extract h_sym - SU3::LatticeAlgebraVector h_sym2(grid); + SU::LatticeAlgebraVector h_sym2(grid); SU_TwoIndex< Nc, Symmetric>::projectOnAlgebra(h_sym2, Ar_sym); - SU3::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2; + SU::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2; std::cout << GridLogMessage << "Projections structure check vector difference (Two Index Symmetric): " << norm2(h_diff_sym) << std::endl; @@ -380,11 +380,11 @@ int main(int argc, char** argv) { // Construct the fundamental matrix in the group - SU3::LatticeMatrix Af_sym(grid); - SU3::FundamentalLieAlgebraMatrix(h_sym,Af_sym); - SU3::LatticeMatrix Ufund2(grid); + SU::LatticeMatrix Af_sym(grid); + SU::FundamentalLieAlgebraMatrix(h_sym,Af_sym); + SU::LatticeMatrix Ufund2(grid); Ufund2 = expMat(Af_sym, 1.0, 16); - SU3::LatticeMatrix UnitCheck2(grid); + SU::LatticeMatrix UnitCheck2(grid); UnitCheck2 = Ufund2 * adj(Ufund2) - uno_f; std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2) << std::endl; @@ -422,14 +422,14 @@ int main(int argc, char** argv) { // Test group structure // (U_f * V_f)_r = U_r * V_r LatticeGaugeField U2A(grid), V2A(grid); - SU3::HotConfiguration(gridRNG, U2A); - SU3::HotConfiguration(gridRNG, V2A); + SU::HotConfiguration(gridRNG, U2A); + SU::HotConfiguration(gridRNG, V2A); LatticeGaugeField UV2A(grid); UV2A = Zero(); for (int mu = 0; mu < Nd; mu++) { - SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu); - SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu); + SU::LatticeMatrix Umu2A = peekLorentz(U2,mu); + SU::LatticeMatrix Vmu2A = peekLorentz(V2,mu); pokeLorentz(UV2A,Umu2A*Vmu2A, mu); } @@ -456,16 +456,16 @@ int main(int argc, char** argv) { // Check correspondence of algebra and group transformations // Create a random vector - SU3::LatticeAlgebraVector h_Asym(grid); + SU::LatticeAlgebraVector h_Asym(grid); typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ar_Asym(grid); random(gridRNG,h_Asym); h_Asym = real(h_Asym); SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym); // Re-extract h_sym - SU3::LatticeAlgebraVector h_Asym2(grid); + SU::LatticeAlgebraVector h_Asym2(grid); SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym); - SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2; + SU::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2; std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl; @@ -490,11 +490,11 @@ int main(int argc, char** argv) { // Construct the fundamental matrix in the group - SU3::LatticeMatrix Af_Asym(grid); - SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym); - SU3::LatticeMatrix Ufund2A(grid); + SU::LatticeMatrix Af_Asym(grid); + SU::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym); + SU::LatticeMatrix Ufund2A(grid); Ufund2A = expMat(Af_Asym, 1.0, 16); - SU3::LatticeMatrix UnitCheck2A(grid); + SU::LatticeMatrix UnitCheck2A(grid); UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f; std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A) << std::endl; @@ -527,8 +527,7 @@ int main(int argc, char** argv) { - Grid_finalize(); } -#endif + diff --git a/tests/core/Test_reunitarise.cc b/tests/core/Test_reunitarise.cc index 0ff8ec7e..e421a2a8 100644 --- a/tests/core/Test_reunitarise.cc +++ b/tests/core/Test_reunitarise.cc @@ -26,7 +26,7 @@ Author: Peter Boyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#if Nc==3 + #include using namespace std; @@ -36,7 +36,7 @@ using namespace Grid; int main (int argc, char ** argv) { Grid_init(&argc,&argv); - + std::vector latt({8,8,8,8}); GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt, GridDefaultSimd(Nd,vComplexD::Nsimd()), @@ -123,14 +123,14 @@ int main (int argc, char ** argv) std::cout << "Determinant defect before projection " <::HotConfiguration(pRNG_f, Umu); RealD checkTolerance = (getPrecision::value == 1) ? 1e-7 : 1e-15; From 371fd123fb5ddf7da5f151e9b2fcb41829b0f7e0 Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Tue, 14 Mar 2023 10:47:07 +0000 Subject: [PATCH 5/9] consequence of iSUnMatrix being no longer a member of the SU class --- Grid/qcd/utils/SUnAdjoint.h | 8 +++++--- Grid/qcd/utils/SUnTwoIndex.h | 8 ++++---- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/Grid/qcd/utils/SUnAdjoint.h b/Grid/qcd/utils/SUnAdjoint.h index c03036ef..3a5954ca 100644 --- a/Grid/qcd/utils/SUnAdjoint.h +++ b/Grid/qcd/utils/SUnAdjoint.h @@ -51,14 +51,16 @@ public: typedef Lattice >, Nd> > LatticeAdjFieldF; typedef Lattice >, Nd> > LatticeAdjFieldD; + template + using iSUnMatrix = iScalar > >; template static void generator(int Index, iSUnAdjointMatrix &iAdjTa) { // returns i(T_Adj)^index necessary for the projectors // see definitions above iAdjTa = Zero(); - Vector::template iSUnMatrix > ta(ncolour * ncolour - 1); - typename SU::template iSUnMatrix tmp; + Vector > ta(ncolour * ncolour - 1); + iSUnMatrix tmp; // FIXME not very efficient to get all the generators everytime for (int a = 0; a < Dimension; a++) SU::generator(a, ta[a]); @@ -66,7 +68,7 @@ public: for (int a = 0; a < Dimension; a++) { tmp = ta[a] * ta[Index] - ta[Index] * ta[a]; for (int b = 0; b < (ncolour * ncolour - 1); b++) { - typename SU::template iSUnMatrix tmp1 = + iSUnMatrix tmp1 = 2.0 * tmp * ta[b]; // 2.0 from the normalization Complex iTr = TensorRemove(timesI(trace(tmp1))); //iAdjTa()()(b, a) = iTr; diff --git a/Grid/qcd/utils/SUnTwoIndex.h b/Grid/qcd/utils/SUnTwoIndex.h index 50ddf5c1..5e283de7 100644 --- a/Grid/qcd/utils/SUnTwoIndex.h +++ b/Grid/qcd/utils/SUnTwoIndex.h @@ -126,10 +126,10 @@ public: template static void generator(int Index, iSUnTwoIndexMatrix &i2indTa) { - Vector::template iSUnMatrix > ta( + Vector > ta( ncolour * ncolour - 1); - Vector::template iSUnMatrix > eij(Dimension); - typename SU::template iSUnMatrix tmp; + Vector > eij(Dimension); + iSUnMatrix tmp; i2indTa = Zero(); for (int a = 0; a < ncolour * ncolour - 1; a++) @@ -140,7 +140,7 @@ public: for (int a = 0; a < Dimension; a++) { tmp = transpose(ta[Index]) * adj(eij[a]) + adj(eij[a]) * ta[Index]; for (int b = 0; b < Dimension; b++) { - typename SU::template iSUnMatrix tmp1 = + iSUnMatrix tmp1 = tmp * eij[b]; Complex iTr = TensorRemove(timesI(trace(tmp1))); i2indTa()()(a, b) = iTr; From ba7f9d7b70321613497f1cbd68197d0f4df08552 Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Wed, 15 Mar 2023 15:55:12 +0000 Subject: [PATCH 6/9] projection on Sp2n algebra, to be used instead of Ta --- Grid/lattice/Lattice_ET.h | 2 + Grid/qcd/utils/Sp2n.impl | 29 +++--- Grid/tensors/Tensor_Ta.h | 84 ++++++++++++++++- tests/sp2n/Test_project_on_Sp.cc | 157 ++++++++++++++++++++++--------- 4 files changed, 214 insertions(+), 58 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 44771a14..0bad3d79 100644 --- a/Grid/qcd/utils/Sp2n.impl +++ b/Grid/qcd/utils/Sp2n.impl @@ -414,13 +414,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.; - // Omega()()(2*ncolour-1-i, i) = -1; - // } for (int i = 0; i < nsp; i++) { Omega()()(i, nsp + i) = 1.; Omega()()(nsp + i, i) = -1; @@ -454,8 +447,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); @@ -473,8 +464,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; @@ -496,3 +485,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..7ceeabba 100644 --- a/Grid/tensors/Tensor_Ta.h +++ b/Grid/tensors/Tensor_Ta.h @@ -66,6 +66,88 @@ template accelerator_inline iMatrix Ta(const iMatrix return ret; } +// for sp2n can't do something as simple as Ta + +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; + // following HiRep's suN_utils.c does a sort of Gram-Schmidt + 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;c2 ProjectOnGroup(const iMatrix &arg) // re-do for sp2n -// Ta cannot be defined here for Sp2n because I need the generators from the Sp class -// It is defined in gauge impl types template accelerator_inline iScalar ProjectOnSpGroup(const iScalar&r) { diff --git a/tests/sp2n/Test_project_on_Sp.cc b/tests/sp2n/Test_project_on_Sp.cc index 36b95cf4..320ab230 100644 --- a/tests/sp2n/Test_project_on_Sp.cc +++ b/tests/sp2n/Test_project_on_Sp.cc @@ -20,6 +20,8 @@ int main (int argc, char **argv) LatticeColourMatrixD aux(&Grid); LatticeColourMatrixD identity(&Grid); + // Will test resimplectification (from ProjectOnGaugeGroup, ProjectOnSpGroup) and projection on the algebra (from ProjectSp2nAlgebra) + const int nsp = Nc / 2; identity = 1.0; @@ -29,15 +31,12 @@ int main (int argc, char **argv) double vol = Umu.Grid()->gSites(); 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); // 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 ); + } + + } + + 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 << "ok" << std::endl; + + std::cout << GridLogMessage << "Now check projection on the algebra" << 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 << "Project on sp2n algebra" << 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 <(in, 1); OmegaInvariance(U); @@ -464,6 +473,8 @@ 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; @@ -485,21 +496,3 @@ 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 7ceeabba..2ee51fb4 100644 --- a/Grid/tensors/Tensor_Ta.h +++ b/Grid/tensors/Tensor_Ta.h @@ -66,88 +66,6 @@ template accelerator_inline iMatrix Ta(const iMatrix return ret; } -// for sp2n can't do something as simple as Ta - -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; - // following HiRep's suN_utils.c does a sort of Gram-Schmidt - 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;c2 ProjectOnGroup(const iMatrix &arg) // re-do for sp2n +// Ta cannot be defined here for Sp2n because I need the generators from the Sp class +// It is defined in gauge impl types template accelerator_inline iScalar ProjectOnSpGroup(const iScalar&r) { diff --git a/tests/sp2n/Test_project_on_Sp.cc b/tests/sp2n/Test_project_on_Sp.cc index 320ab230..36b95cf4 100644 --- a/tests/sp2n/Test_project_on_Sp.cc +++ b/tests/sp2n/Test_project_on_Sp.cc @@ -20,8 +20,6 @@ int main (int argc, char **argv) LatticeColourMatrixD aux(&Grid); LatticeColourMatrixD identity(&Grid); - // Will test resimplectification (from ProjectOnGaugeGroup, ProjectOnSpGroup) and projection on the algebra (from ProjectSp2nAlgebra) - const int nsp = Nc / 2; identity = 1.0; @@ -31,12 +29,15 @@ int main (int argc, char **argv) double vol = Umu.Grid()->gSites(); 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,0); + U = PeekIndex(Umu,2); aux = U*adj(U) - identity; + std::cout <::OmegaInvariance(U); // 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 ); - } - - } + assert(Nc==2); + ColourMatrix A; + A = Zero(); - 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 << "ok" << std::endl; - - std::cout << GridLogMessage << "Now check projection on the algebra" << 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 << "Project on sp2n algebra" << 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 ); - } - } - + Complex a(25041994., 12.); + Complex b(39., 0.22); + Complex d(10000., -2222.3333); + A()()(0,0) = a; + A()()(0,1) = b; + A()()(1,0) = i; + A()()(1,1) = d; + std::cout <