portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-11 22:50:45 +01:00
Code Releases Activity

With function overloading (still dirty).

Browse Source
This commit is contained in:
Julian Lenz 2022-11-30 12:52:33 +00:00
parent 7c2ad4f8c8
commit 2507606bd0
4 changed files with 563 additions and 451 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
Grid/qcd/utils/SUn.h
View File

@ -32,15 +32,17 @@ directory
#ifndef QCD_UTIL_SUN_H #ifndef QCD_UTIL_SUN_H
#define QCD_UTIL_SUN_H #define QCD_UTIL_SUN_H
#define ONLY_IF_SU \ #define ONLY_IF_SU \
typename \ typename dummy_name = group_name, \
typename std::enable_if_t<is_su<group_name>::value, dummy_name> ,\ typename = std::enable_if_t < \
dummy_name = group_name, std::is_same<dummy_name, group_name>::value && \
is_su<dummy_name>::value >
#define ONLY_IF_Sp \ #define ONLY_IF_Sp \
typename \ typename dummy_name = group_name, \
typename std::enable_if_t<is_sp<group_name>::value, dummy_name> ,\ typename = std::enable_if_t < \
dummy_name =group_name, std::is_same<dummy_name, group_name>::value && \
is_sp<dummy_name>::value >
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
namespace GroupName { namespace GroupName {
@ -68,15 +70,18 @@ struct is_sp<GroupName::Sp> {
static const bool value = true; static const bool value = true;
}; };
template<typename group_name> template <typename group_name>
constexpr int compute_adjoint_dimension(int ncolour); constexpr int compute_adjoint_dimension(int ncolour);
template<> template <>
constexpr int compute_adjoint_dimension<GroupName::SU>(int ncolour) { return ncolour * ncolour - 1;} 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 <>
constexpr int compute_adjoint_dimension<GroupName::Sp>(int ncolour) {
return ncolour / 2 * (ncolour + 1);
}
template <int ncolour, class group_name = GroupName::SU> template <int ncolour, class group_name = GroupName::SU>
class GaugeGroup; class GaugeGroup;
@ -91,11 +96,12 @@ template <int ncolour, class group_name>
class GaugeGroup { class GaugeGroup {
public: public:
static const int Dimension = ncolour; static const int Dimension = ncolour;
static const int AdjointDimension = compute_adjoint_dimension<group_name>(ncolour); static const int AdjointDimension =
static const int AlgebraDimension = compute_adjoint_dimension<group_name>(ncolour); 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: // Don't know how to only enable this for Sp:
static const int nsp = ncolour/2; static const int nsp = ncolour / 2;
template <typename vtype> template <typename vtype>
using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >; using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >;
@ -103,6 +109,15 @@ class GaugeGroup {
using iGroupMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >; using iGroupMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >;
template <typename vtype> template <typename vtype>
using iAlgebraVector = iScalar<iScalar<iVector<vtype, AdjointDimension> > >; using iAlgebraVector = iScalar<iScalar<iVector<vtype, AdjointDimension> > >;
static int su2subgroups(void) { return su2subgroups(group_name()); }
template <typename vtype>
using iSUnMatrix = iGroupMatrix<vtype>;
template <typename vtype>
using iSUnAlgebraVector = iAlgebraVector<vtype>;
template <typename vtype>
using iSp2nMatrix = iGroupMatrix<vtype>;
template <typename vtype>
using iSp2nAlgebraVector = iAlgebraVector<vtype>;
////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////
// Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix, // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix,
@ -150,6 +165,17 @@ class GaugeGroup {
#include "Grid/qcd/utils/SUn_impl.h" #include "Grid/qcd/utils/SUn_impl.h"
#include "Grid/qcd/utils/Sp2n_impl.h" #include "Grid/qcd/utils/Sp2n_impl.h"
template <class cplx>
static void generator(int lieIndex, iSp2nMatrix<cplx> &ta) {
return generator(lieIndex, ta, group_name());
}
static void su2SubGroupIndex(int &i1, int &i2, int su2_index) {
return su2SubGroupIndex(i1, i2, su2_index, group_name());
}
static void testGenerators(void) { testGenerators(group_name()); }
static void printGenerators(void) { static void printGenerators(void) {
for (int gen = 0; gen < AdjointDimension; gen++) { for (int gen = 0; gen < AdjointDimension; gen++) {
Matrix ta; Matrix ta;

16
Grid/qcd/utils/SUn_impl.h
View File

@ -6,12 +6,7 @@
// //
// around it. // around it.
template<ONLY_IF_SU> static int su2subgroups(GroupName::SU) { return (ncolour * (ncolour - 1)) / 2; }
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>
using iSUnAlgebraVector = iAlgebraVector<vtype>;
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
// There are N^2-1 generators for SU(N). // There are N^2-1 generators for SU(N).
// //
@ -60,8 +55,8 @@ using iSUnAlgebraVector = iAlgebraVector<vtype>;
// ( -2) // ( -2)
// //
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
template <class cplx, ONLY_IF_SU> template <class cplx>
static void generator(int lieIndex, iSUnMatrix<cplx> &ta) { static void generator(int lieIndex, iSUnMatrix<cplx> &ta, GroupName::SU) {
// map lie index to which type of generator // map lie index to which type of generator
int diagIndex; int diagIndex;
int su2Index; int su2Index;
@ -117,8 +112,7 @@ static void generatorDiagonal(int diagIndex, iSUnMatrix<cplx> &ta) {
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
// Map a su2 subgroup number to the pair of rows that are non zero // Map a su2 subgroup number to the pair of rows that are non zero
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
template <ONLY_IF_SU> static void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::SU) {
static void su2SubGroupIndex(int &i1, int &i2, int su2_index) {
assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2)); assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2));
int spare = su2_index; int spare = su2_index;
@ -471,7 +465,7 @@ static void SubGroupHeatBath(
} }
template <ONLY_IF_SU> template <ONLY_IF_SU>
static void testGenerators(void) { static void testGenerators(GroupName::SU) {
Matrix ta; Matrix ta;
Matrix tb; Matrix tb;
std::cout << GridLogMessage std::cout << GridLogMessage

400
Grid/qcd/utils/Sp2n.h
View File

@ -1,14 +1,14 @@
#ifndef QCD_UTIL_Sp2n_H
#define QCD_UTIL_Sp2n_H
#include "Grid/qcd/utils/SUn.h" #include "Grid/qcd/utils/SUn.h"
// NAMESPACE_BEGIN(Grid);
// #ifndef QCD_UTIL_Sp2n_H //
// #define QCD_UTIL_Sp2n_H
//
// NAMESPACE_BEGIN(Grid);
//
// // Sp(2N) // // Sp(2N)
// // ncolour = 2N // // ncolour = 2N
// //
// //
// template <int ncolour> // template <int ncolour>
// class Sp { // class Sp {
// public: // public:
@ -16,41 +16,42 @@
// static const int Dimension = ncolour; // static const int Dimension = ncolour;
// static const int AlgebraDimension = nsp*(2*nsp +1); // static const int AlgebraDimension = nsp*(2*nsp +1);
// static int su2subgroups(void) { return (nsp * (nsp - 1)) / 2; } // static int su2subgroups(void) { return (nsp * (nsp - 1)) / 2; }
// //
// //
// template <typename vtype> // template <typename vtype>
// using iSp2nMatrix = iScalar<iScalar<iMatrix<vtype, Dimension> > >; // using iSp2nMatrix = iScalar<iScalar<iMatrix<vtype, Dimension> > >;
// template <typename vtype> // template <typename vtype>
// using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >; // using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >;
// template <typename vtype> // template <typename vtype>
// using iSp2nAlgebraVector = iScalar<iScalar<iVector<vtype, AlgebraDimension> > >; // using iSp2nAlgebraVector = iScalar<iScalar<iVector<vtype,
// // AlgebraDimension> > >;
//
// typedef iSp2nMatrix<Complex> Matrix; // typedef iSp2nMatrix<Complex> Matrix;
// typedef iSp2nMatrix<ComplexF> MatrixF; // typedef iSp2nMatrix<ComplexF> MatrixF;
// typedef iSp2nMatrix<ComplexD> MatrixD; // typedef iSp2nMatrix<ComplexD> MatrixD;
// //
// typedef iSp2nMatrix<vComplex> vMatrix; // typedef iSp2nMatrix<vComplex> vMatrix;
// typedef iSp2nMatrix<vComplexF> vMatrixF; // typedef iSp2nMatrix<vComplexF> vMatrixF;
// typedef iSp2nMatrix<vComplexD> vMatrixD; // typedef iSp2nMatrix<vComplexD> vMatrixD;
// //
// typedef iSp2nAlgebraVector<Complex> AlgebraVector; // typedef iSp2nAlgebraVector<Complex> AlgebraVector;
// typedef iSp2nAlgebraVector<ComplexF> AlgebraVectorF; // typedef iSp2nAlgebraVector<ComplexF> AlgebraVectorF;
// typedef iSp2nAlgebraVector<ComplexD> AlgebraVectorD; // typedef iSp2nAlgebraVector<ComplexD> AlgebraVectorD;
// //
// typedef iSp2nAlgebraVector<vComplex> vAlgebraVector; // typedef iSp2nAlgebraVector<vComplex> vAlgebraVector;
// typedef iSp2nAlgebraVector<vComplexF> vAlgebraVectorF; // typedef iSp2nAlgebraVector<vComplexF> vAlgebraVectorF;
// typedef iSp2nAlgebraVector<vComplexD> vAlgebraVectorD; // typedef iSp2nAlgebraVector<vComplexD> vAlgebraVectorD;
// //
// typedef Lattice<vMatrix> LatticeMatrix; // typedef Lattice<vMatrix> LatticeMatrix;
// typedef Lattice<vMatrixF> LatticeMatrixF; // typedef Lattice<vMatrixF> LatticeMatrixF;
// typedef Lattice<vMatrixD> LatticeMatrixD; // typedef Lattice<vMatrixD> LatticeMatrixD;
// //
// typedef Lattice<vAlgebraVector> LatticeAlgebraVector; // typedef Lattice<vAlgebraVector> LatticeAlgebraVector;
// typedef Lattice<vAlgebraVectorF> LatticeAlgebraVectorF; // typedef Lattice<vAlgebraVectorF> LatticeAlgebraVectorF;
// typedef Lattice<vAlgebraVectorD> LatticeAlgebraVectorD; // typedef Lattice<vAlgebraVectorD> LatticeAlgebraVectorD;
// //
// //
// //
// // Sp(2N) has N(2N+1) = 2N^2+N generators // // Sp(2N) has N(2N+1) = 2N^2+N generators
// // // //
// // normalise the generators such that // // normalise the generators such that
@ -60,19 +61,19 @@
// // off diagonal: // // off diagonal:
// // there are 6 types named a,b,c,d and w,z // // there are 6 types named a,b,c,d and w,z
// // abcd are N(N-1)/2 each while wz are N each // // abcd are N(N-1)/2 each while wz are N each
// //
// template <class cplx> // template <class cplx>
// static void generator(int lieIndex, iSp2nMatrix<cplx> &ta) { // static void generator(int lieIndex, iSp2nMatrix<cplx> &ta) {
// // map lie index into type of generators: diagonal, abcd type, wz type // // map lie index into type of generators: diagonal, abcd type, wz
// // type
//
// int diagIndex; // int diagIndex;
// int aIndex, bIndex, cIndex, dIndex; // int aIndex, bIndex, cIndex, dIndex;
// int wIndex, zIndex; // a,b,c,d are N(N-1)/2 and w,z are N // 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 mod = nsp * (nsp-1) * 0.5;
// int offdiag = 2*nsp*nsp; // number of generators not in the cartan subalgebra // int offdiag = 2*nsp*nsp; // number of generators not in the cartan
// int wmod = 4*mod; // subalgebra int wmod = 4*mod; int zmod = wmod+nsp; if (lieIndex >=
// int zmod = wmod+nsp; // offdiag) {
// if (lieIndex >= offdiag) {
// diagIndex = lieIndex - offdiag; // 0, ... ,N-1 // diagIndex = lieIndex - offdiag; // 0, ... ,N-1
// //std::cout << GridLogMessage << "diag type " << std::endl; // //std::cout << GridLogMessage << "diag type " << std::endl;
// generatorDiagtype(diagIndex, ta); // generatorDiagtype(diagIndex, ta);
@ -86,7 +87,8 @@
// } // }
// if ( (lieIndex >= zmod) && (lieIndex < offdiag) ) { // if ( (lieIndex >= zmod) && (lieIndex < offdiag) ) {
// //std::cout << GridLogMessage << "z type " << std::endl; // //std::cout << GridLogMessage << "z type " << std::endl;
// //std::cout << GridLogMessage << "lie index " << lieIndex << std::endl; // //std::cout << GridLogMessage << "lie index " << lieIndex <<
// std::endl;
// //std::cout << GridLogMessage << "z mod " << zmod << std::endl; // //std::cout << GridLogMessage << "z mod " << zmod << std::endl;
// zIndex = lieIndex - zmod; // 0, ... ,N-1 // zIndex = lieIndex - zmod; // 0, ... ,N-1
// generatorZtype(zIndex,ta); // generatorZtype(zIndex,ta);
@ -95,168 +97,171 @@
// if (lieIndex < mod) { // atype 0, ... , N(N-1)/2=mod // if (lieIndex < mod) { // atype 0, ... , N(N-1)/2=mod
// //std::cout << GridLogMessage << "a type " << std::endl; // //std::cout << GridLogMessage << "a type " << std::endl;
// aIndex = lieIndex; // aIndex = lieIndex;
// //std::cout << GridLogMessage << "a indx " << aIndex << std::endl; // //std::cout << GridLogMessage << "a indx " << aIndex <<
// generatorAtype(aIndex, ta); // std::endl; generatorAtype(aIndex, ta); return;
// return;
// } // }
// if ( (lieIndex >= mod) && lieIndex < 2*mod) { // btype mod, ... , 2mod-1 // if ( (lieIndex >= mod) && lieIndex < 2*mod) { // btype mod, ... ,
// 2mod-1
// //std::cout << GridLogMessage << "b type " << std::endl; // //std::cout << GridLogMessage << "b type " << std::endl;
// bIndex = lieIndex - mod; // bIndex = lieIndex - mod;
// generatorBtype(bIndex, ta); // generatorBtype(bIndex, ta);
// return; // return;
// } // }
// if ( (lieIndex >= 2*mod) && lieIndex < 3*mod) { // ctype 2mod, ... , 3mod-1 // if ( (lieIndex >= 2*mod) && lieIndex < 3*mod) { // ctype 2mod, ... ,
// 3mod-1
// //std::cout << GridLogMessage << "c type " << std::endl; // //std::cout << GridLogMessage << "c type " << std::endl;
// cIndex = lieIndex - 2*mod; // cIndex = lieIndex - 2*mod;
// generatorCtype(cIndex, ta); // generatorCtype(cIndex, ta);
// return; // return;
// } // }
// if ( (lieIndex >= 3*mod) && lieIndex < wmod) { // ctype 3mod, ... , 4mod-1 = wmod-1 // if ( (lieIndex >= 3*mod) && lieIndex < wmod) { // ctype 3mod, ... ,
// 4mod-1 = wmod-1
// //std::cout << GridLogMessage << "d type " << std::endl; // //std::cout << GridLogMessage << "d type " << std::endl;
// dIndex = lieIndex - 3*mod; // dIndex = lieIndex - 3*mod;
// generatorDtype(dIndex, ta); // generatorDtype(dIndex, ta);
// return; // return;
// } // }
// //
// } //end of generator // } //end of generator
// //
// template <class cplx> // template <class cplx>
// static void generatorDiagtype(int diagIndex, iSp2nMatrix<cplx> &ta) { // 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(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan
// // subalgebra
//
// ta = Zero(); // ta = Zero();
// RealD nrm = 1.0 / 2; // RealD nrm = 1.0 / 2;
// //
// ta()()(diagIndex,diagIndex) = nrm; // ta()()(diagIndex,diagIndex) = nrm;
// ta()()(diagIndex+nsp,diagIndex+nsp) = -nrm; // ta()()(diagIndex+nsp,diagIndex+nsp) = -nrm;
// } // }
// //
// template <class cplx> // template <class cplx>
// static void generatorAtype(int aIndex, iSp2nMatrix<cplx> &ta) { // 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) // // 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 // // with i<j and i=0,...,N-2
// // follows that j=i+1, ... , N // // follows that j=i+1, ... , N
// int i1, i2; // int i1, i2;
// ta = Zero(); // ta = Zero();
// RealD nrm = 1 / (2 * std::sqrt(2) ); // RealD nrm = 1 / (2 * std::sqrt(2) );
// //
// su2SubGroupIndex(i1, i2, aIndex); // su2SubGroupIndex(i1, i2, aIndex);
// ta()()(i1,i2) = 1; // ta()()(i1,i2) = 1;
// ta()()(i2,i1) = 1; // ta()()(i2,i1) = 1;
// ta()()(i1+nsp,i2+nsp) = -1; // ta()()(i1+nsp,i2+nsp) = -1;
// ta()()(i2+nsp,i1+nsp) = -1; // ta()()(i2+nsp,i1+nsp) = -1;
// //
// ta = ta * nrm; // ta = ta * nrm;
// } // }
// //
// template <class cplx> // template <class cplx>
// static void generatorBtype(int bIndex, iSp2nMatrix<cplx> &ta) { // 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) // // 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 // // with i<j and i=0,...,N-2
// // follows that j=i+1, ... , N-1 // // follows that j=i+1, ... , N-1
// //
// int i1, i2; // int i1, i2;
// ta = Zero(); // ta = Zero();
// cplx i(0.0, 1.0); // cplx i(0.0, 1.0);
// RealD nrm = 1 / (2 * std::sqrt(2)); // RealD nrm = 1 / (2 * std::sqrt(2));
// su2SubGroupIndex(i1, i2, bIndex); // su2SubGroupIndex(i1, i2, bIndex);
// //
// //
// ta()()(i1,i2) = i; // ta()()(i1,i2) = i;
// ta()()(i2,i1) = -i; // ta()()(i2,i1) = -i;
// ta()()(i1+nsp,i2+nsp) = i; // ta()()(i1+nsp,i2+nsp) = i;
// ta()()(i2+nsp,i1+nsp) = -i; // ta()()(i2+nsp,i1+nsp) = -i;
// //
// ta = ta * nrm; // ta = ta * nrm;
// } // }
// //
// template <class cplx> // template <class cplx>
// static void generatorCtype(int cIndex, iSp2nMatrix<cplx> &ta) { // 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) // // ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2)
// //
// //
// int i1, i2; // int i1, i2;
// ta = Zero(); // ta = Zero();
// RealD nrm = 1 / (2 * std::sqrt(2) ); // RealD nrm = 1 / (2 * std::sqrt(2) );
// su2SubGroupIndex(i1, i2, cIndex); // su2SubGroupIndex(i1, i2, cIndex);
// //
// ta()()(i1,i2+nsp) = 1; // ta()()(i1,i2+nsp) = 1;
// ta()()(i2,i1+nsp) = 1; // ta()()(i2,i1+nsp) = 1;
// ta()()(i1+nsp,i2) = 1; // ta()()(i1+nsp,i2) = 1;
// ta()()(i2+nsp,i1) = 1; // ta()()(i2+nsp,i1) = 1;
// //
// ta = ta * nrm; // ta = ta * nrm;
// } // }
// //
// template <class cplx> // template <class cplx>
// static void generatorDtype(int dIndex, iSp2nMatrix<cplx> &ta) { // 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) // // ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2)
// //
// int i1, i2; // int i1, i2;
// ta = Zero(); // ta = Zero();
// cplx i(0.0, 1.0); // cplx i(0.0, 1.0);
// RealD nrm = 1 / (2 * std::sqrt(2) ); // RealD nrm = 1 / (2 * std::sqrt(2) );
// su2SubGroupIndex(i1, i2, dIndex); // su2SubGroupIndex(i1, i2, dIndex);
// //
// ta()()(i1,i2+nsp) = i; // ta()()(i1,i2+nsp) = i;
// ta()()(i2,i1+nsp) = i; // ta()()(i2,i1+nsp) = i;
// ta()()(i1+nsp,i2) = -i; // ta()()(i1+nsp,i2) = -i;
// ta()()(i2+nsp,i1) = -i; // ta()()(i2+nsp,i1) = -i;
// //
// ta = ta * nrm; // ta = ta * nrm;
// } // }
// //
// template <class cplx> // template <class cplx>
// static void generatorWtype(int wIndex, iSp2nMatrix<cplx> &ta) { // static void generatorWtype(int wIndex, iSp2nMatrix<cplx> &ta) {
// //
// // ta(i,i+N) = ta(i+N,i) = 1/2 // // ta(i,i+N) = ta(i+N,i) = 1/2
// //
// ta = Zero(); // ta = Zero();
// RealD nrm = 1.0 / 2; //check // RealD nrm = 1.0 / 2; //check
// //
// ta()()(wIndex,wIndex+nsp) = 1; // ta()()(wIndex,wIndex+nsp) = 1;
// ta()()(wIndex+nsp,wIndex) = 1; // ta()()(wIndex+nsp,wIndex) = 1;
// //
// ta = ta * nrm; // ta = ta * nrm;
// } // }
// //
// template <class cplx> // template <class cplx>
// static void generatorZtype(int zIndex, iSp2nMatrix<cplx> &ta) { // static void generatorZtype(int zIndex, iSp2nMatrix<cplx> &ta) {
// //
// // ta(i,i+N) = - ta(i+N,i) = i/2 // // ta(i,i+N) = - ta(i+N,i) = i/2
// //
// ta = Zero(); // ta = Zero();
// RealD nrm = 1.0 / 2; //check // RealD nrm = 1.0 / 2; //check
// cplx i(0.0, 1.0); // cplx i(0.0, 1.0);
// ta()()(zIndex,zIndex+nsp) = i; // ta()()(zIndex,zIndex+nsp) = i;
// ta()()(zIndex+nsp,zIndex) = -i; // ta()()(zIndex+nsp,zIndex) = -i;
// //
// ta = ta * nrm; // ta = ta * nrm;
// } // }
// //
// //
// //////////////////////////////////////////////////////////////////////// // ////////////////////////////////////////////////////////////////////////
// // Map a su2 subgroup number to the pair of rows that are non zero // // Map a su2 subgroup number to the pair of rows that are non zero
// //////////////////////////////////////////////////////////////////////// // ////////////////////////////////////////////////////////////////////////
// static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { // static void su2SubGroupIndex(int &i1, int &i2, int su2_index) {
// assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2)); // assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2));
// //
// int spare = su2_index; // int spare = su2_index;
// for (i1 = 0; spare >= (nsp - 1 - i1); i1++) { // for (i1 = 0; spare >= (nsp - 1 - i1); i1++) {
// spare = spare - (nsp - 1 - i1); // remove the Nc-1-i1 terms // spare = spare - (nsp - 1 - i1); // remove the Nc-1-i1 terms
// } // }
// i2 = i1 + 1 + spare; // i2 = i1 + 1 + spare;
// } // }
// //
// //
// //
// //
// //
// static void printGenerators(void) { // static void printGenerators(void) {
// for (int gen = 0; gen < AlgebraDimension; gen++) { // for (int gen = 0; gen < AlgebraDimension; gen++) {
// Matrix ta; // Matrix ta;
@ -266,84 +271,89 @@
// std::cout << GridLogMessage << ta << std::endl; // std::cout << GridLogMessage << ta << std::endl;
// } // }
// } // }
// //
// //
// //
// static void testGenerators(void) { // static void testGenerators(void) {
// Matrix ta; // Matrix ta;
// Matrix tb; // Matrix tb;
// std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is 0.5 delta_ab " << std::endl; // std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is
// for (int a = 0; a < AlgebraDimension; a++) { // 0.5 delta_ab " << std::endl; for (int a = 0; a < AlgebraDimension;
// a++) {
// for (int b = 0; b < AlgebraDimension; b++) { // for (int b = 0; b < AlgebraDimension; b++) {
// generator(a,ta); // generator(a,ta);
// generator(b,tb); // generator(b,tb);
// Complex tr = TensorRemove(trace( ta * tb) ); // Complex tr = TensorRemove(trace( ta * tb) );
// std::cout << GridLogMessage << "(" << a << "," << b << ") = " << tr // std::cout << GridLogMessage << "(" << a << "," << b << ") =
// " << tr
// << std::endl; // << std::endl;
// if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6); // if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6);
// if (a != b) assert(abs(tr) < 1.0e-6); // if (a != b) assert(abs(tr) < 1.0e-6);
// //
// } // }
// } // }
// std::cout << GridLogMessage << std::endl; // std::cout << GridLogMessage << std::endl;
// std::cout << GridLogMessage << "Fundamental - Checking if hermitian" << std::endl; // std::cout << GridLogMessage << "Fundamental - Checking if hermitian"
// for (int a = 0; a < AlgebraDimension; a++) { // << std::endl; for (int a = 0; a < AlgebraDimension; a++) {
// generator(a,ta); // generator(a,ta);
// std::cout << GridLogMessage << a << std::endl; // std::cout << GridLogMessage << a << std::endl;
// assert(norm2(ta - adj(ta)) < 1.0e-6); // assert(norm2(ta - adj(ta)) < 1.0e-6);
// } // }
// std::cout << GridLogMessage << std::endl; // std::cout << GridLogMessage << std::endl;
// std::cout << GridLogMessage << "Fundamental - Checking if traceless" << std::endl; // std::cout << GridLogMessage << "Fundamental - Checking if traceless"
// for (int a = 0; a < AlgebraDimension; a++) { // << std::endl; for (int a = 0; a < AlgebraDimension; a++) {
// generator(a, ta); // generator(a, ta);
// Complex tr = TensorRemove(trace(ta)); // Complex tr = TensorRemove(trace(ta));
// std::cout << GridLogMessage << a << std::endl; // std::cout << GridLogMessage << a << std::endl;
// assert(abs(tr) < 1.0e-6); // assert(abs(tr) < 1.0e-6);
// } // }
// //
// } // }
// //
// //
// template <typename LatticeMatrixType> // template <typename LatticeMatrixType>
// static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, double scale = 1.0) // static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out,
// double scale = 1.0)
// { // {
// GridBase *grid = out.Grid(); // GridBase *grid = out.Grid();
// //
// typedef typename LatticeMatrixType::vector_type vector_type; // typedef typename LatticeMatrixType::vector_type vector_type;
// typedef typename LatticeMatrixType::scalar_type scalar_type; // typedef typename LatticeMatrixType::scalar_type scalar_type;
// //
// typedef iSinglet<vector_type> vTComplexType; // typedef iSinglet<vector_type> vTComplexType;
// //
// typedef Lattice<vTComplexType> LatticeComplexType; // typedef Lattice<vTComplexType> LatticeComplexType;
// typedef typename GridTypeMapper<typename LatticeMatrixType::vector_object>::scalar_object MatrixType; // typedef typename GridTypeMapper<typename
// // LatticeMatrixType::vector_object>::scalar_object MatrixType;
//
// LatticeComplexType ca(grid); // LatticeComplexType ca(grid);
// LatticeMatrixType lie(grid); // LatticeMatrixType lie(grid);
// LatticeMatrixType la(grid); // LatticeMatrixType la(grid);
// ComplexD ci(0.0, scale); // ComplexD ci(0.0, scale);
// // ComplexD cone(1.0, 0.0); // // ComplexD cone(1.0, 0.0);
// MatrixType ta; // MatrixType ta;
// //
// lie = Zero(); // lie = Zero();
// //
// for (int a = 0; a < AlgebraDimension; a++) { // for (int a = 0; a < AlgebraDimension; a++) {
// random(pRNG, ca); // random(pRNG, ca);
// //
// ca = (ca + conjugate(ca)) * 0.5; // ca = (ca + conjugate(ca)) * 0.5;
// ca = ca - 0.5; // ca = ca - 0.5;
// //
// generator(a, ta); // generator(a, ta);
// //
// la = ci * ca * ta; // la = ci * ca * ta;
// //
// lie = lie + la; // e^{i la ta} // lie = lie + la; // e^{i la ta}
// //
// } // }
// taExp(lie, out); // taExp(lie, out);
// } // }
// //
// //
// static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG, //same as sun // static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG,
// //same as sun
// LatticeMatrix &out, // LatticeMatrix &out,
// Real scale = 1.0) { // Real scale = 1.0) {
// GridBase *grid = out.Grid(); // GridBase *grid = out.Grid();
@ -351,7 +361,7 @@
// LatticeMatrix la(grid); // LatticeMatrix la(grid);
// Complex ci(0.0, scale); // Complex ci(0.0, scale);
// Matrix ta; // Matrix ta;
// //
// out = Zero(); // out = Zero();
// for (int a = 0; a < AlgebraDimension; a++) { // for (int a = 0; a < AlgebraDimension; a++) {
// gaussian(pRNG, ca); // gaussian(pRNG, ca);
@ -361,7 +371,7 @@
// } // }
// out *= ci; // out *= ci;
// } // }
// //
// static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h, // static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h,
// LatticeMatrix &out, // LatticeMatrix &out,
// Real scale = 1.0) { // Real scale = 1.0) {
@ -369,7 +379,7 @@
// GridBase *grid = out.Grid(); // GridBase *grid = out.Grid();
// LatticeMatrix la(grid); // LatticeMatrix la(grid);
// Matrix ta; // Matrix ta;
// //
// out = Zero(); // out = Zero();
// for (int a = 0; a < AlgebraDimension; a++) { // for (int a = 0; a < AlgebraDimension; a++) {
// generator(a, ta); // generator(a, ta);
@ -377,18 +387,19 @@
// out += la; // out += la;
// } // }
// } // }
// //
// //
// template <typename LatticeMatrixType> // template <typename LatticeMatrixType>
// static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { // same as sun // static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { //
// same as sun
// typedef typename LatticeMatrixType::scalar_type ComplexType; // typedef typename LatticeMatrixType::scalar_type ComplexType;
// //
// LatticeMatrixType xn(x.Grid()); // LatticeMatrixType xn(x.Grid());
// RealD nfac = 1.0; // RealD nfac = 1.0;
// //
// xn = x; // xn = x;
// ex = xn + ComplexType(1.0); // 1+x // ex = xn + ComplexType(1.0); // 1+x
// //
// // Do a 12th order exponentiation // // Do a 12th order exponentiation
// for (int i = 2; i <= 12; ++i) { // for (int i = 2; i <= 12; ++i) {
// nfac = nfac / RealD(i); // 1/2, 1/2.3 ... // nfac = nfac / RealD(i); // 1/2, 1/2.3 ...
@ -396,25 +407,26 @@
// ex = ex + xn * nfac; // x2/2!, x3/3!.... // ex = ex + xn * nfac; // x2/2!, x3/3!....
// } // }
// } // }
// //
// static void projectOnAlgebra(LatticeAlgebraVector &h_out, const LatticeMatrix &in, Real scale = 1.0) { // static void projectOnAlgebra(LatticeAlgebraVector &h_out, const
// LatticeMatrix &in, Real scale = 1.0) {
// conformable(h_out, in); // conformable(h_out, in);
// h_out = Zero(); // h_out = Zero();
// Matrix Ta; // Matrix Ta;
// //
// for (int a = 0; a < AlgebraDimension; a++) { // for (int a = 0; a < AlgebraDimension; a++) {
// generator(a, Ta); // generator(a, Ta);
// pokeColour(h_out, - 2.0 * (trace(timesI(Ta) * in)) * scale, a); // pokeColour(h_out, - 2.0 * (trace(timesI(Ta) * in)) * scale, a);
// } // }
// } // }
// //
// //
// template <typename GaugeField> // template <typename GaugeField>
// static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) { // static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
// typedef typename GaugeField::vector_type vector_type; // typedef typename GaugeField::vector_type vector_type;
// typedef iSp2nMatrix<vector_type> vMatrixType; // typedef iSp2nMatrix<vector_type> vMatrixType;
// typedef Lattice<vMatrixType> LatticeMatrixType; // typedef Lattice<vMatrixType> LatticeMatrixType;
// //
// LatticeMatrixType Umu(out.Grid()); // LatticeMatrixType Umu(out.Grid());
// for (int mu = 0; mu < Nd; mu++) { // for (int mu = 0; mu < Nd; mu++) {
// LieRandomize(pRNG, Umu, 1.0); //def // LieRandomize(pRNG, Umu, 1.0); //def
@ -426,20 +438,20 @@
// typedef typename GaugeField::vector_type vector_type; // typedef typename GaugeField::vector_type vector_type;
// typedef iSp2nMatrix<vector_type> vMatrixType; // typedef iSp2nMatrix<vector_type> vMatrixType;
// typedef Lattice<vMatrixType> LatticeMatrixType; // typedef Lattice<vMatrixType> LatticeMatrixType;
// //
// LatticeMatrixType Umu(out.Grid()); // LatticeMatrixType Umu(out.Grid());
// for(int mu=0;mu<Nd;mu++){ // for(int mu=0;mu<Nd;mu++){
// LieRandomize(pRNG,Umu,0.01); //def // LieRandomize(pRNG,Umu,0.01); //def
// PokeIndex<LorentzIndex>(out,Umu,mu); // PokeIndex<LorentzIndex>(out,Umu,mu);
// } // }
// } // }
// //
// template<typename GaugeField> // template<typename GaugeField>
// static void ColdConfiguration(GaugeField &out){ // static void ColdConfiguration(GaugeField &out){
// typedef typename GaugeField::vector_type vector_type; // typedef typename GaugeField::vector_type vector_type;
// typedef iSp2nMatrix<vector_type> vMatrixType; // typedef iSp2nMatrix<vector_type> vMatrixType;
// typedef Lattice<vMatrixType> LatticeMatrixType; // typedef Lattice<vMatrixType> LatticeMatrixType;
// //
// LatticeMatrixType Umu(out.Grid()); // LatticeMatrixType Umu(out.Grid());
// Umu=1.0; // Umu=1.0;
// for(int mu=0;mu<Nd;mu++){ // for(int mu=0;mu<Nd;mu++){
@ -450,15 +462,15 @@
// static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){ // static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){
// ColdConfiguration(out); // ColdConfiguration(out);
// } // }
// //
// static void OmegaInvariance(ColourMatrix &in) // static void OmegaInvariance(ColourMatrix &in)
// { // {
// //
// ColourMatrix Omega; // ColourMatrix Omega;
// Omega = Zero(); // Omega = Zero();
// //
// std::cout << GridLogMessage << "I am a ColourMatrix" << std::endl; // std::cout << GridLogMessage << "I am a ColourMatrix" << std::endl;
// //
// //for (int i = 0; i < ncolour; i++) wrong?! // //for (int i = 0; i < ncolour; i++) wrong?!
// //{ // //{
// // Omega()()(i, 2*ncolour-1-i) = 1.; // // Omega()()(i, 2*ncolour-1-i) = 1.;
@ -469,113 +481,117 @@
// Omega()()(i, nsp+i) = 1.; // Omega()()(i, nsp+i) = 1.;
// Omega()()(nsp+i, i) = -1; // Omega()()(nsp+i, i) = -1;
// } // }
// //
// auto diff = Omega - (in * Omega * transpose(in) ); // auto diff = Omega - (in * Omega * transpose(in) );
// auto sdiff = norm2(diff); // auto sdiff = norm2(diff);
// if (norm2(sdiff) < 1e-8) // if (norm2(sdiff) < 1e-8)
// { // {
// std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl; // std::cout << GridLogMessage << "Symplectic condition satisfied:
// Omega invariant" << std::endl;
// } else { // } else {
// std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT left invariant by " << norm2(sdiff) << std::endl; // std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT
// left invariant by " << norm2(sdiff) << std::endl;
// } // }
// //
// } // }
// //
// template<typename GaugeField> // template<typename GaugeField>
// static void OmegaInvariance(GaugeField &in) // static void OmegaInvariance(GaugeField &in)
// { // {
// typedef typename GaugeField::vector_type vector_type; // typedef typename GaugeField::vector_type vector_type;
// typedef iSp2nMatrix<vector_type> vMatrixType; // typedef iSp2nMatrix<vector_type> vMatrixType;
// typedef Lattice<vMatrixType> LatticeMatrixType; // typedef Lattice<vMatrixType> LatticeMatrixType;
// //
// LatticeMatrixType U(in.Grid()); // LatticeMatrixType U(in.Grid());
// LatticeMatrixType OmegaLatt(in.Grid()); // LatticeMatrixType OmegaLatt(in.Grid());
// LatticeMatrixType identity(in.Grid()); // LatticeMatrixType identity(in.Grid());
// RealD vol = in.Grid()->gSites(); // RealD vol = in.Grid()->gSites();
// ColourMatrix Omega; // ColourMatrix Omega;
// //
// OmegaLatt = Zero(); // OmegaLatt = Zero();
// Omega = Zero(); // Omega = Zero();
// identity = 1.; // identity = 1.;
// //
// std::cout << GridLogMessage << "I am a GaugeField " << std::endl; // std::cout << GridLogMessage << "I am a GaugeField " << std::endl;
// //
// U = PeekIndex<LorentzIndex>(in,1); // U = PeekIndex<LorentzIndex>(in,1);
// //
// OmegaInvariance(U); // OmegaInvariance(U);
// //
// } // }
// //
// static void OmegaInvariance(LatticeColourMatrixD &in) // static void OmegaInvariance(LatticeColourMatrixD &in)
// { // {
// //
// LatticeColourMatrixD OmegaLatt(in.Grid()); // LatticeColourMatrixD OmegaLatt(in.Grid());
// LatticeColourMatrixD identity(in.Grid()); // LatticeColourMatrixD identity(in.Grid());
// RealD vol = in.Grid()->gSites(); // RealD vol = in.Grid()->gSites();
// ColourMatrix Omega; // ColourMatrix Omega;
// //
// OmegaLatt = Zero(); // OmegaLatt = Zero();
// Omega = Zero(); // Omega = Zero();
// identity = 1.; // identity = 1.;
// //
// std::cout << GridLogMessage << "I am a LatticeColourMatrix " << std::endl; // std::cout << GridLogMessage << "I am a LatticeColourMatrix " <<
// // std::endl;
//
// for (int i = 0; i < nsp; i++) // for (int i = 0; i < nsp; i++)
// { // {
// Omega()()(i, nsp+i) = 1.; // Omega()()(i, nsp+i) = 1.;
// Omega()()(nsp+i, i) = -1; // Omega()()(nsp+i, i) = -1;
// } // }
// //
// std::cout << GridLogMessage << "Omega = " << Omega()() << std::endl; // std::cout << GridLogMessage << "Omega = " << Omega()() <<
// OmegaLatt = OmegaLatt + (identity*Omega); // std::endl; OmegaLatt = OmegaLatt + (identity*Omega);
// //
// auto diff = OmegaLatt - (in*OmegaLatt*transpose(in)); // auto diff = OmegaLatt - (in*OmegaLatt*transpose(in));
// auto sdiff = sum(diff); // auto sdiff = sum(diff);
// //assert( norm2(sdiff) < 1e-8 ); // //assert( norm2(sdiff) < 1e-8 );
// if (norm2(sdiff) < 1e-8) // if (norm2(sdiff) < 1e-8)
// { // {
// std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl; // std::cout << GridLogMessage << "Symplectic condition
// satisfied: Omega invariant" << std::endl;
// } else { // } else {
// std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT left invariant by " << norm2(sdiff) << std::endl; // std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT
// left invariant by " << norm2(sdiff) << std::endl;
// } // }
// //
// //
// } // }
// //
// //
// }; // end of class Sp // }; // end of class Sp
//
// template <int N>
// template<int N> static void ProjectSp2n(
// static void ProjectSp2n(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
// { Umu = ProjectOnSpGroup(Umu);
// Umu = ProjectOnSpGroup(Umu); auto det = Determinant(Umu); // ok ?
// auto det = Determinant(Umu); // ok ?
// det = conjugate(det);
// det = conjugate(det);
// for (int i = 0; i < N; i++) {
// for(int i=0;i<N;i++){ auto element = PeekIndex<ColourIndex>(Umu, N - 1, i);
// auto element = PeekIndex<ColourIndex>(Umu,N-1,i); element = element * det;
// element = element * det; PokeIndex<ColourIndex>(Umu, element, Nc - 1, i);
// PokeIndex<ColourIndex>(Umu,element,Nc-1,i); }
// } }
// } template <int N>
// template<int N> static void ProjectSp2n(
// static void ProjectSp2n(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >,Nd> > &U) Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
// { GridBase *grid = U.Grid();
// GridBase *grid=U.Grid(); for (int mu = 0; mu < Nd; mu++) {
// for(int mu=0;mu<Nd;mu++){ auto Umu = PeekIndex<LorentzIndex>(U, mu);
// auto Umu = PeekIndex<LorentzIndex>(U,mu); Umu = ProjectOnSpGroup(Umu);
// Umu = ProjectOnSpGroup(Umu); ProjectSp2n(Umu);
// ProjectSp2n(Umu); PokeIndex<LorentzIndex>(U, Umu, mu);
// PokeIndex<LorentzIndex>(U,Umu,mu); }
// } }
// }
// typedef Sp<2> Sp2;
// typedef Sp<2> Sp2; typedef Sp<4> Sp4;
// typedef Sp<4> Sp4; typedef Sp<6> Sp6;
// typedef Sp<6> Sp6; typedef Sp<8> Sp8;
// typedef Sp<8> Sp8;
// NAMESPACE_END(Grid);
// NAMESPACE_END(Grid); #endif
// #endif

536
Grid/qcd/utils/Sp2n_impl.h
View File

@ -1,245 +1,321 @@
template<ONLY_IF_Sp> static int su2subgroups(GroupName::Sp) { return (nsp * (nsp - 1)) / 2; }
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 // Sp(2N) has N(2N+1) = 2N^2+N generators
// //
// normalise the generators such that // normalise the generators such that
// Trace ( Ta Tb) = 1/2 delta_ab // Trace ( Ta Tb) = 1/2 delta_ab
// //
// N generators in the cartan, 2N^2 off // N generators in the cartan, 2N^2 off
// off diagonal: // off diagonal:
// there are 6 types named a,b,c,d and w,z // there are 6 types named a,b,c,d and w,z
// abcd are N(N-1)/2 each while wz are N each // 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; template <class cplx, ONLY_IF_Sp>
int aIndex, bIndex, cIndex, dIndex; static void generator(int lieIndex, iSp2nMatrix<cplx> &ta, GroupName::Sp) {
int wIndex, zIndex; // a,b,c,d are N(N-1)/2 and w,z are N // map lie index into type of generators: diagonal, abcd type, wz type
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; int diagIndex;
ta()()(diagIndex+nsp,diagIndex+nsp) = -nrm; 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;
template <class cplx, ONLY_IF_Sp> int offdiag =
static void generatorAtype(int aIndex, iSp2nMatrix<cplx> &ta) { 2 * nsp * nsp; // number of generators not in the cartan subalgebra
int wmod = 4 * mod;
// ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2) int zmod = wmod + nsp;
// with i<j and i=0,...,N-2 if (lieIndex >= offdiag) {
// follows that j=i+1, ... , N diagIndex = lieIndex - offdiag; // 0, ... ,N-1
int i1, i2; // std::cout << GridLogMessage << "diag type " << std::endl;
ta = Zero(); generatorDiagtype(diagIndex, ta);
RealD nrm = 1 / (2 * std::sqrt(2) ); return;
}
su2SubGroupIndex(i1, i2, aIndex); if ((lieIndex >= wmod) && (lieIndex < zmod)) {
ta()()(i1,i2) = 1; // std::cout << GridLogMessage << "w type " << std::endl;
ta()()(i2,i1) = 1; wIndex = lieIndex - wmod; // 0, ... ,N-1
ta()()(i1+nsp,i2+nsp) = -1; generatorWtype(wIndex, ta);
ta()()(i2+nsp,i1+nsp) = -1; 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;
}
ta = ta * nrm; } // end of generator
}
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; template <class cplx, ONLY_IF_Sp>
ta()()(i2,i1) = -i; static void generatorDiagtype(int diagIndex, iSp2nMatrix<cplx> &ta) {
ta()()(i1+nsp,i2+nsp) = i; // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra
ta()()(i2+nsp,i1+nsp) = -i;
ta = ta * nrm; ta = Zero();
} RealD nrm = 1.0 / 2;
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; ta()()(diagIndex, diagIndex) = nrm;
} ta()()(diagIndex + nsp, diagIndex + nsp) = -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; template <class cplx, ONLY_IF_Sp>
ta()()(i2,i1+nsp) = i; static void generatorAtype(int aIndex, iSp2nMatrix<cplx> &ta) {
ta()()(i1+nsp,i2) = -i; // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2)
ta()()(i2+nsp,i1) = -i; // 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));
ta = ta * nrm; su2SubGroupIndex(i1, i2, aIndex);
} ta()()(i1, i2) = 1;
ta()()(i2, i1) = 1;
template <class cplx, ONLY_IF_Sp> ta()()(i1 + nsp, i2 + nsp) = -1;
static void generatorWtype(int wIndex, iSp2nMatrix<cplx> &ta) { ta()()(i2 + nsp, i1 + nsp) = -1;
// ta(i,i+N) = ta(i+N,i) = 1/2
ta = Zero();
RealD nrm = 1.0 / 2; //check
ta()()(wIndex,wIndex+nsp) = 1; ta = ta * nrm;
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; template <class cplx, ONLY_IF_Sp>
for (i1 = 0; spare >= (nsp - 1 - i1); i1++) { static void generatorBtype(int bIndex, iSp2nMatrix<cplx> &ta) {
spare = spare - (nsp - 1 - i1); // remove the Nc-1-i1 terms // 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
i2 = i1 + 1 + spare; // follows that j=i+1, ... , N-1
}
int i1, i2;
template<ONLY_IF_Sp> ta = Zero();
static void testGenerators(void) { cplx i(0.0, 1.0);
Matrix ta; RealD nrm = 1 / (2 * std::sqrt(2));
Matrix tb; su2SubGroupIndex(i1, i2, bIndex);
std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is 0.5 delta_ab " << std::endl;
for (int a = 0; a < AlgebraDimension; a++) { ta()()(i1, i2) = i;
for (int b = 0; b < AlgebraDimension; b++) { ta()()(i2, i1) = -i;
generator(a,ta); ta()()(i1 + nsp, i2 + nsp) = i;
generator(b,tb); ta()()(i2 + nsp, i1 + nsp) = -i;
Complex tr = TensorRemove(trace( ta * tb) );
std::cout << GridLogMessage << "(" << a << "," << b << ") = " << tr 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, GroupName::Sp) {
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;
}
static void testGenerators(GroupName::Sp) {
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; << std::endl;
if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6); if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6);
if (a != b) assert(abs(tr) < 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);
}
} }
}
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);
}
}
template <ONLY_IF_Sp>
static void OmegaInvariance(ColourMatrix &in) {
ColourMatrix Omega;
Omega = Zero();
std::cout << GridLogMessage << "I am a ColourMatrix" << std::endl;
// for (int i = 0; i < ncolour; i++) wrong?!
//{
// Omega()()(i, 2*ncolour-1-i) = 1.;
// Omega()()(2*ncolour-1-i, i) = -1;
// }
for (int i = 0; i < nsp; i++) {
Omega()()(i, nsp + i) = 1.;
Omega()()(nsp + i, i) = -1;
}
auto diff = Omega - (in * Omega * transpose(in));
auto sdiff = norm2(diff);
if (norm2(sdiff) < 1e-8) {
std::cout << GridLogMessage
<< "Symplectic condition satisfied: Omega invariant" << std::endl;
} else {
std::cout << GridLogMessage
<< "WARNING!!!!!! Matrix Omega NOT left invariant by "
<< norm2(sdiff) << std::endl;
}
}
template <typename GaugeField, ONLY_IF_Sp>
static void OmegaInvariance(GaugeField &in) {
typedef typename GaugeField::vector_type vector_type;
typedef iSp2nMatrix<vector_type> vMatrixType;
typedef Lattice<vMatrixType> LatticeMatrixType;
LatticeMatrixType U(in.Grid());
LatticeMatrixType OmegaLatt(in.Grid());
LatticeMatrixType identity(in.Grid());
RealD vol = in.Grid()->gSites();
ColourMatrix Omega;
OmegaLatt = Zero();
Omega = Zero();
identity = 1.;
std::cout << GridLogMessage << "I am a GaugeField " << std::endl;
U = PeekIndex<LorentzIndex>(in, 1);
OmegaInvariance(U);
}
template <ONLY_IF_Sp>
static void OmegaInvariance(LatticeColourMatrixD &in) {
LatticeColourMatrixD OmegaLatt(in.Grid());
LatticeColourMatrixD identity(in.Grid());
RealD vol = in.Grid()->gSites();
ColourMatrix Omega;
OmegaLatt = Zero();
Omega = Zero();
identity = 1.;
std::cout << GridLogMessage << "I am a LatticeColourMatrix " << std::endl;
for (int i = 0; i < nsp; i++) {
Omega()()(i, nsp + i) = 1.;
Omega()()(nsp + i, i) = -1;
}
std::cout << GridLogMessage << "Omega = " << Omega()() << std::endl;
OmegaLatt = OmegaLatt + (identity * Omega);
auto diff = OmegaLatt - (in * OmegaLatt * transpose(in));
auto sdiff = sum(diff);
// assert( norm2(sdiff) < 1e-8 );
if (norm2(sdiff) < 1e-8) {
std::cout << GridLogMessage
<< "Symplectic condition satisfied: Omega invariant" << std::endl;
} else {
std::cout << GridLogMessage
<< "WARNING!!!!!! Matrix Omega NOT left invariant by "
<< norm2(sdiff) << std::endl;
}
}