portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-04 05:54:32 +00:00
Code Releases Activity

Projection on algebra can now be templated. Fix #12

Browse Source
This commit is contained in:
Alessandro Lupo
2023-04-03 16:31:19 +01:00
parent 4275b3f431
commit 026e736dfa
4 changed files with 261 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files

211
tests/sp2n/Test_project_on_Sp.cc
View File

@@ -20,6 +20,12 @@ int main (int argc, char **argv)
LatticeColourMatrixD aux(&Grid);
LatticeColourMatrixD identity(&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
const int nsp = Nc / 2;
identity = 1.0;
@@ -29,15 +35,12 @@ int main (int argc, char **argv)
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);
U = PeekIndex<LorentzIndex>(Umu,0);
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;
@@ -55,35 +58,23 @@ int main (int argc, char **argv)
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;
std::cout << GridLogMessage << "Projecting on Sp2n " << std::endl;
// Testing ProjectOnSpGroup
std::cout << GridLogMessage << "Testing ProjectOnSpGroup" << std::endl;
std::cout << GridLogMessage << "Apply ProjectOnSpGroup to deformed matrix" << std::endl;
U = ProjectOnSpGroup(U);
//U = ProjectOnGroup(U);
aux = U*adj(U) - identity;
std::cout <<GridLogMessage << std::endl;
std::cout <<GridLogMessage << std::endl;
std::cout << GridLogMessage << "Unitary check after Sp(2n) projection " << std::endl;
std::cout << GridLogMessage << "Unitary check after ProjectOnSpGroup " << std::endl;
std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
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);
std::cout << GridLogMessage << "Checking Omega invariance after ProjectOnSpGroup" << std::endl;
Sp<Nc>::OmegaInvariance(U); // no assertion here, but the next check will kill us if we are not simplectic
// 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;
@@ -119,38 +110,158 @@ int main (int argc, char **argv)
}
}
std::cout << GridLogMessage << "ok" << std::endl;
// an explicit check for sp2
/*
if (Nc == 2)
std::cout << GridLogMessage << "Testing ProjectOnGaugeGroup" << std::endl;
U = U + 2932.111*identity;
std::cout << GridLogMessage << "Apply ProjectOnGaugeGroup to deformed matrix" << std::endl;
Sp<Nc>::ProjectOnGaugeGroup(U);
aux = U*adj(U) - identity;
std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
assert( norm2(aux) < 1e-8);
Sp<Nc>::OmegaInvariance(U);
std::cout << GridLogMessage << "Checking the structure is " << std::endl;
std::cout << GridLogMessage << "U = ( W X ) " << std::endl;
std::cout << GridLogMessage << " ( -X^* W^* ) " << std::endl;
std::cout <<GridLogMessage << std::endl;
for (int c1 = 0; c1 < nsp; c1++) //check on W
{
assert(Nc==2);
ColourMatrix A;
A = Zero();
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 );
}
}
Complex a(25041994., 12.);
Complex b(39., 0.22);
Complex d(10000., -2222.3333);
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 );
}
}
std::cout << GridLogMessage << "Testing ProjectGn" << std::endl;
U = U + 2932.111*identity;
std::cout << GridLogMessage << "Apply ProjectGn to deformed matrix" << std::endl;
Sp<Nc>::ProjectGn(U);
aux = U*adj(U) - identity;
std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
assert( norm2(aux) < 1e-8);
std::cout << GridLogMessage << "Det after ProjectGn = " << norm2( Determinant(U) ) / vol << std::endl;
assert( norm2(aux) - 1 < 1e-8);
Sp<Nc>::OmegaInvariance(U);
std::cout << GridLogMessage << "Checking the structure is " << std::endl;
std::cout << GridLogMessage << "U = ( W X ) " << std::endl;
std::cout << GridLogMessage << " ( -X^* W^* ) " << std::endl;
std::cout <<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 );
}
}
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 );
}
}
std::cout << GridLogMessage << "Testing ProjectSp2nAlgebra" << std::endl;
U = PeekIndex<LorentzIndex>(Umu,1);
U = U + 666.*identity;
aux = U*adj(U) - identity;
std::cout << GridLogMessage << "Matrix deformed " << std::endl;
std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
std::cout << GridLogMessage << "Apply ProjectSp2nAlgebra to deformed matrix" << std::endl;
U = ProjectSp2nAlgebra(U);
aux = U*adj(U) - identity;
std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
assert( norm2(aux) < 1e-8);
std::cout << GridLogMessage << "Check that Omega U Omega = conj(U)" << std::endl;
LatticeColourMatrixD Omega(&Grid);
Sp<Nc>::Omega(Omega);
aux = Omega*U*Omega - conjugate(U);
std::cout << GridLogMessage << "Omega U Omega - conj(U) = " << norm2(aux) << std::endl;
assert( norm2(aux) < 1e-8);
std::cout << GridLogMessage << "Checking the structure is " << std::endl;
std::cout << GridLogMessage << "U = ( W X ) " << std::endl;
std::cout << GridLogMessage << " ( X^* -W^* ) " << std::endl;
std::cout <<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 );
}
}
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 );
}
}
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();
}