mirror of
https://github.com/paboyle/Grid.git
synced 2025-10-31 03:54:33 +00:00
Revert "projection on Sp2n algebra, to be used instead of Ta"
This reverts commit ba7f9d7b70.
This commit is contained in:
Alessandro Lupo
@@ -345,7 +345,6 @@ GridUnopClass(UnaryNot, Not(a));
|
||||
GridUnopClass(UnaryTrace, trace(a));
|
||||
GridUnopClass(UnaryTranspose, transpose(a));
|
||||
GridUnopClass(UnaryTa, Ta(a));
|
||||
GridUnopClass(UnaryProjectSp2nAlgebra, ProjectSp2nAlgebra(a));
|
||||
GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a));
|
||||
GridUnopClass(UnaryProjectOnSpGroup, ProjectOnSpGroup(a));
|
||||
GridUnopClass(UnaryTimesI, timesI(a));
|
||||
@@ -458,7 +457,6 @@ GRID_DEF_UNOP(operator!, UnaryNot);
|
||||
GRID_DEF_UNOP(trace, UnaryTrace);
|
||||
GRID_DEF_UNOP(transpose, UnaryTranspose);
|
||||
GRID_DEF_UNOP(Ta, UnaryTa);
|
||||
GRID_DEF_UNOP(ProjectSp2nAlgebra, UnaryProjectSp2nAlgebra);
|
||||
GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup);
|
||||
GRID_DEF_UNOP(ProjectOnSpGroup, UnaryProjectOnSpGroup);
|
||||
GRID_DEF_UNOP(timesI, UnaryTimesI);
|
||||
|
||||
@@ -414,6 +414,13 @@ static void OmegaInvariance(ColourMatrix &in) {
|
||||
Omega = Zero();
|
||||
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++) {
|
||||
Omega()()(i, nsp + i) = 1.;
|
||||
Omega()()(nsp + i, i) = -1;
|
||||
@@ -447,6 +454,8 @@ static void OmegaInvariance(GaugeField &in) {
|
||||
Omega = Zero();
|
||||
identity = 1.;
|
||||
|
||||
std::cout << GridLogMessage << "I am a GaugeField " << std::endl;
|
||||
|
||||
U = PeekIndex<LorentzIndex>(in, 1);
|
||||
|
||||
OmegaInvariance(U);
|
||||
@@ -464,6 +473,8 @@ static void OmegaInvariance(LatticeColourMatrixD &in) {
|
||||
Omega = Zero();
|
||||
identity = 1.;
|
||||
|
||||
std::cout << GridLogMessage << "I am a LatticeColourMatrix " << std::endl;
|
||||
|
||||
for (int i = 0; i < nsp; i++) {
|
||||
Omega()()(i, nsp + i) = 1.;
|
||||
Omega()()(nsp + i, i) = -1;
|
||||
@@ -485,21 +496,3 @@ 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;
|
||||
}
|
||||
|
||||
@@ -66,88 +66,6 @@ template<class vtype,int N> accelerator_inline iMatrix<vtype,N> Ta(const iMatrix
|
||||
return ret;
|
||||
}
|
||||
|
||||
// for sp2n can't do something as simple as Ta
|
||||
|
||||
template<class vtype> accelerator_inline iScalar<vtype> ProjectSp2nAlgebra(const iScalar<vtype>&r)
|
||||
{
|
||||
iScalar<vtype> ret;
|
||||
ret._internal = ProjectSp2nAlgebra(r._internal);
|
||||
return ret;
|
||||
}
|
||||
template<class vtype,int N> accelerator_inline iVector<vtype,N> ProjectSp2nAlgebra(const iVector<vtype,N>&r)
|
||||
{
|
||||
iVector<vtype,N> ret;
|
||||
for(int i=0;i<N;i++){
|
||||
ret._internal[i] = ProjectSp2nAlgebra(r._internal[i]);
|
||||
}
|
||||
return ret;
|
||||
}
|
||||
template<class vtype,int N> accelerator_inline iMatrix<vtype,N> ProjectSp2nAlgebra(const iMatrix<vtype,N> &arg)
|
||||
{
|
||||
iMatrix<vtype,N> ret;
|
||||
// following HiRep's suN_utils.c does a sort of Gram-Schmidt
|
||||
vtype nrm;
|
||||
vtype inner;
|
||||
vtype tmp;
|
||||
|
||||
for(int c1=0;c1<N/2;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]) prn;
|
||||
zeroit(pr);
|
||||
zeroit(prn);
|
||||
|
||||
for(int c=0; c<N; c++)
|
||||
{
|
||||
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} >
|
||||
}
|
||||
|
||||
|
||||
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} )
|
||||
}
|
||||
}
|
||||
|
||||
zeroit(inner);
|
||||
for(int c2=0;c2<N;c2++)
|
||||
{
|
||||
inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
|
||||
}
|
||||
|
||||
nrm = sqrt(inner);
|
||||
nrm = 1.0/nrm;
|
||||
for(int c2=0;c2<N;c2++)
|
||||
{
|
||||
ret._internal[c1][c2]*= nrm;
|
||||
}
|
||||
|
||||
|
||||
for(int c2=0;c2<N/2;c2++)
|
||||
{
|
||||
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
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
///////////////////////////////////////////////
|
||||
// ProjectOnGroup function for scalar, vector, matrix
|
||||
@@ -219,6 +137,8 @@ accelerator_inline iMatrix<vtype,N> ProjectOnGroup(const iMatrix<vtype,N> &arg)
|
||||
|
||||
// re-do for sp2n
|
||||
|
||||
// Ta cannot be defined here for Sp2n because I need the generators from the Sp class
|
||||
// It is defined in gauge impl types
|
||||
|
||||
template<class vtype> accelerator_inline iScalar<vtype> ProjectOnSpGroup(const iScalar<vtype>&r)
|
||||
{
|
||||
|
||||
@@ -20,8 +20,6 @@ int main (int argc, char **argv)
|
||||
LatticeColourMatrixD aux(&Grid);
|
||||
LatticeColourMatrixD identity(&Grid);
|
||||
|
||||
// Will test resimplectification (from ProjectOnGaugeGroup, ProjectOnSpGroup) and projection on the algebra (from ProjectSp2nAlgebra)
|
||||
|
||||
const int nsp = Nc / 2;
|
||||
|
||||
identity = 1.0;
|
||||
@@ -31,12 +29,15 @@ 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,0);
|
||||
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;
|
||||
@@ -54,27 +55,35 @@ 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 << "Simplectify" << std::endl;
|
||||
std::cout <<GridLogMessage << std::endl;
|
||||
|
||||
std::cout <<GridLogMessage << std::endl;
|
||||
std::cout << GridLogMessage << "Projecting on Sp2n " << 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 simplectification " << std::endl;
|
||||
std::cout << GridLogMessage << "Unitary check after Sp(2n) projection " << 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 simplectification = " << norm2( Determinant(U) ) / vol << std::endl;
|
||||
assert( norm2(aux) - 1 < 1e-8);
|
||||
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 simplectification. Will not kill if it fails, but the next check will " << std::endl;
|
||||
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;
|
||||
@@ -112,110 +121,34 @@ int main (int argc, char **argv)
|
||||
|
||||
std::cout << GridLogMessage << "ok" << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "Check the ProjectOnGaugeGroup function to simplectify" << std::endl;
|
||||
U = U + 2932.111*identity;
|
||||
std::cout << GridLogMessage << "Resimplectify 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
|
||||
// an explicit check for sp2
|
||||
/*
|
||||
if (Nc == 2)
|
||||
{
|
||||
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 );
|
||||
}
|
||||
|
||||
}
|
||||
assert(Nc==2);
|
||||
ColourMatrix A;
|
||||
A = Zero();
|
||||
|
||||
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 << "ok" << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "Now check projection on the algebra" << 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 << "Project on sp2n algebra" << 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 );
|
||||
}
|
||||
}
|
||||
|
||||
Complex a(25041994., 12.);
|
||||
Complex b(39., 0.22);
|
||||
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();
|
||||
|
||||
|
||||
Reference in New Issue
Block a user