portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-09-20 09:15:38 +01:00
Code Releases Activity

Merge pull request #19 from LupoA/refactoring_sp2n

Browse Source 
refactoring sp2n
This commit is contained in:
chillenzer 2023-04-05 10:50:58 +00:00 committed by GitHub
commit 6a0eb466ee
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 143 additions and 68 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -170,7 +170,6 @@ class GaugeGroup {
} }
} }
// reunitarise??
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) {
@ -282,6 +281,7 @@ class GaugeGroup {
PokeIndex<LorentzIndex>(out, Umu, mu); PokeIndex<LorentzIndex>(out, Umu, mu);
} }
} }
template <typename GaugeField> template <typename GaugeField>
static void ColdConfiguration(GaugeField &out) { static void ColdConfiguration(GaugeField &out) {
typedef typename GaugeField::vector_type vector_type; typedef typename GaugeField::vector_type vector_type;
@ -294,6 +294,7 @@ class GaugeGroup {
PokeIndex<LorentzIndex>(out, Umu, mu); PokeIndex<LorentzIndex>(out, Umu, mu);
} }
} }
template <typename GaugeField> template <typename GaugeField>
static void ColdConfiguration(GridParallelRNG &pRNG, GaugeField &out) { static void ColdConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
ColdConfiguration(out); ColdConfiguration(out);
@ -303,6 +304,7 @@ class GaugeGroup {
static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) { static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) {
out = Ta(in); out = Ta(in);
} }
template <typename LatticeMatrixType> template <typename LatticeMatrixType>
static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) {
typedef typename LatticeMatrixType::scalar_type ComplexType; typedef typename LatticeMatrixType::scalar_type ComplexType;
@ -320,6 +322,42 @@ class GaugeGroup {
ex = ex + xn * nfac; // x2/2!, x3/3!.... ex = ex + xn * nfac; // x2/2!, x3/3!....
} }
} }
template <int N> // reunitarise, resimplectify...
static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
ProjectOnGaugeGroup(U, group_name());
}
template <int N> // reunitarise, resimplectify...
static void ProjectOnGaugeGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
ProjectOnGaugeGroup(Umu, group_name());
}
template <int N> // reunitarise, resimplectify... previously ProjectSUn
static void ProjectGn(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
ProjectOnGaugeGroup(Umu);
auto det = Determinant(Umu);
det = conjugate(det);
for (int i = 0; i < N; i++) {
auto element = PeekIndex<ColourIndex>(Umu, N - 1, i);
element = element * det;
PokeIndex<ColourIndex>(Umu, element, Nc - 1, i);
}
}
template <int N> // reunitarise, resimplectify... previously ProjectSUn
static void ProjectGn(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U, GroupName::SU) {
// Reunitarise
for (int mu = 0; mu < Nd; mu++) {
auto Umu = PeekIndex<LorentzIndex>(U, mu);
ProjectGn(Umu);
PokeIndex<LorentzIndex>(U, Umu, mu);
}
}
}; };
template <int N> template <int N>
@ -347,6 +385,8 @@ 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) {
@ -361,6 +401,7 @@ static void ProjectSUn(
PokeIndex<ColourIndex>(Umu, element, Nc - 1, i); PokeIndex<ColourIndex>(Umu, element, Nc - 1, i);
} }
} }
template <int N> template <int N>
static void ProjectSUn( static void ProjectSUn(
Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) { Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
@ -373,6 +414,7 @@ static void ProjectSUn(
PokeIndex<LorentzIndex>(U, Umu, mu); PokeIndex<LorentzIndex>(U, Umu, mu);
} }
} }
// Explicit specialisation for SU(3). // Explicit specialisation for SU(3).
// Explicit specialisation for SU(3). // Explicit specialisation for SU(3).
static void ProjectSU3( static void ProjectSU3(

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

@ -508,6 +508,20 @@ static void testGenerators(GroupName::SU) {
std::cout << GridLogMessage << std::endl; std::cout << GridLogMessage << std::endl;
} }
template <int N>
static void ProjectOnGaugeGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::SU) {
Umu = ProjectOnGroup(Umu);
}
template <int N>
static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U, GroupName::SU) {
// Reunitarise
for (int mu = 0; mu < Nd; mu++) {
auto Umu = PeekIndex<LorentzIndex>(U, mu);
Umu = ProjectOnGroup(Umu);
}
}
/* /*
* Fundamental rep gauge xform * Fundamental rep gauge xform
*/ */

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

