portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-13 01:05:36 +00:00
Code Releases Activity

Merge in Alessandro's changes [test fails]

Browse Source
This commit is contained in:
Julian Lenz 2023-05-03 02:53:03 +01:00
commit ac6c7cb8d6
8 changed files with 86 additions and 144 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
Grid/qcd/utils/GaugeGroup.h
View File

@ -303,7 +303,6 @@ class GaugeGroup {
template <typename LatticeMatrixType> template <typename LatticeMatrixType>
static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) { static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) {
taProj(in, out, group_name()); taProj(in, out, group_name());
//out = Ta(in);
} }
template <typename LatticeMatrixType> template <typename LatticeMatrixType>
@ -337,32 +336,20 @@ class GaugeGroup {
ProjectOnGaugeGroup(Umu, group_name()); ProjectOnGaugeGroup(Umu, group_name());
} }
template <class vtype> template <class vtype>
static void ProjectOnGaugeGroup(iScalar<vtype> &r) { static void ProjectOnGaugeGroup(iScalar<vtype> &r) {
ProjectOnGaugeGroup(r, group_name()); ProjectOnGaugeGroup(r, group_name());
}
template <class vtype, int N>
static void ProjectOnGaugeGroup(iVector<vtype,N> &r) {
r = ProjectOnGaugeGroup(r, group_name());
}
template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
static void ProjectOnGaugeGroup(iMatrix<vtype,N> &arg) {
arg = ProjectOnGaugeGroup(arg, group_name());
}
/*
template <int N> // reunitarise, resimplectify...
static void ProjectOnGaugeGroup(iVector<iScalar<iMatrix<vComplexD, N> >, Nd> &U) {
ProjectOnGaugeGroup(U, group_name());
} }
template <int N> // reunitarise, resimplectify... template <class vtype, int N>
static void ProjectOnGaugeGroup(iScalar<iScalar<iMatrix<vComplexD, N> > > &Umu) { static void ProjectOnGaugeGroup(iVector<vtype,N> &r) {
ProjectOnGaugeGroup(Umu, group_name()); ProjectOnGaugeGroup(r, group_name());
}*/ }
template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
static void ProjectOnGaugeGroup(iMatrix<vtype,N> &arg) {
ProjectOnGaugeGroup(arg, group_name());
}
template <int N> // reunitarise, resimplectify... previously ProjectSUn template <int N> // reunitarise, resimplectify... previously ProjectSUn
static void ProjectGn(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) { static void ProjectGn(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
@ -416,7 +403,6 @@ LatticeComplexD Determinant(
return ret; return ret;
} }
template <int N> template <int N>
static void ProjectSUn( static void ProjectSUn(
Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) { Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
@ -445,7 +431,6 @@ static void ProjectSUn(
} }
} }
// Explicit specialisation for SU(3).
// Explicit specialisation for SU(3). // Explicit specialisation for SU(3).
static void ProjectSU3( static void ProjectSU3(
Lattice<iScalar<iScalar<iMatrix<vComplexD, 3> > > > &Umu) { Lattice<iScalar<iScalar<iMatrix<vComplexD, 3> > > > &Umu) {

3
Grid/qcd/utils/SUnAdjoint.h
View File

@ -68,8 +68,7 @@ public:
for (int a = 0; a < Dimension; a++) { for (int a = 0; a < Dimension; a++) {
tmp = ta[a] * ta[Index] - ta[Index] * ta[a]; tmp = ta[a] * ta[Index] - ta[Index] * ta[a];
for (int b = 0; b < (ncolour * ncolour - 1); b++) { for (int b = 0; b < (ncolour * ncolour - 1); b++) {
iSUnMatrix<cplx> tmp1 = iSUnMatrix<cplx> tmp1 = 2.0 * tmp * ta[b]; // 2.0 from the normalization
2.0 * tmp * ta[b]; // 2.0 from the normalization
Complex iTr = TensorRemove(timesI(trace(tmp1))); Complex iTr = TensorRemove(timesI(trace(tmp1)));
//iAdjTa()()(b, a) = iTr; //iAdjTa()()(b, a) = iTr;
iAdjTa()()(a, b) = iTr; iAdjTa()()(a, b) = iTr;

3
Grid/qcd/utils/SUnTwoIndex.h
View File

@ -126,8 +126,7 @@ public:
template <class cplx> template <class cplx>
static void generator(int Index, iSUnTwoIndexMatrix<cplx> &i2indTa) { static void generator(int Index, iSUnTwoIndexMatrix<cplx> &i2indTa) {
Vector<iSUnMatrix<cplx> > ta( Vector<iSUnMatrix<cplx> > ta(ncolour * ncolour - 1);
ncolour * ncolour - 1);
Vector<iSUnMatrix<cplx> > eij(Dimension); Vector<iSUnMatrix<cplx> > eij(Dimension);
iSUnMatrix<cplx> tmp; iSUnMatrix<cplx> tmp;
i2indTa = Zero(); i2indTa = Zero();

168
Grid/tensors/Tensor_Ta.h
View File

@ -66,8 +66,6 @@ template<class vtype,int N> accelerator_inline iMatrix<vtype,N> Ta(const iMatrix
return ret; return ret;
} }
// for sp2n can't be as simple as Ta. We do a Gram-Schmidt
template<class vtype> accelerator_inline iScalar<vtype> SpTa(const iScalar<vtype>&r) template<class vtype> accelerator_inline iScalar<vtype> SpTa(const iScalar<vtype>&r)
{ {
iScalar<vtype> ret; iScalar<vtype> ret;
@ -82,77 +80,47 @@ template<class vtype,int N> accelerator_inline iVector<vtype,N> SpTa(const iVect
} }
return ret; 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)
{ {
iMatrix<vtype,N> ret(arg); // Generalises Ta to Sp2n
vtype nrm; // Applies the following projections
vtype inner; // P_{antihermitian} P_{antihermitian-Sp-algebra} P_{traceless}
vtype tmp; // 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;
nrm = 0.5;
for(int c1=0;c1<N/2;c1++) ret = arg - (trace(arg)*factor);
{
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
}
}
for(int c1=0;c1<N/2;c1++)
return Ta(ret); {
for(int c2=0;c2<N/2;c2++)
{
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++)
{
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
}
}
ret = (ret - adj(ret))*0.5;
return ret;
} }
/////////////////////////////////////////////// ///////////////////////////////////////////////
// ProjectOnGroup function for scalar, vector, matrix // ProjectOnGroup function for scalar, vector, matrix
// Projects on orthogonal, unitary group // Projects on orthogonal, unitary group
/////////////////////////////////////////////// ///////////////////////////////////////////////
template<class vtype> accelerator_inline iScalar<vtype> ProjectOnGroup(const iScalar<vtype>&r) template<class vtype> accelerator_inline iScalar<vtype> ProjectOnGroup(const iScalar<vtype>&r)
{ {
iScalar<vtype> ret; iScalar<vtype> ret;
@ -244,58 +212,52 @@ accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg
iMatrix<vtype,N> ret(arg); iMatrix<vtype,N> ret(arg);
vtype nrm; vtype nrm;
vtype inner; vtype inner;
vtype tmp;
for(int c1=0;c1<N/2;c1++) for(int c1=0;c1<N/2;c1++)
{ {
for (int b=0; b<c1; b++) // remove the b-rows from U_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]) pr;
decltype(ret._internal[b][b]*ret._internal[b][b]) prn; decltype(ret._internal[b][b]*ret._internal[b][b]) prn;
zeroit(pr); zeroit(pr);
zeroit(prn); zeroit(prn);
for(int c=0; c<N; c++) for(int c=0; c<N; c++)
{ {
pr += conjugate(ret._internal[c1][c])*ret._internal[b][c]; // <U_c1 | U_b > 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} > prn += conjugate(ret._internal[c1][c])*ret._internal[b+N/2][c]; // <U_c1 | U_{b+N} >
} }
for(int c=0; c<N; c++) 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} ) 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); zeroit(inner);
for(int c2=0;c2<N;c2++) for(int c2=0;c2<N;c2++)
{ {
inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]); inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
} }
nrm = sqrt(inner); nrm = sqrt(inner);
nrm = 1.0/nrm; nrm = 1.0/nrm;
for(int c2=0;c2<N;c2++) for(int c2=0;c2<N;c2++)
{ {
ret._internal[c1][c2]*= nrm; ret._internal[c1][c2]*= nrm;
} }
for(int c2=0;c2<N/2;c2++)
for(int c2=0;c2<N/2;c2++) {
{ ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]); // down right in the new matrix = (up-left)* of the old matrix
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
}
for(int c2=N/2;c2<N;c2++)
{
ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);; // down left in the new matrix = -(up-right)* of the old
}
} }
return ret; return ret;
} }

