portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-17 15:27:06 +01:00
Code Releases Activity

Correct implementation of SpTa

Browse Source
This commit is contained in:
Alessandro Lupo
2023-04-27 18:17:06 +01:00
parent 5aabe074fe
commit 881b08a465
2 changed files with 30 additions and 120 deletions
Download Patch File Download Diff File Expand all files Collapse all files

66
Grid/tensors/Tensor_Ta.h
View File

@ -82,67 +82,39 @@ template<class vtype,int N> accelerator_inline iVector<vtype,N> SpTa(const iVect
}
return ret;
}
template<class vtype,int N> accelerator_inline iMatrix<vtype,N> SpTa(const iMatrix<vtype,N> &arg)
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;
vtype inner;
vtype tmp;
nrm = 0.5;
ret = arg - (trace(arg)*factor);
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
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++)
{
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
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
}
}
return Ta(ret);
ret = (ret - adj(ret))*0.5;
return ret;
}