From 03235d63686dac3b369fd5208c2bb1ef50c0dbb1 Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Fri, 25 Nov 2022 16:57:40 +0000 Subject: [PATCH 01/17] Fixed type in configure.ac --- configure.ac | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/configure.ac b/configure.ac index 3a4ed469..04a3cdf4 100644 --- a/configure.ac +++ b/configure.ac @@ -838,7 +838,7 @@ AC_CONFIG_FILES(tests/lanczos/Makefile) AC_CONFIG_FILES(tests/smearing/Makefile) AC_CONFIG_FILES(tests/qdpxx/Makefile) AC_CONFIG_FILES(tests/testu01/Makefile) -AC_CONFIG_FILES(tests/Sp2n/Makefile) +AC_CONFIG_FILES(tests/sp2n/Makefile) AC_CONFIG_FILES(benchmarks/Makefile) AC_CONFIG_FILES(examples/Makefile) AC_OUTPUT From 629cb2987ad4293f4cc79aac713a0a6fd8c50fe2 Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Fri, 25 Nov 2022 17:40:21 +0000 Subject: [PATCH 02/17] Fix typo in Makefile.am --- tests/sp2n/Makefile.am | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/sp2n/Makefile.am b/tests/sp2n/Makefile.am index 0ee7421c..164fce52 100644 --- a/tests/sp2n/Makefile.am +++ b/tests/sp2n/Makefile.am @@ -4,4 +4,4 @@ include Make.inc check: tests ./Test_project_on_Sp - ./test_sp2n_lie_gen + ./Test_sp2n_lie_gen From 1aa28b47ae31c32e8a885c2d8c183ab231dc0a6f Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Fri, 25 Nov 2022 17:40:40 +0000 Subject: [PATCH 03/17] Add existing test to check --- tests/sp2n/Makefile.am | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/sp2n/Makefile.am b/tests/sp2n/Makefile.am index 164fce52..346c8ad0 100644 --- a/tests/sp2n/Makefile.am +++ b/tests/sp2n/Makefile.am @@ -5,3 +5,4 @@ include Make.inc check: tests ./Test_project_on_Sp ./Test_sp2n_lie_gen + ./Test_Sp_start From 9273f2937c3322c2fd33dfe06088a2d1566939ca Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Fri, 25 Nov 2022 17:44:08 +0000 Subject: [PATCH 04/17] Autoformat google style --- Grid/qcd/utils/SUn.h | 281 +++++++++++++++++++++---------------------- 1 file changed, 138 insertions(+), 143 deletions(-) diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index 675493b3..0f7482ae 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -28,7 +28,7 @@ with this program; if not, write to the Free Software Foundation, Inc., See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ - /* END LEGAL */ +/* END LEGAL */ #ifndef QCD_UTIL_SUN_H #define QCD_UTIL_SUN_H @@ -36,7 +36,7 @@ NAMESPACE_BEGIN(Grid); template class SU { -public: + public: static const int Dimension = ncolour; static const int AdjointDimension = ncolour * ncolour - 1; static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; } @@ -47,7 +47,7 @@ public: using iSU2Matrix = iScalar > >; template using iSUnAlgebraVector = - iScalar > >; + iScalar > >; ////////////////////////////////////////////////////////////////////////////////////////////////// // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix, @@ -159,7 +159,7 @@ public: else generatorSigmaX(su2Index, ta); } - + template static void generatorSigmaY(int su2Index, iSUnMatrix &ta) { ta = Zero(); @@ -169,7 +169,7 @@ public: ta()()(i2, i1) = 1.0; ta = ta * 0.5; } - + template static void generatorSigmaX(int su2Index, iSUnMatrix &ta) { ta = Zero(); @@ -194,8 +194,6 @@ public: ta = ta * nrm; } - - //////////////////////////////////////////////////////////////////////// // Map a su2 subgroup number to the pair of rows that are non zero //////////////////////////////////////////////////////////////////////// @@ -223,11 +221,10 @@ public: int i0, i1; su2SubGroupIndex(i0, i1, su2_index); - autoView( subgroup_v , subgroup,AcceleratorWrite); - autoView( source_v , source,AcceleratorRead); - autoView( Determinant_v , Determinant,AcceleratorWrite); + autoView(subgroup_v, subgroup, AcceleratorWrite); + autoView(source_v, source, AcceleratorRead); + autoView(Determinant_v, Determinant, AcceleratorWrite); accelerator_for(ss, grid->oSites(), 1, { - subgroup_v[ss]()()(0, 0) = source_v[ss]()()(i0, i0); subgroup_v[ss]()()(0, 1) = source_v[ss]()()(i0, i1); subgroup_v[ss]()()(1, 0) = source_v[ss]()()(i1, i0); @@ -241,7 +238,7 @@ public: // this should be purely real Determinant_v[ss] = - Sigma()()(0, 0) * Sigma()()(1, 1) - Sigma()()(0, 1) * Sigma()()(1, 0); + Sigma()()(0, 0) * Sigma()()(1, 1) - Sigma()()(0, 1) * Sigma()()(1, 0); }); } @@ -257,16 +254,14 @@ public: su2SubGroupIndex(i0, i1, su2_index); dest = 1.0; // start out with identity - autoView( dest_v , dest, AcceleratorWrite); - autoView( subgroup_v, subgroup, AcceleratorRead); - accelerator_for(ss, grid->oSites(),1, - { + autoView(dest_v, dest, AcceleratorWrite); + autoView(subgroup_v, subgroup, AcceleratorRead); + accelerator_for(ss, grid->oSites(), 1, { dest_v[ss]()()(i0, i0) = subgroup_v[ss]()()(0, 0); dest_v[ss]()()(i0, i1) = subgroup_v[ss]()()(0, 1); dest_v[ss]()()(i1, i0) = subgroup_v[ss]()()(1, 0); dest_v[ss]()()(i1, i1) = subgroup_v[ss]()()(1, 1); }); - } /////////////////////////////////////////////// @@ -279,12 +274,12 @@ public: // in action. // /////////////////////////////////////////////// - static void SubGroupHeatBath(GridSerialRNG &sRNG, GridParallelRNG &pRNG, - RealD beta, // coeff multiplying staple in action (with no 1/Nc) - LatticeMatrix &link, - const LatticeMatrix &barestaple, // multiplied by action coeffs so th - int su2_subgroup, int nheatbath, LatticeInteger &wheremask) - { + static void SubGroupHeatBath( + GridSerialRNG &sRNG, GridParallelRNG &pRNG, + RealD beta, // coeff multiplying staple in action (with no 1/Nc) + LatticeMatrix &link, + const LatticeMatrix &barestaple, // multiplied by action coeffs so th + int su2_subgroup, int nheatbath, LatticeInteger &wheremask) { GridBase *grid = link.Grid(); const RealD twopi = 2.0 * M_PI; @@ -297,7 +292,8 @@ public: V = link * staple; // Subgroup manipulation in the lie algebra space - LatticeSU2Matrix u(grid); // Kennedy pendleton "u" real projected normalised Sigma + LatticeSU2Matrix u( + grid); // Kennedy pendleton "u" real projected normalised Sigma LatticeSU2Matrix uinv(grid); LatticeSU2Matrix ua(grid); // a in pauli form LatticeSU2Matrix b(grid); // rotated matrix after hb @@ -390,7 +386,7 @@ public: xi = 0.5 * sqrt(udet); // 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] u = 0.5 * u * - pow(xi, -1.0); // u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] + pow(xi, -1.0); // u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] // Debug test for sanity uinv = adj(u); @@ -405,29 +401,24 @@ public: r) ) = da0 r/2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r) ) - Action factor Q(h) dh = e^-S[h] dh = e^{ xi Tr uh} dh // beta enters - through xi - = e^{2 xi (h.u)} dh - = e^{2 xi h0u0}.e^{2 xi h1u1}.e^{2 xi - h2u2}.e^{2 xi h3u3} dh + Action factor Q(h) dh = e^-S[h] dh = e^{ xi Tr uh} dh // beta + enters through xi = e^{2 xi (h.u)} dh = e^{2 xi h0u0}.e^{2 xi h1u1}.e^{2 + xi h2u2}.e^{2 xi h3u3} dh Therefore for each site, take xi for that site i) generate |a0|<1 with dist (1-a0^2)^0.5 e^{2 xi a0 } da0 - Take alpha = 2 xi = 2 xi [ recall 2 beta/Nc unmod staple norm]; hence 2.0/Nc - factor in Chroma ] - A. Generate two uniformly distributed pseudo-random numbers R and R', R'', - R''' in the unit interval; - B. Set X = -(ln R)/alpha, X' =-(ln R')/alpha; - C. Set C = cos^2(2pi R"), with R" another uniform random number in [0,1] ; - D. Set A = XC; - E. Let d = X'+A; + Take alpha = 2 xi = 2 xi [ recall 2 beta/Nc unmod staple norm]; + hence 2.0/Nc factor in Chroma ] A. Generate two uniformly distributed + pseudo-random numbers R and R', R'', R''' in the unit interval; B. Set X = + -(ln R)/alpha, X' =-(ln R')/alpha; C. Set C = cos^2(2pi R"), with R" + another uniform random number in [0,1] ; D. Set A = XC; E. Let d = X'+A; F. If R'''^2 :> 1 - 0.5 d, go back to A; G. Set a0 = 1 - d; - Note that in step D setting B ~ X - A and using B in place of A in step E will - generate a second independent a 0 value. + Note that in step D setting B ~ X - A and using B in place of A in step E + will generate a second independent a 0 value. */ ///////////////////////////////////////////////////////// @@ -449,7 +440,7 @@ public: LatticeReal alpha(grid); // std::cout< - static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, double scale = 1.0) - { + static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, + double scale = 1.0) { GridBase *grid = out.Grid(); typedef typename LatticeMatrixType::vector_type vector_type; @@ -620,7 +609,8 @@ public: typedef iSinglet vTComplexType; typedef Lattice LatticeComplexType; - typedef typename GridTypeMapper::scalar_object MatrixType; + typedef typename GridTypeMapper< + typename LatticeMatrixType::vector_object>::scalar_object MatrixType; LatticeComplexType ca(grid); LatticeMatrixType lie(grid); @@ -642,7 +632,6 @@ public: la = ci * ca * ta; lie = lie + la; // e^{i la ta} - } taExp(lie, out); } @@ -681,57 +670,61 @@ public: out += la; } } -/* - * Fundamental rep gauge xform - */ - template - static void GaugeTransformFundamental( Fundamental &ferm, GaugeMat &g){ + /* + * Fundamental rep gauge xform + */ + template + static void GaugeTransformFundamental(Fundamental &ferm, GaugeMat &g) { GridBase *grid = ferm._grid; - conformable(grid,g._grid); - ferm = g*ferm; + conformable(grid, g._grid); + ferm = g * ferm; } -/* - * Adjoint rep gauge xform - */ + /* + * Adjoint rep gauge xform + */ - template - static void GaugeTransform( GaugeField &Umu, GaugeMat &g){ + template + static void GaugeTransform(GaugeField &Umu, GaugeMat &g) { GridBase *grid = Umu.Grid(); - conformable(grid,g.Grid()); + conformable(grid, g.Grid()); GaugeMat U(grid); - GaugeMat ag(grid); ag = adj(g); + GaugeMat ag(grid); + ag = adj(g); - for(int mu=0;mu(Umu,mu); - U = g*U*Cshift(ag, mu, 1); - PokeIndex(Umu,U,mu); + for (int mu = 0; mu < Nd; mu++) { + U = PeekIndex(Umu, mu); + U = g * U * Cshift(ag, mu, 1); + PokeIndex(Umu, U, mu); } } - template - static void GaugeTransform( std::vector &U, GaugeMat &g){ + template + static void GaugeTransform(std::vector &U, GaugeMat &g) { GridBase *grid = g.Grid(); - GaugeMat ag(grid); ag = adj(g); - for(int mu=0;mu - static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, GaugeMat &g){ - LieRandomize(pRNG,g,1.0); - GaugeTransform(Umu,g); + template + static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, + GaugeMat &g) { + LieRandomize(pRNG, g, 1.0); + GaugeTransform(Umu, g); } - // 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) { + // 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); + pokeColour(h_out, -2.0 * (trace(timesI(Ta) * in)) * scale, a); } } @@ -747,37 +740,37 @@ public: PokeIndex(out, Umu, mu); } } - template - static void TepidConfiguration(GridParallelRNG &pRNG,GaugeField &out){ + template + static void TepidConfiguration(GridParallelRNG &pRNG, GaugeField &out) { typedef typename GaugeField::vector_type vector_type; typedef iSUnMatrix vMatrixType; typedef Lattice LatticeMatrixType; LatticeMatrixType Umu(out.Grid()); - for(int mu=0;mu(out,Umu,mu); + for (int mu = 0; mu < Nd; mu++) { + LieRandomize(pRNG, Umu, 0.01); + PokeIndex(out, Umu, mu); } } - template - static void ColdConfiguration(GaugeField &out){ + template + static void ColdConfiguration(GaugeField &out) { typedef typename GaugeField::vector_type vector_type; typedef iSUnMatrix vMatrixType; typedef Lattice LatticeMatrixType; LatticeMatrixType Umu(out.Grid()); - Umu=1.0; - for(int mu=0;mu(out,Umu,mu); + Umu = 1.0; + for (int mu = 0; mu < Nd; mu++) { + PokeIndex(out, Umu, mu); } } - template - static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){ + template + static void ColdConfiguration(GridParallelRNG &pRNG, GaugeField &out) { ColdConfiguration(out); } - template - static void taProj( const LatticeMatrixType &in, LatticeMatrixType &out){ + template + static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) { out = Ta(in); } template @@ -799,85 +792,88 @@ public: } }; -template -LatticeComplexD Determinant(const Lattice > > > &Umu) -{ - GridBase *grid=Umu.Grid(); +template +LatticeComplexD Determinant( + const Lattice > > > &Umu) { + GridBase *grid = Umu.Grid(); auto lvol = grid->lSites(); LatticeComplexD ret(grid); - autoView(Umu_v,Umu,CpuRead); - autoView(ret_v,ret,CpuWrite); - thread_for(site,lvol,{ - Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N,N); + autoView(Umu_v, Umu, CpuRead); + autoView(ret_v, ret, CpuWrite); + thread_for(site, lvol, { + Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N, N); Coordinate lcoor; grid->LocalIndexToLocalCoor(site, lcoor); iScalar > > Us; peekLocalSite(Us, Umu_v, lcoor); - for(int i=0;i -static void ProjectSUn(Lattice > > > &Umu) -{ - Umu = ProjectOnGroup(Umu); +template +static void ProjectSUn( + Lattice > > > &Umu) { + Umu = ProjectOnGroup(Umu); auto det = Determinant(Umu); det = conjugate(det); - for(int i=0;i(Umu,N-1,i); + for (int i = 0; i < N; i++) { + auto element = PeekIndex(Umu, N - 1, i); element = element * det; - PokeIndex(Umu,element,Nc-1,i); + PokeIndex(Umu, element, Nc - 1, i); } } -template -static void ProjectSUn(Lattice >,Nd> > &U) -{ - GridBase *grid=U.Grid(); +template +static void ProjectSUn( + Lattice >, Nd> > &U) { + GridBase *grid = U.Grid(); // Reunitarise - for(int mu=0;mu(U,mu); - Umu = ProjectOnGroup(Umu); + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnGroup(Umu); ProjectSUn(Umu); - PokeIndex(U,Umu,mu); + PokeIndex(U, Umu, mu); } } // Explicit specialisation for SU(3). // Explicit specialisation for SU(3). -static void -ProjectSU3 (Lattice > > > &Umu) -{ - GridBase *grid=Umu.Grid(); - const int x=0; - const int y=1; - const int z=2; +static void ProjectSU3( + Lattice > > > &Umu) { + GridBase *grid = Umu.Grid(); + const int x = 0; + const int y = 1; + const int z = 2; // Reunitarise Umu = ProjectOnGroup(Umu); - autoView(Umu_v,Umu,CpuWrite); - thread_for(ss,grid->oSites(),{ - auto cm = Umu_v[ss]; - cm()()(2,x) = adj(cm()()(0,y)*cm()()(1,z)-cm()()(0,z)*cm()()(1,y)); //x= yz-zy - cm()()(2,y) = adj(cm()()(0,z)*cm()()(1,x)-cm()()(0,x)*cm()()(1,z)); //y= zx-xz - cm()()(2,z) = adj(cm()()(0,x)*cm()()(1,y)-cm()()(0,y)*cm()()(1,x)); //z= xy-yx - Umu_v[ss]=cm; + autoView(Umu_v, Umu, CpuWrite); + thread_for(ss, grid->oSites(), { + auto cm = Umu_v[ss]; + cm()()(2, x) = adj(cm()()(0, y) * cm()()(1, z) - + cm()()(0, z) * cm()()(1, y)); // x= yz-zy + cm()()(2, y) = adj(cm()()(0, z) * cm()()(1, x) - + cm()()(0, x) * cm()()(1, z)); // y= zx-xz + cm()()(2, z) = adj(cm()()(0, x) * cm()()(1, y) - + cm()()(0, y) * cm()()(1, x)); // z= xy-yx + Umu_v[ss] = cm; }); } -static void ProjectSU3(Lattice >,Nd> > &U) -{ - GridBase *grid=U.Grid(); +static void ProjectSU3( + Lattice >, Nd> > &U) { + GridBase *grid = U.Grid(); // Reunitarise - for(int mu=0;mu(U,mu); - Umu = ProjectOnGroup(Umu); + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnGroup(Umu); ProjectSU3(Umu); - PokeIndex(U,Umu,mu); + PokeIndex(U, Umu, mu); } } @@ -886,7 +882,6 @@ typedef SU<3> SU3; typedef SU<4> SU4; typedef SU<5> SU5; - typedef SU FundamentalMatrices; NAMESPACE_END(Grid); From b8f1f5d2a3e164007b6901aa40b0cdc1d42bcf56 Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Fri, 25 Nov 2022 17:45:32 +0000 Subject: [PATCH 05/17] Introduce GaugeGroup --- Grid/qcd/utils/SUn.h | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index 0f7482ae..50a035d4 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -33,9 +33,18 @@ directory #define QCD_UTIL_SUN_H NAMESPACE_BEGIN(Grid); +namespace GroupName { +class SU {}; +} // namespace GroupName + +template +class GaugeGroup; template -class SU { +using SU = GaugeGroup; + +template +class GaugeGroup { public: static const int Dimension = ncolour; static const int AdjointDimension = ncolour * ncolour - 1; From 6e750ecb0edc01e633be65f8ddc4bb7cf4932a1d Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Mon, 28 Nov 2022 16:33:46 +0000 Subject: [PATCH 06/17] Remove apparently forgotten file --- tests/core/test_sp2n_lie_gen.cc | 55 --------------------------------- 1 file changed, 55 deletions(-) delete mode 100644 tests/core/test_sp2n_lie_gen.cc diff --git a/tests/core/test_sp2n_lie_gen.cc b/tests/core/test_sp2n_lie_gen.cc deleted file mode 100644 index 6c3766e6..00000000 --- a/tests/core/test_sp2n_lie_gen.cc +++ /dev/null @@ -1,55 +0,0 @@ - - - -#include -#include - -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()); - //GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid); - - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "* Generators for Sp(2)" << std::endl; - std::cout << GridLogMessage << "*********************************************" - << std::endl; - - Sp2::printGenerators(); - Sp2::testGenerators(); - - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "* Generators for Sp(4)" << std::endl; - std::cout << GridLogMessage << "*********************************************" - << std::endl; - - Sp4::printGenerators(); - Sp4::testGenerators(); - - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "* Generators for Sp(6)" << std::endl; - std::cout << GridLogMessage << "*********************************************" - << std::endl; - - Sp6::printGenerators(); - Sp6::testGenerators(); - - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "* Generators for Sp(8)" << std::endl; - std::cout << GridLogMessage << "*********************************************" - << std::endl; - - Sp8::printGenerators(); - Sp8::testGenerators(); - - - Grid_finalize(); -} From 921e23e83ce252044932e7d0455e4dbf904b716c Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Mon, 28 Nov 2022 17:47:50 +0000 Subject: [PATCH 07/17] Separated out everything SU specific --- Grid/qcd/utils/SUn.h | 591 +++----------------------------------- Grid/qcd/utils/SUn_impl.h | 552 +++++++++++++++++++++++++++++++++++ 2 files changed, 586 insertions(+), 557 deletions(-) create mode 100644 Grid/qcd/utils/SUn_impl.h diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index 50a035d4..4a7148e9 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -32,11 +32,25 @@ directory #ifndef QCD_UTIL_SUN_H #define QCD_UTIL_SUN_H +#define ONLY_IF_SU \ + typename dummy_name = group_name, \ + typename = std::enable_if_t::value> + NAMESPACE_BEGIN(Grid); namespace GroupName { class SU {}; } // namespace GroupName +template +struct is_su { + static const bool value = false; +}; + +template <> +struct is_su { + static const bool value = true; +}; + template class GaugeGroup; @@ -50,36 +64,35 @@ class GaugeGroup { static const int AdjointDimension = ncolour * ncolour - 1; static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; } - template - using iSUnMatrix = iScalar > >; template using iSU2Matrix = iScalar > >; template - using iSUnAlgebraVector = - iScalar > >; + using iGroupMatrix = iScalar > >; + template + using iAlgebraVector = iScalar > >; ////////////////////////////////////////////////////////////////////////////////////////////////// // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix, // SU<2>::LatticeMatrix etc... ////////////////////////////////////////////////////////////////////////////////////////////////// - typedef iSUnMatrix Matrix; - typedef iSUnMatrix MatrixF; - typedef iSUnMatrix MatrixD; + typedef iGroupMatrix Matrix; + typedef iGroupMatrix MatrixF; + typedef iGroupMatrix MatrixD; - typedef iSUnMatrix vMatrix; - typedef iSUnMatrix vMatrixF; - typedef iSUnMatrix vMatrixD; + typedef iGroupMatrix vMatrix; + typedef iGroupMatrix vMatrixF; + typedef iGroupMatrix vMatrixD; // For the projectors to the algebra // these should be real... // keeping complex for consistency with the SIMD vector types - typedef iSUnAlgebraVector AlgebraVector; - typedef iSUnAlgebraVector AlgebraVectorF; - typedef iSUnAlgebraVector AlgebraVectorD; + typedef iAlgebraVector AlgebraVector; + typedef iAlgebraVector AlgebraVectorF; + typedef iAlgebraVector AlgebraVectorD; - typedef iSUnAlgebraVector vAlgebraVector; - typedef iSUnAlgebraVector vAlgebraVectorF; - typedef iSUnAlgebraVector vAlgebraVectorD; + typedef iAlgebraVector vAlgebraVector; + typedef iAlgebraVector vAlgebraVectorF; + typedef iAlgebraVector vAlgebraVectorD; typedef Lattice LatticeMatrix; typedef Lattice LatticeMatrixF; @@ -101,462 +114,7 @@ class GaugeGroup { typedef Lattice LatticeSU2MatrixF; typedef Lattice LatticeSU2MatrixD; - //////////////////////////////////////////////////////////////////////// - // There are N^2-1 generators for SU(N). - // - // We take a traceless hermitian generator basis as follows - // - // * Normalisation: trace ta tb = 1/2 delta_ab = T_F delta_ab - // T_F = 1/2 for SU(N) groups - // - // * Off diagonal - // - pairs of rows i1,i2 behaving like pauli matrices signma_x, sigma_y - // - // - there are (Nc-1-i1) slots for i2 on each row [ x 0 x ] - // direct count off each row - // - // - Sum of all pairs is Nc(Nc-1)/2: proof arithmetic series - // - // (Nc-1) + (Nc-2)+... 1 ==> Nc*(Nc-1)/2 - // 1+ 2+ + + Nc-1 - // - // - There are 2 x Nc (Nc-1)/ 2 of these = Nc^2 - Nc - // - // - We enumerate the row-col pairs. - // - for each row col pair there is a (sigma_x) and a (sigma_y) like - // generator - // - // - // t^a_ij = { in 0.. Nc(Nc-1)/2 -1} => 1/2(delta_{i,i1} delta_{j,i2} + - // delta_{i,i1} delta_{j,i2}) - // t^a_ij = { in Nc(Nc-1)/2 ... Nc(Nc-1) - 1} => i/2( delta_{i,i1} - // delta_{j,i2} - i delta_{i,i1} delta_{j,i2}) - // - // * Diagonal; must be traceless and normalised - // - Sequence is - // N (1,-1,0,0...) - // N (1, 1,-2,0...) - // N (1, 1, 1,-3,0...) - // N (1, 1, 1, 1,-4,0...) - // - // where 1/2 = N^2 (1+.. m^2)etc.... for the m-th diagonal generator - // NB this gives the famous SU3 result for su2 index 8 - // - // N= sqrt(1/2 . 1/6 ) = 1/2 . 1/sqrt(3) - // - // ( 1 ) - // ( 1 ) / sqrt(3) /2 = 1/2 lambda_8 - // ( -2) - // - //////////////////////////////////////////////////////////////////////// - template - static void generator(int lieIndex, iSUnMatrix &ta) { - // map lie index to which type of generator - int diagIndex; - int su2Index; - int sigxy; - int NNm1 = ncolour * (ncolour - 1); - if (lieIndex >= NNm1) { - diagIndex = lieIndex - NNm1; - generatorDiagonal(diagIndex, ta); - return; - } - sigxy = lieIndex & 0x1; // even or odd - su2Index = lieIndex >> 1; - if (sigxy) - generatorSigmaY(su2Index, ta); - else - generatorSigmaX(su2Index, ta); - } - - template - static void generatorSigmaY(int su2Index, iSUnMatrix &ta) { - ta = Zero(); - int i1, i2; - su2SubGroupIndex(i1, i2, su2Index); - ta()()(i1, i2) = 1.0; - ta()()(i2, i1) = 1.0; - ta = ta * 0.5; - } - - template - static void generatorSigmaX(int su2Index, iSUnMatrix &ta) { - ta = Zero(); - cplx i(0.0, 1.0); - int i1, i2; - su2SubGroupIndex(i1, i2, su2Index); - ta()()(i1, i2) = i; - ta()()(i2, i1) = -i; - ta = ta * 0.5; - } - - template - static void generatorDiagonal(int diagIndex, iSUnMatrix &ta) { - // diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...) - ta = Zero(); - int k = diagIndex + 1; // diagIndex starts from 0 - for (int i = 0; i <= diagIndex; i++) { // k iterations - ta()()(i, i) = 1.0; - } - ta()()(k, k) = -k; // indexing starts from 0 - RealD nrm = 1.0 / std::sqrt(2.0 * k * (k + 1)); - ta = ta * nrm; - } - - //////////////////////////////////////////////////////////////////////// - // Map a su2 subgroup number to the pair of rows that are non zero - //////////////////////////////////////////////////////////////////////// - static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { - assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2)); - - int spare = su2_index; - for (i1 = 0; spare >= (ncolour - 1 - i1); i1++) { - spare = spare - (ncolour - 1 - i1); // remove the Nc-1-i1 terms - } - i2 = i1 + 1 + spare; - } - - ////////////////////////////////////////////////////////////////////////////////////////// - // Pull out a subgroup and project on to real coeffs x pauli basis - ////////////////////////////////////////////////////////////////////////////////////////// - template - static void su2Extract(Lattice > &Determinant, - Lattice > &subgroup, - const Lattice > &source, - int su2_index) { - GridBase *grid(source.Grid()); - conformable(subgroup, source); - conformable(subgroup, Determinant); - int i0, i1; - su2SubGroupIndex(i0, i1, su2_index); - - autoView(subgroup_v, subgroup, AcceleratorWrite); - autoView(source_v, source, AcceleratorRead); - autoView(Determinant_v, Determinant, AcceleratorWrite); - accelerator_for(ss, grid->oSites(), 1, { - subgroup_v[ss]()()(0, 0) = source_v[ss]()()(i0, i0); - subgroup_v[ss]()()(0, 1) = source_v[ss]()()(i0, i1); - subgroup_v[ss]()()(1, 0) = source_v[ss]()()(i1, i0); - subgroup_v[ss]()()(1, 1) = source_v[ss]()()(i1, i1); - - iSU2Matrix Sigma = subgroup_v[ss]; - - Sigma = Sigma - adj(Sigma) + trace(adj(Sigma)); - - subgroup_v[ss] = Sigma; - - // this should be purely real - Determinant_v[ss] = - Sigma()()(0, 0) * Sigma()()(1, 1) - Sigma()()(0, 1) * Sigma()()(1, 0); - }); - } - - ////////////////////////////////////////////////////////////////////////////////////////// - // Set matrix to one and insert a pauli subgroup - ////////////////////////////////////////////////////////////////////////////////////////// - template - static void su2Insert(const Lattice > &subgroup, - Lattice > &dest, int su2_index) { - GridBase *grid(dest.Grid()); - conformable(subgroup, dest); - int i0, i1; - su2SubGroupIndex(i0, i1, su2_index); - - dest = 1.0; // start out with identity - autoView(dest_v, dest, AcceleratorWrite); - autoView(subgroup_v, subgroup, AcceleratorRead); - accelerator_for(ss, grid->oSites(), 1, { - dest_v[ss]()()(i0, i0) = subgroup_v[ss]()()(0, 0); - dest_v[ss]()()(i0, i1) = subgroup_v[ss]()()(0, 1); - dest_v[ss]()()(i1, i0) = subgroup_v[ss]()()(1, 0); - dest_v[ss]()()(i1, i1) = subgroup_v[ss]()()(1, 1); - }); - } - - /////////////////////////////////////////////// - // Generate e^{ Re Tr Staple Link} dlink - // - // *** Note Staple should be appropriate linear compbination between all - // staples. - // *** If already by beta pass coefficient 1.0. - // *** This routine applies the additional 1/Nc factor that comes after trace - // in action. - // - /////////////////////////////////////////////// - static void SubGroupHeatBath( - GridSerialRNG &sRNG, GridParallelRNG &pRNG, - RealD beta, // coeff multiplying staple in action (with no 1/Nc) - LatticeMatrix &link, - const LatticeMatrix &barestaple, // multiplied by action coeffs so th - int su2_subgroup, int nheatbath, LatticeInteger &wheremask) { - GridBase *grid = link.Grid(); - - const RealD twopi = 2.0 * M_PI; - - LatticeMatrix staple(grid); - - staple = barestaple * (beta / ncolour); - - LatticeMatrix V(grid); - V = link * staple; - - // Subgroup manipulation in the lie algebra space - LatticeSU2Matrix u( - grid); // Kennedy pendleton "u" real projected normalised Sigma - LatticeSU2Matrix uinv(grid); - LatticeSU2Matrix ua(grid); // a in pauli form - LatticeSU2Matrix b(grid); // rotated matrix after hb - - // Some handy constant fields - LatticeComplex ones(grid); - ones = 1.0; - LatticeComplex zeros(grid); - zeros = Zero(); - LatticeReal rones(grid); - rones = 1.0; - LatticeReal rzeros(grid); - rzeros = Zero(); - LatticeComplex udet(grid); // determinant of real(staple) - LatticeInteger mask_true(grid); - mask_true = 1; - LatticeInteger mask_false(grid); - mask_false = 0; - - /* - PLB 156 P393 (1985) (Kennedy and Pendleton) - - Note: absorb "beta" into the def of sigma compared to KP paper; staple - passed to this routine has "beta" already multiplied in - - Action linear in links h and of form: - - beta S = beta Sum_p (1 - 1/Nc Re Tr Plaq ) - - Writing Sigma = 1/Nc (beta Sigma') where sum over staples is "Sigma' " - - beta S = const - beta/Nc Re Tr h Sigma' - = const - Re Tr h Sigma - - Decompose h and Sigma into (1, sigma_j) ; h_i real, h^2=1, Sigma_i complex - arbitrary. - - Tr h Sigma = h_i Sigma_j Tr (sigma_i sigma_j) = h_i Sigma_j 2 delta_ij - Re Tr h Sigma = 2 h_j Re Sigma_j - - Normalised re Sigma_j = xi u_j - - With u_j a unit vector and U can be in SU(2); - - Re Tr h Sigma = 2 h_j Re Sigma_j = 2 xi (h.u) - - 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] - u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] - - xi = sqrt(Det)/2; - - Write a= u h in SU(2); a has pauli decomp a_j; - - Note: Product b' xi is unvariant because scaling Sigma leaves - normalised vector "u" fixed; Can rescale Sigma so b' = 1. - */ - - //////////////////////////////////////////////////////// - // Real part of Pauli decomposition - // Note a subgroup can project to zero in cold start - //////////////////////////////////////////////////////// - su2Extract(udet, u, V, su2_subgroup); - - ////////////////////////////////////////////////////// - // Normalising this vector if possible; else identity - ////////////////////////////////////////////////////// - LatticeComplex xi(grid); - - LatticeSU2Matrix lident(grid); - - SU2Matrix ident = Complex(1.0); - SU2Matrix pauli1; - SU<2>::generator(0, pauli1); - SU2Matrix pauli2; - SU<2>::generator(1, pauli2); - SU2Matrix pauli3; - SU<2>::generator(2, pauli3); - pauli1 = timesI(pauli1) * 2.0; - pauli2 = timesI(pauli2) * 2.0; - pauli3 = timesI(pauli3) * 2.0; - - LatticeComplex cone(grid); - LatticeReal adet(grid); - adet = abs(toReal(udet)); - lident = Complex(1.0); - cone = Complex(1.0); - Real machine_epsilon = 1.0e-7; - u = where(adet > machine_epsilon, u, lident); - udet = where(adet > machine_epsilon, udet, cone); - - xi = 0.5 * sqrt(udet); // 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] - u = 0.5 * u * - pow(xi, -1.0); // u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] - - // Debug test for sanity - uinv = adj(u); - b = u * uinv - 1.0; - assert(norm2(b) < 1.0e-4); - - /* - Measure: Haar measure dh has d^4a delta(1-|a^2|) - In polars: - da = da0 r^2 sin theta dr dtheta dphi delta( 1 - r^2 -a0^2) - = da0 r^2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r)(sqrt(1-a0^) + - r) ) - = da0 r/2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r) ) - - Action factor Q(h) dh = e^-S[h] dh = e^{ xi Tr uh} dh // beta - enters through xi = e^{2 xi (h.u)} dh = e^{2 xi h0u0}.e^{2 xi h1u1}.e^{2 - xi h2u2}.e^{2 xi h3u3} dh - - Therefore for each site, take xi for that site - i) generate |a0|<1 with dist - (1-a0^2)^0.5 e^{2 xi a0 } da0 - - Take alpha = 2 xi = 2 xi [ recall 2 beta/Nc unmod staple norm]; - hence 2.0/Nc factor in Chroma ] A. Generate two uniformly distributed - pseudo-random numbers R and R', R'', R''' in the unit interval; B. Set X = - -(ln R)/alpha, X' =-(ln R')/alpha; C. Set C = cos^2(2pi R"), with R" - another uniform random number in [0,1] ; D. Set A = XC; E. Let d = X'+A; - F. If R'''^2 :> 1 - 0.5 d, go back to A; - G. Set a0 = 1 - d; - - Note that in step D setting B ~ X - A and using B in place of A in step E - will generate a second independent a 0 value. - */ - - ///////////////////////////////////////////////////////// - // count the number of sites by picking "1"'s out of hat - ///////////////////////////////////////////////////////// - Integer hit = 0; - LatticeReal rtmp(grid); - rtmp = where(wheremask, rones, rzeros); - RealD numSites = sum(rtmp); - RealD numAccepted; - LatticeInteger Accepted(grid); - Accepted = Zero(); - LatticeInteger newlyAccepted(grid); - - std::vector xr(4, grid); - std::vector a(4, grid); - LatticeReal d(grid); - d = Zero(); - LatticeReal alpha(grid); - - // std::cout< 1 - 0.5 d, go back to A; - LatticeReal thresh(grid); - thresh = 1.0 - d * 0.5; - xrsq = xr[0] * xr[0]; - LatticeInteger ione(grid); - ione = 1; - LatticeInteger izero(grid); - izero = Zero(); - - newlyAccepted = where(xrsq < thresh, ione, izero); - Accepted = where(newlyAccepted, newlyAccepted, Accepted); - Accepted = where(wheremask, Accepted, izero); - - // FIXME need an iSum for integer to avoid overload on return type?? - rtmp = where(Accepted, rones, rzeros); - numAccepted = sum(rtmp); - - hit++; - - } while ((numAccepted < numSites) && (hit < nheatbath)); - - // G. Set a0 = 1 - d; - a[0] = Zero(); - a[0] = where(wheremask, 1.0 - d, a[0]); - - ////////////////////////////////////////// - // ii) generate a_i uniform on two sphere radius (1-a0^2)^0.5 - ////////////////////////////////////////// - - LatticeReal a123mag(grid); - a123mag = sqrt(abs(1.0 - a[0] * a[0])); - - LatticeReal cos_theta(grid); - LatticeReal sin_theta(grid); - LatticeReal phi(grid); - - random(pRNG, phi); - phi = phi * twopi; // uniform in [0,2pi] - random(pRNG, cos_theta); - cos_theta = (cos_theta * 2.0) - 1.0; // uniform in [-1,1] - sin_theta = sqrt(abs(1.0 - cos_theta * cos_theta)); - - a[1] = a123mag * sin_theta * cos(phi); - a[2] = a123mag * sin_theta * sin(phi); - a[3] = a123mag * cos_theta; - - ua = toComplex(a[0]) * ident + toComplex(a[1]) * pauli1 + - toComplex(a[2]) * pauli2 + toComplex(a[3]) * pauli3; - - b = 1.0; - b = where(wheremask, uinv * ua, b); - su2Insert(b, V, su2_subgroup); - - // mask the assignment back based on Accptance - link = where(Accepted, V * link, link); - - ////////////////////////////// - // Debug Checks - // SU2 check - LatticeSU2Matrix check(grid); // rotated matrix after hb - u = Zero(); - check = ua * adj(ua) - 1.0; - check = where(Accepted, check, u); - assert(norm2(check) < 1.0e-4); - - check = b * adj(b) - 1.0; - check = where(Accepted, check, u); - assert(norm2(check) < 1.0e-4); - - LatticeMatrix Vcheck(grid); - Vcheck = Zero(); - Vcheck = where(Accepted, V * adj(V) - 1.0, Vcheck); - // std::cout< static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, @@ -679,49 +199,6 @@ class GaugeGroup { out += la; } } - /* - * Fundamental rep gauge xform - */ - template - static void GaugeTransformFundamental(Fundamental &ferm, GaugeMat &g) { - GridBase *grid = ferm._grid; - conformable(grid, g._grid); - ferm = g * ferm; - } - /* - * Adjoint rep gauge xform - */ - - template - static void GaugeTransform(GaugeField &Umu, GaugeMat &g) { - GridBase *grid = Umu.Grid(); - conformable(grid, g.Grid()); - - GaugeMat U(grid); - GaugeMat ag(grid); - ag = adj(g); - - for (int mu = 0; mu < Nd; mu++) { - U = PeekIndex(Umu, mu); - U = g * U * Cshift(ag, mu, 1); - PokeIndex(Umu, U, mu); - } - } - template - static void GaugeTransform(std::vector &U, GaugeMat &g) { - GridBase *grid = g.Grid(); - GaugeMat ag(grid); - ag = adj(g); - for (int mu = 0; mu < Nd; mu++) { - U[mu] = g * U[mu] * Cshift(ag, mu, 1); - } - } - template - static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, - GaugeMat &g) { - LieRandomize(pRNG, g, 1.0); - GaugeTransform(Umu, g); - } // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 // ) inverse operation: FundamentalLieAlgebraMatrix @@ -740,7 +217,7 @@ class GaugeGroup { template static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) { typedef typename GaugeField::vector_type vector_type; - typedef iSUnMatrix vMatrixType; + typedef iGroupMatrix vMatrixType; typedef Lattice LatticeMatrixType; LatticeMatrixType Umu(out.Grid()); @@ -752,7 +229,7 @@ class GaugeGroup { template static void TepidConfiguration(GridParallelRNG &pRNG, GaugeField &out) { typedef typename GaugeField::vector_type vector_type; - typedef iSUnMatrix vMatrixType; + typedef iGroupMatrix vMatrixType; typedef Lattice LatticeMatrixType; LatticeMatrixType Umu(out.Grid()); @@ -764,7 +241,7 @@ class GaugeGroup { template static void ColdConfiguration(GaugeField &out) { typedef typename GaugeField::vector_type vector_type; - typedef iSUnMatrix vMatrixType; + typedef iGroupMatrix vMatrixType; typedef Lattice LatticeMatrixType; LatticeMatrixType Umu(out.Grid()); @@ -778,7 +255,7 @@ class GaugeGroup { ColdConfiguration(out); } - template + template static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) { out = Ta(in); } diff --git a/Grid/qcd/utils/SUn_impl.h b/Grid/qcd/utils/SUn_impl.h new file mode 100644 index 00000000..e7b972ef --- /dev/null +++ b/Grid/qcd/utils/SUn_impl.h @@ -0,0 +1,552 @@ +// This file is #included into the body of the class template definition of +// GaugeGroup. So, image there to be +// +// template +// class GaugeGroup { +// +// around it. + +template +using iSUnMatrix = iGroupMatrix; +template +using iSUnAlgebraVector = iAlgebraVector; +//////////////////////////////////////////////////////////////////////// +// There are N^2-1 generators for SU(N). +// +// We take a traceless hermitian generator basis as follows +// +// * Normalisation: trace ta tb = 1/2 delta_ab = T_F delta_ab +// T_F = 1/2 for SU(N) groups +// +// * Off diagonal +// - pairs of rows i1,i2 behaving like pauli matrices signma_x, sigma_y +// +// - there are (Nc-1-i1) slots for i2 on each row [ x 0 x ] +// direct count off each row +// +// - Sum of all pairs is Nc(Nc-1)/2: proof arithmetic series +// +// (Nc-1) + (Nc-2)+... 1 ==> Nc*(Nc-1)/2 +// 1+ 2+ + + Nc-1 +// +// - There are 2 x Nc (Nc-1)/ 2 of these = Nc^2 - Nc +// +// - We enumerate the row-col pairs. +// - for each row col pair there is a (sigma_x) and a (sigma_y) like +// generator +// +// +// t^a_ij = { in 0.. Nc(Nc-1)/2 -1} => 1/2(delta_{i,i1} delta_{j,i2} + +// delta_{i,i1} delta_{j,i2}) +// t^a_ij = { in Nc(Nc-1)/2 ... Nc(Nc-1) - 1} => i/2( delta_{i,i1} +// delta_{j,i2} - i delta_{i,i1} delta_{j,i2}) +// +// * Diagonal; must be traceless and normalised +// - Sequence is +// N (1,-1,0,0...) +// N (1, 1,-2,0...) +// N (1, 1, 1,-3,0...) +// N (1, 1, 1, 1,-4,0...) +// +// where 1/2 = N^2 (1+.. m^2)etc.... for the m-th diagonal generator +// NB this gives the famous SU3 result for su2 index 8 +// +// N= sqrt(1/2 . 1/6 ) = 1/2 . 1/sqrt(3) +// +// ( 1 ) +// ( 1 ) / sqrt(3) /2 = 1/2 lambda_8 +// ( -2) +// +//////////////////////////////////////////////////////////////////////// +template +static void generator(int lieIndex, iSUnMatrix &ta) { + // map lie index to which type of generator + int diagIndex; + int su2Index; + int sigxy; + int NNm1 = ncolour * (ncolour - 1); + if (lieIndex >= NNm1) { + diagIndex = lieIndex - NNm1; + generatorDiagonal(diagIndex, ta); + return; + } + sigxy = lieIndex & 0x1; // even or odd + su2Index = lieIndex >> 1; + if (sigxy) + generatorSigmaY(su2Index, ta); + else + generatorSigmaX(su2Index, ta); +} + +template +static void generatorSigmaY(int su2Index, iSUnMatrix &ta) { + ta = Zero(); + int i1, i2; + su2SubGroupIndex(i1, i2, su2Index); + ta()()(i1, i2) = 1.0; + ta()()(i2, i1) = 1.0; + ta = ta * 0.5; +} + +template +static void generatorSigmaX(int su2Index, iSUnMatrix &ta) { + ta = Zero(); + cplx i(0.0, 1.0); + int i1, i2; + su2SubGroupIndex(i1, i2, su2Index); + ta()()(i1, i2) = i; + ta()()(i2, i1) = -i; + ta = ta * 0.5; +} + +template +static void generatorDiagonal(int diagIndex, iSUnMatrix &ta) { + // diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...) + ta = Zero(); + int k = diagIndex + 1; // diagIndex starts from 0 + for (int i = 0; i <= diagIndex; i++) { // k iterations + ta()()(i, i) = 1.0; + } + ta()()(k, k) = -k; // indexing starts from 0 + RealD nrm = 1.0 / std::sqrt(2.0 * k * (k + 1)); + ta = ta * nrm; +} + +//////////////////////////////////////////////////////////////////////// +// Map a su2 subgroup number to the pair of rows that are non zero +//////////////////////////////////////////////////////////////////////// +template +static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { + assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2)); + + int spare = su2_index; + for (i1 = 0; spare >= (ncolour - 1 - i1); i1++) { + spare = spare - (ncolour - 1 - i1); // remove the Nc-1-i1 terms + } + i2 = i1 + 1 + spare; +} + +////////////////////////////////////////////////////////////////////////////////////////// +// Pull out a subgroup and project on to real coeffs x pauli basis +////////////////////////////////////////////////////////////////////////////////////////// +template +static void su2Extract(Lattice > &Determinant, + Lattice > &subgroup, + const Lattice > &source, + int su2_index) { + GridBase *grid(source.Grid()); + conformable(subgroup, source); + conformable(subgroup, Determinant); + int i0, i1; + su2SubGroupIndex(i0, i1, su2_index); + + autoView(subgroup_v, subgroup, AcceleratorWrite); + autoView(source_v, source, AcceleratorRead); + autoView(Determinant_v, Determinant, AcceleratorWrite); + accelerator_for(ss, grid->oSites(), 1, { + subgroup_v[ss]()()(0, 0) = source_v[ss]()()(i0, i0); + subgroup_v[ss]()()(0, 1) = source_v[ss]()()(i0, i1); + subgroup_v[ss]()()(1, 0) = source_v[ss]()()(i1, i0); + subgroup_v[ss]()()(1, 1) = source_v[ss]()()(i1, i1); + + iSU2Matrix Sigma = subgroup_v[ss]; + + Sigma = Sigma - adj(Sigma) + trace(adj(Sigma)); + + subgroup_v[ss] = Sigma; + + // this should be purely real + Determinant_v[ss] = + Sigma()()(0, 0) * Sigma()()(1, 1) - Sigma()()(0, 1) * Sigma()()(1, 0); + }); +} + +////////////////////////////////////////////////////////////////////////////////////////// +// Set matrix to one and insert a pauli subgroup +////////////////////////////////////////////////////////////////////////////////////////// +template +static void su2Insert(const Lattice > &subgroup, + Lattice > &dest, int su2_index) { + GridBase *grid(dest.Grid()); + conformable(subgroup, dest); + int i0, i1; + su2SubGroupIndex(i0, i1, su2_index); + + dest = 1.0; // start out with identity + autoView(dest_v, dest, AcceleratorWrite); + autoView(subgroup_v, subgroup, AcceleratorRead); + accelerator_for(ss, grid->oSites(), 1, { + dest_v[ss]()()(i0, i0) = subgroup_v[ss]()()(0, 0); + dest_v[ss]()()(i0, i1) = subgroup_v[ss]()()(0, 1); + dest_v[ss]()()(i1, i0) = subgroup_v[ss]()()(1, 0); + dest_v[ss]()()(i1, i1) = subgroup_v[ss]()()(1, 1); + }); +} + +/////////////////////////////////////////////// +// Generate e^{ Re Tr Staple Link} dlink +// +// *** Note Staple should be appropriate linear compbination between all +// staples. +// *** If already by beta pass coefficient 1.0. +// *** This routine applies the additional 1/Nc factor that comes after trace +// in action. +// +/////////////////////////////////////////////// +template +static void SubGroupHeatBath( + GridSerialRNG &sRNG, GridParallelRNG &pRNG, + RealD beta, // coeff multiplying staple in action (with no 1/Nc) + LatticeMatrix &link, + const LatticeMatrix &barestaple, // multiplied by action coeffs so th + int su2_subgroup, int nheatbath, LatticeInteger &wheremask) { + GridBase *grid = link.Grid(); + + const RealD twopi = 2.0 * M_PI; + + LatticeMatrix staple(grid); + + staple = barestaple * (beta / ncolour); + + LatticeMatrix V(grid); + V = link * staple; + + // Subgroup manipulation in the lie algebra space + LatticeSU2Matrix u( + grid); // Kennedy pendleton "u" real projected normalised Sigma + LatticeSU2Matrix uinv(grid); + LatticeSU2Matrix ua(grid); // a in pauli form + LatticeSU2Matrix b(grid); // rotated matrix after hb + + // Some handy constant fields + LatticeComplex ones(grid); + ones = 1.0; + LatticeComplex zeros(grid); + zeros = Zero(); + LatticeReal rones(grid); + rones = 1.0; + LatticeReal rzeros(grid); + rzeros = Zero(); + LatticeComplex udet(grid); // determinant of real(staple) + LatticeInteger mask_true(grid); + mask_true = 1; + LatticeInteger mask_false(grid); + mask_false = 0; + + /* + PLB 156 P393 (1985) (Kennedy and Pendleton) + + Note: absorb "beta" into the def of sigma compared to KP paper; staple + passed to this routine has "beta" already multiplied in + + Action linear in links h and of form: + + beta S = beta Sum_p (1 - 1/Nc Re Tr Plaq ) + + Writing Sigma = 1/Nc (beta Sigma') where sum over staples is "Sigma' " + + beta S = const - beta/Nc Re Tr h Sigma' + = const - Re Tr h Sigma + + Decompose h and Sigma into (1, sigma_j) ; h_i real, h^2=1, Sigma_i complex + arbitrary. + + Tr h Sigma = h_i Sigma_j Tr (sigma_i sigma_j) = h_i Sigma_j 2 delta_ij + Re Tr h Sigma = 2 h_j Re Sigma_j + + Normalised re Sigma_j = xi u_j + + With u_j a unit vector and U can be in SU(2); + + Re Tr h Sigma = 2 h_j Re Sigma_j = 2 xi (h.u) + + 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] + u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] + + xi = sqrt(Det)/2; + + Write a= u h in SU(2); a has pauli decomp a_j; + + Note: Product b' xi is unvariant because scaling Sigma leaves + normalised vector "u" fixed; Can rescale Sigma so b' = 1. + */ + + //////////////////////////////////////////////////////// + // Real part of Pauli decomposition + // Note a subgroup can project to zero in cold start + //////////////////////////////////////////////////////// + su2Extract(udet, u, V, su2_subgroup); + + ////////////////////////////////////////////////////// + // Normalising this vector if possible; else identity + ////////////////////////////////////////////////////// + LatticeComplex xi(grid); + + LatticeSU2Matrix lident(grid); + + SU2Matrix ident = Complex(1.0); + SU2Matrix pauli1; + SU<2>::generator(0, pauli1); + SU2Matrix pauli2; + SU<2>::generator(1, pauli2); + SU2Matrix pauli3; + SU<2>::generator(2, pauli3); + pauli1 = timesI(pauli1) * 2.0; + pauli2 = timesI(pauli2) * 2.0; + pauli3 = timesI(pauli3) * 2.0; + + LatticeComplex cone(grid); + LatticeReal adet(grid); + adet = abs(toReal(udet)); + lident = Complex(1.0); + cone = Complex(1.0); + Real machine_epsilon = 1.0e-7; + u = where(adet > machine_epsilon, u, lident); + udet = where(adet > machine_epsilon, udet, cone); + + xi = 0.5 * sqrt(udet); // 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] + u = 0.5 * u * pow(xi, -1.0); // u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] + + // Debug test for sanity + uinv = adj(u); + b = u * uinv - 1.0; + assert(norm2(b) < 1.0e-4); + + /* + Measure: Haar measure dh has d^4a delta(1-|a^2|) + In polars: + da = da0 r^2 sin theta dr dtheta dphi delta( 1 - r^2 -a0^2) + = da0 r^2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r)(sqrt(1-a0^) + + r) ) + = da0 r/2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r) ) + + Action factor Q(h) dh = e^-S[h] dh = e^{ xi Tr uh} dh // beta + enters through xi = e^{2 xi (h.u)} dh = e^{2 xi h0u0}.e^{2 xi h1u1}.e^{2 + xi h2u2}.e^{2 xi h3u3} dh + + Therefore for each site, take xi for that site + i) generate |a0|<1 with dist + (1-a0^2)^0.5 e^{2 xi a0 } da0 + + Take alpha = 2 xi = 2 xi [ recall 2 beta/Nc unmod staple norm]; + hence 2.0/Nc factor in Chroma ] A. Generate two uniformly distributed + pseudo-random numbers R and R', R'', R''' in the unit interval; B. Set X = + -(ln R)/alpha, X' =-(ln R')/alpha; C. Set C = cos^2(2pi R"), with R" + another uniform random number in [0,1] ; D. Set A = XC; E. Let d = X'+A; + F. If R'''^2 :> 1 - 0.5 d, go back to A; + G. Set a0 = 1 - d; + + Note that in step D setting B ~ X - A and using B in place of A in step E + will generate a second independent a 0 value. + */ + + ///////////////////////////////////////////////////////// + // count the number of sites by picking "1"'s out of hat + ///////////////////////////////////////////////////////// + Integer hit = 0; + LatticeReal rtmp(grid); + rtmp = where(wheremask, rones, rzeros); + RealD numSites = sum(rtmp); + RealD numAccepted; + LatticeInteger Accepted(grid); + Accepted = Zero(); + LatticeInteger newlyAccepted(grid); + + std::vector xr(4, grid); + std::vector a(4, grid); + LatticeReal d(grid); + d = Zero(); + LatticeReal alpha(grid); + + // std::cout< 1 - 0.5 d, go back to A; + LatticeReal thresh(grid); + thresh = 1.0 - d * 0.5; + xrsq = xr[0] * xr[0]; + LatticeInteger ione(grid); + ione = 1; + LatticeInteger izero(grid); + izero = Zero(); + + newlyAccepted = where(xrsq < thresh, ione, izero); + Accepted = where(newlyAccepted, newlyAccepted, Accepted); + Accepted = where(wheremask, Accepted, izero); + + // FIXME need an iSum for integer to avoid overload on return type?? + rtmp = where(Accepted, rones, rzeros); + numAccepted = sum(rtmp); + + hit++; + + } while ((numAccepted < numSites) && (hit < nheatbath)); + + // G. Set a0 = 1 - d; + a[0] = Zero(); + a[0] = where(wheremask, 1.0 - d, a[0]); + + ////////////////////////////////////////// + // ii) generate a_i uniform on two sphere radius (1-a0^2)^0.5 + ////////////////////////////////////////// + + LatticeReal a123mag(grid); + a123mag = sqrt(abs(1.0 - a[0] * a[0])); + + LatticeReal cos_theta(grid); + LatticeReal sin_theta(grid); + LatticeReal phi(grid); + + random(pRNG, phi); + phi = phi * twopi; // uniform in [0,2pi] + random(pRNG, cos_theta); + cos_theta = (cos_theta * 2.0) - 1.0; // uniform in [-1,1] + sin_theta = sqrt(abs(1.0 - cos_theta * cos_theta)); + + a[1] = a123mag * sin_theta * cos(phi); + a[2] = a123mag * sin_theta * sin(phi); + a[3] = a123mag * cos_theta; + + ua = toComplex(a[0]) * ident + toComplex(a[1]) * pauli1 + + toComplex(a[2]) * pauli2 + toComplex(a[3]) * pauli3; + + b = 1.0; + b = where(wheremask, uinv * ua, b); + su2Insert(b, V, su2_subgroup); + + // mask the assignment back based on Accptance + link = where(Accepted, V * link, link); + + ////////////////////////////// + // Debug Checks + // SU2 check + LatticeSU2Matrix check(grid); // rotated matrix after hb + u = Zero(); + check = ua * adj(ua) - 1.0; + check = where(Accepted, check, u); + assert(norm2(check) < 1.0e-4); + + check = b * adj(b) - 1.0; + check = where(Accepted, check, u); + assert(norm2(check) < 1.0e-4); + + LatticeMatrix Vcheck(grid); + Vcheck = Zero(); + Vcheck = where(Accepted, V * adj(V) - 1.0, Vcheck); + // std::cout< +static void testGenerators(void) { + Matrix ta; + Matrix tb; + std::cout << GridLogMessage + << "Fundamental - Checking trace ta tb is 0.5 delta_ab" + << std::endl; + for (int a = 0; a < AdjointDimension; a++) { + for (int b = 0; b < AdjointDimension; b++) { + generator(a, ta); + generator(b, tb); + Complex tr = TensorRemove(trace(ta * tb)); + std::cout << GridLogMessage << "(" << a << "," << b << ") = " << tr + << std::endl; + if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6); + if (a != b) assert(abs(tr) < 1.0e-6); + } + std::cout << GridLogMessage << std::endl; + } + std::cout << GridLogMessage << "Fundamental - Checking if hermitian" + << std::endl; + for (int a = 0; a < AdjointDimension; a++) { + generator(a, ta); + std::cout << GridLogMessage << a << std::endl; + assert(norm2(ta - adj(ta)) < 1.0e-6); + } + std::cout << GridLogMessage << std::endl; + + std::cout << GridLogMessage << "Fundamental - Checking if traceless" + << std::endl; + for (int a = 0; a < AdjointDimension; a++) { + generator(a, ta); + Complex tr = TensorRemove(trace(ta)); + std::cout << GridLogMessage << a << " " << std::endl; + assert(abs(tr) < 1.0e-6); + } + std::cout << GridLogMessage << std::endl; +} + +/* + * Fundamental rep gauge xform + */ +template +static void GaugeTransformFundamental(Fundamental &ferm, GaugeMat &g) { + GridBase *grid = ferm._grid; + conformable(grid, g._grid); + ferm = g * ferm; +} +/* + * Adjoint rep gauge xform + */ + +template +static void GaugeTransform(GaugeField &Umu, GaugeMat &g) { + GridBase *grid = Umu.Grid(); + conformable(grid, g.Grid()); + + GaugeMat U(grid); + GaugeMat ag(grid); + ag = adj(g); + + for (int mu = 0; mu < Nd; mu++) { + U = PeekIndex(Umu, mu); + U = g * U * Cshift(ag, mu, 1); + PokeIndex(Umu, U, mu); + } +} +template +static void GaugeTransform(std::vector &U, GaugeMat &g) { + GridBase *grid = g.Grid(); + GaugeMat ag(grid); + ag = adj(g); + for (int mu = 0; mu < Nd; mu++) { + U[mu] = g * U[mu] * Cshift(ag, mu, 1); + } +} +template +static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, + GaugeMat &g) { + LieRandomize(pRNG, g, 1.0); + GaugeTransform(Umu, g); +} From 54c8025aad0263b8258d8e6b2a683689f6ec9260 Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Mon, 28 Nov 2022 17:50:38 +0000 Subject: [PATCH 08/17] Remove unnecessary pwd in scripts/filelist --- scripts/filelist | 1 - 1 file changed, 1 deletion(-) diff --git a/scripts/filelist b/scripts/filelist index 4e64cd62..d91da5f1 100755 --- a/scripts/filelist +++ b/scripts/filelist @@ -37,7 +37,6 @@ cd $home/tests dirs=`find . -type d -not -path '*/\.*'` for subdir in $dirs; do cd $home/tests/$subdir - pwd TESTS=`ls T*.cc` TESTLIST=`echo ${TESTS} | sed s/.cc//g ` PREF=`[ $subdir = '.' ] && echo noinst || echo EXTRA` From 7c2ad4f8c8c42aee5fc8aba0dd3be0d970776b3c Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Wed, 30 Nov 2022 11:57:39 +0000 Subject: [PATCH 09/17] Attempt with SFINAE (failed) --- Grid/qcd/utils/SUn.h | 42 +- Grid/qcd/utils/SUn_impl.h | 2 + Grid/qcd/utils/Sp2n.h | 1161 ++++++++++++++++++------------------ Grid/qcd/utils/Sp2n_impl.h | 245 ++++++++ 4 files changed, 866 insertions(+), 584 deletions(-) create mode 100644 Grid/qcd/utils/Sp2n_impl.h diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index 4a7148e9..0c961ab4 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -33,12 +33,19 @@ directory #define QCD_UTIL_SUN_H #define ONLY_IF_SU \ - typename dummy_name = group_name, \ - typename = std::enable_if_t::value> + typename \ + typename std::enable_if_t::value, dummy_name> ,\ + dummy_name = group_name, + +#define ONLY_IF_Sp \ + typename \ + typename std::enable_if_t::value, dummy_name> ,\ + dummy_name =group_name, NAMESPACE_BEGIN(Grid); namespace GroupName { class SU {}; +class Sp {}; } // namespace GroupName template @@ -51,18 +58,44 @@ struct is_su { static const bool value = true; }; +template +struct is_sp { + static const bool value = false; +}; + +template <> +struct is_sp { + static const bool value = true; +}; + +template +constexpr int compute_adjoint_dimension(int ncolour); + +template<> +constexpr int compute_adjoint_dimension(int ncolour) { return ncolour * ncolour - 1;} + +template<> +constexpr int compute_adjoint_dimension(int ncolour) { return ncolour/2*(ncolour+1);} + + template class GaugeGroup; template using SU = GaugeGroup; +template +using Sp = GaugeGroup; + template class GaugeGroup { public: static const int Dimension = ncolour; - static const int AdjointDimension = ncolour * ncolour - 1; - static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; } + static const int AdjointDimension = compute_adjoint_dimension(ncolour); + static const int AlgebraDimension = compute_adjoint_dimension(ncolour); + // Don't know how to only enable this for Sp: + static const int nsp = ncolour/2; + template using iSU2Matrix = iScalar > >; @@ -115,6 +148,7 @@ class GaugeGroup { typedef Lattice LatticeSU2MatrixD; #include "Grid/qcd/utils/SUn_impl.h" +#include "Grid/qcd/utils/Sp2n_impl.h" static void printGenerators(void) { for (int gen = 0; gen < AdjointDimension; gen++) { diff --git a/Grid/qcd/utils/SUn_impl.h b/Grid/qcd/utils/SUn_impl.h index e7b972ef..30d8ff14 100644 --- a/Grid/qcd/utils/SUn_impl.h +++ b/Grid/qcd/utils/SUn_impl.h @@ -6,6 +6,8 @@ // // around it. +template + static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; } template using iSUnMatrix = iGroupMatrix; template diff --git a/Grid/qcd/utils/Sp2n.h b/Grid/qcd/utils/Sp2n.h index 19e30568..1f7009e3 100644 --- a/Grid/qcd/utils/Sp2n.h +++ b/Grid/qcd/utils/Sp2n.h @@ -1,580 +1,581 @@ - -#ifndef QCD_UTIL_Sp2n_H -#define QCD_UTIL_Sp2n_H - -NAMESPACE_BEGIN(Grid); - -// Sp(2N) -// ncolour = 2N - - -template -class Sp { -public: - static const int nsp = ncolour/2; - static const int Dimension = ncolour; - static const int AlgebraDimension = nsp*(2*nsp +1); - static int su2subgroups(void) { return (nsp * (nsp - 1)) / 2; } - - - template - using iSp2nMatrix = iScalar > >; - template - using iSU2Matrix = iScalar > >; - template - using iSp2nAlgebraVector = iScalar > >; - - typedef iSp2nMatrix Matrix; - typedef iSp2nMatrix MatrixF; - typedef iSp2nMatrix MatrixD; - - typedef iSp2nMatrix vMatrix; - typedef iSp2nMatrix vMatrixF; - typedef iSp2nMatrix vMatrixD; - - typedef iSp2nAlgebraVector AlgebraVector; - typedef iSp2nAlgebraVector AlgebraVectorF; - typedef iSp2nAlgebraVector AlgebraVectorD; - - typedef iSp2nAlgebraVector vAlgebraVector; - typedef iSp2nAlgebraVector vAlgebraVectorF; - typedef iSp2nAlgebraVector vAlgebraVectorD; - - typedef Lattice LatticeMatrix; - typedef Lattice LatticeMatrixF; - typedef Lattice LatticeMatrixD; - - typedef Lattice LatticeAlgebraVector; - typedef Lattice LatticeAlgebraVectorF; - typedef Lattice LatticeAlgebraVectorD; - - - - // Sp(2N) has N(2N+1) = 2N^2+N generators - // - // normalise the generators such that - // Trace ( Ta Tb) = 1/2 delta_ab - // - // N generators in the cartan, 2N^2 off - // off diagonal: - // there are 6 types named a,b,c,d and w,z - // abcd are N(N-1)/2 each while wz are N each - - template - static void generator(int lieIndex, iSp2nMatrix &ta) { - // map lie index into type of generators: diagonal, abcd type, wz type - - int diagIndex; - int aIndex, bIndex, cIndex, dIndex; - int wIndex, zIndex; // a,b,c,d are N(N-1)/2 and w,z are N - int mod = nsp * (nsp-1) * 0.5; - int offdiag = 2*nsp*nsp; // number of generators not in the cartan subalgebra - int wmod = 4*mod; - int zmod = wmod+nsp; - if (lieIndex >= offdiag) { - diagIndex = lieIndex - offdiag; // 0, ... ,N-1 - //std::cout << GridLogMessage << "diag type " << std::endl; - generatorDiagtype(diagIndex, ta); - return; - } - if ( (lieIndex >= wmod) && (lieIndex < zmod) ) { - //std::cout << GridLogMessage << "w type " << std::endl; - wIndex = lieIndex- wmod; // 0, ... ,N-1 - generatorWtype(wIndex,ta); - return; - } - if ( (lieIndex >= zmod) && (lieIndex < offdiag) ) { - //std::cout << GridLogMessage << "z type " << std::endl; - //std::cout << GridLogMessage << "lie index " << lieIndex << std::endl; - //std::cout << GridLogMessage << "z mod " << zmod << std::endl; - zIndex = lieIndex - zmod; // 0, ... ,N-1 - generatorZtype(zIndex,ta); - return; - } - if (lieIndex < mod) { // atype 0, ... , N(N-1)/2=mod - //std::cout << GridLogMessage << "a type " << std::endl; - aIndex = lieIndex; - //std::cout << GridLogMessage << "a indx " << aIndex << std::endl; - generatorAtype(aIndex, ta); - return; - } - if ( (lieIndex >= mod) && lieIndex < 2*mod) { // btype mod, ... , 2mod-1 - //std::cout << GridLogMessage << "b type " << std::endl; - bIndex = lieIndex - mod; - generatorBtype(bIndex, ta); - return; - } - if ( (lieIndex >= 2*mod) && lieIndex < 3*mod) { // ctype 2mod, ... , 3mod-1 - //std::cout << GridLogMessage << "c type " << std::endl; - cIndex = lieIndex - 2*mod; - generatorCtype(cIndex, ta); - return; - } - if ( (lieIndex >= 3*mod) && lieIndex < wmod) { // ctype 3mod, ... , 4mod-1 = wmod-1 - //std::cout << GridLogMessage << "d type " << std::endl; - dIndex = lieIndex - 3*mod; - generatorDtype(dIndex, ta); - return; - } - - } //end of generator - - template - static void generatorDiagtype(int diagIndex, iSp2nMatrix &ta) { - - // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra - - ta = Zero(); - RealD nrm = 1.0 / 2; - - ta()()(diagIndex,diagIndex) = nrm; - ta()()(diagIndex+nsp,diagIndex+nsp) = -nrm; - } - - template - static void generatorAtype(int aIndex, iSp2nMatrix &ta) { - - // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2) - // with i - static void generatorBtype(int bIndex, iSp2nMatrix &ta) { - - // ta(i,j) = -ta(j,i) = ta(i+N,j+N) = -ta(j+N,i+N) = i / 1/ 2 sqrt(2) - // with i - static void generatorCtype(int cIndex, iSp2nMatrix &ta) { - - // ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2) - - - int i1, i2; - ta = Zero(); - RealD nrm = 1 / (2 * std::sqrt(2) ); - su2SubGroupIndex(i1, i2, cIndex); - - ta()()(i1,i2+nsp) = 1; - ta()()(i2,i1+nsp) = 1; - ta()()(i1+nsp,i2) = 1; - ta()()(i2+nsp,i1) = 1; - - ta = ta * nrm; - } - - template - static void generatorDtype(int dIndex, iSp2nMatrix &ta) { - - // ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2) - - int i1, i2; - ta = Zero(); - cplx i(0.0, 1.0); - RealD nrm = 1 / (2 * std::sqrt(2) ); - su2SubGroupIndex(i1, i2, dIndex); - - ta()()(i1,i2+nsp) = i; - ta()()(i2,i1+nsp) = i; - ta()()(i1+nsp,i2) = -i; - ta()()(i2+nsp,i1) = -i; - - ta = ta * nrm; - } - - template - static void generatorWtype(int wIndex, iSp2nMatrix &ta) { - - // ta(i,i+N) = ta(i+N,i) = 1/2 - - ta = Zero(); - RealD nrm = 1.0 / 2; //check - - ta()()(wIndex,wIndex+nsp) = 1; - ta()()(wIndex+nsp,wIndex) = 1; - - ta = ta * nrm; - } - - template - static void generatorZtype(int zIndex, iSp2nMatrix &ta) { - - // ta(i,i+N) = - ta(i+N,i) = i/2 - - ta = Zero(); - RealD nrm = 1.0 / 2; //check - cplx i(0.0, 1.0); - ta()()(zIndex,zIndex+nsp) = i; - ta()()(zIndex+nsp,zIndex) = -i; - - ta = ta * nrm; - } - - - //////////////////////////////////////////////////////////////////////// - // Map a su2 subgroup number to the pair of rows that are non zero - //////////////////////////////////////////////////////////////////////// - static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { - assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2)); - - int spare = su2_index; - for (i1 = 0; spare >= (nsp - 1 - i1); i1++) { - spare = spare - (nsp - 1 - i1); // remove the Nc-1-i1 terms - } - i2 = i1 + 1 + spare; - } - - - - - - static void printGenerators(void) { - for (int gen = 0; gen < AlgebraDimension; gen++) { - Matrix ta; - generator(gen, ta); - std::cout << GridLogMessage << "Nc = " << ncolour << std::endl; - std::cout << GridLogMessage << " t_" << gen << std::endl; - std::cout << GridLogMessage << ta << std::endl; - } - } - - - - static void testGenerators(void) { - Matrix ta; - Matrix tb; - std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is 0.5 delta_ab " << std::endl; - for (int a = 0; a < AlgebraDimension; a++) { - for (int b = 0; b < AlgebraDimension; b++) { - generator(a,ta); - generator(b,tb); - Complex tr = TensorRemove(trace( ta * tb) ); - std::cout << GridLogMessage << "(" << a << "," << b << ") = " << tr - << std::endl; - if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6); - if (a != b) assert(abs(tr) < 1.0e-6); - - } - } - std::cout << GridLogMessage << std::endl; - std::cout << GridLogMessage << "Fundamental - Checking if hermitian" << std::endl; - for (int a = 0; a < AlgebraDimension; a++) { - generator(a,ta); - std::cout << GridLogMessage << a << std::endl; - assert(norm2(ta - adj(ta)) < 1.0e-6); - } - std::cout << GridLogMessage << std::endl; - std::cout << GridLogMessage << "Fundamental - Checking if traceless" << std::endl; - for (int a = 0; a < AlgebraDimension; a++) { - generator(a, ta); - Complex tr = TensorRemove(trace(ta)); - std::cout << GridLogMessage << a << std::endl; - assert(abs(tr) < 1.0e-6); - } - - } - - - 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::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 < AlgebraDimension; 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, //same as sun - 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 < AlgebraDimension; a++) { - gaussian(pRNG, ca); - generator(a, ta); - la = toComplex(ca) * ta; - out += la; - } - out *= ci; - } - - 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 < AlgebraDimension; a++) { - generator(a, ta); - la = peekColour(h, a) * timesI(ta) * scale; - out += la; - } - } - - - template - static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { // same as sun - typedef typename LatticeMatrixType::scalar_type ComplexType; - - LatticeMatrixType xn(x.Grid()); - RealD nfac = 1.0; - - xn = x; - ex = xn + ComplexType(1.0); // 1+x - - // Do a 12th order exponentiation - for (int i = 2; i <= 12; ++i) { - nfac = nfac / RealD(i); // 1/2, 1/2.3 ... - xn = xn * x; // x2, x3,x4.... - ex = ex + xn * nfac; // x2/2!, x3/3!.... - } - } - - 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 < AlgebraDimension; a++) { - generator(a, Ta); - pokeColour(h_out, - 2.0 * (trace(timesI(Ta) * in)) * scale, a); - } - } - - - template - static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) { - typedef typename GaugeField::vector_type vector_type; - typedef iSp2nMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType Umu(out.Grid()); - for (int mu = 0; mu < Nd; mu++) { - LieRandomize(pRNG, Umu, 1.0); //def - PokeIndex(out, Umu, mu); - } - } - template - static void TepidConfiguration(GridParallelRNG &pRNG,GaugeField &out){ - typedef typename GaugeField::vector_type vector_type; - typedef iSp2nMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType Umu(out.Grid()); - for(int mu=0;mu(out,Umu,mu); - } - } - - template - static void ColdConfiguration(GaugeField &out){ - typedef typename GaugeField::vector_type vector_type; - typedef iSp2nMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType Umu(out.Grid()); - Umu=1.0; - for(int mu=0;mu(out,Umu,mu); - } - } - template - static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){ - ColdConfiguration(out); - } - - static void OmegaInvariance(ColourMatrix &in) - { - - ColourMatrix Omega; - Omega = Zero(); - - 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; - } - - auto diff = Omega - (in * Omega * transpose(in) ); - auto sdiff = norm2(diff); - if (norm2(sdiff) < 1e-8) - { - std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl; - } else { - std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT left invariant by " << norm2(sdiff) << std::endl; - } - - } - - template - static void OmegaInvariance(GaugeField &in) - { - typedef typename GaugeField::vector_type vector_type; - typedef iSp2nMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType U(in.Grid()); - LatticeMatrixType OmegaLatt(in.Grid()); - LatticeMatrixType identity(in.Grid()); - RealD vol = in.Grid()->gSites(); - ColourMatrix Omega; - - OmegaLatt = Zero(); - Omega = Zero(); - identity = 1.; - - std::cout << GridLogMessage << "I am a GaugeField " << std::endl; - - U = PeekIndex(in,1); - - OmegaInvariance(U); - - } - - static void OmegaInvariance(LatticeColourMatrixD &in) - { - - LatticeColourMatrixD OmegaLatt(in.Grid()); - LatticeColourMatrixD identity(in.Grid()); - RealD vol = in.Grid()->gSites(); - ColourMatrix Omega; - - OmegaLatt = Zero(); - 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; - } - - std::cout << GridLogMessage << "Omega = " << Omega()() << std::endl; - OmegaLatt = OmegaLatt + (identity*Omega); - - auto diff = OmegaLatt - (in*OmegaLatt*transpose(in)); - auto sdiff = sum(diff); - //assert( norm2(sdiff) < 1e-8 ); - if (norm2(sdiff) < 1e-8) - { - std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl; - } else { - std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT left invariant by " << norm2(sdiff) << std::endl; - } - - - } - - -}; // end of class Sp - - - template - static void ProjectSp2n(Lattice > > > &Umu) - { - Umu = ProjectOnSpGroup(Umu); - auto det = Determinant(Umu); // ok ? - - det = conjugate(det); - - for(int i=0;i(Umu,N-1,i); - element = element * det; - PokeIndex(Umu,element,Nc-1,i); - } - } - template - static void ProjectSp2n(Lattice >,Nd> > &U) - { - GridBase *grid=U.Grid(); - for(int mu=0;mu(U,mu); - Umu = ProjectOnSpGroup(Umu); - ProjectSp2n(Umu); - PokeIndex(U,Umu,mu); - } - } - -typedef Sp<2> Sp2; -typedef Sp<4> Sp4; -typedef Sp<6> Sp6; -typedef Sp<8> Sp8; - -NAMESPACE_END(Grid); -#endif +#include "Grid/qcd/utils/SUn.h" +// +// #ifndef QCD_UTIL_Sp2n_H +// #define QCD_UTIL_Sp2n_H +// +// NAMESPACE_BEGIN(Grid); +// +// // Sp(2N) +// // ncolour = 2N +// +// +// template +// class Sp { +// public: +// static const int nsp = ncolour/2; +// static const int Dimension = ncolour; +// static const int AlgebraDimension = nsp*(2*nsp +1); +// static int su2subgroups(void) { return (nsp * (nsp - 1)) / 2; } +// +// +// template +// using iSp2nMatrix = iScalar > >; +// template +// using iSU2Matrix = iScalar > >; +// template +// using iSp2nAlgebraVector = iScalar > >; +// +// typedef iSp2nMatrix Matrix; +// typedef iSp2nMatrix MatrixF; +// typedef iSp2nMatrix MatrixD; +// +// typedef iSp2nMatrix vMatrix; +// typedef iSp2nMatrix vMatrixF; +// typedef iSp2nMatrix vMatrixD; +// +// typedef iSp2nAlgebraVector AlgebraVector; +// typedef iSp2nAlgebraVector AlgebraVectorF; +// typedef iSp2nAlgebraVector AlgebraVectorD; +// +// typedef iSp2nAlgebraVector vAlgebraVector; +// typedef iSp2nAlgebraVector vAlgebraVectorF; +// typedef iSp2nAlgebraVector vAlgebraVectorD; +// +// typedef Lattice LatticeMatrix; +// typedef Lattice LatticeMatrixF; +// typedef Lattice LatticeMatrixD; +// +// typedef Lattice LatticeAlgebraVector; +// typedef Lattice LatticeAlgebraVectorF; +// typedef Lattice LatticeAlgebraVectorD; +// +// +// +// // Sp(2N) has N(2N+1) = 2N^2+N generators +// // +// // normalise the generators such that +// // Trace ( Ta Tb) = 1/2 delta_ab +// // +// // N generators in the cartan, 2N^2 off +// // off diagonal: +// // there are 6 types named a,b,c,d and w,z +// // abcd are N(N-1)/2 each while wz are N each +// +// template +// static void generator(int lieIndex, iSp2nMatrix &ta) { +// // map lie index into type of generators: diagonal, abcd type, wz type +// +// int diagIndex; +// int aIndex, bIndex, cIndex, dIndex; +// int wIndex, zIndex; // a,b,c,d are N(N-1)/2 and w,z are N +// int mod = nsp * (nsp-1) * 0.5; +// int offdiag = 2*nsp*nsp; // number of generators not in the cartan subalgebra +// int wmod = 4*mod; +// int zmod = wmod+nsp; +// if (lieIndex >= offdiag) { +// diagIndex = lieIndex - offdiag; // 0, ... ,N-1 +// //std::cout << GridLogMessage << "diag type " << std::endl; +// generatorDiagtype(diagIndex, ta); +// return; +// } +// if ( (lieIndex >= wmod) && (lieIndex < zmod) ) { +// //std::cout << GridLogMessage << "w type " << std::endl; +// wIndex = lieIndex- wmod; // 0, ... ,N-1 +// generatorWtype(wIndex,ta); +// return; +// } +// if ( (lieIndex >= zmod) && (lieIndex < offdiag) ) { +// //std::cout << GridLogMessage << "z type " << std::endl; +// //std::cout << GridLogMessage << "lie index " << lieIndex << std::endl; +// //std::cout << GridLogMessage << "z mod " << zmod << std::endl; +// zIndex = lieIndex - zmod; // 0, ... ,N-1 +// generatorZtype(zIndex,ta); +// return; +// } +// if (lieIndex < mod) { // atype 0, ... , N(N-1)/2=mod +// //std::cout << GridLogMessage << "a type " << std::endl; +// aIndex = lieIndex; +// //std::cout << GridLogMessage << "a indx " << aIndex << std::endl; +// generatorAtype(aIndex, ta); +// return; +// } +// if ( (lieIndex >= mod) && lieIndex < 2*mod) { // btype mod, ... , 2mod-1 +// //std::cout << GridLogMessage << "b type " << std::endl; +// bIndex = lieIndex - mod; +// generatorBtype(bIndex, ta); +// return; +// } +// if ( (lieIndex >= 2*mod) && lieIndex < 3*mod) { // ctype 2mod, ... , 3mod-1 +// //std::cout << GridLogMessage << "c type " << std::endl; +// cIndex = lieIndex - 2*mod; +// generatorCtype(cIndex, ta); +// return; +// } +// if ( (lieIndex >= 3*mod) && lieIndex < wmod) { // ctype 3mod, ... , 4mod-1 = wmod-1 +// //std::cout << GridLogMessage << "d type " << std::endl; +// dIndex = lieIndex - 3*mod; +// generatorDtype(dIndex, ta); +// return; +// } +// +// } //end of generator +// +// template +// static void generatorDiagtype(int diagIndex, iSp2nMatrix &ta) { +// +// // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra +// +// ta = Zero(); +// RealD nrm = 1.0 / 2; +// +// ta()()(diagIndex,diagIndex) = nrm; +// ta()()(diagIndex+nsp,diagIndex+nsp) = -nrm; +// } +// +// template +// static void generatorAtype(int aIndex, iSp2nMatrix &ta) { +// +// // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2) +// // with i +// static void generatorBtype(int bIndex, iSp2nMatrix &ta) { +// +// // ta(i,j) = -ta(j,i) = ta(i+N,j+N) = -ta(j+N,i+N) = i / 1/ 2 sqrt(2) +// // with i +// static void generatorCtype(int cIndex, iSp2nMatrix &ta) { +// +// // ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2) +// +// +// int i1, i2; +// ta = Zero(); +// RealD nrm = 1 / (2 * std::sqrt(2) ); +// su2SubGroupIndex(i1, i2, cIndex); +// +// ta()()(i1,i2+nsp) = 1; +// ta()()(i2,i1+nsp) = 1; +// ta()()(i1+nsp,i2) = 1; +// ta()()(i2+nsp,i1) = 1; +// +// ta = ta * nrm; +// } +// +// template +// static void generatorDtype(int dIndex, iSp2nMatrix &ta) { +// +// // ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2) +// +// int i1, i2; +// ta = Zero(); +// cplx i(0.0, 1.0); +// RealD nrm = 1 / (2 * std::sqrt(2) ); +// su2SubGroupIndex(i1, i2, dIndex); +// +// ta()()(i1,i2+nsp) = i; +// ta()()(i2,i1+nsp) = i; +// ta()()(i1+nsp,i2) = -i; +// ta()()(i2+nsp,i1) = -i; +// +// ta = ta * nrm; +// } +// +// template +// static void generatorWtype(int wIndex, iSp2nMatrix &ta) { +// +// // ta(i,i+N) = ta(i+N,i) = 1/2 +// +// ta = Zero(); +// RealD nrm = 1.0 / 2; //check +// +// ta()()(wIndex,wIndex+nsp) = 1; +// ta()()(wIndex+nsp,wIndex) = 1; +// +// ta = ta * nrm; +// } +// +// template +// static void generatorZtype(int zIndex, iSp2nMatrix &ta) { +// +// // ta(i,i+N) = - ta(i+N,i) = i/2 +// +// ta = Zero(); +// RealD nrm = 1.0 / 2; //check +// cplx i(0.0, 1.0); +// ta()()(zIndex,zIndex+nsp) = i; +// ta()()(zIndex+nsp,zIndex) = -i; +// +// ta = ta * nrm; +// } +// +// +// //////////////////////////////////////////////////////////////////////// +// // Map a su2 subgroup number to the pair of rows that are non zero +// //////////////////////////////////////////////////////////////////////// +// static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { +// assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2)); +// +// int spare = su2_index; +// for (i1 = 0; spare >= (nsp - 1 - i1); i1++) { +// spare = spare - (nsp - 1 - i1); // remove the Nc-1-i1 terms +// } +// i2 = i1 + 1 + spare; +// } +// +// +// +// +// +// static void printGenerators(void) { +// for (int gen = 0; gen < AlgebraDimension; gen++) { +// Matrix ta; +// generator(gen, ta); +// std::cout << GridLogMessage << "Nc = " << ncolour << std::endl; +// std::cout << GridLogMessage << " t_" << gen << std::endl; +// std::cout << GridLogMessage << ta << std::endl; +// } +// } +// +// +// +// static void testGenerators(void) { +// Matrix ta; +// Matrix tb; +// std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is 0.5 delta_ab " << std::endl; +// for (int a = 0; a < AlgebraDimension; a++) { +// for (int b = 0; b < AlgebraDimension; b++) { +// generator(a,ta); +// generator(b,tb); +// Complex tr = TensorRemove(trace( ta * tb) ); +// std::cout << GridLogMessage << "(" << a << "," << b << ") = " << tr +// << std::endl; +// if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6); +// if (a != b) assert(abs(tr) < 1.0e-6); +// +// } +// } +// std::cout << GridLogMessage << std::endl; +// std::cout << GridLogMessage << "Fundamental - Checking if hermitian" << std::endl; +// for (int a = 0; a < AlgebraDimension; a++) { +// generator(a,ta); +// std::cout << GridLogMessage << a << std::endl; +// assert(norm2(ta - adj(ta)) < 1.0e-6); +// } +// std::cout << GridLogMessage << std::endl; +// std::cout << GridLogMessage << "Fundamental - Checking if traceless" << std::endl; +// for (int a = 0; a < AlgebraDimension; a++) { +// generator(a, ta); +// Complex tr = TensorRemove(trace(ta)); +// std::cout << GridLogMessage << a << std::endl; +// assert(abs(tr) < 1.0e-6); +// } +// +// } +// +// +// 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::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 < AlgebraDimension; 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, //same as sun +// 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 < AlgebraDimension; a++) { +// gaussian(pRNG, ca); +// generator(a, ta); +// la = toComplex(ca) * ta; +// out += la; +// } +// out *= ci; +// } +// +// 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 < AlgebraDimension; a++) { +// generator(a, ta); +// la = peekColour(h, a) * timesI(ta) * scale; +// out += la; +// } +// } +// +// +// template +// static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { // same as sun +// typedef typename LatticeMatrixType::scalar_type ComplexType; +// +// LatticeMatrixType xn(x.Grid()); +// RealD nfac = 1.0; +// +// xn = x; +// ex = xn + ComplexType(1.0); // 1+x +// +// // Do a 12th order exponentiation +// for (int i = 2; i <= 12; ++i) { +// nfac = nfac / RealD(i); // 1/2, 1/2.3 ... +// xn = xn * x; // x2, x3,x4.... +// ex = ex + xn * nfac; // x2/2!, x3/3!.... +// } +// } +// +// 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 < AlgebraDimension; a++) { +// generator(a, Ta); +// pokeColour(h_out, - 2.0 * (trace(timesI(Ta) * in)) * scale, a); +// } +// } +// +// +// template +// static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) { +// typedef typename GaugeField::vector_type vector_type; +// typedef iSp2nMatrix vMatrixType; +// typedef Lattice LatticeMatrixType; +// +// LatticeMatrixType Umu(out.Grid()); +// for (int mu = 0; mu < Nd; mu++) { +// LieRandomize(pRNG, Umu, 1.0); //def +// PokeIndex(out, Umu, mu); +// } +// } +// template +// static void TepidConfiguration(GridParallelRNG &pRNG,GaugeField &out){ +// typedef typename GaugeField::vector_type vector_type; +// typedef iSp2nMatrix vMatrixType; +// typedef Lattice LatticeMatrixType; +// +// LatticeMatrixType Umu(out.Grid()); +// for(int mu=0;mu(out,Umu,mu); +// } +// } +// +// template +// static void ColdConfiguration(GaugeField &out){ +// typedef typename GaugeField::vector_type vector_type; +// typedef iSp2nMatrix vMatrixType; +// typedef Lattice LatticeMatrixType; +// +// LatticeMatrixType Umu(out.Grid()); +// Umu=1.0; +// for(int mu=0;mu(out,Umu,mu); +// } +// } +// template +// static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){ +// ColdConfiguration(out); +// } +// +// static void OmegaInvariance(ColourMatrix &in) +// { +// +// ColourMatrix Omega; +// Omega = Zero(); +// +// 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; +// } +// +// auto diff = Omega - (in * Omega * transpose(in) ); +// auto sdiff = norm2(diff); +// if (norm2(sdiff) < 1e-8) +// { +// std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl; +// } else { +// std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT left invariant by " << norm2(sdiff) << std::endl; +// } +// +// } +// +// template +// static void OmegaInvariance(GaugeField &in) +// { +// typedef typename GaugeField::vector_type vector_type; +// typedef iSp2nMatrix vMatrixType; +// typedef Lattice LatticeMatrixType; +// +// LatticeMatrixType U(in.Grid()); +// LatticeMatrixType OmegaLatt(in.Grid()); +// LatticeMatrixType identity(in.Grid()); +// RealD vol = in.Grid()->gSites(); +// ColourMatrix Omega; +// +// OmegaLatt = Zero(); +// Omega = Zero(); +// identity = 1.; +// +// std::cout << GridLogMessage << "I am a GaugeField " << std::endl; +// +// U = PeekIndex(in,1); +// +// OmegaInvariance(U); +// +// } +// +// static void OmegaInvariance(LatticeColourMatrixD &in) +// { +// +// LatticeColourMatrixD OmegaLatt(in.Grid()); +// LatticeColourMatrixD identity(in.Grid()); +// RealD vol = in.Grid()->gSites(); +// ColourMatrix Omega; +// +// OmegaLatt = Zero(); +// 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; +// } +// +// std::cout << GridLogMessage << "Omega = " << Omega()() << std::endl; +// OmegaLatt = OmegaLatt + (identity*Omega); +// +// auto diff = OmegaLatt - (in*OmegaLatt*transpose(in)); +// auto sdiff = sum(diff); +// //assert( norm2(sdiff) < 1e-8 ); +// if (norm2(sdiff) < 1e-8) +// { +// std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl; +// } else { +// std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT left invariant by " << norm2(sdiff) << std::endl; +// } +// +// +// } +// +// +// }; // end of class Sp +// +// +// template +// static void ProjectSp2n(Lattice > > > &Umu) +// { +// Umu = ProjectOnSpGroup(Umu); +// auto det = Determinant(Umu); // ok ? +// +// det = conjugate(det); +// +// for(int i=0;i(Umu,N-1,i); +// element = element * det; +// PokeIndex(Umu,element,Nc-1,i); +// } +// } +// template +// static void ProjectSp2n(Lattice >,Nd> > &U) +// { +// GridBase *grid=U.Grid(); +// for(int mu=0;mu(U,mu); +// Umu = ProjectOnSpGroup(Umu); +// ProjectSp2n(Umu); +// PokeIndex(U,Umu,mu); +// } +// } +// +// typedef Sp<2> Sp2; +// typedef Sp<4> Sp4; +// typedef Sp<6> Sp6; +// typedef Sp<8> Sp8; +// +// NAMESPACE_END(Grid); +// #endif diff --git a/Grid/qcd/utils/Sp2n_impl.h b/Grid/qcd/utils/Sp2n_impl.h new file mode 100644 index 00000000..03a42d70 --- /dev/null +++ b/Grid/qcd/utils/Sp2n_impl.h @@ -0,0 +1,245 @@ + +template + static int su2subgroups(void) { return (nsp * (nsp - 1)) / 2; } +template +using iSp2nMatrix = iGroupMatrix; +template +using iSp2nAlgebraVector = iAlgebraVector; + + // Sp(2N) has N(2N+1) = 2N^2+N generators + // + // normalise the generators such that + // Trace ( Ta Tb) = 1/2 delta_ab + // + // N generators in the cartan, 2N^2 off + // off diagonal: + // there are 6 types named a,b,c,d and w,z + // abcd are N(N-1)/2 each while wz are N each + + template + static void generator(int lieIndex, iSp2nMatrix &ta) { + // map lie index into type of generators: diagonal, abcd type, wz type + + int diagIndex; + int aIndex, bIndex, cIndex, dIndex; + int wIndex, zIndex; // a,b,c,d are N(N-1)/2 and w,z are N + int mod = nsp * (nsp-1) * 0.5; + int offdiag = 2*nsp*nsp; // number of generators not in the cartan subalgebra + int wmod = 4*mod; + int zmod = wmod+nsp; + if (lieIndex >= offdiag) { + diagIndex = lieIndex - offdiag; // 0, ... ,N-1 + //std::cout << GridLogMessage << "diag type " << std::endl; + generatorDiagtype(diagIndex, ta); + return; + } + if ( (lieIndex >= wmod) && (lieIndex < zmod) ) { + //std::cout << GridLogMessage << "w type " << std::endl; + wIndex = lieIndex- wmod; // 0, ... ,N-1 + generatorWtype(wIndex,ta); + return; + } + if ( (lieIndex >= zmod) && (lieIndex < offdiag) ) { + //std::cout << GridLogMessage << "z type " << std::endl; + //std::cout << GridLogMessage << "lie index " << lieIndex << std::endl; + //std::cout << GridLogMessage << "z mod " << zmod << std::endl; + zIndex = lieIndex - zmod; // 0, ... ,N-1 + generatorZtype(zIndex,ta); + return; + } + if (lieIndex < mod) { // atype 0, ... , N(N-1)/2=mod + //std::cout << GridLogMessage << "a type " << std::endl; + aIndex = lieIndex; + //std::cout << GridLogMessage << "a indx " << aIndex << std::endl; + generatorAtype(aIndex, ta); + return; + } + if ( (lieIndex >= mod) && lieIndex < 2*mod) { // btype mod, ... , 2mod-1 + //std::cout << GridLogMessage << "b type " << std::endl; + bIndex = lieIndex - mod; + generatorBtype(bIndex, ta); + return; + } + if ( (lieIndex >= 2*mod) && lieIndex < 3*mod) { // ctype 2mod, ... , 3mod-1 + //std::cout << GridLogMessage << "c type " << std::endl; + cIndex = lieIndex - 2*mod; + generatorCtype(cIndex, ta); + return; + } + if ( (lieIndex >= 3*mod) && lieIndex < wmod) { // ctype 3mod, ... , 4mod-1 = wmod-1 + //std::cout << GridLogMessage << "d type " << std::endl; + dIndex = lieIndex - 3*mod; + generatorDtype(dIndex, ta); + return; + } + + } //end of generator + + template + static void generatorDiagtype(int diagIndex, iSp2nMatrix &ta) { + + // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra + + ta = Zero(); + RealD nrm = 1.0 / 2; + + ta()()(diagIndex,diagIndex) = nrm; + ta()()(diagIndex+nsp,diagIndex+nsp) = -nrm; + } + + template + static void generatorAtype(int aIndex, iSp2nMatrix &ta) { + + // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2) + // with i + static void generatorBtype(int bIndex, iSp2nMatrix &ta) { + + // ta(i,j) = -ta(j,i) = ta(i+N,j+N) = -ta(j+N,i+N) = i / 1/ 2 sqrt(2) + // with i + static void generatorCtype(int cIndex, iSp2nMatrix &ta) { + + // ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2) + + + int i1, i2; + ta = Zero(); + RealD nrm = 1 / (2 * std::sqrt(2) ); + su2SubGroupIndex(i1, i2, cIndex); + + ta()()(i1,i2+nsp) = 1; + ta()()(i2,i1+nsp) = 1; + ta()()(i1+nsp,i2) = 1; + ta()()(i2+nsp,i1) = 1; + + ta = ta * nrm; + } + + template + static void generatorDtype(int dIndex, iSp2nMatrix &ta) { + + // ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2) + + int i1, i2; + ta = Zero(); + cplx i(0.0, 1.0); + RealD nrm = 1 / (2 * std::sqrt(2) ); + su2SubGroupIndex(i1, i2, dIndex); + + ta()()(i1,i2+nsp) = i; + ta()()(i2,i1+nsp) = i; + ta()()(i1+nsp,i2) = -i; + ta()()(i2+nsp,i1) = -i; + + ta = ta * nrm; + } + + template + static void generatorWtype(int wIndex, iSp2nMatrix &ta) { + + // ta(i,i+N) = ta(i+N,i) = 1/2 + + ta = Zero(); + RealD nrm = 1.0 / 2; //check + + ta()()(wIndex,wIndex+nsp) = 1; + ta()()(wIndex+nsp,wIndex) = 1; + + ta = ta * nrm; + } + + template + static void generatorZtype(int zIndex, iSp2nMatrix &ta) { + + // ta(i,i+N) = - ta(i+N,i) = i/2 + + ta = Zero(); + RealD nrm = 1.0 / 2; //check + cplx i(0.0, 1.0); + ta()()(zIndex,zIndex+nsp) = i; + ta()()(zIndex+nsp,zIndex) = -i; + + ta = ta * nrm; + } + + //////////////////////////////////////////////////////////////////////// + // Map a su2 subgroup number to the pair of rows that are non zero + //////////////////////////////////////////////////////////////////////// + template + static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { + assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2)); + + int spare = su2_index; + for (i1 = 0; spare >= (nsp - 1 - i1); i1++) { + spare = spare - (nsp - 1 - i1); // remove the Nc-1-i1 terms + } + i2 = i1 + 1 + spare; + } + + template + static void testGenerators(void) { + Matrix ta; + Matrix tb; + std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is 0.5 delta_ab " << std::endl; + for (int a = 0; a < AlgebraDimension; a++) { + for (int b = 0; b < AlgebraDimension; b++) { + generator(a,ta); + generator(b,tb); + Complex tr = TensorRemove(trace( ta * tb) ); + std::cout << GridLogMessage << "(" << a << "," << b << ") = " << tr + << std::endl; + if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6); + if (a != b) assert(abs(tr) < 1.0e-6); + + } + } + std::cout << GridLogMessage << std::endl; + std::cout << GridLogMessage << "Fundamental - Checking if hermitian" << std::endl; + for (int a = 0; a < AlgebraDimension; a++) { + generator(a,ta); + std::cout << GridLogMessage << a << std::endl; + assert(norm2(ta - adj(ta)) < 1.0e-6); + } + std::cout << GridLogMessage << std::endl; + std::cout << GridLogMessage << "Fundamental - Checking if traceless" << std::endl; + for (int a = 0; a < AlgebraDimension; a++) { + generator(a, ta); + Complex tr = TensorRemove(trace(ta)); + std::cout << GridLogMessage << a << std::endl; + assert(abs(tr) < 1.0e-6); + } + + } + From 2507606bd0771fc14e8b859615d6ecce9af494b7 Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Wed, 30 Nov 2022 12:52:33 +0000 Subject: [PATCH 10/17] With function overloading (still dirty). --- Grid/qcd/utils/SUn.h | 62 +++-- Grid/qcd/utils/SUn_impl.h | 16 +- Grid/qcd/utils/Sp2n.h | 400 ++++++++++++++------------- Grid/qcd/utils/Sp2n_impl.h | 536 +++++++++++++++++++++---------------- 4 files changed, 563 insertions(+), 451 deletions(-) diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index 0c961ab4..6b7e9079 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -32,15 +32,17 @@ directory #ifndef QCD_UTIL_SUN_H #define QCD_UTIL_SUN_H -#define ONLY_IF_SU \ - typename \ - typename std::enable_if_t::value, dummy_name> ,\ - dummy_name = group_name, +#define ONLY_IF_SU \ + typename dummy_name = group_name, \ + typename = std::enable_if_t < \ + std::is_same::value && \ + is_su::value > -#define ONLY_IF_Sp \ - typename \ - typename std::enable_if_t::value, dummy_name> ,\ - dummy_name =group_name, +#define ONLY_IF_Sp \ + typename dummy_name = group_name, \ + typename = std::enable_if_t < \ + std::is_same::value && \ + is_sp::value > NAMESPACE_BEGIN(Grid); namespace GroupName { @@ -68,15 +70,18 @@ struct is_sp { static const bool value = true; }; -template +template constexpr int compute_adjoint_dimension(int ncolour); -template<> -constexpr int compute_adjoint_dimension(int ncolour) { return ncolour * ncolour - 1;} - -template<> -constexpr int compute_adjoint_dimension(int ncolour) { return ncolour/2*(ncolour+1);} +template <> +constexpr int compute_adjoint_dimension(int ncolour) { + return ncolour * ncolour - 1; +} +template <> +constexpr int compute_adjoint_dimension(int ncolour) { + return ncolour / 2 * (ncolour + 1); +} template class GaugeGroup; @@ -91,11 +96,12 @@ template class GaugeGroup { public: static const int Dimension = ncolour; - static const int AdjointDimension = compute_adjoint_dimension(ncolour); - static const int AlgebraDimension = compute_adjoint_dimension(ncolour); + static const int AdjointDimension = + compute_adjoint_dimension(ncolour); + static const int AlgebraDimension = + compute_adjoint_dimension(ncolour); // Don't know how to only enable this for Sp: - static const int nsp = ncolour/2; - + static const int nsp = ncolour / 2; template using iSU2Matrix = iScalar > >; @@ -103,6 +109,15 @@ class GaugeGroup { using iGroupMatrix = iScalar > >; template using iAlgebraVector = iScalar > >; + static int su2subgroups(void) { return su2subgroups(group_name()); } + template + using iSUnMatrix = iGroupMatrix; + template + using iSUnAlgebraVector = iAlgebraVector; + template + using iSp2nMatrix = iGroupMatrix; + template + using iSp2nAlgebraVector = iAlgebraVector; ////////////////////////////////////////////////////////////////////////////////////////////////// // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix, @@ -150,6 +165,17 @@ class GaugeGroup { #include "Grid/qcd/utils/SUn_impl.h" #include "Grid/qcd/utils/Sp2n_impl.h" + template + static void generator(int lieIndex, iSp2nMatrix &ta) { + return generator(lieIndex, ta, group_name()); + } + + static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { + return su2SubGroupIndex(i1, i2, su2_index, group_name()); + } + + static void testGenerators(void) { testGenerators(group_name()); } + static void printGenerators(void) { for (int gen = 0; gen < AdjointDimension; gen++) { Matrix ta; diff --git a/Grid/qcd/utils/SUn_impl.h b/Grid/qcd/utils/SUn_impl.h index 30d8ff14..2d97dc17 100644 --- a/Grid/qcd/utils/SUn_impl.h +++ b/Grid/qcd/utils/SUn_impl.h @@ -6,12 +6,7 @@ // // around it. -template - static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; } -template -using iSUnMatrix = iGroupMatrix; -template -using iSUnAlgebraVector = iAlgebraVector; +static int su2subgroups(GroupName::SU) { return (ncolour * (ncolour - 1)) / 2; } //////////////////////////////////////////////////////////////////////// // There are N^2-1 generators for SU(N). // @@ -60,8 +55,8 @@ using iSUnAlgebraVector = iAlgebraVector; // ( -2) // //////////////////////////////////////////////////////////////////////// -template -static void generator(int lieIndex, iSUnMatrix &ta) { +template +static void generator(int lieIndex, iSUnMatrix &ta, GroupName::SU) { // map lie index to which type of generator int diagIndex; int su2Index; @@ -117,8 +112,7 @@ static void generatorDiagonal(int diagIndex, iSUnMatrix &ta) { //////////////////////////////////////////////////////////////////////// // Map a su2 subgroup number to the pair of rows that are non zero //////////////////////////////////////////////////////////////////////// -template -static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { +static void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::SU) { assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2)); int spare = su2_index; @@ -471,7 +465,7 @@ static void SubGroupHeatBath( } template -static void testGenerators(void) { +static void testGenerators(GroupName::SU) { Matrix ta; Matrix tb; std::cout << GridLogMessage diff --git a/Grid/qcd/utils/Sp2n.h b/Grid/qcd/utils/Sp2n.h index 1f7009e3..201e1e7c 100644 --- a/Grid/qcd/utils/Sp2n.h +++ b/Grid/qcd/utils/Sp2n.h @@ -1,14 +1,14 @@ + +#ifndef QCD_UTIL_Sp2n_H +#define QCD_UTIL_Sp2n_H + #include "Grid/qcd/utils/SUn.h" -// -// #ifndef QCD_UTIL_Sp2n_H -// #define QCD_UTIL_Sp2n_H -// -// NAMESPACE_BEGIN(Grid); -// +NAMESPACE_BEGIN(Grid); +// // // Sp(2N) // // ncolour = 2N -// -// +// +// // template // class Sp { // public: @@ -16,41 +16,42 @@ // static const int Dimension = ncolour; // static const int AlgebraDimension = nsp*(2*nsp +1); // static int su2subgroups(void) { return (nsp * (nsp - 1)) / 2; } -// -// +// +// // template // using iSp2nMatrix = iScalar > >; // template // using iSU2Matrix = iScalar > >; // template -// using iSp2nAlgebraVector = iScalar > >; -// +// using iSp2nAlgebraVector = iScalar > >; +// // typedef iSp2nMatrix Matrix; // typedef iSp2nMatrix MatrixF; // typedef iSp2nMatrix MatrixD; -// +// // typedef iSp2nMatrix vMatrix; // typedef iSp2nMatrix vMatrixF; // typedef iSp2nMatrix vMatrixD; -// +// // typedef iSp2nAlgebraVector AlgebraVector; // typedef iSp2nAlgebraVector AlgebraVectorF; // typedef iSp2nAlgebraVector AlgebraVectorD; -// +// // typedef iSp2nAlgebraVector vAlgebraVector; // typedef iSp2nAlgebraVector vAlgebraVectorF; // typedef iSp2nAlgebraVector vAlgebraVectorD; -// +// // typedef Lattice LatticeMatrix; // typedef Lattice LatticeMatrixF; // typedef Lattice LatticeMatrixD; -// +// // typedef Lattice LatticeAlgebraVector; // typedef Lattice LatticeAlgebraVectorF; // typedef Lattice LatticeAlgebraVectorD; -// -// -// +// +// +// // // Sp(2N) has N(2N+1) = 2N^2+N generators // // // // normalise the generators such that @@ -60,19 +61,19 @@ // // off diagonal: // // there are 6 types named a,b,c,d and w,z // // abcd are N(N-1)/2 each while wz are N each -// +// // template // static void generator(int lieIndex, iSp2nMatrix &ta) { -// // map lie index into type of generators: diagonal, abcd type, wz type -// +// // map lie index into type of generators: diagonal, abcd type, wz +// type +// // int diagIndex; // int aIndex, bIndex, cIndex, dIndex; // int wIndex, zIndex; // a,b,c,d are N(N-1)/2 and w,z are N // int mod = nsp * (nsp-1) * 0.5; -// int offdiag = 2*nsp*nsp; // number of generators not in the cartan subalgebra -// int wmod = 4*mod; -// int zmod = wmod+nsp; -// if (lieIndex >= offdiag) { +// int offdiag = 2*nsp*nsp; // number of generators not in the cartan +// subalgebra int wmod = 4*mod; int zmod = wmod+nsp; if (lieIndex >= +// offdiag) { // diagIndex = lieIndex - offdiag; // 0, ... ,N-1 // //std::cout << GridLogMessage << "diag type " << std::endl; // generatorDiagtype(diagIndex, ta); @@ -86,7 +87,8 @@ // } // if ( (lieIndex >= zmod) && (lieIndex < offdiag) ) { // //std::cout << GridLogMessage << "z type " << std::endl; -// //std::cout << GridLogMessage << "lie index " << lieIndex << std::endl; +// //std::cout << GridLogMessage << "lie index " << lieIndex << +// std::endl; // //std::cout << GridLogMessage << "z mod " << zmod << std::endl; // zIndex = lieIndex - zmod; // 0, ... ,N-1 // generatorZtype(zIndex,ta); @@ -95,168 +97,171 @@ // if (lieIndex < mod) { // atype 0, ... , N(N-1)/2=mod // //std::cout << GridLogMessage << "a type " << std::endl; // aIndex = lieIndex; -// //std::cout << GridLogMessage << "a indx " << aIndex << std::endl; -// generatorAtype(aIndex, ta); -// return; +// //std::cout << GridLogMessage << "a indx " << aIndex << +// std::endl; generatorAtype(aIndex, ta); return; // } -// if ( (lieIndex >= mod) && lieIndex < 2*mod) { // btype mod, ... , 2mod-1 +// if ( (lieIndex >= mod) && lieIndex < 2*mod) { // btype mod, ... , +// 2mod-1 // //std::cout << GridLogMessage << "b type " << std::endl; // bIndex = lieIndex - mod; // generatorBtype(bIndex, ta); // return; // } -// if ( (lieIndex >= 2*mod) && lieIndex < 3*mod) { // ctype 2mod, ... , 3mod-1 +// if ( (lieIndex >= 2*mod) && lieIndex < 3*mod) { // ctype 2mod, ... , +// 3mod-1 // //std::cout << GridLogMessage << "c type " << std::endl; // cIndex = lieIndex - 2*mod; // generatorCtype(cIndex, ta); // return; // } -// if ( (lieIndex >= 3*mod) && lieIndex < wmod) { // ctype 3mod, ... , 4mod-1 = wmod-1 +// if ( (lieIndex >= 3*mod) && lieIndex < wmod) { // ctype 3mod, ... , +// 4mod-1 = wmod-1 // //std::cout << GridLogMessage << "d type " << std::endl; // dIndex = lieIndex - 3*mod; // generatorDtype(dIndex, ta); // return; // } -// +// // } //end of generator -// +// // template // static void generatorDiagtype(int diagIndex, iSp2nMatrix &ta) { -// -// // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra -// +// +// // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan +// subalgebra +// // ta = Zero(); // RealD nrm = 1.0 / 2; -// +// // ta()()(diagIndex,diagIndex) = nrm; // ta()()(diagIndex+nsp,diagIndex+nsp) = -nrm; // } -// +// // template // static void generatorAtype(int aIndex, iSp2nMatrix &ta) { -// +// // // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2) // // with i // static void generatorBtype(int bIndex, iSp2nMatrix &ta) { -// +// // // ta(i,j) = -ta(j,i) = ta(i+N,j+N) = -ta(j+N,i+N) = i / 1/ 2 sqrt(2) // // with i // static void generatorCtype(int cIndex, iSp2nMatrix &ta) { -// +// // // ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2) -// -// +// +// // int i1, i2; // ta = Zero(); // RealD nrm = 1 / (2 * std::sqrt(2) ); // su2SubGroupIndex(i1, i2, cIndex); -// +// // ta()()(i1,i2+nsp) = 1; // ta()()(i2,i1+nsp) = 1; // ta()()(i1+nsp,i2) = 1; // ta()()(i2+nsp,i1) = 1; -// +// // ta = ta * nrm; // } -// +// // template // static void generatorDtype(int dIndex, iSp2nMatrix &ta) { -// +// // // ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2) -// +// // int i1, i2; // ta = Zero(); // cplx i(0.0, 1.0); // RealD nrm = 1 / (2 * std::sqrt(2) ); // su2SubGroupIndex(i1, i2, dIndex); -// +// // ta()()(i1,i2+nsp) = i; // ta()()(i2,i1+nsp) = i; // ta()()(i1+nsp,i2) = -i; // ta()()(i2+nsp,i1) = -i; -// +// // ta = ta * nrm; // } -// +// // template // static void generatorWtype(int wIndex, iSp2nMatrix &ta) { -// +// // // ta(i,i+N) = ta(i+N,i) = 1/2 -// +// // ta = Zero(); // RealD nrm = 1.0 / 2; //check -// +// // ta()()(wIndex,wIndex+nsp) = 1; // ta()()(wIndex+nsp,wIndex) = 1; -// +// // ta = ta * nrm; // } -// +// // template // static void generatorZtype(int zIndex, iSp2nMatrix &ta) { -// +// // // ta(i,i+N) = - ta(i+N,i) = i/2 -// +// // ta = Zero(); // RealD nrm = 1.0 / 2; //check // cplx i(0.0, 1.0); // ta()()(zIndex,zIndex+nsp) = i; // ta()()(zIndex+nsp,zIndex) = -i; -// +// // ta = ta * nrm; // } -// -// +// +// // //////////////////////////////////////////////////////////////////////// // // Map a su2 subgroup number to the pair of rows that are non zero // //////////////////////////////////////////////////////////////////////// // static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { // assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2)); -// +// // int spare = su2_index; // for (i1 = 0; spare >= (nsp - 1 - i1); i1++) { // spare = spare - (nsp - 1 - i1); // remove the Nc-1-i1 terms // } // i2 = i1 + 1 + spare; // } -// -// -// -// -// +// +// +// +// +// // static void printGenerators(void) { // for (int gen = 0; gen < AlgebraDimension; gen++) { // Matrix ta; @@ -266,84 +271,89 @@ // std::cout << GridLogMessage << ta << std::endl; // } // } -// -// -// +// +// +// // static void testGenerators(void) { // Matrix ta; // Matrix tb; -// std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is 0.5 delta_ab " << std::endl; -// for (int a = 0; a < AlgebraDimension; a++) { +// std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is +// 0.5 delta_ab " << std::endl; for (int a = 0; a < AlgebraDimension; +// a++) { // for (int b = 0; b < AlgebraDimension; b++) { // generator(a,ta); // generator(b,tb); // Complex tr = TensorRemove(trace( ta * tb) ); -// std::cout << GridLogMessage << "(" << a << "," << b << ") = " << tr +// std::cout << GridLogMessage << "(" << a << "," << b << ") = +// " << tr // << std::endl; // if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6); // if (a != b) assert(abs(tr) < 1.0e-6); -// +// // } // } // std::cout << GridLogMessage << std::endl; -// std::cout << GridLogMessage << "Fundamental - Checking if hermitian" << std::endl; -// for (int a = 0; a < AlgebraDimension; a++) { +// std::cout << GridLogMessage << "Fundamental - Checking if hermitian" +// << std::endl; for (int a = 0; a < AlgebraDimension; a++) { // generator(a,ta); // std::cout << GridLogMessage << a << std::endl; // assert(norm2(ta - adj(ta)) < 1.0e-6); // } // std::cout << GridLogMessage << std::endl; -// std::cout << GridLogMessage << "Fundamental - Checking if traceless" << std::endl; -// for (int a = 0; a < AlgebraDimension; a++) { +// std::cout << GridLogMessage << "Fundamental - Checking if traceless" +// << std::endl; for (int a = 0; a < AlgebraDimension; a++) { // generator(a, ta); // Complex tr = TensorRemove(trace(ta)); // std::cout << GridLogMessage << a << std::endl; // assert(abs(tr) < 1.0e-6); // } -// +// // } -// -// +// +// // template -// static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, double scale = 1.0) +// 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::scalar_object MatrixType; -// +// typedef typename GridTypeMapper::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 < AlgebraDimension; 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, //same as sun +// +// +// static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG, +// //same as sun // LatticeMatrix &out, // Real scale = 1.0) { // GridBase *grid = out.Grid(); @@ -351,7 +361,7 @@ // LatticeMatrix la(grid); // Complex ci(0.0, scale); // Matrix ta; -// +// // out = Zero(); // for (int a = 0; a < AlgebraDimension; a++) { // gaussian(pRNG, ca); @@ -361,7 +371,7 @@ // } // out *= ci; // } -// +// // static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h, // LatticeMatrix &out, // Real scale = 1.0) { @@ -369,7 +379,7 @@ // GridBase *grid = out.Grid(); // LatticeMatrix la(grid); // Matrix ta; -// +// // out = Zero(); // for (int a = 0; a < AlgebraDimension; a++) { // generator(a, ta); @@ -377,18 +387,19 @@ // out += la; // } // } -// -// +// +// // template -// static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { // same as sun +// static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { // +// same as sun // typedef typename LatticeMatrixType::scalar_type ComplexType; -// +// // LatticeMatrixType xn(x.Grid()); // RealD nfac = 1.0; -// +// // xn = x; // ex = xn + ComplexType(1.0); // 1+x -// +// // // Do a 12th order exponentiation // for (int i = 2; i <= 12; ++i) { // nfac = nfac / RealD(i); // 1/2, 1/2.3 ... @@ -396,25 +407,26 @@ // ex = ex + xn * nfac; // x2/2!, x3/3!.... // } // } -// -// static void projectOnAlgebra(LatticeAlgebraVector &h_out, const LatticeMatrix &in, Real scale = 1.0) { +// +// 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 < AlgebraDimension; a++) { // generator(a, Ta); // pokeColour(h_out, - 2.0 * (trace(timesI(Ta) * in)) * scale, a); // } // } -// -// +// +// // template // static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) { // typedef typename GaugeField::vector_type vector_type; // typedef iSp2nMatrix vMatrixType; // typedef Lattice LatticeMatrixType; -// +// // LatticeMatrixType Umu(out.Grid()); // for (int mu = 0; mu < Nd; mu++) { // LieRandomize(pRNG, Umu, 1.0); //def @@ -426,20 +438,20 @@ // typedef typename GaugeField::vector_type vector_type; // typedef iSp2nMatrix vMatrixType; // typedef Lattice LatticeMatrixType; -// +// // LatticeMatrixType Umu(out.Grid()); // for(int mu=0;mu(out,Umu,mu); // } // } -// +// // template // static void ColdConfiguration(GaugeField &out){ // typedef typename GaugeField::vector_type vector_type; // typedef iSp2nMatrix vMatrixType; // typedef Lattice LatticeMatrixType; -// +// // LatticeMatrixType Umu(out.Grid()); // Umu=1.0; // for(int mu=0;mu // static void OmegaInvariance(GaugeField &in) // { // typedef typename GaugeField::vector_type vector_type; // typedef iSp2nMatrix vMatrixType; // typedef Lattice LatticeMatrixType; -// +// // LatticeMatrixType U(in.Grid()); // LatticeMatrixType OmegaLatt(in.Grid()); // LatticeMatrixType identity(in.Grid()); // RealD vol = in.Grid()->gSites(); // ColourMatrix Omega; -// +// // OmegaLatt = Zero(); // Omega = Zero(); // identity = 1.; -// +// // std::cout << GridLogMessage << "I am a GaugeField " << std::endl; -// +// // U = PeekIndex(in,1); -// +// // OmegaInvariance(U); -// +// // } -// +// // static void OmegaInvariance(LatticeColourMatrixD &in) // { -// +// // LatticeColourMatrixD OmegaLatt(in.Grid()); // LatticeColourMatrixD identity(in.Grid()); // RealD vol = in.Grid()->gSites(); // ColourMatrix Omega; -// +// // OmegaLatt = Zero(); // Omega = Zero(); // identity = 1.; -// -// std::cout << GridLogMessage << "I am a LatticeColourMatrix " << std::endl; -// +// +// 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; // } -// -// std::cout << GridLogMessage << "Omega = " << Omega()() << std::endl; -// OmegaLatt = OmegaLatt + (identity*Omega); -// +// +// std::cout << GridLogMessage << "Omega = " << Omega()() << +// std::endl; OmegaLatt = OmegaLatt + (identity*Omega); +// // auto diff = OmegaLatt - (in*OmegaLatt*transpose(in)); // auto sdiff = sum(diff); // //assert( norm2(sdiff) < 1e-8 ); // if (norm2(sdiff) < 1e-8) // { -// std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl; +// std::cout << GridLogMessage << "Symplectic condition +// satisfied: Omega invariant" << std::endl; // } else { -// std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT left invariant by " << norm2(sdiff) << std::endl; +// std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT +// left invariant by " << norm2(sdiff) << std::endl; // } -// -// +// +// // } -// -// +// +// // }; // end of class Sp -// -// -// template -// static void ProjectSp2n(Lattice > > > &Umu) -// { -// Umu = ProjectOnSpGroup(Umu); -// auto det = Determinant(Umu); // ok ? -// -// det = conjugate(det); -// -// for(int i=0;i(Umu,N-1,i); -// element = element * det; -// PokeIndex(Umu,element,Nc-1,i); -// } -// } -// template -// static void ProjectSp2n(Lattice >,Nd> > &U) -// { -// GridBase *grid=U.Grid(); -// for(int mu=0;mu(U,mu); -// Umu = ProjectOnSpGroup(Umu); -// ProjectSp2n(Umu); -// PokeIndex(U,Umu,mu); -// } -// } -// -// typedef Sp<2> Sp2; -// typedef Sp<4> Sp4; -// typedef Sp<6> Sp6; -// typedef Sp<8> Sp8; -// -// NAMESPACE_END(Grid); -// #endif + +template +static void ProjectSp2n( + Lattice > > > &Umu) { + Umu = ProjectOnSpGroup(Umu); + auto det = Determinant(Umu); // ok ? + + 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 ProjectSp2n( + Lattice >, Nd> > &U) { + GridBase *grid = U.Grid(); + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnSpGroup(Umu); + ProjectSp2n(Umu); + PokeIndex(U, Umu, mu); + } +} + +typedef Sp<2> Sp2; +typedef Sp<4> Sp4; +typedef Sp<6> Sp6; +typedef Sp<8> Sp8; + +NAMESPACE_END(Grid); +#endif diff --git a/Grid/qcd/utils/Sp2n_impl.h b/Grid/qcd/utils/Sp2n_impl.h index 03a42d70..9399d5b2 100644 --- a/Grid/qcd/utils/Sp2n_impl.h +++ b/Grid/qcd/utils/Sp2n_impl.h @@ -1,245 +1,321 @@ -template - static int su2subgroups(void) { return (nsp * (nsp - 1)) / 2; } -template -using iSp2nMatrix = iGroupMatrix; -template -using iSp2nAlgebraVector = iAlgebraVector; +static int su2subgroups(GroupName::Sp) { return (nsp * (nsp - 1)) / 2; } - // Sp(2N) has N(2N+1) = 2N^2+N generators - // - // normalise the generators such that - // Trace ( Ta Tb) = 1/2 delta_ab - // - // N generators in the cartan, 2N^2 off - // off diagonal: - // there are 6 types named a,b,c,d and w,z - // abcd are N(N-1)/2 each while wz are N each - - template - static void generator(int lieIndex, iSp2nMatrix &ta) { - // map lie index into type of generators: diagonal, abcd type, wz type +// Sp(2N) has N(2N+1) = 2N^2+N generators +// +// normalise the generators such that +// Trace ( Ta Tb) = 1/2 delta_ab +// +// N generators in the cartan, 2N^2 off +// off diagonal: +// there are 6 types named a,b,c,d and w,z +// abcd are N(N-1)/2 each while wz are N each - int diagIndex; - int aIndex, bIndex, cIndex, dIndex; - int wIndex, zIndex; // a,b,c,d are N(N-1)/2 and w,z are N - int mod = nsp * (nsp-1) * 0.5; - int offdiag = 2*nsp*nsp; // number of generators not in the cartan subalgebra - int wmod = 4*mod; - int zmod = wmod+nsp; - if (lieIndex >= offdiag) { - diagIndex = lieIndex - offdiag; // 0, ... ,N-1 - //std::cout << GridLogMessage << "diag type " << std::endl; - generatorDiagtype(diagIndex, ta); - return; - } - if ( (lieIndex >= wmod) && (lieIndex < zmod) ) { - //std::cout << GridLogMessage << "w type " << std::endl; - wIndex = lieIndex- wmod; // 0, ... ,N-1 - generatorWtype(wIndex,ta); - return; - } - if ( (lieIndex >= zmod) && (lieIndex < offdiag) ) { - //std::cout << GridLogMessage << "z type " << std::endl; - //std::cout << GridLogMessage << "lie index " << lieIndex << std::endl; - //std::cout << GridLogMessage << "z mod " << zmod << std::endl; - zIndex = lieIndex - zmod; // 0, ... ,N-1 - generatorZtype(zIndex,ta); - return; - } - if (lieIndex < mod) { // atype 0, ... , N(N-1)/2=mod - //std::cout << GridLogMessage << "a type " << std::endl; - aIndex = lieIndex; - //std::cout << GridLogMessage << "a indx " << aIndex << std::endl; - generatorAtype(aIndex, ta); - return; - } - if ( (lieIndex >= mod) && lieIndex < 2*mod) { // btype mod, ... , 2mod-1 - //std::cout << GridLogMessage << "b type " << std::endl; - bIndex = lieIndex - mod; - generatorBtype(bIndex, ta); - return; - } - if ( (lieIndex >= 2*mod) && lieIndex < 3*mod) { // ctype 2mod, ... , 3mod-1 - //std::cout << GridLogMessage << "c type " << std::endl; - cIndex = lieIndex - 2*mod; - generatorCtype(cIndex, ta); - return; - } - if ( (lieIndex >= 3*mod) && lieIndex < wmod) { // ctype 3mod, ... , 4mod-1 = wmod-1 - //std::cout << GridLogMessage << "d type " << std::endl; - dIndex = lieIndex - 3*mod; - generatorDtype(dIndex, ta); - return; - } - - } //end of generator - - template - static void generatorDiagtype(int diagIndex, iSp2nMatrix &ta) { - - // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra - - ta = Zero(); - RealD nrm = 1.0 / 2; +template +static void generator(int lieIndex, iSp2nMatrix &ta, GroupName::Sp) { + // map lie index into type of generators: diagonal, abcd type, wz type - ta()()(diagIndex,diagIndex) = nrm; - ta()()(diagIndex+nsp,diagIndex+nsp) = -nrm; - } - - template - static void generatorAtype(int aIndex, iSp2nMatrix &ta) { - - // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2) - // with i= offdiag) { + diagIndex = lieIndex - offdiag; // 0, ... ,N-1 + // std::cout << GridLogMessage << "diag type " << std::endl; + generatorDiagtype(diagIndex, ta); + return; + } + if ((lieIndex >= wmod) && (lieIndex < zmod)) { + // std::cout << GridLogMessage << "w type " << std::endl; + wIndex = lieIndex - wmod; // 0, ... ,N-1 + generatorWtype(wIndex, ta); + return; + } + if ((lieIndex >= zmod) && (lieIndex < offdiag)) { + // std::cout << GridLogMessage << "z type " << std::endl; + // std::cout << GridLogMessage << "lie index " << lieIndex << std::endl; + // std::cout << GridLogMessage << "z mod " << zmod << std::endl; + zIndex = lieIndex - zmod; // 0, ... ,N-1 + generatorZtype(zIndex, ta); + return; + } + if (lieIndex < mod) { // atype 0, ... , N(N-1)/2=mod + // std::cout << GridLogMessage << "a type " << std::endl; + aIndex = lieIndex; + // std::cout << GridLogMessage << "a indx " << aIndex << std::endl; + generatorAtype(aIndex, ta); + return; + } + if ((lieIndex >= mod) && lieIndex < 2 * mod) { // btype mod, ... , 2mod-1 + // std::cout << GridLogMessage << "b type " << std::endl; + bIndex = lieIndex - mod; + generatorBtype(bIndex, ta); + return; + } + if ((lieIndex >= 2 * mod) && + lieIndex < 3 * mod) { // ctype 2mod, ... , 3mod-1 + // std::cout << GridLogMessage << "c type " << std::endl; + cIndex = lieIndex - 2 * mod; + generatorCtype(cIndex, ta); + return; + } + if ((lieIndex >= 3 * mod) && + lieIndex < wmod) { // ctype 3mod, ... , 4mod-1 = wmod-1 + // std::cout << GridLogMessage << "d type " << std::endl; + dIndex = lieIndex - 3 * mod; + generatorDtype(dIndex, ta); + return; + } - ta = ta * nrm; - } - - template - static void generatorBtype(int bIndex, iSp2nMatrix &ta) { - - // ta(i,j) = -ta(j,i) = ta(i+N,j+N) = -ta(j+N,i+N) = i / 1/ 2 sqrt(2) - // with i +static void generatorDiagtype(int diagIndex, iSp2nMatrix &ta) { + // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra - ta = ta * nrm; - } - - template - static void generatorCtype(int cIndex, iSp2nMatrix &ta) { - - // ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2) - - - int i1, i2; - ta = Zero(); - RealD nrm = 1 / (2 * std::sqrt(2) ); - su2SubGroupIndex(i1, i2, cIndex); - - ta()()(i1,i2+nsp) = 1; - ta()()(i2,i1+nsp) = 1; - ta()()(i1+nsp,i2) = 1; - ta()()(i2+nsp,i1) = 1; + ta = Zero(); + RealD nrm = 1.0 / 2; - ta = ta * nrm; - } - - template - static void generatorDtype(int dIndex, iSp2nMatrix &ta) { - - // ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2) - - int i1, i2; - ta = Zero(); - cplx i(0.0, 1.0); - RealD nrm = 1 / (2 * std::sqrt(2) ); - su2SubGroupIndex(i1, i2, dIndex); + ta()()(diagIndex, diagIndex) = nrm; + ta()()(diagIndex + nsp, diagIndex + nsp) = -nrm; +} - ta()()(i1,i2+nsp) = i; - ta()()(i2,i1+nsp) = i; - ta()()(i1+nsp,i2) = -i; - ta()()(i2+nsp,i1) = -i; +template +static void generatorAtype(int aIndex, iSp2nMatrix &ta) { + // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2) + // with i - static void generatorWtype(int wIndex, iSp2nMatrix &ta) { - - // ta(i,i+N) = ta(i+N,i) = 1/2 - - ta = Zero(); - RealD nrm = 1.0 / 2; //check + su2SubGroupIndex(i1, i2, aIndex); + ta()()(i1, i2) = 1; + ta()()(i2, i1) = 1; + ta()()(i1 + nsp, i2 + nsp) = -1; + ta()()(i2 + nsp, i1 + nsp) = -1; - ta()()(wIndex,wIndex+nsp) = 1; - ta()()(wIndex+nsp,wIndex) = 1; - - ta = ta * nrm; - } - - template - static void generatorZtype(int zIndex, iSp2nMatrix &ta) { - - // ta(i,i+N) = - ta(i+N,i) = i/2 - - ta = Zero(); - RealD nrm = 1.0 / 2; //check - cplx i(0.0, 1.0); - ta()()(zIndex,zIndex+nsp) = i; - ta()()(zIndex+nsp,zIndex) = -i; - - ta = ta * nrm; - } - - //////////////////////////////////////////////////////////////////////// - // Map a su2 subgroup number to the pair of rows that are non zero - //////////////////////////////////////////////////////////////////////// - template - static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { - assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2)); + ta = ta * nrm; +} - int spare = su2_index; - for (i1 = 0; spare >= (nsp - 1 - i1); i1++) { - spare = spare - (nsp - 1 - i1); // remove the Nc-1-i1 terms - } - i2 = i1 + 1 + spare; - } - - template - static void testGenerators(void) { - Matrix ta; - Matrix tb; - std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is 0.5 delta_ab " << std::endl; - for (int a = 0; a < AlgebraDimension; a++) { - for (int b = 0; b < AlgebraDimension; b++) { - generator(a,ta); - generator(b,tb); - Complex tr = TensorRemove(trace( ta * tb) ); - std::cout << GridLogMessage << "(" << a << "," << b << ") = " << tr +template +static void generatorBtype(int bIndex, iSp2nMatrix &ta) { + // ta(i,j) = -ta(j,i) = ta(i+N,j+N) = -ta(j+N,i+N) = i / 1/ 2 sqrt(2) + // with i +static void generatorCtype(int cIndex, iSp2nMatrix &ta) { + // ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2) + + int i1, i2; + ta = Zero(); + RealD nrm = 1 / (2 * std::sqrt(2)); + su2SubGroupIndex(i1, i2, cIndex); + + ta()()(i1, i2 + nsp) = 1; + ta()()(i2, i1 + nsp) = 1; + ta()()(i1 + nsp, i2) = 1; + ta()()(i2 + nsp, i1) = 1; + + ta = ta * nrm; +} + +template +static void generatorDtype(int dIndex, iSp2nMatrix &ta) { + // ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2) + + int i1, i2; + ta = Zero(); + cplx i(0.0, 1.0); + RealD nrm = 1 / (2 * std::sqrt(2)); + su2SubGroupIndex(i1, i2, dIndex); + + ta()()(i1, i2 + nsp) = i; + ta()()(i2, i1 + nsp) = i; + ta()()(i1 + nsp, i2) = -i; + ta()()(i2 + nsp, i1) = -i; + + ta = ta * nrm; +} + +template +static void generatorWtype(int wIndex, iSp2nMatrix &ta) { + // ta(i,i+N) = ta(i+N,i) = 1/2 + + ta = Zero(); + RealD nrm = 1.0 / 2; // check + + ta()()(wIndex, wIndex + nsp) = 1; + ta()()(wIndex + nsp, wIndex) = 1; + + ta = ta * nrm; +} + +template +static void generatorZtype(int zIndex, iSp2nMatrix &ta) { + // ta(i,i+N) = - ta(i+N,i) = i/2 + + ta = Zero(); + RealD nrm = 1.0 / 2; // check + cplx i(0.0, 1.0); + ta()()(zIndex, zIndex + nsp) = i; + ta()()(zIndex + nsp, zIndex) = -i; + + ta = ta * nrm; +} + +//////////////////////////////////////////////////////////////////////// +// Map a su2 subgroup number to the pair of rows that are non zero +//////////////////////////////////////////////////////////////////////// +template +static void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::Sp) { + assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2)); + + int spare = su2_index; + for (i1 = 0; spare >= (nsp - 1 - i1); i1++) { + spare = spare - (nsp - 1 - i1); // remove the Nc-1-i1 terms + } + i2 = i1 + 1 + spare; +} + +static void testGenerators(GroupName::Sp) { + Matrix ta; + Matrix tb; + std::cout << GridLogMessage + << "Fundamental - Checking trace ta tb is 0.5 delta_ab " + << std::endl; + for (int a = 0; a < AlgebraDimension; a++) { + for (int b = 0; b < AlgebraDimension; b++) { + generator(a, ta); + generator(b, tb); + Complex tr = TensorRemove(trace(ta * tb)); + std::cout << GridLogMessage << "(" << a << "," << b << ") = " << tr << std::endl; - if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6); - if (a != b) assert(abs(tr) < 1.0e-6); - - } - } - std::cout << GridLogMessage << std::endl; - std::cout << GridLogMessage << "Fundamental - Checking if hermitian" << std::endl; - for (int a = 0; a < AlgebraDimension; a++) { - generator(a,ta); - std::cout << GridLogMessage << a << std::endl; - assert(norm2(ta - adj(ta)) < 1.0e-6); - } - std::cout << GridLogMessage << std::endl; - std::cout << GridLogMessage << "Fundamental - Checking if traceless" << std::endl; - for (int a = 0; a < AlgebraDimension; a++) { - generator(a, ta); - Complex tr = TensorRemove(trace(ta)); - std::cout << GridLogMessage << a << std::endl; - assert(abs(tr) < 1.0e-6); - } - + if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6); + if (a != b) assert(abs(tr) < 1.0e-6); } + } + std::cout << GridLogMessage << std::endl; + std::cout << GridLogMessage << "Fundamental - Checking if hermitian" + << std::endl; + for (int a = 0; a < AlgebraDimension; a++) { + generator(a, ta); + std::cout << GridLogMessage << a << std::endl; + assert(norm2(ta - adj(ta)) < 1.0e-6); + } + std::cout << GridLogMessage << std::endl; + std::cout << GridLogMessage << "Fundamental - Checking if traceless" + << std::endl; + for (int a = 0; a < AlgebraDimension; a++) { + generator(a, ta); + Complex tr = TensorRemove(trace(ta)); + std::cout << GridLogMessage << a << std::endl; + assert(abs(tr) < 1.0e-6); + } +} + +template +static void OmegaInvariance(ColourMatrix &in) { + ColourMatrix Omega; + Omega = Zero(); + + 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; + } + + auto diff = Omega - (in * Omega * transpose(in)); + auto sdiff = norm2(diff); + if (norm2(sdiff) < 1e-8) { + std::cout << GridLogMessage + << "Symplectic condition satisfied: Omega invariant" << std::endl; + } else { + std::cout << GridLogMessage + << "WARNING!!!!!! Matrix Omega NOT left invariant by " + << norm2(sdiff) << std::endl; + } +} + +template +static void OmegaInvariance(GaugeField &in) { + typedef typename GaugeField::vector_type vector_type; + typedef iSp2nMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType U(in.Grid()); + LatticeMatrixType OmegaLatt(in.Grid()); + LatticeMatrixType identity(in.Grid()); + RealD vol = in.Grid()->gSites(); + ColourMatrix Omega; + + OmegaLatt = Zero(); + Omega = Zero(); + identity = 1.; + + std::cout << GridLogMessage << "I am a GaugeField " << std::endl; + + U = PeekIndex(in, 1); + + OmegaInvariance(U); +} + +template +static void OmegaInvariance(LatticeColourMatrixD &in) { + LatticeColourMatrixD OmegaLatt(in.Grid()); + LatticeColourMatrixD identity(in.Grid()); + RealD vol = in.Grid()->gSites(); + ColourMatrix Omega; + + OmegaLatt = Zero(); + 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; + } + + std::cout << GridLogMessage << "Omega = " << Omega()() << std::endl; + OmegaLatt = OmegaLatt + (identity * Omega); + + auto diff = OmegaLatt - (in * OmegaLatt * transpose(in)); + auto sdiff = sum(diff); + // assert( norm2(sdiff) < 1e-8 ); + if (norm2(sdiff) < 1e-8) { + std::cout << GridLogMessage + << "Symplectic condition satisfied: Omega invariant" << std::endl; + } else { + std::cout << GridLogMessage + << "WARNING!!!!!! Matrix Omega NOT left invariant by " + << norm2(sdiff) << std::endl; + } +} From 55c008da21be79c5d1d41c69be0aed1f0570be0e Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Wed, 30 Nov 2022 13:12:21 +0000 Subject: [PATCH 11/17] Removed forward declaration --- Grid/qcd/utils/SUn.h | 15 ++++++--------- Grid/qcd/utils/SUn_impl.h | 6 +++--- 2 files changed, 9 insertions(+), 12 deletions(-) diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index 6b7e9079..30ac36a3 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -83,15 +83,6 @@ constexpr int compute_adjoint_dimension(int ncolour) { return ncolour / 2 * (ncolour + 1); } -template -class GaugeGroup; - -template -using SU = GaugeGroup; - -template -using Sp = GaugeGroup; - template class GaugeGroup { public: @@ -423,6 +414,12 @@ static void ProjectSU3( } } +template +using SU = GaugeGroup; + +template +using Sp = GaugeGroup; + typedef SU<2> SU2; typedef SU<3> SU3; typedef SU<4> SU4; diff --git a/Grid/qcd/utils/SUn_impl.h b/Grid/qcd/utils/SUn_impl.h index 2d97dc17..be0eb42d 100644 --- a/Grid/qcd/utils/SUn_impl.h +++ b/Grid/qcd/utils/SUn_impl.h @@ -282,11 +282,11 @@ static void SubGroupHeatBath( SU2Matrix ident = Complex(1.0); SU2Matrix pauli1; - SU<2>::generator(0, pauli1); + GaugeGroup<2, GroupName::SU>::generator(0, pauli1); SU2Matrix pauli2; - SU<2>::generator(1, pauli2); + GaugeGroup<2, GroupName::SU>::generator(1, pauli2); SU2Matrix pauli3; - SU<2>::generator(2, pauli3); + GaugeGroup<2, GroupName::SU>::generator(2, pauli3); pauli1 = timesI(pauli1) * 2.0; pauli2 = timesI(pauli2) * 2.0; pauli3 = timesI(pauli3) * 2.0; From b8b3ae6ac14bad667a02078c205a38cb83144026 Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Wed, 30 Nov 2022 13:28:09 +0000 Subject: [PATCH 12/17] Make helper functions private --- Grid/qcd/utils/SUn.h | 1 + Grid/qcd/utils/SUn_impl.h | 2 ++ Grid/qcd/utils/Sp2n_impl.h | 2 ++ 3 files changed, 5 insertions(+) diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index 30ac36a3..a6e6c36c 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -156,6 +156,7 @@ class GaugeGroup { #include "Grid/qcd/utils/SUn_impl.h" #include "Grid/qcd/utils/Sp2n_impl.h" + public: template static void generator(int lieIndex, iSp2nMatrix &ta) { return generator(lieIndex, ta, group_name()); diff --git a/Grid/qcd/utils/SUn_impl.h b/Grid/qcd/utils/SUn_impl.h index be0eb42d..6081f79c 100644 --- a/Grid/qcd/utils/SUn_impl.h +++ b/Grid/qcd/utils/SUn_impl.h @@ -6,6 +6,7 @@ // // around it. +private: static int su2subgroups(GroupName::SU) { return (ncolour * (ncolour - 1)) / 2; } //////////////////////////////////////////////////////////////////////// // There are N^2-1 generators for SU(N). @@ -122,6 +123,7 @@ static void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::SU) { i2 = i1 + 1 + spare; } +public: ////////////////////////////////////////////////////////////////////////////////////////// // Pull out a subgroup and project on to real coeffs x pauli basis ////////////////////////////////////////////////////////////////////////////////////////// diff --git a/Grid/qcd/utils/Sp2n_impl.h b/Grid/qcd/utils/Sp2n_impl.h index 9399d5b2..addbd118 100644 --- a/Grid/qcd/utils/Sp2n_impl.h +++ b/Grid/qcd/utils/Sp2n_impl.h @@ -1,4 +1,5 @@ +private: static int su2subgroups(GroupName::Sp) { return (nsp * (nsp - 1)) / 2; } // Sp(2N) has N(2N+1) = 2N^2+N generators @@ -233,6 +234,7 @@ static void testGenerators(GroupName::Sp) { } } +public: template static void OmegaInvariance(ColourMatrix &in) { ColourMatrix Omega; From fa71b46a41c4cc6455c223aab12b9a2c324c084a Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Wed, 30 Nov 2022 14:27:19 +0000 Subject: [PATCH 13/17] Hide nsp --- Grid/qcd/utils/SUn.h | 2 -- Grid/qcd/utils/Sp2n_impl.h | 21 ++++++++++++++++----- tests/sp2n/Test_Sp_start.cc | 2 +- tests/sp2n/Test_project_on_Sp.cc | 3 +-- 4 files changed, 18 insertions(+), 10 deletions(-) diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index a6e6c36c..2cedf67d 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -91,8 +91,6 @@ class GaugeGroup { compute_adjoint_dimension(ncolour); static const int AlgebraDimension = compute_adjoint_dimension(ncolour); - // Don't know how to only enable this for Sp: - static const int nsp = ncolour / 2; template using iSU2Matrix = iScalar > >; diff --git a/Grid/qcd/utils/Sp2n_impl.h b/Grid/qcd/utils/Sp2n_impl.h index addbd118..28dac5f8 100644 --- a/Grid/qcd/utils/Sp2n_impl.h +++ b/Grid/qcd/utils/Sp2n_impl.h @@ -1,6 +1,6 @@ private: -static int su2subgroups(GroupName::Sp) { return (nsp * (nsp - 1)) / 2; } +static int su2subgroups(GroupName::Sp) { return (ncolour/2 * (ncolour/2 - 1)) / 2; } // Sp(2N) has N(2N+1) = 2N^2+N generators // @@ -16,14 +16,15 @@ template static void generator(int lieIndex, iSp2nMatrix &ta, GroupName::Sp) { // map lie index into type of generators: diagonal, abcd type, wz type + const int nsp = ncolour/2; int diagIndex; int aIndex, bIndex, cIndex, dIndex; int wIndex, zIndex; // a,b,c,d are N(N-1)/2 and w,z are N - int mod = nsp * (nsp - 1) * 0.5; - int offdiag = + const int mod = nsp * (nsp - 1) * 0.5; + const int offdiag = 2 * nsp * nsp; // number of generators not in the cartan subalgebra - int wmod = 4 * mod; - int zmod = wmod + nsp; + const int wmod = 4 * mod; + const int zmod = wmod + nsp; if (lieIndex >= offdiag) { diagIndex = lieIndex - offdiag; // 0, ... ,N-1 // std::cout << GridLogMessage << "diag type " << std::endl; @@ -78,6 +79,7 @@ template static void generatorDiagtype(int diagIndex, iSp2nMatrix &ta) { // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra + const int nsp=ncolour/2; ta = Zero(); RealD nrm = 1.0 / 2; @@ -91,6 +93,7 @@ static void generatorAtype(int aIndex, iSp2nMatrix &ta) { // with i &ta) { // with i static void generatorCtype(int cIndex, iSp2nMatrix &ta) { // ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2) + const int nsp=ncolour/2; int i1, i2; ta = Zero(); RealD nrm = 1 / (2 * std::sqrt(2)); @@ -144,6 +149,7 @@ template static void generatorDtype(int dIndex, iSp2nMatrix &ta) { // ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2) + const int nsp=ncolour/2; int i1, i2; ta = Zero(); cplx i(0.0, 1.0); @@ -162,6 +168,7 @@ template static void generatorWtype(int wIndex, iSp2nMatrix &ta) { // ta(i,i+N) = ta(i+N,i) = 1/2 + const int nsp=ncolour/2; ta = Zero(); RealD nrm = 1.0 / 2; // check @@ -175,6 +182,7 @@ template static void generatorZtype(int zIndex, iSp2nMatrix &ta) { // ta(i,i+N) = - ta(i+N,i) = i/2 + const int nsp=ncolour/2; ta = Zero(); RealD nrm = 1.0 / 2; // check cplx i(0.0, 1.0); @@ -189,6 +197,7 @@ static void generatorZtype(int zIndex, iSp2nMatrix &ta) { //////////////////////////////////////////////////////////////////////// template static void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::Sp) { + const int nsp=ncolour/2; assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2)); int spare = su2_index; @@ -239,6 +248,7 @@ template static void OmegaInvariance(ColourMatrix &in) { ColourMatrix Omega; Omega = Zero(); + const int nsp=ncolour/2; std::cout << GridLogMessage << "I am a ColourMatrix" << std::endl; @@ -289,6 +299,7 @@ static void OmegaInvariance(GaugeField &in) { template static void OmegaInvariance(LatticeColourMatrixD &in) { + const int nsp=ncolour/2; LatticeColourMatrixD OmegaLatt(in.Grid()); LatticeColourMatrixD identity(in.Grid()); RealD vol = in.Grid()->gSites(); diff --git a/tests/sp2n/Test_Sp_start.cc b/tests/sp2n/Test_Sp_start.cc index 16a66367..061acbcc 100644 --- a/tests/sp2n/Test_Sp_start.cc +++ b/tests/sp2n/Test_Sp_start.cc @@ -21,7 +21,7 @@ int main (int argc, char **argv) double vol = Umu.Grid()->gSites(); - const int nsp = Sp::nsp; + const int nsp = Nc/2; identity = 1.; Cidentity = 1.; diff --git a/tests/sp2n/Test_project_on_Sp.cc b/tests/sp2n/Test_project_on_Sp.cc index d1cf33bd..36b95cf4 100644 --- a/tests/sp2n/Test_project_on_Sp.cc +++ b/tests/sp2n/Test_project_on_Sp.cc @@ -20,8 +20,7 @@ int main (int argc, char **argv) LatticeColourMatrixD aux(&Grid); LatticeColourMatrixD identity(&Grid); - //const int nsp = Nc / 2; - const int nsp = Sp::nsp; + const int nsp = Nc / 2; identity = 1.0; RealD epsilon = 0.01; From a13820656a03cca25f74cc0300d0cc8c66cdf6a2 Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Wed, 30 Nov 2022 15:09:03 +0000 Subject: [PATCH 14/17] Removed iSUnMatrix, etc. --- Grid/qcd/action/gauge/GaugeImplTypes.h | 2 +- Grid/qcd/utils/SUn.h | 10 +--------- Grid/qcd/utils/SUn_impl.h | 12 ++++++------ Grid/qcd/utils/Sp2nTwoIndex.h | 8 ++++---- Grid/qcd/utils/Sp2n_impl.h | 18 +++++++++--------- 5 files changed, 21 insertions(+), 29 deletions(-) diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h index dc39c30b..5a8887f2 100644 --- a/Grid/qcd/action/gauge/GaugeImplTypes.h +++ b/Grid/qcd/action/gauge/GaugeImplTypes.h @@ -150,7 +150,7 @@ public: P = Ta(P); //const int nsp = Nc / 2; - Sp::iSp2nMatrix gen; + Sp::iGroupMatrix gen; auto Psum = P; diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index 2cedf67d..7a096df2 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -99,14 +99,6 @@ class GaugeGroup { template using iAlgebraVector = iScalar > >; static int su2subgroups(void) { return su2subgroups(group_name()); } - template - using iSUnMatrix = iGroupMatrix; - template - using iSUnAlgebraVector = iAlgebraVector; - template - using iSp2nMatrix = iGroupMatrix; - template - using iSp2nAlgebraVector = iAlgebraVector; ////////////////////////////////////////////////////////////////////////////////////////////////// // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix, @@ -156,7 +148,7 @@ class GaugeGroup { public: template - static void generator(int lieIndex, iSp2nMatrix &ta) { + static void generator(int lieIndex, iGroupMatrix &ta) { return generator(lieIndex, ta, group_name()); } diff --git a/Grid/qcd/utils/SUn_impl.h b/Grid/qcd/utils/SUn_impl.h index 6081f79c..61b19fbc 100644 --- a/Grid/qcd/utils/SUn_impl.h +++ b/Grid/qcd/utils/SUn_impl.h @@ -57,7 +57,7 @@ static int su2subgroups(GroupName::SU) { return (ncolour * (ncolour - 1)) / 2; } // //////////////////////////////////////////////////////////////////////// template -static void generator(int lieIndex, iSUnMatrix &ta, GroupName::SU) { +static void generator(int lieIndex, iGroupMatrix &ta, GroupName::SU) { // map lie index to which type of generator int diagIndex; int su2Index; @@ -77,7 +77,7 @@ static void generator(int lieIndex, iSUnMatrix &ta, GroupName::SU) { } template -static void generatorSigmaY(int su2Index, iSUnMatrix &ta) { +static void generatorSigmaY(int su2Index, iGroupMatrix &ta) { ta = Zero(); int i1, i2; su2SubGroupIndex(i1, i2, su2Index); @@ -87,7 +87,7 @@ static void generatorSigmaY(int su2Index, iSUnMatrix &ta) { } template -static void generatorSigmaX(int su2Index, iSUnMatrix &ta) { +static void generatorSigmaX(int su2Index, iGroupMatrix &ta) { ta = Zero(); cplx i(0.0, 1.0); int i1, i2; @@ -98,7 +98,7 @@ static void generatorSigmaX(int su2Index, iSUnMatrix &ta) { } template -static void generatorDiagonal(int diagIndex, iSUnMatrix &ta) { +static void generatorDiagonal(int diagIndex, iGroupMatrix &ta) { // diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...) ta = Zero(); int k = diagIndex + 1; // diagIndex starts from 0 @@ -130,7 +130,7 @@ public: template static void su2Extract(Lattice > &Determinant, Lattice > &subgroup, - const Lattice > &source, + const Lattice > &source, int su2_index) { GridBase *grid(source.Grid()); conformable(subgroup, source); @@ -164,7 +164,7 @@ static void su2Extract(Lattice > &Determinant, ////////////////////////////////////////////////////////////////////////////////////////// template static void su2Insert(const Lattice > &subgroup, - Lattice > &dest, int su2_index) { + Lattice > &dest, int su2_index) { GridBase *grid(dest.Grid()); conformable(subgroup, dest); int i0, i1; diff --git a/Grid/qcd/utils/Sp2nTwoIndex.h b/Grid/qcd/utils/Sp2nTwoIndex.h index 3b685802..98306970 100644 --- a/Grid/qcd/utils/Sp2nTwoIndex.h +++ b/Grid/qcd/utils/Sp2nTwoIndex.h @@ -175,10 +175,10 @@ public: template static void generator(int Index, iSp2nTwoIndexMatrix &i2indTa) { - Vector::template iSp2nMatrix > ta( + Vector > ta( NumGenerators); - Vector::template iSp2nMatrix > eij(Dimension); - typename Sp::template iSp2nMatrix tmp; + Vector > eij(Dimension); + iSp2nMatrix tmp; i2indTa = Zero(); for (int a = 0; a < NumGenerators; a++) @@ -189,7 +189,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 Sp::template iSp2nMatrix tmp1 = + iSp2nMatrix tmp1 = tmp * eij[b]; Complex iTr = TensorRemove(timesI(trace(tmp1))); i2indTa()()(a, b) = iTr; diff --git a/Grid/qcd/utils/Sp2n_impl.h b/Grid/qcd/utils/Sp2n_impl.h index 28dac5f8..6d559c0a 100644 --- a/Grid/qcd/utils/Sp2n_impl.h +++ b/Grid/qcd/utils/Sp2n_impl.h @@ -13,7 +13,7 @@ static int su2subgroups(GroupName::Sp) { return (ncolour/2 * (ncolour/2 - 1)) / // abcd are N(N-1)/2 each while wz are N each template -static void generator(int lieIndex, iSp2nMatrix &ta, GroupName::Sp) { +static void generator(int lieIndex, iGroupMatrix &ta, GroupName::Sp) { // map lie index into type of generators: diagonal, abcd type, wz type const int nsp = ncolour/2; @@ -76,7 +76,7 @@ static void generator(int lieIndex, iSp2nMatrix &ta, GroupName::Sp) { } // end of generator template -static void generatorDiagtype(int diagIndex, iSp2nMatrix &ta) { +static void generatorDiagtype(int diagIndex, iGroupMatrix &ta) { // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra const int nsp=ncolour/2; @@ -88,7 +88,7 @@ static void generatorDiagtype(int diagIndex, iSp2nMatrix &ta) { } template -static void generatorAtype(int aIndex, iSp2nMatrix &ta) { +static void generatorAtype(int aIndex, iGroupMatrix &ta) { // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2) // with i &ta) { } template -static void generatorBtype(int bIndex, iSp2nMatrix &ta) { +static void generatorBtype(int bIndex, iGroupMatrix &ta) { // ta(i,j) = -ta(j,i) = ta(i+N,j+N) = -ta(j+N,i+N) = i / 1/ 2 sqrt(2) // with i &ta) { } template -static void generatorCtype(int cIndex, iSp2nMatrix &ta) { +static void generatorCtype(int cIndex, iGroupMatrix &ta) { // ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2) const int nsp=ncolour/2; @@ -146,7 +146,7 @@ static void generatorCtype(int cIndex, iSp2nMatrix &ta) { } template -static void generatorDtype(int dIndex, iSp2nMatrix &ta) { +static void generatorDtype(int dIndex, iGroupMatrix &ta) { // ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2) const int nsp=ncolour/2; @@ -165,7 +165,7 @@ static void generatorDtype(int dIndex, iSp2nMatrix &ta) { } template -static void generatorWtype(int wIndex, iSp2nMatrix &ta) { +static void generatorWtype(int wIndex, iGroupMatrix &ta) { // ta(i,i+N) = ta(i+N,i) = 1/2 const int nsp=ncolour/2; @@ -179,7 +179,7 @@ static void generatorWtype(int wIndex, iSp2nMatrix &ta) { } template -static void generatorZtype(int zIndex, iSp2nMatrix &ta) { +static void generatorZtype(int zIndex, iGroupMatrix &ta) { // ta(i,i+N) = - ta(i+N,i) = i/2 const int nsp=ncolour/2; @@ -277,7 +277,7 @@ static void OmegaInvariance(ColourMatrix &in) { template static void OmegaInvariance(GaugeField &in) { typedef typename GaugeField::vector_type vector_type; - typedef iSp2nMatrix vMatrixType; + typedef iGroupMatrix vMatrixType; typedef Lattice LatticeMatrixType; LatticeMatrixType U(in.Grid()); From 7bcf33def98cced5bcfe5ee70527c1671617539a Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Wed, 30 Nov 2022 16:59:46 +0000 Subject: [PATCH 15/17] Removed Sp2n.h --- Grid/qcd/utils/SUn.h | 32 ++ Grid/qcd/utils/Sp2n.h | 597 -------------------------------- Grid/qcd/utils/Utils.h | 1 - tests/sp2n/Test_sp2n_lie_gen.cc | 4 - 4 files changed, 32 insertions(+), 602 deletions(-) delete mode 100644 Grid/qcd/utils/Sp2n.h diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index 7a096df2..88476f47 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -418,5 +418,37 @@ typedef SU<5> SU5; typedef SU FundamentalMatrices; +template +static void ProjectSp2n( + Lattice > > > &Umu) { + Umu = ProjectOnSpGroup(Umu); + auto det = Determinant(Umu); // ok ? + + 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 ProjectSp2n( + Lattice >, Nd> > &U) { + GridBase *grid = U.Grid(); + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnSpGroup(Umu); + ProjectSp2n(Umu); + PokeIndex(U, Umu, mu); + } +} + +typedef Sp<2> Sp2; +typedef Sp<4> Sp4; +typedef Sp<6> Sp6; +typedef Sp<8> Sp8; + + NAMESPACE_END(Grid); #endif diff --git a/Grid/qcd/utils/Sp2n.h b/Grid/qcd/utils/Sp2n.h deleted file mode 100644 index 201e1e7c..00000000 --- a/Grid/qcd/utils/Sp2n.h +++ /dev/null @@ -1,597 +0,0 @@ - -#ifndef QCD_UTIL_Sp2n_H -#define QCD_UTIL_Sp2n_H - -#include "Grid/qcd/utils/SUn.h" -NAMESPACE_BEGIN(Grid); -// -// // Sp(2N) -// // ncolour = 2N -// -// -// template -// class Sp { -// public: -// static const int nsp = ncolour/2; -// static const int Dimension = ncolour; -// static const int AlgebraDimension = nsp*(2*nsp +1); -// static int su2subgroups(void) { return (nsp * (nsp - 1)) / 2; } -// -// -// template -// using iSp2nMatrix = iScalar > >; -// template -// using iSU2Matrix = iScalar > >; -// template -// using iSp2nAlgebraVector = iScalar > >; -// -// typedef iSp2nMatrix Matrix; -// typedef iSp2nMatrix MatrixF; -// typedef iSp2nMatrix MatrixD; -// -// typedef iSp2nMatrix vMatrix; -// typedef iSp2nMatrix vMatrixF; -// typedef iSp2nMatrix vMatrixD; -// -// typedef iSp2nAlgebraVector AlgebraVector; -// typedef iSp2nAlgebraVector AlgebraVectorF; -// typedef iSp2nAlgebraVector AlgebraVectorD; -// -// typedef iSp2nAlgebraVector vAlgebraVector; -// typedef iSp2nAlgebraVector vAlgebraVectorF; -// typedef iSp2nAlgebraVector vAlgebraVectorD; -// -// typedef Lattice LatticeMatrix; -// typedef Lattice LatticeMatrixF; -// typedef Lattice LatticeMatrixD; -// -// typedef Lattice LatticeAlgebraVector; -// typedef Lattice LatticeAlgebraVectorF; -// typedef Lattice LatticeAlgebraVectorD; -// -// -// -// // Sp(2N) has N(2N+1) = 2N^2+N generators -// // -// // normalise the generators such that -// // Trace ( Ta Tb) = 1/2 delta_ab -// // -// // N generators in the cartan, 2N^2 off -// // off diagonal: -// // there are 6 types named a,b,c,d and w,z -// // abcd are N(N-1)/2 each while wz are N each -// -// template -// static void generator(int lieIndex, iSp2nMatrix &ta) { -// // map lie index into type of generators: diagonal, abcd type, wz -// type -// -// int diagIndex; -// int aIndex, bIndex, cIndex, dIndex; -// int wIndex, zIndex; // a,b,c,d are N(N-1)/2 and w,z are N -// int mod = nsp * (nsp-1) * 0.5; -// int offdiag = 2*nsp*nsp; // number of generators not in the cartan -// subalgebra int wmod = 4*mod; int zmod = wmod+nsp; if (lieIndex >= -// offdiag) { -// diagIndex = lieIndex - offdiag; // 0, ... ,N-1 -// //std::cout << GridLogMessage << "diag type " << std::endl; -// generatorDiagtype(diagIndex, ta); -// return; -// } -// if ( (lieIndex >= wmod) && (lieIndex < zmod) ) { -// //std::cout << GridLogMessage << "w type " << std::endl; -// wIndex = lieIndex- wmod; // 0, ... ,N-1 -// generatorWtype(wIndex,ta); -// return; -// } -// if ( (lieIndex >= zmod) && (lieIndex < offdiag) ) { -// //std::cout << GridLogMessage << "z type " << std::endl; -// //std::cout << GridLogMessage << "lie index " << lieIndex << -// std::endl; -// //std::cout << GridLogMessage << "z mod " << zmod << std::endl; -// zIndex = lieIndex - zmod; // 0, ... ,N-1 -// generatorZtype(zIndex,ta); -// return; -// } -// if (lieIndex < mod) { // atype 0, ... , N(N-1)/2=mod -// //std::cout << GridLogMessage << "a type " << std::endl; -// aIndex = lieIndex; -// //std::cout << GridLogMessage << "a indx " << aIndex << -// std::endl; generatorAtype(aIndex, ta); return; -// } -// if ( (lieIndex >= mod) && lieIndex < 2*mod) { // btype mod, ... , -// 2mod-1 -// //std::cout << GridLogMessage << "b type " << std::endl; -// bIndex = lieIndex - mod; -// generatorBtype(bIndex, ta); -// return; -// } -// if ( (lieIndex >= 2*mod) && lieIndex < 3*mod) { // ctype 2mod, ... , -// 3mod-1 -// //std::cout << GridLogMessage << "c type " << std::endl; -// cIndex = lieIndex - 2*mod; -// generatorCtype(cIndex, ta); -// return; -// } -// if ( (lieIndex >= 3*mod) && lieIndex < wmod) { // ctype 3mod, ... , -// 4mod-1 = wmod-1 -// //std::cout << GridLogMessage << "d type " << std::endl; -// dIndex = lieIndex - 3*mod; -// generatorDtype(dIndex, ta); -// return; -// } -// -// } //end of generator -// -// template -// static void generatorDiagtype(int diagIndex, iSp2nMatrix &ta) { -// -// // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan -// subalgebra -// -// ta = Zero(); -// RealD nrm = 1.0 / 2; -// -// ta()()(diagIndex,diagIndex) = nrm; -// ta()()(diagIndex+nsp,diagIndex+nsp) = -nrm; -// } -// -// template -// static void generatorAtype(int aIndex, iSp2nMatrix &ta) { -// -// // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2) -// // with i -// static void generatorBtype(int bIndex, iSp2nMatrix &ta) { -// -// // ta(i,j) = -ta(j,i) = ta(i+N,j+N) = -ta(j+N,i+N) = i / 1/ 2 sqrt(2) -// // with i -// static void generatorCtype(int cIndex, iSp2nMatrix &ta) { -// -// // ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2) -// -// -// int i1, i2; -// ta = Zero(); -// RealD nrm = 1 / (2 * std::sqrt(2) ); -// su2SubGroupIndex(i1, i2, cIndex); -// -// ta()()(i1,i2+nsp) = 1; -// ta()()(i2,i1+nsp) = 1; -// ta()()(i1+nsp,i2) = 1; -// ta()()(i2+nsp,i1) = 1; -// -// ta = ta * nrm; -// } -// -// template -// static void generatorDtype(int dIndex, iSp2nMatrix &ta) { -// -// // ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2) -// -// int i1, i2; -// ta = Zero(); -// cplx i(0.0, 1.0); -// RealD nrm = 1 / (2 * std::sqrt(2) ); -// su2SubGroupIndex(i1, i2, dIndex); -// -// ta()()(i1,i2+nsp) = i; -// ta()()(i2,i1+nsp) = i; -// ta()()(i1+nsp,i2) = -i; -// ta()()(i2+nsp,i1) = -i; -// -// ta = ta * nrm; -// } -// -// template -// static void generatorWtype(int wIndex, iSp2nMatrix &ta) { -// -// // ta(i,i+N) = ta(i+N,i) = 1/2 -// -// ta = Zero(); -// RealD nrm = 1.0 / 2; //check -// -// ta()()(wIndex,wIndex+nsp) = 1; -// ta()()(wIndex+nsp,wIndex) = 1; -// -// ta = ta * nrm; -// } -// -// template -// static void generatorZtype(int zIndex, iSp2nMatrix &ta) { -// -// // ta(i,i+N) = - ta(i+N,i) = i/2 -// -// ta = Zero(); -// RealD nrm = 1.0 / 2; //check -// cplx i(0.0, 1.0); -// ta()()(zIndex,zIndex+nsp) = i; -// ta()()(zIndex+nsp,zIndex) = -i; -// -// ta = ta * nrm; -// } -// -// -// //////////////////////////////////////////////////////////////////////// -// // Map a su2 subgroup number to the pair of rows that are non zero -// //////////////////////////////////////////////////////////////////////// -// static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { -// assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2)); -// -// int spare = su2_index; -// for (i1 = 0; spare >= (nsp - 1 - i1); i1++) { -// spare = spare - (nsp - 1 - i1); // remove the Nc-1-i1 terms -// } -// i2 = i1 + 1 + spare; -// } -// -// -// -// -// -// static void printGenerators(void) { -// for (int gen = 0; gen < AlgebraDimension; gen++) { -// Matrix ta; -// generator(gen, ta); -// std::cout << GridLogMessage << "Nc = " << ncolour << std::endl; -// std::cout << GridLogMessage << " t_" << gen << std::endl; -// std::cout << GridLogMessage << ta << std::endl; -// } -// } -// -// -// -// static void testGenerators(void) { -// Matrix ta; -// Matrix tb; -// std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is -// 0.5 delta_ab " << std::endl; for (int a = 0; a < AlgebraDimension; -// a++) { -// for (int b = 0; b < AlgebraDimension; b++) { -// generator(a,ta); -// generator(b,tb); -// Complex tr = TensorRemove(trace( ta * tb) ); -// std::cout << GridLogMessage << "(" << a << "," << b << ") = -// " << tr -// << std::endl; -// if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6); -// if (a != b) assert(abs(tr) < 1.0e-6); -// -// } -// } -// std::cout << GridLogMessage << std::endl; -// std::cout << GridLogMessage << "Fundamental - Checking if hermitian" -// << std::endl; for (int a = 0; a < AlgebraDimension; a++) { -// generator(a,ta); -// std::cout << GridLogMessage << a << std::endl; -// assert(norm2(ta - adj(ta)) < 1.0e-6); -// } -// std::cout << GridLogMessage << std::endl; -// std::cout << GridLogMessage << "Fundamental - Checking if traceless" -// << std::endl; for (int a = 0; a < AlgebraDimension; a++) { -// generator(a, ta); -// Complex tr = TensorRemove(trace(ta)); -// std::cout << GridLogMessage << a << std::endl; -// assert(abs(tr) < 1.0e-6); -// } -// -// } -// -// -// 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::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 < AlgebraDimension; 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, -// //same as sun -// 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 < AlgebraDimension; a++) { -// gaussian(pRNG, ca); -// generator(a, ta); -// la = toComplex(ca) * ta; -// out += la; -// } -// out *= ci; -// } -// -// 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 < AlgebraDimension; a++) { -// generator(a, ta); -// la = peekColour(h, a) * timesI(ta) * scale; -// out += la; -// } -// } -// -// -// template -// static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { // -// same as sun -// typedef typename LatticeMatrixType::scalar_type ComplexType; -// -// LatticeMatrixType xn(x.Grid()); -// RealD nfac = 1.0; -// -// xn = x; -// ex = xn + ComplexType(1.0); // 1+x -// -// // Do a 12th order exponentiation -// for (int i = 2; i <= 12; ++i) { -// nfac = nfac / RealD(i); // 1/2, 1/2.3 ... -// xn = xn * x; // x2, x3,x4.... -// ex = ex + xn * nfac; // x2/2!, x3/3!.... -// } -// } -// -// 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 < AlgebraDimension; a++) { -// generator(a, Ta); -// pokeColour(h_out, - 2.0 * (trace(timesI(Ta) * in)) * scale, a); -// } -// } -// -// -// template -// static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) { -// typedef typename GaugeField::vector_type vector_type; -// typedef iSp2nMatrix vMatrixType; -// typedef Lattice LatticeMatrixType; -// -// LatticeMatrixType Umu(out.Grid()); -// for (int mu = 0; mu < Nd; mu++) { -// LieRandomize(pRNG, Umu, 1.0); //def -// PokeIndex(out, Umu, mu); -// } -// } -// template -// static void TepidConfiguration(GridParallelRNG &pRNG,GaugeField &out){ -// typedef typename GaugeField::vector_type vector_type; -// typedef iSp2nMatrix vMatrixType; -// typedef Lattice LatticeMatrixType; -// -// LatticeMatrixType Umu(out.Grid()); -// for(int mu=0;mu(out,Umu,mu); -// } -// } -// -// template -// static void ColdConfiguration(GaugeField &out){ -// typedef typename GaugeField::vector_type vector_type; -// typedef iSp2nMatrix vMatrixType; -// typedef Lattice LatticeMatrixType; -// -// LatticeMatrixType Umu(out.Grid()); -// Umu=1.0; -// for(int mu=0;mu(out,Umu,mu); -// } -// } -// template -// static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){ -// ColdConfiguration(out); -// } -// -// static void OmegaInvariance(ColourMatrix &in) -// { -// -// ColourMatrix Omega; -// Omega = Zero(); -// -// 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; -// } -// -// auto diff = Omega - (in * Omega * transpose(in) ); -// auto sdiff = norm2(diff); -// if (norm2(sdiff) < 1e-8) -// { -// std::cout << GridLogMessage << "Symplectic condition satisfied: -// Omega invariant" << std::endl; -// } else { -// std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT -// left invariant by " << norm2(sdiff) << std::endl; -// } -// -// } -// -// template -// static void OmegaInvariance(GaugeField &in) -// { -// typedef typename GaugeField::vector_type vector_type; -// typedef iSp2nMatrix vMatrixType; -// typedef Lattice LatticeMatrixType; -// -// LatticeMatrixType U(in.Grid()); -// LatticeMatrixType OmegaLatt(in.Grid()); -// LatticeMatrixType identity(in.Grid()); -// RealD vol = in.Grid()->gSites(); -// ColourMatrix Omega; -// -// OmegaLatt = Zero(); -// Omega = Zero(); -// identity = 1.; -// -// std::cout << GridLogMessage << "I am a GaugeField " << std::endl; -// -// U = PeekIndex(in,1); -// -// OmegaInvariance(U); -// -// } -// -// static void OmegaInvariance(LatticeColourMatrixD &in) -// { -// -// LatticeColourMatrixD OmegaLatt(in.Grid()); -// LatticeColourMatrixD identity(in.Grid()); -// RealD vol = in.Grid()->gSites(); -// ColourMatrix Omega; -// -// OmegaLatt = Zero(); -// 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; -// } -// -// std::cout << GridLogMessage << "Omega = " << Omega()() << -// std::endl; OmegaLatt = OmegaLatt + (identity*Omega); -// -// auto diff = OmegaLatt - (in*OmegaLatt*transpose(in)); -// auto sdiff = sum(diff); -// //assert( norm2(sdiff) < 1e-8 ); -// if (norm2(sdiff) < 1e-8) -// { -// std::cout << GridLogMessage << "Symplectic condition -// satisfied: Omega invariant" << std::endl; -// } else { -// std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT -// left invariant by " << norm2(sdiff) << std::endl; -// } -// -// -// } -// -// -// }; // end of class Sp - -template -static void ProjectSp2n( - Lattice > > > &Umu) { - Umu = ProjectOnSpGroup(Umu); - auto det = Determinant(Umu); // ok ? - - 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 ProjectSp2n( - Lattice >, Nd> > &U) { - GridBase *grid = U.Grid(); - for (int mu = 0; mu < Nd; mu++) { - auto Umu = PeekIndex(U, mu); - Umu = ProjectOnSpGroup(Umu); - ProjectSp2n(Umu); - PokeIndex(U, Umu, mu); - } -} - -typedef Sp<2> Sp2; -typedef Sp<4> Sp4; -typedef Sp<6> Sp6; -typedef Sp<8> Sp8; - -NAMESPACE_END(Grid); -#endif diff --git a/Grid/qcd/utils/Utils.h b/Grid/qcd/utils/Utils.h index c654ce18..256fd18e 100644 --- a/Grid/qcd/utils/Utils.h +++ b/Grid/qcd/utils/Utils.h @@ -9,7 +9,6 @@ // Include representations #include -#include #include #include #include diff --git a/tests/sp2n/Test_sp2n_lie_gen.cc b/tests/sp2n/Test_sp2n_lie_gen.cc index 6c3766e6..5c6430a3 100644 --- a/tests/sp2n/Test_sp2n_lie_gen.cc +++ b/tests/sp2n/Test_sp2n_lie_gen.cc @@ -1,8 +1,4 @@ - - - #include -#include using namespace Grid; From 505fa49983c6745d1819a3f49da6ebd9e63e9511 Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Wed, 30 Nov 2022 17:09:48 +0000 Subject: [PATCH 16/17] Renamed SUn.h -> GaugeGroup.h --- Grid/qcd/utils/{SUn.h => GaugeGroup.h} | 0 Grid/qcd/utils/Utils.h | 2 +- tests/core/Test_lie_generators.cc | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) rename Grid/qcd/utils/{SUn.h => GaugeGroup.h} (100%) diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/GaugeGroup.h similarity index 100% rename from Grid/qcd/utils/SUn.h rename to Grid/qcd/utils/GaugeGroup.h diff --git a/Grid/qcd/utils/Utils.h b/Grid/qcd/utils/Utils.h index 256fd18e..87196cdc 100644 --- a/Grid/qcd/utils/Utils.h +++ b/Grid/qcd/utils/Utils.h @@ -8,7 +8,7 @@ #include // Include representations -#include +#include #include #include #include diff --git a/tests/core/Test_lie_generators.cc b/tests/core/Test_lie_generators.cc index e044378c..372e1408 100644 --- a/tests/core/Test_lie_generators.cc +++ b/tests/core/Test_lie_generators.cc @@ -32,7 +32,7 @@ directory #include -#include +#include #include #include From 0af7d5a793f7052a95e716ecbc5e8f125ac12101 Mon Sep 17 00:00:00 2001 From: Julian Lenz Date: Wed, 30 Nov 2022 17:12:00 +0000 Subject: [PATCH 17/17] Rename Grid/qcd/utils/_impl.h -> Grid/qcd/utils/.h --- Grid/qcd/utils/GaugeGroup.h | 4 ++-- Grid/qcd/utils/{SUn_impl.h => SUn.h} | 0 Grid/qcd/utils/{Sp2n_impl.h => Sp2n.h} | 0 3 files changed, 2 insertions(+), 2 deletions(-) rename Grid/qcd/utils/{SUn_impl.h => SUn.h} (100%) rename Grid/qcd/utils/{Sp2n_impl.h => Sp2n.h} (100%) diff --git a/Grid/qcd/utils/GaugeGroup.h b/Grid/qcd/utils/GaugeGroup.h index 88476f47..24b89945 100644 --- a/Grid/qcd/utils/GaugeGroup.h +++ b/Grid/qcd/utils/GaugeGroup.h @@ -143,8 +143,8 @@ class GaugeGroup { typedef Lattice LatticeSU2MatrixF; typedef Lattice LatticeSU2MatrixD; -#include "Grid/qcd/utils/SUn_impl.h" -#include "Grid/qcd/utils/Sp2n_impl.h" +#include "Grid/qcd/utils/SUn.h" +#include "Grid/qcd/utils/Sp2n.h" public: template diff --git a/Grid/qcd/utils/SUn_impl.h b/Grid/qcd/utils/SUn.h similarity index 100% rename from Grid/qcd/utils/SUn_impl.h rename to Grid/qcd/utils/SUn.h diff --git a/Grid/qcd/utils/Sp2n_impl.h b/Grid/qcd/utils/Sp2n.h similarity index 100% rename from Grid/qcd/utils/Sp2n_impl.h rename to Grid/qcd/utils/Sp2n.h