@ -51,14 +51,16 @@ public:
typedef Lattice<iVector<iScalar<iMatrix<vComplexF, Dimension> >, Nd> > LatticeAdjFieldF; typedef Lattice<iVector<iScalar<iMatrix<vComplexF, Dimension> >, Nd> > LatticeAdjFieldF;
typedef Lattice<iVector<iScalar<iMatrix<vComplexD, Dimension> >, Nd> > LatticeAdjFieldD; typedef Lattice<iVector<iScalar<iMatrix<vComplexD, Dimension> >, Nd> > LatticeAdjFieldD;
template <typename vtype>
using iSUnMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >;
template <class cplx> template <class cplx>
static void generator(int Index, iSUnAdjointMatrix<cplx> &iAdjTa) { static void generator(int Index, iSUnAdjointMatrix<cplx> &iAdjTa) {
// returns i(T_Adj)^index necessary for the projectors // returns i(T_Adj)^index necessary for the projectors
// see definitions above // see definitions above
iAdjTa = Zero(); iAdjTa = Zero();
Vector<typename SU<ncolour>::template iSUnMatrix<cplx> > ta(ncolour * ncolour - 1); Vector<iSUnMatrix<cplx> > ta(ncolour * ncolour - 1);
typename SU<ncolour>::template iSUnMatrix<cplx> tmp; iSUnMatrix<cplx> tmp;
// FIXME not very efficient to get all the generators everytime // FIXME not very efficient to get all the generators everytime
for (int a = 0; a < Dimension; a++) SU<ncolour>::generator(a, ta[a]); for (int a = 0; a < Dimension; a++) SU<ncolour>::generator(a, ta[a]);
@ -66,7 +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++) {
typename SU<ncolour>::template 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;

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

@ -126,10 +126,10 @@ public:
template <class cplx> template <class cplx>
static void generator(int Index, iSUnTwoIndexMatrix<cplx> &i2indTa) { static void generator(int Index, iSUnTwoIndexMatrix<cplx> &i2indTa) {
Vector<typename SU<ncolour>::template iSUnMatrix<cplx> > ta( Vector<iSUnMatrix<cplx> > ta(
ncolour * ncolour - 1); ncolour * ncolour - 1);
Vector<typename SU<ncolour>::template iSUnMatrix<cplx> > eij(Dimension); Vector<iSUnMatrix<cplx> > eij(Dimension);
typename SU<ncolour>::template iSUnMatrix<cplx> tmp; iSUnMatrix<cplx> tmp;
i2indTa = Zero(); i2indTa = Zero();
for (int a = 0; a < ncolour * ncolour - 1; a++) for (int a = 0; a < ncolour * ncolour - 1; a++)
@ -140,7 +140,7 @@ public:
for (int a = 0; a < Dimension; a++) { for (int a = 0; a < Dimension; a++) {
tmp = transpose(ta[Index]) * adj(eij[a]) + adj(eij[a]) * ta[Index]; tmp = transpose(ta[Index]) * adj(eij[a]) + adj(eij[a]) * ta[Index];
for (int b = 0; b < Dimension; b++) { for (int b = 0; b < Dimension; b++) {
typename SU<ncolour>::template iSUnMatrix<cplx> tmp1 = iSUnMatrix<cplx> tmp1 =
tmp * eij[b]; tmp * eij[b];
Complex iTr = TensorRemove(timesI(trace(tmp1))); Complex iTr = TensorRemove(timesI(trace(tmp1)));
i2indTa()()(a, b) = iTr; i2indTa()()(a, b) = iTr;

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

@ -253,6 +253,20 @@ static void testGenerators(GroupName::Sp) {
} }
} }
template <int N>
static void ProjectOnGaugeGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::Sp) {
Umu = ProjectOnSpGroup(Umu);
}
template <int N>
static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U, GroupName::Sp) {
// Reunitarise
for (int mu = 0; mu < Nd; mu++) {
auto Umu = PeekIndex<LorentzIndex>(U, mu);
Umu = ProjectOnSpGroup(Umu);
}
}
public: public:
template <ONLY_IF_Sp> template <ONLY_IF_Sp>
static void OmegaInvariance(ColourMatrix &in) { static void OmegaInvariance(ColourMatrix &in) {

109
tests/core/Test_lie_generators.cc
View File

@ -28,6 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#include <Grid/Grid.h> #include <Grid/Grid.h>
#include <Grid/qcd/utils/CovariantCshift.h> #include <Grid/qcd/utils/CovariantCshift.h>
@ -42,7 +43,7 @@ directory
using namespace std; using namespace std;
using namespace Grid; using namespace Grid;
;
int main(int argc, char** argv) { int main(int argc, char** argv) {
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
@ -66,7 +67,7 @@ int main(int argc, char** argv) {
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
std::cout << GridLogMessage << "* Generators for SU(Nc" << std::endl; std::cout << GridLogMessage << "* Generators for SU(3)" << std::endl;
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
SU3::printGenerators(); SU3::printGenerators();
@ -87,22 +88,22 @@ int main(int argc, char** argv) {
// Projectors // Projectors
GridParallelRNG gridRNG(grid); GridParallelRNG gridRNG(grid);
gridRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9})); gridRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
SU3Adjoint::LatticeAdjMatrix Gauss(grid); SU_Adjoint<Nc>::LatticeAdjMatrix Gauss(grid);
SU3::LatticeAlgebraVector ha(grid); SU<Nc>::LatticeAlgebraVector ha(grid);
SU3::LatticeAlgebraVector hb(grid); SU<Nc>::LatticeAlgebraVector hb(grid);
random(gridRNG,Gauss); random(gridRNG,Gauss);
std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl;
SU3Adjoint::projectOnAlgebra(ha, Gauss); SU_Adjoint<Nc>::projectOnAlgebra(ha, Gauss);
std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl;
std::cout << GridLogMessage << "Start projector" << std::endl; std::cout << GridLogMessage << "Start projector" << std::endl;
SU3Adjoint::projector(hb, Gauss); SU_Adjoint<Nc>::projector(hb, Gauss);
std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "end projector" << std::endl;
std::cout << GridLogMessage << "ReStart projector" << std::endl; std::cout << GridLogMessage << "ReStart projector" << std::endl;
SU3Adjoint::projector(hb, Gauss); SU_Adjoint<Nc>::projector(hb, Gauss);
std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "end projector" << std::endl;
SU3::LatticeAlgebraVector diff = ha -hb; SU<Nc>::LatticeAlgebraVector diff = ha -hb;
std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl; std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl;
@ -114,8 +115,8 @@ int main(int argc, char** argv) {
LatticeGaugeField U(grid), V(grid); LatticeGaugeField U(grid), V(grid);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U); SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V); SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V);
// Adjoint representation // Adjoint representation
// Test group structure // Test group structure
@ -123,8 +124,8 @@ int main(int argc, char** argv) {
LatticeGaugeField UV(grid); LatticeGaugeField UV(grid);
UV = Zero(); UV = Zero();
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
SU3::LatticeMatrix Umu = peekLorentz(U,mu); SU<Nc>::LatticeMatrix Umu = peekLorentz(U,mu);
SU3::LatticeMatrix Vmu = peekLorentz(V,mu); SU<Nc>::LatticeMatrix Vmu = peekLorentz(V,mu);
pokeLorentz(UV,Umu*Vmu, mu); pokeLorentz(UV,Umu*Vmu, mu);
} }
@ -151,16 +152,16 @@ int main(int argc, char** argv) {
// Check correspondence of algebra and group transformations // Check correspondence of algebra and group transformations
// Create a random vector // Create a random vector
SU3::LatticeAlgebraVector h_adj(grid); SU<Nc>::LatticeAlgebraVector h_adj(grid);
typename AdjointRep<Nc>::LatticeMatrix Ar(grid); typename AdjointRep<Nc>::LatticeMatrix Ar(grid);
random(gridRNG,h_adj); random(gridRNG,h_adj);
h_adj = real(h_adj); h_adj = real(h_adj);
SU_Adjoint<Nc>::AdjointLieAlgebraMatrix(h_adj,Ar); SU_Adjoint<Nc>::AdjointLieAlgebraMatrix(h_adj,Ar);
// Re-extract h_adj // Re-extract h_adj
SU3::LatticeAlgebraVector h_adj2(grid); SU<Nc>::LatticeAlgebraVector h_adj2(grid);
SU_Adjoint<Nc>::projectOnAlgebra(h_adj2, Ar); SU_Adjoint<Nc>::projectOnAlgebra(h_adj2, Ar);
SU3::LatticeAlgebraVector h_diff = h_adj - h_adj2; SU<Nc>::LatticeAlgebraVector h_diff = h_adj - h_adj2;
std::cout << GridLogMessage << "Projections structure check vector difference (Adjoint representation) : " << norm2(h_diff) << std::endl; std::cout << GridLogMessage << "Projections structure check vector difference (Adjoint representation) : " << norm2(h_diff) << std::endl;
// Exponentiate // Exponentiate
@ -183,14 +184,14 @@ int main(int argc, char** argv) {
// Construct the fundamental matrix in the group // Construct the fundamental matrix in the group
SU3::LatticeMatrix Af(grid); SU<Nc>::LatticeMatrix Af(grid);
SU3::FundamentalLieAlgebraMatrix(h_adj,Af); SU<Nc>::FundamentalLieAlgebraMatrix(h_adj,Af);
SU3::LatticeMatrix Ufund(grid); SU<Nc>::LatticeMatrix Ufund(grid);
Ufund = expMat(Af, 1.0, 16); Ufund = expMat(Af, 1.0, 16);
// Check unitarity // Check unitarity
SU3::LatticeMatrix uno_f(grid); SU<Nc>::LatticeMatrix uno_f(grid);
uno_f = 1.0; uno_f = 1.0;
SU3::LatticeMatrix UnitCheck(grid); SU<Nc>::LatticeMatrix UnitCheck(grid);
UnitCheck = Ufund * adj(Ufund) - uno_f; UnitCheck = Ufund * adj(Ufund) - uno_f;
std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck) std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck)
<< std::endl; << std::endl;
@ -260,20 +261,20 @@ int main(int argc, char** argv) {
std::cout << GridLogMessage << "Test for the Two Index Symmetric projectors" std::cout << GridLogMessage << "Test for the Two Index Symmetric projectors"
<< std::endl; << std::endl;
// Projectors // Projectors
SU3TwoIndexSymm::LatticeTwoIndexMatrix Gauss2(grid); SU_TwoIndex<Nc, Symmetric>::LatticeTwoIndexMatrix Gauss2(grid);
random(gridRNG,Gauss2); random(gridRNG,Gauss2);
std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl;
SU3TwoIndexSymm::projectOnAlgebra(ha, Gauss2); SU_TwoIndex<Nc, Symmetric>::projectOnAlgebra(ha, Gauss2);
std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl;
std::cout << GridLogMessage << "Start projector" << std::endl; std::cout << GridLogMessage << "Start projector" << std::endl;
SU3TwoIndexSymm::projector(hb, Gauss2); SU_TwoIndex<Nc, Symmetric>::projector(hb, Gauss2);
std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "end projector" << std::endl;
std::cout << GridLogMessage << "ReStart projector" << std::endl; std::cout << GridLogMessage << "ReStart projector" << std::endl;
SU3TwoIndexSymm::projector(hb, Gauss2); SU_TwoIndex<Nc, Symmetric>::projector(hb, Gauss2);
std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "end projector" << std::endl;
SU3::LatticeAlgebraVector diff2 = ha - hb; SU<Nc>::LatticeAlgebraVector diff2 = ha - hb;
std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl; std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl;
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
@ -284,20 +285,20 @@ int main(int argc, char** argv) {
std::cout << GridLogMessage << "Test for the Two index anti-Symmetric projectors" std::cout << GridLogMessage << "Test for the Two index anti-Symmetric projectors"
<< std::endl; << std::endl;
// Projectors // Projectors
SU3TwoIndexAntiSymm::LatticeTwoIndexMatrix Gauss2a(grid); SU_TwoIndex<Nc, AntiSymmetric>::LatticeTwoIndexMatrix Gauss2a(grid);
random(gridRNG,Gauss2a); random(gridRNG,Gauss2a);
std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl;
SU3TwoIndexAntiSymm::projectOnAlgebra(ha, Gauss2a); SU_TwoIndex<Nc, AntiSymmetric>::projectOnAlgebra(ha, Gauss2a);
std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl;
std::cout << GridLogMessage << "Start projector" << std::endl; std::cout << GridLogMessage << "Start projector" << std::endl;
SU3TwoIndexAntiSymm::projector(hb, Gauss2a); SU_TwoIndex<Nc, AntiSymmetric>::projector(hb, Gauss2a);
std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "end projector" << std::endl;
std::cout << GridLogMessage << "ReStart projector" << std::endl; std::cout << GridLogMessage << "ReStart projector" << std::endl;
SU3TwoIndexAntiSymm::projector(hb, Gauss2a); SU_TwoIndex<Nc, AntiSymmetric>::projector(hb, Gauss2a);
std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "end projector" << std::endl;
SU3::LatticeAlgebraVector diff2a = ha - hb; SU<Nc>::LatticeAlgebraVector diff2a = ha - hb;
std::cout << GridLogMessage << "Difference: " << norm2(diff2a) << std::endl; std::cout << GridLogMessage << "Difference: " << norm2(diff2a) << std::endl;
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
@ -311,14 +312,14 @@ int main(int argc, char** argv) {
// Test group structure // Test group structure
// (U_f * V_f)_r = U_r * V_r // (U_f * V_f)_r = U_r * V_r
LatticeGaugeField U2(grid), V2(grid); LatticeGaugeField U2(grid), V2(grid);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2); SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U2);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2); SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V2);
LatticeGaugeField UV2(grid); LatticeGaugeField UV2(grid);
UV2 = Zero(); UV2 = Zero();
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
SU3::LatticeMatrix Umu2 = peekLorentz(U2,mu); SU<Nc>::LatticeMatrix Umu2 = peekLorentz(U2,mu);
SU3::LatticeMatrix Vmu2 = peekLorentz(V2,mu); SU<Nc>::LatticeMatrix Vmu2 = peekLorentz(V2,mu);
pokeLorentz(UV2,Umu2*Vmu2, mu); pokeLorentz(UV2,Umu2*Vmu2, mu);
} }
@ -345,16 +346,16 @@ int main(int argc, char** argv) {
// Check correspondence of algebra and group transformations // Check correspondence of algebra and group transformations
// Create a random vector // Create a random vector
SU3::LatticeAlgebraVector h_sym(grid); SU<Nc>::LatticeAlgebraVector h_sym(grid);
typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ar_sym(grid); typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ar_sym(grid);
random(gridRNG,h_sym); random(gridRNG,h_sym);
h_sym = real(h_sym); h_sym = real(h_sym);
SU_TwoIndex<Nc,Symmetric>::TwoIndexLieAlgebraMatrix(h_sym,Ar_sym); SU_TwoIndex<Nc,Symmetric>::TwoIndexLieAlgebraMatrix(h_sym,Ar_sym);
// Re-extract h_sym // Re-extract h_sym
SU3::LatticeAlgebraVector h_sym2(grid); SU<Nc>::LatticeAlgebraVector h_sym2(grid);
SU_TwoIndex< Nc, Symmetric>::projectOnAlgebra(h_sym2, Ar_sym); SU_TwoIndex< Nc, Symmetric>::projectOnAlgebra(h_sym2, Ar_sym);
SU3::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2; SU<Nc>::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2;
std::cout << GridLogMessage << "Projections structure check vector difference (Two Index Symmetric): " << norm2(h_diff_sym) << std::endl; std::cout << GridLogMessage << "Projections structure check vector difference (Two Index Symmetric): " << norm2(h_diff_sym) << std::endl;
@ -379,11 +380,11 @@ int main(int argc, char** argv) {
// Construct the fundamental matrix in the group // Construct the fundamental matrix in the group
SU3::LatticeMatrix Af_sym(grid); SU<Nc>::LatticeMatrix Af_sym(grid);
SU3::FundamentalLieAlgebraMatrix(h_sym,Af_sym); SU<Nc>::FundamentalLieAlgebraMatrix(h_sym,Af_sym);
SU3::LatticeMatrix Ufund2(grid); SU<Nc>::LatticeMatrix Ufund2(grid);
Ufund2 = expMat(Af_sym, 1.0, 16); Ufund2 = expMat(Af_sym, 1.0, 16);
SU3::LatticeMatrix UnitCheck2(grid); SU<Nc>::LatticeMatrix UnitCheck2(grid);
UnitCheck2 = Ufund2 * adj(Ufund2) - uno_f; UnitCheck2 = Ufund2 * adj(Ufund2) - uno_f;
std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2) std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2)
<< std::endl; << std::endl;
@ -421,14 +422,14 @@ int main(int argc, char** argv) {
// Test group structure // Test group structure
// (U_f * V_f)_r = U_r * V_r // (U_f * V_f)_r = U_r * V_r
LatticeGaugeField U2A(grid), V2A(grid); LatticeGaugeField U2A(grid), V2A(grid);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2A); SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U2A);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2A); SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V2A);
LatticeGaugeField UV2A(grid); LatticeGaugeField UV2A(grid);
UV2A = Zero(); UV2A = Zero();
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu); SU<Nc>::LatticeMatrix Umu2A = peekLorentz(U2,mu);
SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu); SU<Nc>::LatticeMatrix Vmu2A = peekLorentz(V2,mu);
pokeLorentz(UV2A,Umu2A*Vmu2A, mu); pokeLorentz(UV2A,Umu2A*Vmu2A, mu);
} }
@ -455,16 +456,16 @@ int main(int argc, char** argv) {
// Check correspondence of algebra and group transformations // Check correspondence of algebra and group transformations
// Create a random vector // Create a random vector
SU3::LatticeAlgebraVector h_Asym(grid); SU<Nc>::LatticeAlgebraVector h_Asym(grid);
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ar_Asym(grid); typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ar_Asym(grid);
random(gridRNG,h_Asym); random(gridRNG,h_Asym);
h_Asym = real(h_Asym); h_Asym = real(h_Asym);
SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym); SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym);
// Re-extract h_sym // Re-extract h_sym
SU3::LatticeAlgebraVector h_Asym2(grid); SU<Nc>::LatticeAlgebraVector h_Asym2(grid);
SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym); SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym);
SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2; SU<Nc>::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2;
std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl; std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl;
@ -489,11 +490,11 @@ int main(int argc, char** argv) {
// Construct the fundamental matrix in the group // Construct the fundamental matrix in the group
SU3::LatticeMatrix Af_Asym(grid); SU<Nc>::LatticeMatrix Af_Asym(grid);
SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym); SU<Nc>::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym);
SU3::LatticeMatrix Ufund2A(grid); SU<Nc>::LatticeMatrix Ufund2A(grid);
Ufund2A = expMat(Af_Asym, 1.0, 16); Ufund2A = expMat(Af_Asym, 1.0, 16);
SU3::LatticeMatrix UnitCheck2A(grid); SU<Nc>::LatticeMatrix UnitCheck2A(grid);
UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f; UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f;
std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A) std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A)
<< std::endl; << std::endl;
@ -527,6 +528,6 @@ int main(int argc, char** argv) {
Grid_finalize(); Grid_finalize();
} }

6
tests/core/Test_reunitarise.cc
View File

@ -26,6 +26,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
See the full license in the file "LICENSE" in the top level distribution directory See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#include <Grid/Grid.h> #include <Grid/Grid.h>
using namespace std; using namespace std;
@ -122,14 +123,14 @@ int main (int argc, char ** argv)
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
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
ProjectSUn(UU); ProjectSUn(UU);
detUU= Determinant(UU); detUU= Determinant(UU);
detUU= detUU -1.0; detUU= detUU -1.0;
@ -143,3 +144,4 @@ int main (int argc, char ** argv)

2
tests/solver/Test_coarse_even_odd.cc
View File

@ -93,7 +93,7 @@ int main(int argc, char** argv) {
// Setup of Dirac Matrix and Operator // // Setup of Dirac Matrix and Operator //
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////
LatticeGaugeField Umu(Grid_f); SU3::HotConfiguration(pRNG_f, Umu); LatticeGaugeField Umu(Grid_f); SU<Nc>::HotConfiguration(pRNG_f, Umu);
RealD checkTolerance = (getPrecision<LatticeFermion>::value == 1) ? 1e-7 : 1e-15; RealD checkTolerance = (getPrecision<LatticeFermion>::value == 1) ? 1e-7 : 1e-15;