portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

project on sp2n

Browse Source
This commit is contained in:
Alessandro Lupo 2021-10-11 16:28:15 +01:00
parent 7ff3e5eed4
commit 7e3b298d3d
2 changed files with 110 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Grid/lattice/Lattice_ET.h
View File

@ -346,6 +346,7 @@ GridUnopClass(UnaryTrace, trace(a));
GridUnopClass(UnaryTranspose, transpose(a));
GridUnopClass(UnaryTa, Ta(a));
GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a));
GridUnopClass(UnaryProjectOnSpGroup, ProjectOnSpGroup(a));
GridUnopClass(UnaryTimesI, timesI(a));
GridUnopClass(UnaryTimesMinusI, timesMinusI(a));
GridUnopClass(UnaryAbs, abs(a));
@ -457,6 +458,7 @@ GRID_DEF_UNOP(trace, UnaryTrace);
GRID_DEF_UNOP(transpose, UnaryTranspose);
GRID_DEF_UNOP(Ta, UnaryTa);
GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup);
GRID_DEF_UNOP(ProjectOnSpGroup, UnaryProjectOnSpGroup);
GRID_DEF_UNOP(timesI, UnaryTimesI);
GRID_DEF_UNOP(timesMinusI, UnaryTimesMinusI);
GRID_DEF_UNOP(abs, UnaryAbs); // abs overloaded in cmath C++98; DON'T do the

108
Grid/tensors/Tensor_Ta.h
View File

@ -135,6 +135,114 @@ accelerator_inline iMatrix<vtype,N> ProjectOnGroup(const iMatrix<vtype,N> &arg)
return ret;
}
// re-do for sp2n
template<class vtype> accelerator_inline iScalar<vtype> ProjectOnSpGroup(const iScalar<vtype>&r)
{
iScalar<vtype> ret;
ret._internal = ProjectOnSpGroup(r._internal);
return ret;
}
template<class vtype,int N> accelerator_inline iVector<vtype,N> ProjectOnSpGroup(const iVector<vtype,N>&r)
{
iVector<vtype,N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = ProjectOnSpGroup(r._internal[i]);
}
return ret;
}
// int N is 2n in Sp(2n)
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;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++)
//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]);
// 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]);
}
}
// assuming the determinant is ok
return ret;
}
NAMESPACE_END(Grid);
#endif