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

add projection on sp2n algebra

Browse Source
This commit is contained in:
Alessandro Lupo
2021-10-26 10:20:44 +01:00
parent 4d8ae6221c
commit 4044536eea
3 changed files with 134 additions and 83 deletions
Download Patch File Download Diff File Expand all files Collapse all files

32
Grid/qcd/action/gauge/GaugeImplTypes.h
View File

@ -140,8 +140,38 @@ public:
PokeIndex<LorentzIndex>(P, Pmu, mu);
}
}
static inline Field projectForce(Field &P) { return Ta(P); }
static inline Field projectForce(Field &P)
{
if (isSp2n == true)
{
P = Ta(P);
const int nsp = Nc / 2;
Sp<nsp>::iSp2nMatrix<Complex> gen;
auto Psum = P;
Psum = Zero();
for (int a = 0; a < Sp<nsp>::AlgebraDimension; a++)
{
Sp<nsp>::generator(a, gen);
auto coeff = 2. * trace(P * gen);
Psum += coeff * gen;
}
return Psum;
} else
{
return Ta(P);
}
}
static inline void update_field(Field& P, Field& U, double ep){
//static std::chrono::duration<double> diff;

85
Grid/tensors/Tensor_Ta.h
View File

@ -137,6 +137,9 @@ 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)
{
iScalar<vtype> ret;
@ -217,88 +220,6 @@ accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg
return ret;
}
/*
template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg)
{
// need a check for the group type?
iMatrix<vtype,N> ret(arg);
vtype nrm;
vtype inner;
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]);
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++)
ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]); // lupo: c1,c2? dont we want c1 c2+N/2
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];
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 c1 0, ... N/2-1 loop
// 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