mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-09 23:45:36 +00:00
fix sp2n projection
This commit is contained in:
Alessandro Lupo
parent
d4866157fe
commit
8c0fbcccae
@ -161,88 +161,75 @@ accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg
|
||||
iMatrix<vtype,N> ret(arg);
|
||||
vtype nrm;
|
||||
vtype inner;
|
||||
for(int c1=0;c1<N;c1++){
|
||||
for(int c1=0;c1<N/2;c1++){
|
||||
|
||||
// Normalises row c1
|
||||
zeroit(inner);
|
||||
for(int c2=0;c2<N;c2++)
|
||||
inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
|
||||
// Normalises row c1
|
||||
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;
|
||||
nrm = sqrt(inner);
|
||||
nrm = 1.0/nrm;
|
||||
for(int c2=0;c2<N;c2++)
|
||||
ret._internal[c1][c2]*= nrm;
|
||||
|
||||
// Compute row c1+N2/2: c1+N/2 = - \Omega c1*
|
||||
for(int c2=0;c2<N/2;c2++)
|
||||
//tmp = conjugate(ret._internal[c1][c2]);
|
||||
//ret._internal[c1+N/2][c2+N/2] = tmp;
|
||||
ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]);
|
||||
// Compute row c1+N2/2: c1+N/2 = - \Omega c1*
|
||||
for(int c2=0;c2<N/2;c2++)
|
||||
ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]);
|
||||
|
||||
for(int c2=N/2;c2<N;c2++)
|
||||
//tmp = conjugate(ret._internal[c1][c2]);
|
||||
// ret._internal[c1+N/2][c2-N/2] = -tmp;
|
||||
ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);
|
||||
for(int c2=N/2;c2<N;c2++)
|
||||
ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);
|
||||
|
||||
// Remove c1 from rows c1+1...N/2-1
|
||||
for (int b=c1+1; b<N/2; ++b){
|
||||
decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
|
||||
zeroit(pr);
|
||||
for(int c=0; c<N; ++c)
|
||||
pr += conjugate(ret._internal[c1][c])*ret._internal[b][c];
|
||||
// Remove c1 from rows c1+1...N/2-1
|
||||
for (int b=c1+1; b<N/2; ++b){
|
||||
decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
|
||||
zeroit(pr);
|
||||
for(int c=0; c<N; ++c)
|
||||
pr += conjugate(ret._internal[c1][c])*ret._internal[b][c];
|
||||
|
||||
for(int c=0; c<N; ++c){
|
||||
ret._internal[b][c] -= pr * ret._internal[c1][c];
|
||||
}
|
||||
}
|
||||
|
||||
// Remove c1+N/2 from rows c1+1...N/2-1
|
||||
for (int b=c1+1; b<N; ++b){
|
||||
decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
|
||||
zeroit(pr);
|
||||
for(int c=0; c<N; ++c)
|
||||
pr += conjugate(ret._internal[c1+N/2][c])*ret._internal[b][c];
|
||||
|
||||
for(int c=0; c<N; ++c){
|
||||
ret._internal[b][c] -= pr * ret._internal[c1+N/2][c];
|
||||
}
|
||||
}
|
||||
} //end for loop over c1
|
||||
|
||||
// Compute the last row
|
||||
{
|
||||
// Normalise last row
|
||||
int c1 = N/2-1;
|
||||
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;
|
||||
}
|
||||
// Compute row cN-1+N2/2: cN-1+N/2 = - \Omega cN-1*
|
||||
for(int c2=0;c2<N/2;c2++)
|
||||
{
|
||||
//tmp = conjugate(ret._internal[c1][c2]);
|
||||
// ret._internal[c1+N/2][c2+N/2] = tmp;
|
||||
ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]);
|
||||
}
|
||||
for(int c2=N/2;c2<N;c2++)
|
||||
{
|
||||
//tmp = conjugate(ret._internal[c1][c2]);
|
||||
//ret._internal[c1+N/2][c2-N/2] = -tmp;
|
||||
ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);
|
||||
}
|
||||
for(int c=0; c<N; ++c){
|
||||
ret._internal[b][c] -= pr * ret._internal[c1][c];
|
||||
}
|
||||
// assuming the determinant is ok
|
||||
return ret;
|
||||
}
|
||||
|
||||
// Remove c1+N/2 from rows c1+1...N/2-1
|
||||
for (int b=c1+1; b<N; ++b){
|
||||
decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
|
||||
zeroit(pr);
|
||||
for(int c=0; c<N; ++c)
|
||||
pr += conjugate(ret._internal[c1+N/2][c])*ret._internal[b][c];
|
||||
|
||||
for(int c=0; c<N; ++c){
|
||||
ret._internal[b][c] -= pr * ret._internal[c1+N/2][c];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Compute the last row
|
||||
{
|
||||
// Normalise last row
|
||||
int c1 = N/2-1;
|
||||
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;
|
||||
|
||||
// Compute row cN-1+N2/2: cN-1+N/2 = - \Omega cN-1*
|
||||
for(int c2=0;c2<N/2;c2++)
|
||||
ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]);
|
||||
|
||||
for(int c2=N/2;c2<N;c2++)
|
||||
ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);
|
||||
}
|
||||
// assuming the determinant is ok
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
|
||||
#endif
|
||||
|
||||
@ -17,6 +17,14 @@ int main (int argc, char **argv)
|
||||
|
||||
LatticeGaugeField Umu(&Grid);
|
||||
LatticeColourMatrixD U(&Grid);
|
||||
LatticeColourMatrixD aux(&Grid);
|
||||
LatticeColourMatrixD identity(&Grid);
|
||||
|
||||
identity = 1.0;
|
||||
RealD epsilon = 0.01;
|
||||
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});
|
||||
@ -24,10 +32,122 @@ int main (int argc, char **argv)
|
||||
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(sseeds);
|
||||
|
||||
SU<Nc>::HotConfiguration(pRNG,Umu);
|
||||
U = PeekIndex<LorentzIndex>(Umu,0);
|
||||
U = PeekIndex<LorentzIndex>(Umu,2);
|
||||
|
||||
// it is unitary
|
||||
aux = U*adj(U) - identity;
|
||||
std::cout <<GridLogMessage << std::endl;
|
||||
std::cout << GridLogMessage << "Unitary check of starting random SUn matrix " << std::endl;
|
||||
std::cout <<GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
|
||||
std::cout <<GridLogMessage << std::endl;
|
||||
std::cout <<GridLogMessage << std::endl;
|
||||
|
||||
U = U + epsilon*identity;
|
||||
aux = U*adj(U) - identity;
|
||||
|
||||
std::cout << GridLogMessage << "Unitary matrix modified " << 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;
|
||||
|
||||
U = ProjectOnSpGroup(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 << "U adjU - 1 = " << norm2(aux) << std::endl;
|
||||
assert( norm2(aux) < 1e-6);
|
||||
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;
|
||||
// checks on elements
|
||||
// an Sp2n matrix is a block matrix of the type
|
||||
//
|
||||
//
|
||||
// U = ( W X )
|
||||
// ( -X^* W^* )
|
||||
//
|
||||
// will check this
|
||||
//
|
||||
std::cout << GridLogMessage << "Checking that U[ i , j ] - conj U[ i+N/2, j+N/2 ] = 0 for i = 0, ... , N/2 " << std::endl;
|
||||
std::cout <<GridLogMessage << std::endl;
|
||||
const int nsp = Nc / 2;
|
||||
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);
|
||||
std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
|
||||
assert( amizeroo.real() < 10e-6 );
|
||||
amizeroo *= i;
|
||||
assert( amizeroo.real() < 10e-6 );
|
||||
std::cout << GridLogMessage << "ok " << std::endl;
|
||||
}
|
||||
|
||||
}
|
||||
std::cout <<GridLogMessage << std::endl;
|
||||
std::cout << GridLogMessage << "Checking that U [ i , j+N/2 ] + conj U [ i+N/2, j+N/2 ] = 0 for i = 0, ... , N/2 " << std::endl;
|
||||
std::cout <<GridLogMessage << std::endl;
|
||||
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);
|
||||
std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
|
||||
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.);
|
||||
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();
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user