portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-04 19:25:56 +01:00
Code Releases Activity

Merge pull request #20 from LupoA/sp2n/unify_gaugeimpltypes

Browse Source 
Sp2n/unify gaugeimpltypes
This commit is contained in:
chillenzer 2023-05-03 15:17:10 +00:00 committed by GitHub
commit 3953312a93
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
19 changed files with 422 additions and 337 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Grid/lattice/Lattice_ET.h
View File

@ -345,6 +345,7 @@ GridUnopClass(UnaryNot, Not(a));
GridUnopClass(UnaryTrace, trace(a)); GridUnopClass(UnaryTrace, trace(a));
GridUnopClass(UnaryTranspose, transpose(a)); GridUnopClass(UnaryTranspose, transpose(a));
GridUnopClass(UnaryTa, Ta(a)); GridUnopClass(UnaryTa, Ta(a));
GridUnopClass(UnarySpTa, SpTa(a));
GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a)); GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a));
GridUnopClass(UnaryProjectOnSpGroup, ProjectOnSpGroup(a)); GridUnopClass(UnaryProjectOnSpGroup, ProjectOnSpGroup(a));
GridUnopClass(UnaryTimesI, timesI(a)); GridUnopClass(UnaryTimesI, timesI(a));
@ -457,6 +458,7 @@ GRID_DEF_UNOP(operator!, UnaryNot);
GRID_DEF_UNOP(trace, UnaryTrace); GRID_DEF_UNOP(trace, UnaryTrace);
GRID_DEF_UNOP(transpose, UnaryTranspose); GRID_DEF_UNOP(transpose, UnaryTranspose);
GRID_DEF_UNOP(Ta, UnaryTa); GRID_DEF_UNOP(Ta, UnaryTa);
GRID_DEF_UNOP(SpTa, UnarySpTa);
GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup); GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup);
GRID_DEF_UNOP(ProjectOnSpGroup, UnaryProjectOnSpGroup); GRID_DEF_UNOP(ProjectOnSpGroup, UnaryProjectOnSpGroup);
GRID_DEF_UNOP(timesI, UnaryTimesI); GRID_DEF_UNOP(timesI, UnaryTimesI);

6
Grid/qcd/action/gauge/Gauge.h
View File

@ -39,9 +39,9 @@ NAMESPACE_BEGIN(Grid);
typedef WilsonGaugeAction<PeriodicGimplR> WilsonGaugeActionR; typedef WilsonGaugeAction<PeriodicGimplR> WilsonGaugeActionR;
typedef WilsonGaugeAction<PeriodicGimplF> WilsonGaugeActionF; typedef WilsonGaugeAction<PeriodicGimplF> WilsonGaugeActionF;
typedef WilsonGaugeAction<PeriodicGimplD> WilsonGaugeActionD; typedef WilsonGaugeAction<PeriodicGimplD> WilsonGaugeActionD;
typedef WilsonGaugeAction<SymplPeriodicGimplR> SymplWilsonGaugeActionR; typedef WilsonGaugeAction<SpPeriodicGimplR> SpWilsonGaugeActionR;
typedef WilsonGaugeAction<SymplPeriodicGimplF> SymplWilsonGaugeActionF; typedef WilsonGaugeAction<SpPeriodicGimplF> SpWilsonGaugeActionF;
typedef WilsonGaugeAction<SymplPeriodicGimplD> SymplWilsonGaugeActionD; typedef WilsonGaugeAction<SpPeriodicGimplD> SpWilsonGaugeActionD;
typedef PlaqPlusRectangleAction<PeriodicGimplR> PlaqPlusRectangleActionR; typedef PlaqPlusRectangleAction<PeriodicGimplR> PlaqPlusRectangleActionR;
typedef PlaqPlusRectangleAction<PeriodicGimplF> PlaqPlusRectangleActionF; typedef PlaqPlusRectangleAction<PeriodicGimplF> PlaqPlusRectangleActionF;
typedef PlaqPlusRectangleAction<PeriodicGimplD> PlaqPlusRectangleActionD; typedef PlaqPlusRectangleAction<PeriodicGimplD> PlaqPlusRectangleActionD;

101
Grid/qcd/action/gauge/GaugeImplTypes.h
View File

