From 046a23121e6c186d13b69ec52781aec0c48a1db0 Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Tue, 5 Oct 2021 15:51:22 +0100 Subject: [PATCH 001/129] sp2n generators --- Grid/qcd/utils/Sp2n.h | 396 ++++++++++++++++++++++++++++++++ Grid/qcd/utils/Utils.h | 1 + tests/core/test_sp2n_lie_gen.cc | 46 ++++ 3 files changed, 443 insertions(+) create mode 100644 Grid/qcd/utils/Sp2n.h create mode 100644 tests/core/test_sp2n_lie_gen.cc diff --git a/Grid/qcd/utils/Sp2n.h b/Grid/qcd/utils/Sp2n.h new file mode 100644 index 00000000..c8960383 --- /dev/null +++ b/Grid/qcd/utils/Sp2n.h @@ -0,0 +1,396 @@ + +#ifndef QCD_UTIL_Sp2n_H +#define QCD_UTIL_Sp2n_H + +NAMESPACE_BEGIN(Grid); + +// Sp(2N) +// ncolour = N +template +class Sp { +public: + static const int Dimension = ncolour*2; + static const int AlgebraDimension = ncolour*(2*ncolour +1); + static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; } + + + template + using iSp2nMatrix = iScalar > >; + template + using iSU2Matrix = iScalar > >; + template + using iSp2nAlgebraVector = iScalar > >; + + typedef iSp2nMatrix Matrix; + + typedef iSp2nMatrix vMatrix; + + + typedef Lattice LatticeMatrix; + // 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 = ncolour * (ncolour-1) * 0.5; + int offdiag = 2*ncolour*ncolour; // number of generators not in the cartan subalgebra + int wmod = 4*mod; + int zmod = wmod+ncolour; + 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+ncolour,diagIndex+ncolour) = -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+ncolour) = 1; + ta()()(i2,i1+ncolour) = 1; + ta()()(i1+ncolour,i2) = 1; + ta()()(i2+ncolour,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+ncolour) = i; + ta()()(i2,i1+ncolour) = i; + ta()()(i1+ncolour,i2) = -i; + ta()()(i2+ncolour,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+ncolour) = 1; + ta()()(wIndex+ncolour,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+ncolour) = i; + ta()()(zIndex+ncolour,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 < (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; + } + + + + + + static void printGenerators(void) { + for (int gen = 0; gen < AlgebraDimension; gen++) { + Matrix ta; + generator(gen, ta); + std::cout << GridLogMessage << "Nc (2n) = " << 2*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); + } + + } + + + 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; + } + + + + 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!.... + } + } + + + + + +}; // end of class Sp + + + +typedef Sp<2> Sp4; +typedef Sp<3> Sp6; +typedef Sp<4> Sp8; + +NAMESPACE_END(Grid); +#endif + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +/* } + sigxy = lieIndex & 0x1; // 1 if odd, 0 if even + su2Index = lieIndex >> 1 ; // where to put the sigma_x(y) + //for the even(odd) lieindex, sigmax(y) + if (sigxy) + //put sigmay at su2index + else + //put sigmax */ diff --git a/Grid/qcd/utils/Utils.h b/Grid/qcd/utils/Utils.h index 8ce3df37..35594656 100644 --- a/Grid/qcd/utils/Utils.h +++ b/Grid/qcd/utils/Utils.h @@ -9,6 +9,7 @@ // Include representations #include +#include #include #include diff --git a/tests/core/test_sp2n_lie_gen.cc b/tests/core/test_sp2n_lie_gen.cc new file mode 100644 index 00000000..5dbba9d3 --- /dev/null +++ b/tests/core/test_sp2n_lie_gen.cc @@ -0,0 +1,46 @@ + + + +#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(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 11fb943b1e756942e52cc24bda111df64b866a00 Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Mon, 11 Oct 2021 16:21:25 +0100 Subject: [PATCH 002/129] gauge and fermion implementation for sp2n --- Grid/qcd/action/fermion/Fermion.h | 5 + Grid/qcd/action/fermion/WilsonFermion.h | 2 + Grid/qcd/action/fermion/WilsonImpl.h | 10 + Grid/qcd/action/gauge/Gauge.h | 3 + Grid/qcd/action/gauge/GaugeImplTypes.h | 102 +++++++++-- Grid/qcd/action/gauge/GaugeImplementations.h | 5 + Grid/qcd/representations/Representations.h | 1 + Grid/qcd/representations/spfundamental.h | 42 +++++ Grid/qcd/utils/Sp2n.h | 181 +++++++++++++------ tests/core/test_sp2n_lie_gen.cc | 11 +- 10 files changed, 292 insertions(+), 70 deletions(-) create mode 100644 Grid/qcd/representations/spfundamental.h diff --git a/Grid/qcd/action/fermion/Fermion.h b/Grid/qcd/action/fermion/Fermion.h index 59fc17d5..b988b306 100644 --- a/Grid/qcd/action/fermion/Fermion.h +++ b/Grid/qcd/action/fermion/Fermion.h @@ -115,6 +115,11 @@ typedef WilsonFermion WilsonFermionR; typedef WilsonFermion WilsonFermionF; typedef WilsonFermion WilsonFermionD; + +typedef WilsonFermion SpWilsonFermionR; +typedef WilsonFermion SpWilsonFermionF; +typedef WilsonFermion SpWilsonFermionD; + //typedef WilsonFermion WilsonFermionRL; //typedef WilsonFermion WilsonFermionFH; //typedef WilsonFermion WilsonFermionDF; diff --git a/Grid/qcd/action/fermion/WilsonFermion.h b/Grid/qcd/action/fermion/WilsonFermion.h index bf8926d0..4a2d1ddd 100644 --- a/Grid/qcd/action/fermion/WilsonFermion.h +++ b/Grid/qcd/action/fermion/WilsonFermion.h @@ -208,5 +208,7 @@ public: typedef WilsonFermion WilsonFermionF; typedef WilsonFermion WilsonFermionD; +//typedef WilsonFermion SpWilsonFermionF; +//typedef WilsonFermion SpWilsonFermionD; NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/WilsonImpl.h b/Grid/qcd/action/fermion/WilsonImpl.h index 2685796d..01545d3b 100644 --- a/Grid/qcd/action/fermion/WilsonImpl.h +++ b/Grid/qcd/action/fermion/WilsonImpl.h @@ -267,6 +267,16 @@ typedef WilsonImpl W typedef WilsonImpl WilsonTwoIndexAntiSymmetricImplF; // Float typedef WilsonImpl WilsonTwoIndexAntiSymmetricImplD; // Double +//sp 2n + +typedef WilsonImpl SpWilsonImplR; // Real.. whichever prec +typedef WilsonImpl SpWilsonImplF; // Float +typedef WilsonImpl SpWilsonImplD; // Double + +//typedef WilsonImpl SpZWilsonImplR; // Real.. whichever prec +//typedef WilsonImpl SpZWilsonImplF; // Float +//typedef WilsonImpl SpZWilsonImplD; // Double + NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/gauge/Gauge.h b/Grid/qcd/action/gauge/Gauge.h index b3b34eba..535f0845 100644 --- a/Grid/qcd/action/gauge/Gauge.h +++ b/Grid/qcd/action/gauge/Gauge.h @@ -39,6 +39,9 @@ NAMESPACE_BEGIN(Grid); typedef WilsonGaugeAction WilsonGaugeActionR; typedef WilsonGaugeAction WilsonGaugeActionF; typedef WilsonGaugeAction WilsonGaugeActionD; +typedef WilsonGaugeAction SymplWilsonGaugeActionR; +typedef WilsonGaugeAction SymplWilsonGaugeActionF; +typedef WilsonGaugeAction SymplWilsonGaugeActionD; typedef PlaqPlusRectangleAction PlaqPlusRectangleActionR; typedef PlaqPlusRectangleAction PlaqPlusRectangleActionF; typedef PlaqPlusRectangleAction PlaqPlusRectangleActionD; diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h index 2499e0e9..2183d2e5 100644 --- a/Grid/qcd/action/gauge/GaugeImplTypes.h +++ b/Grid/qcd/action/gauge/GaugeImplTypes.h @@ -61,7 +61,8 @@ NAMESPACE_BEGIN(Grid); typedef typename Impl::Field Field; // hardcodes the exponential approximation in the template -template class GaugeImplTypes { +//template class GaugeImplTypes { +template class GaugeImplTypes { public: typedef S Simd; typedef typename Simd::scalar_type scalar_type; @@ -70,6 +71,7 @@ public: template using iImplGaugeLink = iScalar > >; template using iImplGaugeField = iVector >, Nd>; + typedef iImplScalar SiteComplex; typedef iImplGaugeLink SiteLink; typedef iImplGaugeField SiteField; @@ -117,8 +119,17 @@ public: // LinkField Pmu(P.Grid()); Pmu = Zero(); - for (int mu = 0; mu < Nd; mu++) { - SU::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); + for (int mu = 0; mu < Nd; mu++) + { + if (isSp2n == true) + { + const int nSp = Nrepresentation/2; + Sp::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); + } else + { + + SU::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); + } RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR) ; Pmu = Pmu*scale; PokeIndex(P, Pmu, mu); @@ -135,14 +146,21 @@ public: autoView(P_v,P,AcceleratorRead); accelerator_for(ss, P.Grid()->oSites(),1,{ for (int mu = 0; mu < Nd; mu++) { - U_v[ss](mu) = ProjectOnGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu)); + if (isSp2n == true) + { + U_v[ss](mu) = ProjectOnSpGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu)); + } else + { + U_v[ss](mu) = ProjectOnGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu)); + } + } }); //auto end = std::chrono::high_resolution_clock::now(); // diff += end - start; // std::cout << "Time to exponentiate matrix " << diff.count() << " s\n"; } - + static inline RealD FieldSquareNorm(Field& U){ LatticeComplex Hloc(U.Grid()); Hloc = Zero(); @@ -154,21 +172,77 @@ public: return Hsum.real(); } - static inline void Project(Field &U) { - ProjectSUn(U); + static inline void Project(Field &U) + { + if (isSp2n == true) + { + ProjectSp2n(U); + } else + { + ProjectSUn(U); + } } - static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) { - SU::HotConfiguration(pRNG, U); + static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) + { + SU::HotConfiguration(pRNG, U); } static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) { SU::TepidConfiguration(pRNG, U); } - static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) { - SU::ColdConfiguration(pRNG, U); + static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) + { + if (isSp2n == true) + { + const int nSp = Nrepresentation/2; + Sp::ColdConfiguration(pRNG, U); + } else + { + SU::ColdConfiguration(pRNG, U); + } } + + //sp2n... see sp2n.h +/* + static inline void generate_sp2n_momenta(Field &P, GridSerialRNG & sRNG, GridParallelRNG &pRNG) + { + + const int nSp = Nrepresentation/2; + LinkField Pmu(P.Grid()); + Pmu = Zero(); + for (int mu = 0; mu < Nd; mu++) { + + Sp::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); + RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR) ; + Pmu = Pmu*scale; + PokeIndex(P, Pmu, mu); + } + } + + static inline void update_sp2n_field(Field& P, Field& U, double ep){ + + autoView(U_v,U,AcceleratorWrite); + autoView(P_v,P,AcceleratorRead); + accelerator_for(ss, P.Grid()->oSites(),1,{ + for (int mu = 0; mu < Nd; mu++) { + U_v[ss](mu) = ProjectOnSpGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu)); + } + }); + } + + static inline void SpProject(Field &U) { + ProjectSp2n(U); + } + + + static inline void ColdSpConfiguration(GridParallelRNG &pRNG, Field &U) { + const int nSp = Nrepresentation/2; + Sp::ColdConfiguration(pRNG, U); + } +*/ + }; @@ -176,10 +250,16 @@ typedef GaugeImplTypes GimplTypesR; typedef GaugeImplTypes GimplTypesF; typedef GaugeImplTypes GimplTypesD; +typedef GaugeImplTypes SymplGimplTypesR; +typedef GaugeImplTypes SymplGimplTypesF; +typedef GaugeImplTypes SymplGimplTypesD; + typedef GaugeImplTypes::AdjointDimension> GimplAdjointTypesR; typedef GaugeImplTypes::AdjointDimension> GimplAdjointTypesF; typedef GaugeImplTypes::AdjointDimension> GimplAdjointTypesD; + + NAMESPACE_END(Grid); #endif // GRID_GAUGE_IMPL_TYPES_H diff --git a/Grid/qcd/action/gauge/GaugeImplementations.h b/Grid/qcd/action/gauge/GaugeImplementations.h index 16147c77..bc056263 100644 --- a/Grid/qcd/action/gauge/GaugeImplementations.h +++ b/Grid/qcd/action/gauge/GaugeImplementations.h @@ -155,6 +155,11 @@ typedef ConjugateGaugeImpl ConjugateGimplR; // Real.. whichever pre typedef ConjugateGaugeImpl ConjugateGimplF; // Float typedef ConjugateGaugeImpl ConjugateGimplD; // Double +typedef PeriodicGaugeImpl SymplPeriodicGimplR; // Real.. whichever prec +typedef PeriodicGaugeImpl SymplPeriodicGimplF; // Float +typedef PeriodicGaugeImpl SymplPeriodicGimplD; // Double + + NAMESPACE_END(Grid); #endif diff --git a/Grid/qcd/representations/Representations.h b/Grid/qcd/representations/Representations.h index 22311be0..d408b22a 100644 --- a/Grid/qcd/representations/Representations.h +++ b/Grid/qcd/representations/Representations.h @@ -4,6 +4,7 @@ #include #include #include +#include #include #endif diff --git a/Grid/qcd/representations/spfundamental.h b/Grid/qcd/representations/spfundamental.h new file mode 100644 index 00000000..bab9a005 --- /dev/null +++ b/Grid/qcd/representations/spfundamental.h @@ -0,0 +1,42 @@ + + +#ifndef SPFUNDAMENTAL_H +#define SPFUNDAMENTAL_H + +NAMESPACE_BEGIN(Grid); + +/* + * This is an helper class for the HMC + * Empty since HMC updates already the fundamental representation + */ + +template +class SpFundamentalRep { +public: + static const int Dimension = ncolour; + static const int nSp = ncolour/2; + static const bool isFundamental = true; + + // typdef to be used by the Representations class in HMC to get the + // types for the higher representation fields + typedef typename Sp::LatticeMatrix LatticeMatrix; + typedef LatticeGaugeField LatticeField; + + explicit SpFundamentalRep(GridBase* grid) {} //do nothing + void update_representation(const LatticeGaugeField& Uin) {} // do nothing + + LatticeField RtoFundamentalProject(const LatticeField& in, Real scale = 1.0) const{ + return (scale * in); + } + +}; + + + + + +typedef SpFundamentalRep SpFundamentalRepresentation; + +NAMESPACE_END(Grid); + +#endif diff --git a/Grid/qcd/utils/Sp2n.h b/Grid/qcd/utils/Sp2n.h index c8960383..51976146 100644 --- a/Grid/qcd/utils/Sp2n.h +++ b/Grid/qcd/utils/Sp2n.h @@ -6,6 +6,13 @@ NAMESPACE_BEGIN(Grid); // Sp(2N) // ncolour = N + +// need to be careful with n and 2n +// I am defining the Sp class for Sp(2n) to be such that the template variable ncolour +// is n inside 2n and the typedef at the end of the file should eliminate possible confusion. + +// the other routines, like projectOnSp2n, N will be the dimension of the actual number of colors for consistency with the sun routines + template class Sp { public: @@ -22,11 +29,31 @@ public: 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 @@ -322,75 +349,113 @@ public: } } + +/* template + static void HotSpConfiguration(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 TepidSpConfiguration(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); + } + }; // end of class Sp + +template //same of su + LatticeComplexD SpDeterminant(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); + Coordinate lcoor; + grid->LocalIndexToLocalCoor(site, lcoor); + iScalar > > Us; + peekLocalSite(Us, Umu_v, lcoor); + for(int i=0;i + static void ProjectSp2n(Lattice > > > &Umu) + { + Umu = ProjectOnSpGroup(Umu); + auto det = SpDeterminant(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<1> Sp2; typedef Sp<2> Sp4; typedef Sp<3> Sp6; typedef Sp<4> Sp8; NAMESPACE_END(Grid); #endif - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -/* } - sigxy = lieIndex & 0x1; // 1 if odd, 0 if even - su2Index = lieIndex >> 1 ; // where to put the sigma_x(y) - //for the even(odd) lieindex, sigmax(y) - if (sigxy) - //put sigmay at su2index - else - //put sigmax */ diff --git a/tests/core/test_sp2n_lie_gen.cc b/tests/core/test_sp2n_lie_gen.cc index 5dbba9d3..6c3766e6 100644 --- a/tests/core/test_sp2n_lie_gen.cc +++ b/tests/core/test_sp2n_lie_gen.cc @@ -2,7 +2,7 @@ #include -#include +#include using namespace Grid; @@ -13,6 +13,15 @@ int main(int argc, char** argv) { //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; From 4f5fe5792069c91c75ac68ff0569e5aa52787085 Mon Sep 17 00:00:00 2001 From: Alessandro Lupo Date: Mon, 11 Oct 2021 16:28:15 +0100 Subject: [PATCH 003/129] project on sp2n --- Grid/lattice/Lattice_ET.h | 2 + Grid/tensors/Tensor_Ta.h | 108 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 110 insertions(+) diff --git a/Grid/lattice/Lattice_ET.h b/Grid/lattice/Lattice_ET.h index 4a8a7423..9042eb8f 100644 --- a/Grid/lattice/Lattice_ET.h +++ b/Grid/lattice/Lattice_ET.h @@ -346,6 +346,7 @@ GridUnopClass(UnaryTrace, trace(a)); GridUnopClass(UnaryTranspose, transpose(a)); GridUnopClass(UnaryTa, Ta(a)); GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a)); +GridUnopClass(UnaryProjectOnSpGroup, ProjectOnSpGroup(a)); GridUnopClass(UnaryTimesI, timesI(a)); GridUnopClass(UnaryTimesMinusI, timesMinusI(a)); GridUnopClass(UnaryAbs, abs(a)); @@ -457,6 +458,7 @@ GRID_DEF_UNOP(trace, UnaryTrace); GRID_DEF_UNOP(transpose, UnaryTranspose); GRID_DEF_UNOP(Ta, UnaryTa); GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup); +GRID_DEF_UNOP(ProjectOnSpGroup, UnaryProjectOnSpGroup); GRID_DEF_UNOP(timesI, UnaryTimesI); GRID_DEF_UNOP(timesMinusI, UnaryTimesMinusI); GRID_DEF_UNOP(abs, UnaryAbs); // abs overloaded in cmath C++98; DON'T do the diff --git a/Grid/tensors/Tensor_Ta.h b/Grid/tensors/Tensor_Ta.h index 90e57b2b..f53296bc 100644 --- a/Grid/tensors/Tensor_Ta.h +++ b/Grid/tensors/Tensor_Ta.h @@ -135,6 +135,114 @@ accelerator_inline iMatrix ProjectOnGroup(const iMatrix &arg) return ret; } +// re-do for sp2n + +template accelerator_inline iScalar ProjectOnSpGroup(const iScalar&r) +{ + iScalar ret; + ret._internal = ProjectOnSpGroup(r._internal); + return ret; +} +template accelerator_inline iVector ProjectOnSpGroup(const iVector&r) +{ + iVector ret; + for(int i=0;i::TensorLevel == 0 >::type * =nullptr> +accelerator_inline iMatrix ProjectOnSpGroup(const iMatrix &arg) +{ + // need a check for the group type? + iMatrix ret(arg); + vtype nrm; + vtype inner; + for(int c1=0;c1 Date: Mon, 11 Oct 2021 16:32:10 +0100 Subject: [PATCH 004/129] sp fermion instantiation --- Grid/Makefile.am | 1 + .../WilsonCloverFermionInstantiationSpWilsonImplD.cc | 1 + .../WilsonFermionInstantiationSpWilsonImplD.cc | 1 + .../WilsonKernelsInstantiationSpWilsonImplD.cc | 1 + .../WilsonTMFermionInstantiationSpWilsonImplD.cc | 1 + .../fermion/instantiation/SpWilsonImplD/impl.h | 1 + .../WilsonCloverFermionInstantiationSpWilsonImplF.cc | 1 + .../WilsonFermionInstantiationSpWilsonImplF.cc | 1 + .../WilsonKernelsInstantiationSpWilsonImplF.cc | 1 + .../WilsonTMFermionInstantiationSpWilsonImplF.cc | 1 + .../fermion/instantiation/SpWilsonImplF/impl.h | 1 + .../fermion/instantiation/generate_instantiations.sh | 2 ++ Grid/qcd/hmc/GenericHMCrunner.h | 12 ++++++++++++ scripts/filelist | 2 ++ 14 files changed, 27 insertions(+) create mode 120000 Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonCloverFermionInstantiationSpWilsonImplD.cc create mode 120000 Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonFermionInstantiationSpWilsonImplD.cc create mode 120000 Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonKernelsInstantiationSpWilsonImplD.cc create mode 120000 Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonTMFermionInstantiationSpWilsonImplD.cc create mode 100644 Grid/qcd/action/fermion/instantiation/SpWilsonImplD/impl.h create mode 120000 Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonCloverFermionInstantiationSpWilsonImplF.cc create mode 120000 Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonFermionInstantiationSpWilsonImplF.cc create mode 120000 Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonKernelsInstantiationSpWilsonImplF.cc create mode 120000 Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonTMFermionInstantiationSpWilsonImplF.cc create mode 100644 Grid/qcd/action/fermion/instantiation/SpWilsonImplF/impl.h diff --git a/Grid/Makefile.am b/Grid/Makefile.am index 7c3c151b..fda92792 100644 --- a/Grid/Makefile.am +++ b/Grid/Makefile.am @@ -63,6 +63,7 @@ if BUILD_GPARITY extra_sources+=$(GP_FERMION_FILES) endif if BUILD_FERMION_REPS + extra_sources+=$(SP_FERMION_FILES) extra_sources+=$(ADJ_FERMION_FILES) extra_sources+=$(TWOIND_FERMION_FILES) endif diff --git a/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonCloverFermionInstantiationSpWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonCloverFermionInstantiationSpWilsonImplD.cc new file mode 120000 index 00000000..9cc05107 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonCloverFermionInstantiationSpWilsonImplD.cc @@ -0,0 +1 @@ +../WilsonCloverFermionInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonFermionInstantiationSpWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonFermionInstantiationSpWilsonImplD.cc new file mode 120000 index 00000000..5f6ab65e --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonFermionInstantiationSpWilsonImplD.cc @@ -0,0 +1 @@ +../WilsonFermionInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonKernelsInstantiationSpWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonKernelsInstantiationSpWilsonImplD.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonKernelsInstantiationSpWilsonImplD.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonTMFermionInstantiationSpWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonTMFermionInstantiationSpWilsonImplD.cc new file mode 120000 index 00000000..d5789bcf --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonTMFermionInstantiationSpWilsonImplD.cc @@ -0,0 +1 @@ +../WilsonTMFermionInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/impl.h b/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/impl.h new file mode 100644 index 00000000..5aa843dd --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/SpWilsonImplD/impl.h @@ -0,0 +1 @@ +#define IMPLEMENTATION SpWilsonImplD diff --git a/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonCloverFermionInstantiationSpWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonCloverFermionInstantiationSpWilsonImplF.cc new file mode 120000 index 00000000..9cc05107 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonCloverFermionInstantiationSpWilsonImplF.cc @@ -0,0 +1 @@ +../WilsonCloverFermionInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonFermionInstantiationSpWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonFermionInstantiationSpWilsonImplF.cc new file mode 120000 index 00000000..5f6ab65e --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonFermionInstantiationSpWilsonImplF.cc @@ -0,0 +1 @@ +../WilsonFermionInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonKernelsInstantiationSpWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonKernelsInstantiationSpWilsonImplF.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonKernelsInstantiationSpWilsonImplF.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonTMFermionInstantiationSpWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonTMFermionInstantiationSpWilsonImplF.cc new file mode 120000 index 00000000..d5789bcf --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonTMFermionInstantiationSpWilsonImplF.cc @@ -0,0 +1 @@ +../WilsonTMFermionInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/impl.h b/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/impl.h new file mode 100644 index 00000000..761b0c82 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/SpWilsonImplF/impl.h @@ -0,0 +1 @@ +#define IMPLEMENTATION SpWilsonImplF diff --git a/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh b/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh index d7553cdb..9c834f6e 100755 --- a/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh +++ b/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh @@ -9,6 +9,8 @@ STAG5_IMPL_LIST="" WILSON_IMPL_LIST=" \ WilsonImplF \ WilsonImplD \ + SpWilsonImplF \ + SpWilsonImplD \ WilsonAdjImplF \ WilsonAdjImplD \ WilsonTwoIndexSymmetricImplF \ diff --git a/Grid/qcd/hmc/GenericHMCrunner.h b/Grid/qcd/hmc/GenericHMCrunner.h index 98e8175a..16143455 100644 --- a/Grid/qcd/hmc/GenericHMCrunner.h +++ b/Grid/qcd/hmc/GenericHMCrunner.h @@ -204,6 +204,18 @@ template ; +// sp2n + +template