2
tests/core/Test_lie_generators.cc
View File

@ -47,7 +47,7 @@ using namespace Grid;
int main(int argc, char** argv) { int main(int argc, char** argv) {
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
std::vector<int> latt({4, 4, 4, 8}); std::vector<int> latt({4, 4, 4, 8});
GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid( GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid(
latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());

2
tests/core/Test_reunitarise.cc
View File

@ -36,7 +36,7 @@ using namespace Grid;
int main (int argc, char ** argv) int main (int argc, char ** argv)
{ {
Grid_init(&argc,&argv); Grid_init(&argc,&argv);
std::vector<int> latt({8,8,8,8}); std::vector<int> latt({8,8,8,8});
GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt, GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt,
GridDefaultSimd(Nd,vComplexD::Nsimd()), GridDefaultSimd(Nd,vComplexD::Nsimd()),

2
tests/sp2n/Test_project_on_Sp.cc
View File

@ -144,10 +144,8 @@ bool is_element_of_sp2n_algebra(T U) {
aux = U - adj(U); aux = U - adj(U);
std::cout << GridLogMessage << "T - Tda = " << norm2(aux) << std::endl; std::cout << GridLogMessage << "T - Tda = " << norm2(aux) << std::endl;
assert(norm2(aux) < 1e-8);
aux = U + adj(U); aux = U + adj(U);
std::cout << GridLogMessage << "T + Tda = " << norm2(aux) << std::endl; std::cout << GridLogMessage << "T + Tda = " << norm2(aux) << std::endl;
assert(norm2(aux) < 1e-8);
std::cout << GridLogMessage << "Check that Omega U Omega = conj(U)" std::cout << GridLogMessage << "Check that Omega U Omega = conj(U)"
<< std::endl; << std::endl;

11
tests/sp2n/Test_sp2n_lie_gen.cc
View File

@ -10,14 +10,14 @@ int main(int argc, char** argv) {
//latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); //latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
//GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid); //GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid);
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
std::cout << GridLogMessage << "* Generators for Sp(2)" << std::endl; std::cout << GridLogMessage << "* Generators for Sp(2)" << std::endl;
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
Sp2::printGenerators(); Sp2::printGenerators();
Sp2::testGenerators(); Sp2::testGenerators();
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
@ -46,6 +46,5 @@ int main(int argc, char** argv) {
Sp8::printGenerators(); Sp8::printGenerators();
Sp8::testGenerators(); Sp8::testGenerators();
Grid_finalize(); Grid_finalize();
} }