@ -62,7 +62,7 @@ NAMESPACE_BEGIN(Grid);
// hardcodes the exponential approximation in the template // hardcodes the exponential approximation in the template
//template <class S, int Nrepresentation = Nc, int Nexp = 12 > class GaugeImplTypes { //template <class S, int Nrepresentation = Nc, int Nexp = 12 > class GaugeImplTypes {
template <class S, int Nrepresentation = Nc, int Nexp = 12, bool isSp2n = false > class GaugeImplTypes { template <class S, int Nrepresentation = Nc, int Nexp = 12, class Group = SU<Nc> > class GaugeImplTypes {
public: public:
typedef S Simd; typedef S Simd;
typedef typename Simd::scalar_type scalar_type; typedef typename Simd::scalar_type scalar_type;
@ -80,8 +80,6 @@ public:
typedef Lattice<SiteLink> LinkField; typedef Lattice<SiteLink> LinkField;
typedef Lattice<SiteField> Field; typedef Lattice<SiteField> Field;
typedef SU<Nrepresentation> Group;
// Guido: we can probably separate the types from the HMC functions // Guido: we can probably separate the types from the HMC functions
// this will create 2 kind of implementations // this will create 2 kind of implementations
// probably confusing the users // probably confusing the users
@ -124,53 +122,18 @@ public:
for (int mu = 0; mu < Nd; mu++) for (int mu = 0; mu < Nd; mu++)
{ {
if (isSp2n == true) Group::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
{
//const int nSp = Nrepresentation/2;
Sp<Nrepresentation>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
} else
{
Group::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
}
RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR) ; RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR) ;
Pmu = Pmu*scale; Pmu = Pmu*scale;
PokeIndex<LorentzIndex>(P, Pmu, mu); PokeIndex<LorentzIndex>(P, Pmu, mu);
} }
} }
static inline Field projectForce(Field &P) static inline Field projectForce(Field &P)
{ {
if (isSp2n == true) Field ret(P.Grid());
{ Group::taProj(P, ret);
P = Ta(P); return ret;
//const int nsp = Nc / 2;
Sp<Nc>::iGroupMatrix<Complex> gen;
auto Psum = P;
Psum = Zero();
for (int a = 0; a < Sp<Nrepresentation>::AlgebraDimension; a++)
{
Sp<Nrepresentation>::generator(a, gen);
auto coeff = 2. * trace(P * gen);
Psum += coeff * gen;
}
return Psum;
} else
{
return Ta(P);
}
} }
static inline void update_field(Field& P, Field& U, double ep){ static inline void update_field(Field& P, Field& U, double ep){
@ -181,14 +144,8 @@ public:
autoView(P_v,P,AcceleratorRead); autoView(P_v,P,AcceleratorRead);
accelerator_for(ss, P.Grid()->oSites(),1,{ accelerator_for(ss, P.Grid()->oSites(),1,{
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
if (isSp2n == true) U_v[ss](mu) = Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu);
{ Group::ProjectOnGaugeGroup(U_v[ss](mu));
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(); //auto end = std::chrono::high_resolution_clock::now();
@ -209,52 +166,22 @@ public:
static inline void Project(Field &U) static inline void Project(Field &U)
{ {
if (isSp2n == true) Group::ProjectGn(U);
{
ProjectSp2n(U);
} else
{
ProjectSUn(U);
}
} }
static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U)
{ {
if (isSp2n == true) Group::HotConfiguration(pRNG, U);
{
//const int nSp = Nrepresentation/2;
Sp<Nrepresentation>::HotConfiguration(pRNG, U);
} else
{
Group::HotConfiguration(pRNG, U);
}
} }
static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U)
{ {
if (isSp2n == true) Group::TepidConfiguration(pRNG, U);
{
//const int nSp = Nrepresentation/2;
Sp<Nrepresentation>::TepidConfiguration(pRNG, U);
} else
{
Group::TepidConfiguration(pRNG, U);
}
} }
static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U)
{ {
if (isSp2n == true) Group::ColdConfiguration(pRNG, U);
{
//const int nSp = Nrepresentation/2;
Sp<Nrepresentation>::ColdConfiguration(pRNG, U);
} else
{
Group::ColdConfiguration(pRNG, U);
}
} }
}; };
@ -264,9 +191,9 @@ typedef GaugeImplTypes<vComplex, Nc> GimplTypesR;
typedef GaugeImplTypes<vComplexF, Nc> GimplTypesF; typedef GaugeImplTypes<vComplexF, Nc> GimplTypesF;
typedef GaugeImplTypes<vComplexD, Nc> GimplTypesD; typedef GaugeImplTypes<vComplexD, Nc> GimplTypesD;
typedef GaugeImplTypes<vComplex, Nc, 12, true> SymplGimplTypesR; typedef GaugeImplTypes<vComplex, Nc, 12, Sp<Nc> > SpGimplTypesR;
typedef GaugeImplTypes<vComplexF, Nc, 12, true> SymplGimplTypesF; typedef GaugeImplTypes<vComplexF, Nc, 12, Sp<Nc> > SpGimplTypesF;
typedef GaugeImplTypes<vComplexD, Nc, 12, true> SymplGimplTypesD; typedef GaugeImplTypes<vComplexD, Nc, 12, Sp<Nc> > SpGimplTypesD;
typedef GaugeImplTypes<vComplex, SU<Nc>::AdjointDimension> GimplAdjointTypesR; typedef GaugeImplTypes<vComplex, SU<Nc>::AdjointDimension> GimplAdjointTypesR;
typedef GaugeImplTypes<vComplexF, SU<Nc>::AdjointDimension> GimplAdjointTypesF; typedef GaugeImplTypes<vComplexF, SU<Nc>::AdjointDimension> GimplAdjointTypesF;

6
Grid/qcd/action/gauge/GaugeImplementations.h
View File

@ -155,9 +155,9 @@ typedef ConjugateGaugeImpl<GimplTypesR> ConjugateGimplR; // Real.. whichever pre
typedef ConjugateGaugeImpl<GimplTypesF> ConjugateGimplF; // Float typedef ConjugateGaugeImpl<GimplTypesF> ConjugateGimplF; // Float
typedef ConjugateGaugeImpl<GimplTypesD> ConjugateGimplD; // Double typedef ConjugateGaugeImpl<GimplTypesD> ConjugateGimplD; // Double
typedef PeriodicGaugeImpl<SymplGimplTypesR> SymplPeriodicGimplR; // Real.. whichever prec typedef PeriodicGaugeImpl<SpGimplTypesR> SpPeriodicGimplR; // Real.. whichever prec
typedef PeriodicGaugeImpl<SymplGimplTypesF> SymplPeriodicGimplF; // Float typedef PeriodicGaugeImpl<SpGimplTypesF> SpPeriodicGimplF; // Float
typedef PeriodicGaugeImpl<SymplGimplTypesD> SymplPeriodicGimplD; // Double typedef PeriodicGaugeImpl<SpGimplTypesD> SpPeriodicGimplD; // Double
NAMESPACE_END(Grid); NAMESPACE_END(Grid);

4
Grid/qcd/hmc/GenericHMCrunner.h
View File

@ -207,12 +207,12 @@ using GenericHMCRunnerHirep =
// sp2n // sp2n
template <template <typename, typename, typename> class Integrator> template <template <typename, typename, typename> class Integrator>
using GenericSp2nHMCRunner = HMCWrapperTemplate<SymplPeriodicGimplR, Integrator>; using GenericSp2nHMCRunner = HMCWrapperTemplate<SpPeriodicGimplR, Integrator>;
template <class RepresentationsPolicy, template <class RepresentationsPolicy,
template <typename, typename, typename> class Integrator> template <typename, typename, typename> class Integrator>
using GenericSp2nHMCRunnerHirep = using GenericSp2nHMCRunnerHirep =
HMCWrapperTemplate<SymplPeriodicGimplR, Integrator, RepresentationsPolicy>; HMCWrapperTemplate<SpPeriodicGimplR, Integrator, RepresentationsPolicy>;

29
Grid/qcd/utils/GaugeGroup.h
View File

@ -170,7 +170,7 @@ class GaugeGroup {
} }
} }
template <typename LatticeMatrixType> template <typename LatticeMatrixType>
static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out,
double scale = 1.0) { double scale = 1.0) {
GridBase *grid = out.Grid(); GridBase *grid = out.Grid();
@ -300,9 +300,9 @@ class GaugeGroup {
ColdConfiguration(out); ColdConfiguration(out);
} }
template <typename LatticeMatrixType, ONLY_IF_SU> template <typename LatticeMatrixType>
static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) { static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) {
out = Ta(in); taProj(in, out, group_name());
} }
template <typename LatticeMatrixType> template <typename LatticeMatrixType>
@ -325,7 +325,10 @@ class GaugeGroup {
template <int N> // reunitarise, resimplectify... template <int N> // reunitarise, resimplectify...
static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) { static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
ProjectOnGaugeGroup(U, group_name()); for (int mu = 0; mu < Nd; mu++) {
auto Umu = PeekIndex<LorentzIndex>(U, mu);
ProjectOnGaugeGroup(Umu);
}
} }
template <int N> // reunitarise, resimplectify... template <int N> // reunitarise, resimplectify...
@ -333,6 +336,20 @@ class GaugeGroup {
ProjectOnGaugeGroup(Umu, group_name()); ProjectOnGaugeGroup(Umu, group_name());
} }
template <class vtype>
static void ProjectOnGaugeGroup(iScalar<vtype> &r) {
ProjectOnGaugeGroup(r, group_name());
}
template <class vtype, int N>
static void ProjectOnGaugeGroup(iVector<vtype,N> &r) {
ProjectOnGaugeGroup(r, group_name());
}
template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
static void ProjectOnGaugeGroup(iMatrix<vtype,N> &arg) {
ProjectOnGaugeGroup(arg, group_name());
}
template <int N> // reunitarise, resimplectify... previously ProjectSUn template <int N> // reunitarise, resimplectify... previously ProjectSUn
static void ProjectGn(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) { static void ProjectGn(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
@ -349,7 +366,7 @@ class GaugeGroup {
} }
template <int N> // reunitarise, resimplectify... previously ProjectSUn template <int N> // reunitarise, resimplectify... previously ProjectSUn
static void ProjectGn(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U, GroupName::SU) { static void ProjectGn(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
// Reunitarise // Reunitarise
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
auto Umu = PeekIndex<LorentzIndex>(U, mu); auto Umu = PeekIndex<LorentzIndex>(U, mu);
@ -386,7 +403,6 @@ LatticeComplexD Determinant(
return ret; return ret;
} }
template <int N> template <int N>
static void ProjectSUn( static void ProjectSUn(
Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) { Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
@ -415,7 +431,6 @@ static void ProjectSUn(
} }
} }
// Explicit specialisation for SU(3).
// Explicit specialisation for SU(3). // Explicit specialisation for SU(3).
static void ProjectSU3( static void ProjectSU3(
Lattice<iScalar<iScalar<iMatrix<vComplexD, 3> > > > &Umu) { Lattice<iScalar<iScalar<iMatrix<vComplexD, 3> > > > &Umu) {

27
Grid/qcd/utils/SUn.impl
View File

@ -511,15 +511,26 @@ static void testGenerators(GroupName::SU) {
template <int N> template <int N>
static void ProjectOnGaugeGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::SU) { static void ProjectOnGaugeGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::SU) {
Umu = ProjectOnGroup(Umu); Umu = ProjectOnGroup(Umu);
} }
template <int N> template <class vtype>
static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U, GroupName::SU) { static void ProjectOnGaugeGroup(iScalar<vtype> &r, GroupName::SU) {
// Reunitarise r = ProjectOnGroup(r);
for (int mu = 0; mu < Nd; mu++) { }
auto Umu = PeekIndex<LorentzIndex>(U, mu);
Umu = ProjectOnGroup(Umu); template <class vtype, int N>
} static void ProjectOnGaugeGroup(iVector<vtype,N> &r, GroupName::SU) {
r = ProjectOnGroup(r);
}
template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
static void ProjectOnGaugeGroup(iMatrix<vtype,N> &arg, GroupName::SU) {
arg = ProjectOnGroup(arg);
}
template <typename LatticeMatrixType>
static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out, GroupName::SU) {
out = Ta(in);
} }
/* /*

3
Grid/qcd/utils/SUnAdjoint.h
View File

@ -68,8 +68,7 @@ public:
for (int a = 0; a < Dimension; a++) { for (int a = 0; a < Dimension; a++) {
tmp = ta[a] * ta[Index] - ta[Index] * ta[a]; tmp = ta[a] * ta[Index] - ta[Index] * ta[a];
for (int b = 0; b < (ncolour * ncolour - 1); b++) { for (int b = 0; b < (ncolour * ncolour - 1); b++) {
iSUnMatrix<cplx> tmp1 = iSUnMatrix<cplx> tmp1 = 2.0 * tmp * ta[b]; // 2.0 from the normalization
2.0 * tmp * ta[b]; // 2.0 from the normalization
Complex iTr = TensorRemove(timesI(trace(tmp1))); Complex iTr = TensorRemove(timesI(trace(tmp1)));
//iAdjTa()()(b, a) = iTr; //iAdjTa()()(b, a) = iTr;
iAdjTa()()(a, b) = iTr; iAdjTa()()(a, b) = iTr;

3
Grid/qcd/utils/SUnTwoIndex.h
View File

@ -126,8 +126,7 @@ public:
template <class cplx> template <class cplx>
static void generator(int Index, iSUnTwoIndexMatrix<cplx> &i2indTa) { static void generator(int Index, iSUnTwoIndexMatrix<cplx> &i2indTa) {
Vector<iSUnMatrix<cplx> > ta( Vector<iSUnMatrix<cplx> > ta(ncolour * ncolour - 1);
ncolour * ncolour - 1);
Vector<iSUnMatrix<cplx> > eij(Dimension); Vector<iSUnMatrix<cplx> > eij(Dimension);
iSUnMatrix<cplx> tmp; iSUnMatrix<cplx> tmp;
i2indTa = Zero(); i2indTa = Zero();

56
Grid/qcd/utils/Sp2n.impl
View File

@ -256,15 +256,26 @@ static void testGenerators(GroupName::Sp) {
template <int N> template <int N>
static void ProjectOnGaugeGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::Sp) { static void ProjectOnGaugeGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::Sp) {
Umu = ProjectOnSpGroup(Umu); Umu = ProjectOnSpGroup(Umu);
} }
template <int N> template <class vtype>
static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U, GroupName::Sp) { static void ProjectOnGaugeGroup(iScalar<vtype> &r, GroupName::Sp) {
// Reunitarise r = ProjectOnSpGroup(r);
for (int mu = 0; mu < Nd; mu++) { }
auto Umu = PeekIndex<LorentzIndex>(U, mu);
Umu = ProjectOnSpGroup(Umu); template <class vtype, int N>
} static void ProjectOnGaugeGroup(iVector<vtype,N> &r, GroupName::Sp) {
r = ProjectOnSpGroup(r);
}
template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
static void ProjectOnGaugeGroup(iMatrix<vtype,N> &arg, GroupName::Sp) {
arg = ProjectOnSpGroup(arg);
}
template <typename LatticeMatrixType>
static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out, GroupName::Sp) {
out = SpTa(in);
} }
public: public:
@ -274,13 +285,6 @@ static void OmegaInvariance(ColourMatrix &in) {
Omega = Zero(); Omega = Zero();
const int nsp=ncolour/2; const int nsp=ncolour/2;
std::cout << GridLogMessage << "I am a ColourMatrix" << std::endl;
// for (int i = 0; i < ncolour; i++) wrong?!
//{
// Omega()()(i, 2*ncolour-1-i) = 1.;
// Omega()()(2*ncolour-1-i, i) = -1;
// }
for (int i = 0; i < nsp; i++) { for (int i = 0; i < nsp; i++) {
Omega()()(i, nsp + i) = 1.; Omega()()(i, nsp + i) = 1.;
Omega()()(nsp + i, i) = -1; Omega()()(nsp + i, i) = -1;
@ -314,8 +318,6 @@ static void OmegaInvariance(GaugeField &in) {
Omega = Zero(); Omega = Zero();
identity = 1.; identity = 1.;
std::cout << GridLogMessage << "I am a GaugeField " << std::endl;
U = PeekIndex<LorentzIndex>(in, 1); U = PeekIndex<LorentzIndex>(in, 1);
OmegaInvariance(U); OmegaInvariance(U);
@ -333,8 +335,6 @@ static void OmegaInvariance(LatticeColourMatrixD &in) {
Omega = Zero(); Omega = Zero();
identity = 1.; identity = 1.;
std::cout << GridLogMessage << "I am a LatticeColourMatrix " << std::endl;
for (int i = 0; i < nsp; i++) { for (int i = 0; i < nsp; i++) {
Omega()()(i, nsp + i) = 1.; Omega()()(i, nsp + i) = 1.;
Omega()()(nsp + i, i) = -1; Omega()()(nsp + i, i) = -1;
@ -356,3 +356,21 @@ static void OmegaInvariance(LatticeColourMatrixD &in) {
} }
} }
template <ONLY_IF_Sp>
static void Omega(LatticeColourMatrixD &in) {
const int nsp=ncolour/2;
LatticeColourMatrixD OmegaLatt(in.Grid());
LatticeColourMatrixD identity(in.Grid());
ColourMatrix Omega;
OmegaLatt = Zero();
Omega = Zero();
identity = 1.;
for (int i = 0; i < nsp; i++) {
Omega()()(i, nsp + i) = 1.;
Omega()()(nsp + i, i) = -1;
}
OmegaLatt = OmegaLatt + (identity * Omega);
in = OmegaLatt;
}

124
Grid/tensors/Tensor_Ta.h
View File

@ -66,13 +66,61 @@ template<class vtype,int N> accelerator_inline iMatrix<vtype,N> Ta(const iMatrix
return ret; return ret;
} }
template<class vtype> accelerator_inline iScalar<vtype> SpTa(const iScalar<vtype>&r)
{
iScalar<vtype> ret;
ret._internal = SpTa(r._internal);
return ret;
}
template<class vtype,int N> accelerator_inline iVector<vtype,N> SpTa(const iVector<vtype,N>&r)
{
iVector<vtype,N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = SpTa(r._internal[i]);
}
return ret;
}
template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
accelerator_inline iMatrix<vtype,N> SpTa(const iMatrix<vtype,N> &arg)
{
// Generalises Ta to Sp2n
// Applies the following projections
// P_{antihermitian} P_{antihermitian-Sp-algebra} P_{traceless}
// where the ordering matters
// P_{traceless} subtracts the trace
// P_{antihermitian-Sp-algebra} provides the block structure of the algebra based on U = exp(T) i.e. anti-hermitian generators
// P_{antihermitian} does in-adj(in) / 2
iMatrix<vtype,N> ret(arg);
double factor = (1.0/(double)N);
vtype nrm;
nrm = 0.5;
ret = arg - (trace(arg)*factor);
for(int c1=0;c1<N/2;c1++)
{
for(int c2=0;c2<N/2;c2++)
{
ret._internal[c1][c2] = nrm*(conjugate(ret._internal[c1+N/2][c2+N/2]) + ret._internal[c1][c2]); // new[up-left] = old[up-left]+old*[down-right]
ret._internal[c1][c2+N/2] = nrm*(ret._internal[c1][c2+N/2] - conjugate(ret._internal[c1+N/2][c2])); // new[up-right] = old[up-right]-old*[down-left]
}
for(int c2=N/2;c2<N;c2++)
{
ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]); // reconstructs lower blocks
ret._internal[c1+N/2][c2] = conjugate(ret._internal[c1][c2-N/2]); // from upper blocks
}
}
ret = (ret - adj(ret))*0.5;
return ret;
}
/////////////////////////////////////////////// ///////////////////////////////////////////////
// ProjectOnGroup function for scalar, vector, matrix // ProjectOnGroup function for scalar, vector, matrix
// Projects on orthogonal, unitary group // Projects on orthogonal, unitary group
/////////////////////////////////////////////// ///////////////////////////////////////////////
template<class vtype> accelerator_inline iScalar<vtype> ProjectOnGroup(const iScalar<vtype>&r) template<class vtype> accelerator_inline iScalar<vtype> ProjectOnGroup(const iScalar<vtype>&r)
{ {
iScalar<vtype> ret; iScalar<vtype> ret;
@ -164,58 +212,52 @@ accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg
iMatrix<vtype,N> ret(arg); iMatrix<vtype,N> ret(arg);
vtype nrm; vtype nrm;
vtype inner; vtype inner;
vtype tmp;
for(int c1=0;c1<N/2;c1++) for(int c1=0;c1<N/2;c1++)
{ {
for (int b=0; b<c1; b++) // remove the b-rows from U_c1 for (int b=0; b<c1; b++) // remove the b-rows from U_c1
{ {
decltype(ret._internal[b][b]*ret._internal[b][b]) pr; decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
decltype(ret._internal[b][b]*ret._internal[b][b]) prn; decltype(ret._internal[b][b]*ret._internal[b][b]) prn;
zeroit(pr); zeroit(pr);
zeroit(prn); zeroit(prn);
for(int c=0; c<N; c++) for(int c=0; c<N; c++)
{ {
pr += conjugate(ret._internal[c1][c])*ret._internal[b][c]; // <U_c1 | U_b > pr += conjugate(ret._internal[c1][c])*ret._internal[b][c]; // <U_c1 | U_b >
prn += conjugate(ret._internal[c1][c])*ret._internal[b+N/2][c]; // <U_c1 | U_{b+N} > prn += conjugate(ret._internal[c1][c])*ret._internal[b+N/2][c]; // <U_c1 | U_{b+N} >
} }
for(int c=0; c<N; c++) for(int c=0; c<N; c++)
{ {
ret._internal[c1][c] -= (conjugate(pr) * ret._internal[b][c] + conjugate(prn) * ret._internal[b+N/2][c] ); // U_c1 -= ( <U_c1 | U_b > U_b + <U_c1 | U_{b+N} > U_{b+N} ) ret._internal[c1][c] -= (conjugate(pr) * ret._internal[b][c] + conjugate(prn) * ret._internal[b+N/2][c] ); // U_c1 -= ( <U_c1 | U_b > U_b + <U_c1 | U_{b+N} > U_{b+N} )
}
} }
}
zeroit(inner); zeroit(inner);
for(int c2=0;c2<N;c2++) for(int c2=0;c2<N;c2++)
{ {
inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]); inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
} }
nrm = sqrt(inner); nrm = sqrt(inner);
nrm = 1.0/nrm; nrm = 1.0/nrm;
for(int c2=0;c2<N;c2++) for(int c2=0;c2<N;c2++)
{ {
ret._internal[c1][c2]*= nrm; ret._internal[c1][c2]*= nrm;
} }
for(int c2=0;c2<N/2;c2++)
for(int c2=0;c2<N/2;c2++) {
{ ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]); // down right in the new matrix = (up-left)* of the old matrix
tmp = conjugate(ret._internal[c1][c2]); // (up-left)* of the old matrix }
ret._internal[c1+N/2][c2+N/2] = tmp; // down right in the new matrix = (up-left)* of the old matrix
}
for(int c2=N/2;c2<N;c2++)
{
tmp = conjugate(ret._internal[c1][c2]); // (up-right)* of the old
ret._internal[c1+N/2][c2-N/2] = -tmp; // down left in the new matrix = -(up-right)* of the old
}
for(int c2=N/2;c2<N;c2++)
{
ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);; // down left in the new matrix = -(up-right)* of the old
}
} }
return ret; return ret;
} }

2
tests/core/Test_lie_generators.cc
View File

@ -47,7 +47,7 @@ using namespace Grid;
int main(int argc, char** argv) { int main(int argc, char** argv) {
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
std::vector<int> latt({4, 4, 4, 8}); std::vector<int> latt({4, 4, 4, 8});
GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid( GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid(
latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());

14
tests/core/Test_reunitarise.cc
View File

@ -36,7 +36,7 @@ using namespace Grid;
int main (int argc, char ** argv) int main (int argc, char ** argv)
{ {
Grid_init(&argc,&argv); Grid_init(&argc,&argv);
std::vector<int> latt({8,8,8,8}); std::vector<int> latt({8,8,8,8});
GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt, GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt,
GridDefaultSimd(Nd,vComplexD::Nsimd()), GridDefaultSimd(Nd,vComplexD::Nsimd()),
@ -122,15 +122,16 @@ int main (int argc, char ** argv)
detU=detU-1.0; detU=detU-1.0;
std::cout << "Determinant defect before projection " <<norm2(detU)<<std::endl; std::cout << "Determinant defect before projection " <<norm2(detU)<<std::endl;
tmp = U*adj(U) - ident; tmp = U*adj(U) - ident;
std::cout << "Unitarity check before projection " << norm2(tmp)<<std::endl; std::cout << "Unitarity check before projection " << norm2(tmp)<<std::endl;
#if Nc==3
#if Nc==3
ProjectSU3(U); ProjectSU3(U);
detU= Determinant(U) ; detU= Determinant(U) ;
detU= detU -1.0; detU= detU -1.0;
std::cout << "Determinant ProjectSU3 defect " <<norm2(detU)<<std::endl; std::cout << "Determinant ProjectSU3 defect " <<norm2(detU)<<std::endl;
tmp = U*adj(U) - ident; tmp = U*adj(U) - ident;
std::cout << "Unitarity check after projection " << norm2(tmp)<<std::endl; std::cout << "Unitarity check after projection " << norm2(tmp)<<std::endl;
#endif #endif
ProjectSUn(UU); ProjectSUn(UU);
detUU= Determinant(UU); detUU= Determinant(UU);
detUU= detUU -1.0; detUU= detUU -1.0;
@ -140,8 +141,3 @@ int main (int argc, char ** argv)
Grid_finalize(); Grid_finalize();
} }

2
tests/sp2n/Test_hmc_Sp_WF_2_Fund_3_2AS.cc
View File

@ -48,7 +48,7 @@ int main(int argc, char **argv) {
RealD beta = 6 ; RealD beta = 6 ;
SymplWilsonGaugeActionR Waction(beta); SpWilsonGaugeActionR Waction(beta);
auto GridPtr = TheHMC.Resources.GetCartesian(); auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); auto GridRBPtr = TheHMC.Resources.GetRBCartesian();

2
tests/sp2n/Test_hmc_Sp_Wilson2ASFermionGauge.cc
View File

@ -38,7 +38,7 @@ int main(int argc, char **argv) {
RealD beta = 6.7 ; RealD beta = 6.7 ;
SymplWilsonGaugeActionR Waction(beta); SpWilsonGaugeActionR Waction(beta);
auto GridPtr = TheHMC.Resources.GetCartesian(); auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); auto GridRBPtr = TheHMC.Resources.GetRBCartesian();

2
tests/sp2n/Test_hmc_Sp_WilsonFundFermionGauge.cc
View File

@ -38,7 +38,7 @@ int main(int argc, char **argv) {
RealD beta = 7.2 ; RealD beta = 7.2 ;
SymplWilsonGaugeActionR Waction(beta); SpWilsonGaugeActionR Waction(beta);
auto GridPtr = TheHMC.Resources.GetCartesian(); auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); auto GridRBPtr = TheHMC.Resources.GetRBCartesian();

2
tests/sp2n/Test_hmc_Sp_pureGaugeWilson.cc
View File

@ -80,7 +80,7 @@ int main(int argc, char **argv)
// that have a complex construction // that have a complex construction
// standard // standard
RealD beta = 8.0 ; RealD beta = 8.0 ;
SymplWilsonGaugeActionR Waction(beta); SpWilsonGaugeActionR Waction(beta);
ActionLevel<HMCWrapper::Field> Level1(1); ActionLevel<HMCWrapper::Field> Level1(1);
Level1.push_back(&Waction); Level1.push_back(&Waction);

363
tests/sp2n/Test_project_on_Sp.cc
View File

@ -2,155 +2,232 @@
using namespace Grid; using namespace Grid;
int main (int argc, char **argv) template <typename T>
{ bool has_correct_group_block_structure(const T& U) {
Grid_init(&argc,&argv); std::cout << GridLogMessage << "Checking the structure is " << std::endl;
std::cout << GridLogMessage << "U = ( W X ) " << std::endl;
std::cout << GridLogMessage << " ( -X^* W^* ) " << std::endl;
Coordinate latt_size = GridDefaultLatt(); std::cout << GridLogMessage << std::endl;
Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
Coordinate mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size,simd_layout,mpi_layout); const int nsp = Nc / 2;
GridRedBlackCartesian RBGrid(&Grid); Complex i(0., 1.);
for (int c1 = 0; c1 < nsp; c1++) // check on W
{
LatticeGaugeField Umu(&Grid); for (int c2 = 0; c2 < nsp; c2++) {
LatticeColourMatrixD U(&Grid); auto W = PeekIndex<ColourIndex>(U, c1, c2);
LatticeColourMatrixD aux(&Grid); auto Wstar = PeekIndex<ColourIndex>(U, c1 + nsp, c2 + nsp);
LatticeColourMatrixD identity(&Grid); auto Ww = conjugate(Wstar);
auto amizero = sum(W - Ww);
const int nsp = Nc / 2; auto amizeroo = TensorRemove(amizero);
assert(amizeroo.real() < 10e-6);
identity = 1.0; amizeroo *= i;
RealD epsilon = 0.01; assert(amizeroo.real() < 10e-6);
Complex i(0., 1.);
RealD u = 0.;
double vol = Umu.Grid()->gSites();
std::vector<int> pseeds({1,2,3,4,5});
std::vector<int> sseeds({6,7,8,9,10});
GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(pseeds);
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(sseeds);
SU<Nc>::HotConfiguration(pRNG,Umu);
U = PeekIndex<LorentzIndex>(Umu,2);
aux = U*adj(U) - identity;
std::cout <<GridLogMessage << std::endl;
std::cout << GridLogMessage << "Starting with random SUn matrix " << std::endl;
std::cout << GridLogMessage << "Unitary check " << std::endl;
std::cout <<GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
assert ( norm2(aux) < 1e-8 );
std::cout <<GridLogMessage << std::endl;
if (Nc != 2)
{
std::cout << GridLogMessage << "This matrix should not leave Omega invariant, expect a warning" << std::endl;
} }
Sp<Nc>::OmegaInvariance(U); }
std::cout <<GridLogMessage << std::endl;
U = U + epsilon*identity;
aux = U*adj(U) - identity;
std::cout << GridLogMessage << "Unitary matrix deformed " << std::endl;
std::cout << GridLogMessage << "now U adjU - 1 = " << norm2(aux) << std::endl;
std::cout <<GridLogMessage << std::endl;
std::cout <<GridLogMessage << std::endl; for (int c1 = 0; c1 < nsp; c1++) {
std::cout << GridLogMessage << "Projecting on Sp2n " << std::endl; for (int c2 = 0; c2 < nsp; c2++) {
auto X = PeekIndex<ColourIndex>(U, c1, c2 + nsp);
U = ProjectOnSpGroup(U); auto minusXstar = PeekIndex<ColourIndex>(U, c1 + nsp, c2);
//U = ProjectOnGroup(U); auto minusXx = conjugate(minusXstar);
aux = U*adj(U) - identity; auto amizero = sum(X + minusXx);
std::cout <<GridLogMessage << std::endl; auto amizeroo = TensorRemove(amizero);
std::cout <<GridLogMessage << std::endl; assert(amizeroo.real() < 10e-6);
std::cout << GridLogMessage << "Unitary check after Sp(2n) projection " << std::endl; amizeroo *= i;
std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl; assert(amizeroo.real() < 10e-6);
assert( norm2(aux) < 1e-8);
std::cout <<GridLogMessage << std::endl;
std::cout <<GridLogMessage << std::endl;
// checks on determinant
std::cout << GridLogMessage << "Det after Projection on Sp2n = " << norm2( Determinant(U) ) / vol << std::endl;
std::cout <<GridLogMessage << std::endl;
std::cout <<GridLogMessage << std::endl;
// actual sp2n check
std::cout << GridLogMessage << "Checking invariance after projection "<< std::endl;
Sp<Nc>::OmegaInvariance(U);
// checks on elements
std::cout <<GridLogMessage << std::endl;
std::cout <<GridLogMessage << std::endl;
std::cout << GridLogMessage << "Checking the structure is " << std::endl;
std::cout << GridLogMessage << "U = ( W X ) " << std::endl;
std::cout << GridLogMessage << " ( -X^* W^* ) " << std::endl;
std::cout <<GridLogMessage << std::endl;
for (int c1 = 0; c1 < nsp; c1++) //check on W
{
for (int c2 = 0; c2 < nsp; c2++)
{
auto W = PeekIndex<ColourIndex>(U,c1,c2);
auto Wstar = PeekIndex<ColourIndex>(U,c1+nsp,c2+nsp);
auto Ww = conjugate( Wstar );
auto amizero = sum(W - Ww);
auto amizeroo = TensorRemove(amizero);
assert( amizeroo.real() < 10e-6 );
amizeroo *= i;
assert( amizeroo.real() < 10e-6 );
}
} }
}
return true;
};
for (int c1 = 0; c1 < nsp ; c1++) template <typename T>
{ bool is_element_of_sp2n_group(const T& U) {
for (int c2 = 0; c2 < nsp; c2++) LatticeColourMatrixD aux(U.Grid());
{ LatticeColourMatrixD identity(U.Grid());
auto X = PeekIndex<ColourIndex>(U,c1,c2+nsp); identity = 1.0;
auto minusXstar = PeekIndex<ColourIndex>(U,c1+nsp,c2); LatticeColourMatrixD Omega(U.Grid());
auto minusXx = conjugate(minusXstar); Sp<Nc>::Omega(Omega);
auto amizero = sum (X + minusXx);
auto amizeroo = TensorRemove(amizero);
assert( amizeroo.real() < 10e-6 );
amizeroo *= i;
assert( amizeroo.real() < 10e-6 );
}
}
std::cout << GridLogMessage << "ok" << std::endl;
// an explicit check for sp2
/*
if (Nc == 2)
{
assert(Nc==2);
ColourMatrix A;
A = Zero();
Complex a(25041994., 12.); std::cout << GridLogMessage << "Check matrix is non-zero " << std::endl;
Complex b(39., 0.22); assert(norm2(U) > 1e-8);
Complex d(10000., -2222.3333);
A()()(0,0) = a;
A()()(0,1) = b;
A()()(1,0) = i;
A()()(1,1) = d;
std::cout <<GridLogMessage << std::endl;
std::cout <<GridLogMessage << std::endl;
std::cout << GridLogMessage << "An explicit check for Sp2" << std::endl;
std::cout <<GridLogMessage << std::endl;
std::cout << GridLogMessage << "Building a non unitary matrix by hand with funny entries " << std::endl;
std::cout << GridLogMessage << "A = " << A << std::endl;
std::cout << GridLogMessage << "Projecting on Sp2 " << std::endl;
A = ProjectOnSpGroup(A);
std::cout << GridLogMessage << "now A = " << A << std::endl;
std::cout << GridLogMessage << "A(0,0) - conjA(1,1) = " << A()()(0,0) - adj ( A()()(1,1) )<< std::endl;
std::cout << GridLogMessage << "A(0,1) + conjA(1,0) = " << A()()(0,1) + adj ( A()()(1,0) )<< std::endl;
}*/
Grid_finalize();
std::cout << GridLogMessage << "Unitary check" << std::endl;
aux = U * adj(U) - identity;
std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
assert(norm2(aux) < 1e-8);
aux = Omega - (U * Omega * transpose(U));
std::cout << GridLogMessage << "Omega - U Omega transpose(U) = " << norm2(aux)
<< std::endl;
assert(norm2(aux) < 1e-8);
std::cout << GridLogMessage
<< "|Det| = " << norm2(Determinant(U)) / U.Grid()->gSites()
<< std::endl;
assert(norm2(Determinant(U)) / U.Grid()->gSites() - 1 < 1e-8);
return has_correct_group_block_structure(U);
}
template <typename T>
void test_group_projections(T U) {
RealD Delta = 666.;
LatticeColourMatrixD identity(U.Grid());
identity = 1.0;
std::cout << GridLogMessage << "# # # #" << std::endl;
std::cout << GridLogMessage << "Group" << std::endl;
std::cout << GridLogMessage << "# # # #" << std::endl;
std::cout << GridLogMessage << std::endl;
std::string name = "ProjectOnSpGroup";
std::cout << GridLogMessage << "Testing " << name << std::endl;
std::cout << GridLogMessage << "Apply to deformed matrix" << std::endl;
U = U + Delta * identity;
U = ProjectOnSpGroup(U);
assert(is_element_of_sp2n_group(U));
name = "ProjectOnGaugeGroup";
std::cout << GridLogMessage << "Testing " << name << std::endl;
std::cout << GridLogMessage << "Apply to deformed matrix" << std::endl;
U = U + Delta * identity;
Sp<Nc>::ProjectOnGaugeGroup(U);
assert(is_element_of_sp2n_group(U));
name = "ProjectGn";
std::cout << GridLogMessage << "Testing " << name << std::endl;
std::cout << GridLogMessage << "Apply to deformed matrix" << std::endl;
U = U + Delta * identity;
Sp<Nc>::ProjectGn(U);
assert(is_element_of_sp2n_group(U));
}
template <typename T>
bool has_correct_algebra_block_structure(const T& U) {
// this only checks for the anti-hermitian part of the algebra
const int nsp = Nc / 2;
Complex i(0., 1.);
std::cout << GridLogMessage << "Checking the structure is " << std::endl;
std::cout << GridLogMessage << "U = ( W X ) " << std::endl;
std::cout << GridLogMessage << " ( -X^* W^* ) " << std::endl;
for (int c1 = 0; c1 < nsp; c1++) // check on W
{
for (int c2 = 0; c2 < nsp; c2++) {
auto W = PeekIndex<ColourIndex>(U, c1, c2);
auto Wstar = PeekIndex<ColourIndex>(U, c1 + nsp, c2 + nsp);
auto Ww = conjugate(Wstar);
auto amizero = sum(W - Ww);
auto amizeroo = TensorRemove(amizero);
assert(amizeroo.real() < 10e-6);
amizeroo *= i;
assert(amizeroo.real() < 10e-6);
}
}
for (int c1 = 0; c1 < nsp; c1++) {
for (int c2 = 0; c2 < nsp; c2++) {
auto X = PeekIndex<ColourIndex>(U, c1, c2 + nsp);
auto minusXstar = PeekIndex<ColourIndex>(U, c1 + nsp, c2);
auto minusXx = conjugate(minusXstar);
auto amizero = sum(X + minusXx);
auto amizeroo = TensorRemove(amizero);
assert(amizeroo.real() < 10e-6);
amizeroo *= i;
assert(amizeroo.real() < 10e-6);
}
}
return true;
}
template <typename T>
bool is_element_of_sp2n_algebra(const T& U) {
LatticeColourMatrixD aux(U.Grid());
LatticeColourMatrixD identity(U.Grid());
identity = 1.0;
LatticeColourMatrixD Omega(U.Grid());
Sp<Nc>::Omega(Omega);
std::cout << GridLogMessage << "Check matrix is non-zero " << std::endl;
assert(norm2(U) > 1e-8);
aux = U - adj(U);
std::cout << GridLogMessage << "T - Tda = " << norm2(aux)
<< " (not supposed to vanish)" << std::endl;
aux = U + adj(U);
std::cout << GridLogMessage << "T + Tda = " << norm2(aux)
<< " (supposed to vanish)" << std::endl;
assert(norm2(aux) - 1 < 1e-8);
std::cout << GridLogMessage << "Check that Omega T Omega + conj(T) = 0 "
<< std::endl;
aux = Omega * U * Omega + conjugate(U);
assert(norm2(aux) < 1e-8);
return has_correct_algebra_block_structure(U);
}
template <typename T>
void test_algebra_projections(T U) {
RealD Delta = 666.;
LatticeColourMatrixD tmp(U.Grid());
LatticeColourMatrixD identity(U.Grid());
identity = 1.0;
std::cout << GridLogMessage << "# # # #" << std::endl;
std::cout << GridLogMessage << "Algebra" << std::endl;
std::cout << GridLogMessage << "# # # #" << std::endl;
std::cout << GridLogMessage << std::endl;
std::string name = "SpTa";
std::cout << GridLogMessage << "Testing " << name << std::endl;
std::cout << GridLogMessage << "Apply to deformed matrix" << std::endl;
U = U + Delta * identity;
U = SpTa(U);
assert(is_element_of_sp2n_algebra(U));
name = "TaProj";
std::cout << GridLogMessage << "Testing " << name << std::endl;
std::cout << GridLogMessage << "Apply to deformed matrix" << std::endl;
U = U + Delta * identity;
Sp<Nc>::taProj(U, tmp);
U = tmp;
assert(is_element_of_sp2n_algebra(U));
}
int main(int argc, char** argv) {
Grid_init(&argc, &argv);
Coordinate latt_size = GridDefaultLatt();
Coordinate simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
Coordinate mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size, simd_layout, mpi_layout);
LatticeGaugeField Umu(&Grid);
LatticeColourMatrixD U(&Grid);
// Will test resimplectification-related functionalities (from
// ProjectOnGaugeGroup, ProjectOnSpGroup, ProjectGn) and projection on the
// algebra (from ProjectSp2nAlgebra) we work with matrices with positive
// determinant so detU = 1 even if in principle ProjectOnGaugeGroup and
// ProjectOnSpGroup allow for detU=-1 so the checks will be the same for the
// three functions NB only ProjectGn is the proper simplectification function
std::vector<int> pseeds({1, 2, 3, 4, 5});
GridParallelRNG pRNG(&Grid);
pRNG.SeedFixedIntegers(pseeds);
SU<Nc>::HotConfiguration(pRNG, Umu);
U = PeekIndex<LorentzIndex>(Umu, 0);
test_group_projections(U);
U = PeekIndex<LorentzIndex>(Umu, 1);
test_algebra_projections(U);
Grid_finalize();
} }

11
tests/sp2n/Test_sp2n_lie_gen.cc
View File

@ -10,14 +10,14 @@ int main(int argc, char** argv) {
//latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); //latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
//GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid); //GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid);
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
std::cout << GridLogMessage << "* Generators for Sp(2)" << std::endl; std::cout << GridLogMessage << "* Generators for Sp(2)" << std::endl;
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
Sp2::printGenerators(); Sp2::printGenerators();
Sp2::testGenerators(); Sp2::testGenerators();
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
@ -46,6 +46,5 @@ int main(int argc, char** argv) {
Sp8::printGenerators(); Sp8::printGenerators();
Sp8::testGenerators(); Sp8::testGenerators();
Grid_finalize(); Grid_finalize();
} }