1
0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00

Attempt with SFINAE (failed)

This commit is contained in:
Julian Lenz 2022-11-30 11:57:39 +00:00
parent 54c8025aad
commit 7c2ad4f8c8
4 changed files with 866 additions and 584 deletions

View File

@ -33,12 +33,19 @@ directory
#define QCD_UTIL_SUN_H
#define ONLY_IF_SU \
typename dummy_name = group_name, \
typename = std::enable_if_t<is_su<dummy_name>::value>
typename \
typename std::enable_if_t<is_su<group_name>::value, dummy_name> ,\
dummy_name = group_name,
#define ONLY_IF_Sp \
typename \
typename std::enable_if_t<is_sp<group_name>::value, dummy_name> ,\
dummy_name =group_name,
NAMESPACE_BEGIN(Grid);
namespace GroupName {
class SU {};
class Sp {};
} // namespace GroupName
template <typename group_name>
@ -51,18 +58,44 @@ struct is_su<GroupName::SU> {
static const bool value = true;
};
template <typename group_name>
struct is_sp {
static const bool value = false;
};
template <>
struct is_sp<GroupName::Sp> {
static const bool value = true;
};
template<typename group_name>
constexpr int compute_adjoint_dimension(int ncolour);
template<>
constexpr int compute_adjoint_dimension<GroupName::SU>(int ncolour) { return ncolour * ncolour - 1;}
template<>
constexpr int compute_adjoint_dimension<GroupName::Sp>(int ncolour) { return ncolour/2*(ncolour+1);}
template <int ncolour, class group_name = GroupName::SU>
class GaugeGroup;
template <int ncolour>
using SU = GaugeGroup<ncolour, GroupName::SU>;
template <int ncolour>
using Sp = GaugeGroup<ncolour, GroupName::Sp>;
template <int ncolour, class group_name>
class GaugeGroup {
public:
static const int Dimension = ncolour;
static const int AdjointDimension = ncolour * ncolour - 1;
static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; }
static const int AdjointDimension = compute_adjoint_dimension<group_name>(ncolour);
static const int AlgebraDimension = compute_adjoint_dimension<group_name>(ncolour);
// Don't know how to only enable this for Sp:
static const int nsp = ncolour/2;
template <typename vtype>
using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >;
@ -115,6 +148,7 @@ class GaugeGroup {
typedef Lattice<vSU2MatrixD> LatticeSU2MatrixD;
#include "Grid/qcd/utils/SUn_impl.h"
#include "Grid/qcd/utils/Sp2n_impl.h"
static void printGenerators(void) {
for (int gen = 0; gen < AdjointDimension; gen++) {

View File

@ -6,6 +6,8 @@
//
// around it.
template<ONLY_IF_SU>
static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; }
template <typename vtype, ONLY_IF_SU>
using iSUnMatrix = iGroupMatrix<vtype>;
template <typename vtype, ONLY_IF_SU>

File diff suppressed because it is too large Load Diff

245
Grid/qcd/utils/Sp2n_impl.h Normal file
View File

@ -0,0 +1,245 @@
template<ONLY_IF_Sp>
static int su2subgroups(void) { return (nsp * (nsp - 1)) / 2; }
template <typename vtype, ONLY_IF_Sp>
using iSp2nMatrix = iGroupMatrix<vtype>;
template <typename vtype, ONLY_IF_Sp>
using iSp2nAlgebraVector = iAlgebraVector<vtype>;
// Sp(2N) has N(2N+1) = 2N^2+N generators
//
// normalise the generators such that
// Trace ( Ta Tb) = 1/2 delta_ab
//
// N generators in the cartan, 2N^2 off
// off diagonal:
// there are 6 types named a,b,c,d and w,z
// abcd are N(N-1)/2 each while wz are N each
template <class cplx, ONLY_IF_Sp>
static void generator(int lieIndex, iSp2nMatrix<cplx> &ta) {
// map lie index into type of generators: diagonal, abcd type, wz type
int diagIndex;
int aIndex, bIndex, cIndex, dIndex;
int wIndex, zIndex; // a,b,c,d are N(N-1)/2 and w,z are N
int mod = nsp * (nsp-1) * 0.5;
int offdiag = 2*nsp*nsp; // number of generators not in the cartan subalgebra
int wmod = 4*mod;
int zmod = wmod+nsp;
if (lieIndex >= offdiag) {
diagIndex = lieIndex - offdiag; // 0, ... ,N-1
//std::cout << GridLogMessage << "diag type " << std::endl;
generatorDiagtype(diagIndex, ta);
return;
}
if ( (lieIndex >= wmod) && (lieIndex < zmod) ) {
//std::cout << GridLogMessage << "w type " << std::endl;
wIndex = lieIndex- wmod; // 0, ... ,N-1
generatorWtype(wIndex,ta);
return;
}
if ( (lieIndex >= zmod) && (lieIndex < offdiag) ) {
//std::cout << GridLogMessage << "z type " << std::endl;
//std::cout << GridLogMessage << "lie index " << lieIndex << std::endl;
//std::cout << GridLogMessage << "z mod " << zmod << std::endl;
zIndex = lieIndex - zmod; // 0, ... ,N-1
generatorZtype(zIndex,ta);
return;
}
if (lieIndex < mod) { // atype 0, ... , N(N-1)/2=mod
//std::cout << GridLogMessage << "a type " << std::endl;
aIndex = lieIndex;
//std::cout << GridLogMessage << "a indx " << aIndex << std::endl;
generatorAtype(aIndex, ta);
return;
}
if ( (lieIndex >= mod) && lieIndex < 2*mod) { // btype mod, ... , 2mod-1
//std::cout << GridLogMessage << "b type " << std::endl;
bIndex = lieIndex - mod;
generatorBtype(bIndex, ta);
return;
}
if ( (lieIndex >= 2*mod) && lieIndex < 3*mod) { // ctype 2mod, ... , 3mod-1
//std::cout << GridLogMessage << "c type " << std::endl;
cIndex = lieIndex - 2*mod;
generatorCtype(cIndex, ta);
return;
}
if ( (lieIndex >= 3*mod) && lieIndex < wmod) { // ctype 3mod, ... , 4mod-1 = wmod-1
//std::cout << GridLogMessage << "d type " << std::endl;
dIndex = lieIndex - 3*mod;
generatorDtype(dIndex, ta);
return;
}
} //end of generator
template <class cplx, ONLY_IF_Sp>
static void generatorDiagtype(int diagIndex, iSp2nMatrix<cplx> &ta) {
// ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra
ta = Zero();
RealD nrm = 1.0 / 2;
ta()()(diagIndex,diagIndex) = nrm;
ta()()(diagIndex+nsp,diagIndex+nsp) = -nrm;
}
template <class cplx, ONLY_IF_Sp>
static void generatorAtype(int aIndex, iSp2nMatrix<cplx> &ta) {
// ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2)
// with i<j and i=0,...,N-2
// follows that j=i+1, ... , N
int i1, i2;
ta = Zero();
RealD nrm = 1 / (2 * std::sqrt(2) );
su2SubGroupIndex(i1, i2, aIndex);
ta()()(i1,i2) = 1;
ta()()(i2,i1) = 1;
ta()()(i1+nsp,i2+nsp) = -1;
ta()()(i2+nsp,i1+nsp) = -1;
ta = ta * nrm;
}
template <class cplx, ONLY_IF_Sp>
static void generatorBtype(int bIndex, iSp2nMatrix<cplx> &ta) {
// ta(i,j) = -ta(j,i) = ta(i+N,j+N) = -ta(j+N,i+N) = i / 1/ 2 sqrt(2)
// with i<j and i=0,...,N-2
// follows that j=i+1, ... , N-1
int i1, i2;
ta = Zero();
cplx i(0.0, 1.0);
RealD nrm = 1 / (2 * std::sqrt(2));
su2SubGroupIndex(i1, i2, bIndex);
ta()()(i1,i2) = i;
ta()()(i2,i1) = -i;
ta()()(i1+nsp,i2+nsp) = i;
ta()()(i2+nsp,i1+nsp) = -i;
ta = ta * nrm;
}
template <class cplx, ONLY_IF_Sp>
static void generatorCtype(int cIndex, iSp2nMatrix<cplx> &ta) {
// ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2)
int i1, i2;
ta = Zero();
RealD nrm = 1 / (2 * std::sqrt(2) );
su2SubGroupIndex(i1, i2, cIndex);
ta()()(i1,i2+nsp) = 1;
ta()()(i2,i1+nsp) = 1;
ta()()(i1+nsp,i2) = 1;
ta()()(i2+nsp,i1) = 1;
ta = ta * nrm;
}
template <class cplx, ONLY_IF_Sp>
static void generatorDtype(int dIndex, iSp2nMatrix<cplx> &ta) {
// ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2)
int i1, i2;
ta = Zero();
cplx i(0.0, 1.0);
RealD nrm = 1 / (2 * std::sqrt(2) );
su2SubGroupIndex(i1, i2, dIndex);
ta()()(i1,i2+nsp) = i;
ta()()(i2,i1+nsp) = i;
ta()()(i1+nsp,i2) = -i;
ta()()(i2+nsp,i1) = -i;
ta = ta * nrm;
}
template <class cplx, ONLY_IF_Sp>
static void generatorWtype(int wIndex, iSp2nMatrix<cplx> &ta) {
// ta(i,i+N) = ta(i+N,i) = 1/2
ta = Zero();
RealD nrm = 1.0 / 2; //check
ta()()(wIndex,wIndex+nsp) = 1;
ta()()(wIndex+nsp,wIndex) = 1;
ta = ta * nrm;
}
template <class cplx, ONLY_IF_Sp>
static void generatorZtype(int zIndex, iSp2nMatrix<cplx> &ta) {
// ta(i,i+N) = - ta(i+N,i) = i/2
ta = Zero();
RealD nrm = 1.0 / 2; //check
cplx i(0.0, 1.0);
ta()()(zIndex,zIndex+nsp) = i;
ta()()(zIndex+nsp,zIndex) = -i;
ta = ta * nrm;
}
////////////////////////////////////////////////////////////////////////
// Map a su2 subgroup number to the pair of rows that are non zero
////////////////////////////////////////////////////////////////////////
template<ONLY_IF_Sp>
static void su2SubGroupIndex(int &i1, int &i2, int su2_index) {
assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2));
int spare = su2_index;
for (i1 = 0; spare >= (nsp - 1 - i1); i1++) {
spare = spare - (nsp - 1 - i1); // remove the Nc-1-i1 terms
}
i2 = i1 + 1 + spare;
}
template<ONLY_IF_Sp>
static void testGenerators(void) {
Matrix ta;
Matrix tb;
std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is 0.5 delta_ab " << std::endl;
for (int a = 0; a < AlgebraDimension; a++) {
for (int b = 0; b < AlgebraDimension; b++) {
generator(a,ta);
generator(b,tb);
Complex tr = TensorRemove(trace( ta * tb) );
std::cout << GridLogMessage << "(" << a << "," << b << ") = " << tr
<< std::endl;
if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6);
if (a != b) assert(abs(tr) < 1.0e-6);
}
}
std::cout << GridLogMessage << std::endl;
std::cout << GridLogMessage << "Fundamental - Checking if hermitian" << std::endl;
for (int a = 0; a < AlgebraDimension; a++) {
generator(a,ta);
std::cout << GridLogMessage << a << std::endl;
assert(norm2(ta - adj(ta)) < 1.0e-6);
}
std::cout << GridLogMessage << std::endl;
std::cout << GridLogMessage << "Fundamental - Checking if traceless" << std::endl;
for (int a = 0; a < AlgebraDimension; a++) {
generator(a, ta);
Complex tr = TensorRemove(trace(ta));
std::cout << GridLogMessage << a << std::endl;
assert(abs(tr) < 1.0e-6);
}
}