diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h index dc39c30b..5a8887f2 100644 --- a/Grid/qcd/action/gauge/GaugeImplTypes.h +++ b/Grid/qcd/action/gauge/GaugeImplTypes.h @@ -150,7 +150,7 @@ public: P = Ta(P); //const int nsp = Nc / 2; - Sp::iSp2nMatrix gen; + Sp::iGroupMatrix gen; auto Psum = P; diff --git a/Grid/qcd/utils/GaugeGroup.h b/Grid/qcd/utils/GaugeGroup.h new file mode 100644 index 00000000..24b89945 --- /dev/null +++ b/Grid/qcd/utils/GaugeGroup.h @@ -0,0 +1,454 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/utils/SUn.h + +Copyright (C) 2015 + +Author: Azusa Yamaguchi +Author: Peter Boyle +Author: neo +Author: paboyle + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ +/* END LEGAL */ +#ifndef QCD_UTIL_SUN_H +#define QCD_UTIL_SUN_H + +#define ONLY_IF_SU \ + typename dummy_name = group_name, \ + typename = std::enable_if_t < \ + std::is_same::value && \ + is_su::value > + +#define ONLY_IF_Sp \ + typename dummy_name = group_name, \ + typename = std::enable_if_t < \ + std::is_same::value && \ + is_sp::value > + +NAMESPACE_BEGIN(Grid); +namespace GroupName { +class SU {}; +class Sp {}; +} // namespace GroupName + +template +struct is_su { + static const bool value = false; +}; + +template <> +struct is_su { + static const bool value = true; +}; + +template +struct is_sp { + static const bool value = false; +}; + +template <> +struct is_sp { + static const bool value = true; +}; + +template +constexpr int compute_adjoint_dimension(int ncolour); + +template <> +constexpr int compute_adjoint_dimension(int ncolour) { + return ncolour * ncolour - 1; +} + +template <> +constexpr int compute_adjoint_dimension(int ncolour) { + return ncolour / 2 * (ncolour + 1); +} + +template +class GaugeGroup { + public: + static const int Dimension = ncolour; + static const int AdjointDimension = + compute_adjoint_dimension(ncolour); + static const int AlgebraDimension = + compute_adjoint_dimension(ncolour); + + template + using iSU2Matrix = iScalar > >; + template + using iGroupMatrix = iScalar > >; + template + using iAlgebraVector = iScalar > >; + static int su2subgroups(void) { return su2subgroups(group_name()); } + + ////////////////////////////////////////////////////////////////////////////////////////////////// + // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix, + // SU<2>::LatticeMatrix etc... + ////////////////////////////////////////////////////////////////////////////////////////////////// + typedef iGroupMatrix Matrix; + typedef iGroupMatrix MatrixF; + typedef iGroupMatrix MatrixD; + + typedef iGroupMatrix vMatrix; + typedef iGroupMatrix vMatrixF; + typedef iGroupMatrix vMatrixD; + + // For the projectors to the algebra + // these should be real... + // keeping complex for consistency with the SIMD vector types + typedef iAlgebraVector AlgebraVector; + typedef iAlgebraVector AlgebraVectorF; + typedef iAlgebraVector AlgebraVectorD; + + typedef iAlgebraVector vAlgebraVector; + typedef iAlgebraVector vAlgebraVectorF; + typedef iAlgebraVector vAlgebraVectorD; + + typedef Lattice LatticeMatrix; + typedef Lattice LatticeMatrixF; + typedef Lattice LatticeMatrixD; + + typedef Lattice LatticeAlgebraVector; + typedef Lattice LatticeAlgebraVectorF; + typedef Lattice LatticeAlgebraVectorD; + + typedef iSU2Matrix SU2Matrix; + typedef iSU2Matrix SU2MatrixF; + typedef iSU2Matrix SU2MatrixD; + + typedef iSU2Matrix vSU2Matrix; + typedef iSU2Matrix vSU2MatrixF; + typedef iSU2Matrix vSU2MatrixD; + + typedef Lattice LatticeSU2Matrix; + typedef Lattice LatticeSU2MatrixF; + typedef Lattice LatticeSU2MatrixD; + +#include "Grid/qcd/utils/SUn.h" +#include "Grid/qcd/utils/Sp2n.h" + + public: + template + static void generator(int lieIndex, iGroupMatrix &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) { + for (int gen = 0; gen < AdjointDimension; gen++) { + Matrix ta; + generator(gen, ta); + std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen + << std::endl; + std::cout << GridLogMessage << ta << std::endl; + } + } + + // reunitarise?? + template + static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, + double scale = 1.0) { + GridBase *grid = out.Grid(); + + typedef typename LatticeMatrixType::vector_type vector_type; + typedef typename LatticeMatrixType::scalar_type scalar_type; + + typedef iSinglet vTComplexType; + + typedef Lattice LatticeComplexType; + typedef typename GridTypeMapper< + typename LatticeMatrixType::vector_object>::scalar_object MatrixType; + + LatticeComplexType ca(grid); + LatticeMatrixType lie(grid); + LatticeMatrixType la(grid); + ComplexD ci(0.0, scale); + // ComplexD cone(1.0, 0.0); + MatrixType ta; + + lie = Zero(); + + for (int a = 0; a < AdjointDimension; a++) { + random(pRNG, ca); + + ca = (ca + conjugate(ca)) * 0.5; + ca = ca - 0.5; + + generator(a, ta); + + la = ci * ca * ta; + + lie = lie + la; // e^{i la ta} + } + taExp(lie, out); + } + + static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG, + LatticeMatrix &out, + Real scale = 1.0) { + GridBase *grid = out.Grid(); + LatticeReal ca(grid); + LatticeMatrix la(grid); + Complex ci(0.0, scale); + Matrix ta; + + out = Zero(); + for (int a = 0; a < AdjointDimension; a++) { + gaussian(pRNG, ca); + generator(a, ta); + la = toComplex(ca) * ta; + out += la; + } + out *= ci; + } + + static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h, + LatticeMatrix &out, + Real scale = 1.0) { + conformable(h, out); + GridBase *grid = out.Grid(); + LatticeMatrix la(grid); + Matrix ta; + + out = Zero(); + for (int a = 0; a < AdjointDimension; a++) { + generator(a, ta); + la = peekColour(h, a) * timesI(ta) * scale; + out += la; + } + } + + // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 + // ) inverse operation: FundamentalLieAlgebraMatrix + static void projectOnAlgebra(LatticeAlgebraVector &h_out, + const LatticeMatrix &in, Real scale = 1.0) { + conformable(h_out, in); + h_out = Zero(); + Matrix Ta; + + for (int a = 0; a < AdjointDimension; a++) { + generator(a, Ta); + pokeColour(h_out, -2.0 * (trace(timesI(Ta) * in)) * scale, a); + } + } + + template + static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) { + typedef typename GaugeField::vector_type vector_type; + typedef iGroupMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType Umu(out.Grid()); + for (int mu = 0; mu < Nd; mu++) { + LieRandomize(pRNG, Umu, 1.0); + PokeIndex(out, Umu, mu); + } + } + template + static void TepidConfiguration(GridParallelRNG &pRNG, GaugeField &out) { + typedef typename GaugeField::vector_type vector_type; + typedef iGroupMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType Umu(out.Grid()); + for (int mu = 0; mu < Nd; mu++) { + LieRandomize(pRNG, Umu, 0.01); + PokeIndex(out, Umu, mu); + } + } + template + static void ColdConfiguration(GaugeField &out) { + typedef typename GaugeField::vector_type vector_type; + typedef iGroupMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType Umu(out.Grid()); + Umu = 1.0; + for (int mu = 0; mu < Nd; mu++) { + PokeIndex(out, Umu, mu); + } + } + template + static void ColdConfiguration(GridParallelRNG &pRNG, GaugeField &out) { + ColdConfiguration(out); + } + + template + static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) { + out = Ta(in); + } + template + static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { + typedef typename LatticeMatrixType::scalar_type ComplexType; + + LatticeMatrixType xn(x.Grid()); + RealD nfac = 1.0; + + xn = x; + ex = xn + ComplexType(1.0); // 1+x + + // Do a 12th order exponentiation + for (int i = 2; i <= 12; ++i) { + nfac = nfac / RealD(i); // 1/2, 1/2.3 ... + xn = xn * x; // x2, x3,x4.... + ex = ex + xn * nfac; // x2/2!, x3/3!.... + } + } +}; + +template +LatticeComplexD Determinant( + const Lattice > > > &Umu) { + GridBase *grid = Umu.Grid(); + auto lvol = grid->lSites(); + LatticeComplexD ret(grid); + + autoView(Umu_v, Umu, CpuRead); + autoView(ret_v, ret, CpuWrite); + thread_for(site, lvol, { + Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N, N); + Coordinate lcoor; + grid->LocalIndexToLocalCoor(site, lcoor); + iScalar > > Us; + peekLocalSite(Us, Umu_v, lcoor); + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + EigenU(i, j) = Us()()(i, j); + } + } + ComplexD det = EigenU.determinant(); + pokeLocalSite(det, ret_v, lcoor); + }); + return ret; +} +template +static void ProjectSUn( + Lattice > > > &Umu) { + Umu = ProjectOnGroup(Umu); + auto det = Determinant(Umu); + + det = conjugate(det); + + for (int i = 0; i < N; i++) { + auto element = PeekIndex(Umu, N - 1, i); + element = element * det; + PokeIndex(Umu, element, Nc - 1, i); + } +} +template +static void ProjectSUn( + Lattice >, Nd> > &U) { + GridBase *grid = U.Grid(); + // Reunitarise + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnGroup(Umu); + ProjectSUn(Umu); + PokeIndex(U, Umu, mu); + } +} +// Explicit specialisation for SU(3). +// Explicit specialisation for SU(3). +static void ProjectSU3( + Lattice > > > &Umu) { + GridBase *grid = Umu.Grid(); + const int x = 0; + const int y = 1; + const int z = 2; + // Reunitarise + Umu = ProjectOnGroup(Umu); + autoView(Umu_v, Umu, CpuWrite); + thread_for(ss, grid->oSites(), { + auto cm = Umu_v[ss]; + cm()()(2, x) = adj(cm()()(0, y) * cm()()(1, z) - + cm()()(0, z) * cm()()(1, y)); // x= yz-zy + cm()()(2, y) = adj(cm()()(0, z) * cm()()(1, x) - + cm()()(0, x) * cm()()(1, z)); // y= zx-xz + cm()()(2, z) = adj(cm()()(0, x) * cm()()(1, y) - + cm()()(0, y) * cm()()(1, x)); // z= xy-yx + Umu_v[ss] = cm; + }); +} +static void ProjectSU3( + Lattice >, Nd> > &U) { + GridBase *grid = U.Grid(); + // Reunitarise + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnGroup(Umu); + ProjectSU3(Umu); + PokeIndex(U, Umu, mu); + } +} + +template +using SU = GaugeGroup; + +template +using Sp = GaugeGroup; + +typedef SU<2> SU2; +typedef SU<3> SU3; +typedef SU<4> SU4; +typedef SU<5> SU5; + +typedef SU FundamentalMatrices; + +template +static void ProjectSp2n( + Lattice > > > &Umu) { + Umu = ProjectOnSpGroup(Umu); + auto det = Determinant(Umu); // ok ? + + det = conjugate(det); + + for (int i = 0; i < N; i++) { + auto element = PeekIndex(Umu, N - 1, i); + element = element * det; + PokeIndex(Umu, element, Nc - 1, i); + } +} +template +static void ProjectSp2n( + Lattice >, Nd> > &U) { + GridBase *grid = U.Grid(); + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnSpGroup(Umu); + ProjectSp2n(Umu); + PokeIndex(U, Umu, mu); + } +} + +typedef Sp<2> Sp2; +typedef Sp<4> Sp4; +typedef Sp<6> Sp6; +typedef Sp<8> Sp8; + + +NAMESPACE_END(Grid); +#endif diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index 675493b3..61b19fbc 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -1,893 +1,550 @@ -/************************************************************************************* +// This file is #included into the body of the class template definition of +// GaugeGroup. So, image there to be +// +// template +// class GaugeGroup { +// +// around it. -Grid physics library, www.github.com/paboyle/Grid +private: +static int su2subgroups(GroupName::SU) { return (ncolour * (ncolour - 1)) / 2; } +//////////////////////////////////////////////////////////////////////// +// There are N^2-1 generators for SU(N). +// +// We take a traceless hermitian generator basis as follows +// +// * Normalisation: trace ta tb = 1/2 delta_ab = T_F delta_ab +// T_F = 1/2 for SU(N) groups +// +// * Off diagonal +// - pairs of rows i1,i2 behaving like pauli matrices signma_x, sigma_y +// +// - there are (Nc-1-i1) slots for i2 on each row [ x 0 x ] +// direct count off each row +// +// - Sum of all pairs is Nc(Nc-1)/2: proof arithmetic series +// +// (Nc-1) + (Nc-2)+... 1 ==> Nc*(Nc-1)/2 +// 1+ 2+ + + Nc-1 +// +// - There are 2 x Nc (Nc-1)/ 2 of these = Nc^2 - Nc +// +// - We enumerate the row-col pairs. +// - for each row col pair there is a (sigma_x) and a (sigma_y) like +// generator +// +// +// t^a_ij = { in 0.. Nc(Nc-1)/2 -1} => 1/2(delta_{i,i1} delta_{j,i2} + +// delta_{i,i1} delta_{j,i2}) +// t^a_ij = { in Nc(Nc-1)/2 ... Nc(Nc-1) - 1} => i/2( delta_{i,i1} +// delta_{j,i2} - i delta_{i,i1} delta_{j,i2}) +// +// * Diagonal; must be traceless and normalised +// - Sequence is +// N (1,-1,0,0...) +// N (1, 1,-2,0...) +// N (1, 1, 1,-3,0...) +// N (1, 1, 1, 1,-4,0...) +// +// where 1/2 = N^2 (1+.. m^2)etc.... for the m-th diagonal generator +// NB this gives the famous SU3 result for su2 index 8 +// +// N= sqrt(1/2 . 1/6 ) = 1/2 . 1/sqrt(3) +// +// ( 1 ) +// ( 1 ) / sqrt(3) /2 = 1/2 lambda_8 +// ( -2) +// +//////////////////////////////////////////////////////////////////////// +template +static void generator(int lieIndex, iGroupMatrix &ta, GroupName::SU) { + // map lie index to which type of generator + int diagIndex; + int su2Index; + int sigxy; + int NNm1 = ncolour * (ncolour - 1); + if (lieIndex >= NNm1) { + diagIndex = lieIndex - NNm1; + generatorDiagonal(diagIndex, ta); + return; + } + sigxy = lieIndex & 0x1; // even or odd + su2Index = lieIndex >> 1; + if (sigxy) + generatorSigmaY(su2Index, ta); + else + generatorSigmaX(su2Index, ta); +} -Source file: ./lib/qcd/utils/SUn.h +template +static void generatorSigmaY(int su2Index, iGroupMatrix &ta) { + ta = Zero(); + int i1, i2; + su2SubGroupIndex(i1, i2, su2Index); + ta()()(i1, i2) = 1.0; + ta()()(i2, i1) = 1.0; + ta = ta * 0.5; +} -Copyright (C) 2015 +template +static void generatorSigmaX(int su2Index, iGroupMatrix &ta) { + ta = Zero(); + cplx i(0.0, 1.0); + int i1, i2; + su2SubGroupIndex(i1, i2, su2Index); + ta()()(i1, i2) = i; + ta()()(i2, i1) = -i; + ta = ta * 0.5; +} -Author: Azusa Yamaguchi -Author: Peter Boyle -Author: neo -Author: paboyle +template +static void generatorDiagonal(int diagIndex, iGroupMatrix &ta) { + // diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...) + ta = Zero(); + int k = diagIndex + 1; // diagIndex starts from 0 + for (int i = 0; i <= diagIndex; i++) { // k iterations + ta()()(i, i) = 1.0; + } + ta()()(k, k) = -k; // indexing starts from 0 + RealD nrm = 1.0 / std::sqrt(2.0 * k * (k + 1)); + ta = ta * nrm; +} -This program is free software; you can redistribute it and/or modify -it under the terms of the GNU General Public License as published by -the Free Software Foundation; either version 2 of the License, or -(at your option) any later version. +//////////////////////////////////////////////////////////////////////// +// Map a su2 subgroup number to the pair of rows that are non zero +//////////////////////////////////////////////////////////////////////// +static void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::SU) { + assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2)); -This program is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. + int spare = su2_index; + for (i1 = 0; spare >= (ncolour - 1 - i1); i1++) { + spare = spare - (ncolour - 1 - i1); // remove the Nc-1-i1 terms + } + i2 = i1 + 1 + spare; +} -You should have received a copy of the GNU General Public License along -with this program; if not, write to the Free Software Foundation, Inc., -51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - -See the full license in the file "LICENSE" in the top level distribution -directory -*************************************************************************************/ - /* END LEGAL */ -#ifndef QCD_UTIL_SUN_H -#define QCD_UTIL_SUN_H - -NAMESPACE_BEGIN(Grid); - -template -class SU { public: - static const int Dimension = ncolour; - static const int AdjointDimension = ncolour * ncolour - 1; - static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; } - - template - using iSUnMatrix = iScalar > >; - template - using iSU2Matrix = iScalar > >; - template - using iSUnAlgebraVector = - iScalar > >; - - ////////////////////////////////////////////////////////////////////////////////////////////////// - // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix, - // SU<2>::LatticeMatrix etc... - ////////////////////////////////////////////////////////////////////////////////////////////////// - typedef iSUnMatrix Matrix; - typedef iSUnMatrix MatrixF; - typedef iSUnMatrix MatrixD; - - typedef iSUnMatrix vMatrix; - typedef iSUnMatrix vMatrixF; - typedef iSUnMatrix vMatrixD; - - // For the projectors to the algebra - // these should be real... - // keeping complex for consistency with the SIMD vector types - typedef iSUnAlgebraVector AlgebraVector; - typedef iSUnAlgebraVector AlgebraVectorF; - typedef iSUnAlgebraVector AlgebraVectorD; - - typedef iSUnAlgebraVector vAlgebraVector; - typedef iSUnAlgebraVector vAlgebraVectorF; - typedef iSUnAlgebraVector vAlgebraVectorD; - - typedef Lattice LatticeMatrix; - typedef Lattice LatticeMatrixF; - typedef Lattice LatticeMatrixD; - - typedef Lattice LatticeAlgebraVector; - typedef Lattice LatticeAlgebraVectorF; - typedef Lattice LatticeAlgebraVectorD; - - typedef iSU2Matrix SU2Matrix; - typedef iSU2Matrix SU2MatrixF; - typedef iSU2Matrix SU2MatrixD; - - typedef iSU2Matrix vSU2Matrix; - typedef iSU2Matrix vSU2MatrixF; - typedef iSU2Matrix vSU2MatrixD; - - typedef Lattice LatticeSU2Matrix; - typedef Lattice LatticeSU2MatrixF; - typedef Lattice LatticeSU2MatrixD; - - //////////////////////////////////////////////////////////////////////// - // There are N^2-1 generators for SU(N). - // - // We take a traceless hermitian generator basis as follows - // - // * Normalisation: trace ta tb = 1/2 delta_ab = T_F delta_ab - // T_F = 1/2 for SU(N) groups - // - // * Off diagonal - // - pairs of rows i1,i2 behaving like pauli matrices signma_x, sigma_y - // - // - there are (Nc-1-i1) slots for i2 on each row [ x 0 x ] - // direct count off each row - // - // - Sum of all pairs is Nc(Nc-1)/2: proof arithmetic series - // - // (Nc-1) + (Nc-2)+... 1 ==> Nc*(Nc-1)/2 - // 1+ 2+ + + Nc-1 - // - // - There are 2 x Nc (Nc-1)/ 2 of these = Nc^2 - Nc - // - // - We enumerate the row-col pairs. - // - for each row col pair there is a (sigma_x) and a (sigma_y) like - // generator - // - // - // t^a_ij = { in 0.. Nc(Nc-1)/2 -1} => 1/2(delta_{i,i1} delta_{j,i2} + - // delta_{i,i1} delta_{j,i2}) - // t^a_ij = { in Nc(Nc-1)/2 ... Nc(Nc-1) - 1} => i/2( delta_{i,i1} - // delta_{j,i2} - i delta_{i,i1} delta_{j,i2}) - // - // * Diagonal; must be traceless and normalised - // - Sequence is - // N (1,-1,0,0...) - // N (1, 1,-2,0...) - // N (1, 1, 1,-3,0...) - // N (1, 1, 1, 1,-4,0...) - // - // where 1/2 = N^2 (1+.. m^2)etc.... for the m-th diagonal generator - // NB this gives the famous SU3 result for su2 index 8 - // - // N= sqrt(1/2 . 1/6 ) = 1/2 . 1/sqrt(3) - // - // ( 1 ) - // ( 1 ) / sqrt(3) /2 = 1/2 lambda_8 - // ( -2) - // - //////////////////////////////////////////////////////////////////////// - template - static void generator(int lieIndex, iSUnMatrix &ta) { - // map lie index to which type of generator - int diagIndex; - int su2Index; - int sigxy; - int NNm1 = ncolour * (ncolour - 1); - if (lieIndex >= NNm1) { - diagIndex = lieIndex - NNm1; - generatorDiagonal(diagIndex, ta); - return; - } - sigxy = lieIndex & 0x1; // even or odd - su2Index = lieIndex >> 1; - if (sigxy) - generatorSigmaY(su2Index, ta); - else - generatorSigmaX(su2Index, ta); - } - - template - static void generatorSigmaY(int su2Index, iSUnMatrix &ta) { - ta = Zero(); - int i1, i2; - su2SubGroupIndex(i1, i2, su2Index); - ta()()(i1, i2) = 1.0; - ta()()(i2, i1) = 1.0; - ta = ta * 0.5; - } - - template - static void generatorSigmaX(int su2Index, iSUnMatrix &ta) { - ta = Zero(); - cplx i(0.0, 1.0); - int i1, i2; - su2SubGroupIndex(i1, i2, su2Index); - ta()()(i1, i2) = i; - ta()()(i2, i1) = -i; - ta = ta * 0.5; - } - - template - static void generatorDiagonal(int diagIndex, iSUnMatrix &ta) { - // diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...) - ta = Zero(); - int k = diagIndex + 1; // diagIndex starts from 0 - for (int i = 0; i <= diagIndex; i++) { // k iterations - ta()()(i, i) = 1.0; - } - ta()()(k, k) = -k; // indexing starts from 0 - RealD nrm = 1.0 / std::sqrt(2.0 * k * (k + 1)); - ta = ta * nrm; - } - - - - //////////////////////////////////////////////////////////////////////// - // Map a su2 subgroup number to the pair of rows that are non zero - //////////////////////////////////////////////////////////////////////// - static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { - assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2)); - - int spare = su2_index; - for (i1 = 0; spare >= (ncolour - 1 - i1); i1++) { - spare = spare - (ncolour - 1 - i1); // remove the Nc-1-i1 terms - } - i2 = i1 + 1 + spare; - } - - ////////////////////////////////////////////////////////////////////////////////////////// - // Pull out a subgroup and project on to real coeffs x pauli basis - ////////////////////////////////////////////////////////////////////////////////////////// - template - static void su2Extract(Lattice > &Determinant, - Lattice > &subgroup, - const Lattice > &source, - int su2_index) { - GridBase *grid(source.Grid()); - conformable(subgroup, source); - conformable(subgroup, Determinant); - int i0, i1; - su2SubGroupIndex(i0, i1, su2_index); - - autoView( subgroup_v , subgroup,AcceleratorWrite); - autoView( source_v , source,AcceleratorRead); - autoView( Determinant_v , Determinant,AcceleratorWrite); - accelerator_for(ss, grid->oSites(), 1, { - - subgroup_v[ss]()()(0, 0) = source_v[ss]()()(i0, i0); - subgroup_v[ss]()()(0, 1) = source_v[ss]()()(i0, i1); - subgroup_v[ss]()()(1, 0) = source_v[ss]()()(i1, i0); - subgroup_v[ss]()()(1, 1) = source_v[ss]()()(i1, i1); - - iSU2Matrix Sigma = subgroup_v[ss]; - - Sigma = Sigma - adj(Sigma) + trace(adj(Sigma)); - - subgroup_v[ss] = Sigma; - - // this should be purely real - Determinant_v[ss] = - Sigma()()(0, 0) * Sigma()()(1, 1) - Sigma()()(0, 1) * Sigma()()(1, 0); - }); - } - - ////////////////////////////////////////////////////////////////////////////////////////// - // Set matrix to one and insert a pauli subgroup - ////////////////////////////////////////////////////////////////////////////////////////// - template - static void su2Insert(const Lattice > &subgroup, - Lattice > &dest, int su2_index) { - GridBase *grid(dest.Grid()); - conformable(subgroup, dest); - int i0, i1; - su2SubGroupIndex(i0, i1, su2_index); - - dest = 1.0; // start out with identity - autoView( dest_v , dest, AcceleratorWrite); - autoView( subgroup_v, subgroup, AcceleratorRead); - accelerator_for(ss, grid->oSites(),1, - { - dest_v[ss]()()(i0, i0) = subgroup_v[ss]()()(0, 0); - dest_v[ss]()()(i0, i1) = subgroup_v[ss]()()(0, 1); - dest_v[ss]()()(i1, i0) = subgroup_v[ss]()()(1, 0); - dest_v[ss]()()(i1, i1) = subgroup_v[ss]()()(1, 1); - }); - - } - - /////////////////////////////////////////////// - // Generate e^{ Re Tr Staple Link} dlink - // - // *** Note Staple should be appropriate linear compbination between all - // staples. - // *** If already by beta pass coefficient 1.0. - // *** This routine applies the additional 1/Nc factor that comes after trace - // in action. - // - /////////////////////////////////////////////// - static void SubGroupHeatBath(GridSerialRNG &sRNG, GridParallelRNG &pRNG, - RealD beta, // coeff multiplying staple in action (with no 1/Nc) - LatticeMatrix &link, - const LatticeMatrix &barestaple, // multiplied by action coeffs so th - int su2_subgroup, int nheatbath, LatticeInteger &wheremask) - { - GridBase *grid = link.Grid(); - - const RealD twopi = 2.0 * M_PI; - - LatticeMatrix staple(grid); - - staple = barestaple * (beta / ncolour); - - LatticeMatrix V(grid); - V = link * staple; - - // Subgroup manipulation in the lie algebra space - LatticeSU2Matrix u(grid); // Kennedy pendleton "u" real projected normalised Sigma - LatticeSU2Matrix uinv(grid); - LatticeSU2Matrix ua(grid); // a in pauli form - LatticeSU2Matrix b(grid); // rotated matrix after hb - - // Some handy constant fields - LatticeComplex ones(grid); - ones = 1.0; - LatticeComplex zeros(grid); - zeros = Zero(); - LatticeReal rones(grid); - rones = 1.0; - LatticeReal rzeros(grid); - rzeros = Zero(); - LatticeComplex udet(grid); // determinant of real(staple) - LatticeInteger mask_true(grid); - mask_true = 1; - LatticeInteger mask_false(grid); - mask_false = 0; - - /* - PLB 156 P393 (1985) (Kennedy and Pendleton) - - Note: absorb "beta" into the def of sigma compared to KP paper; staple - passed to this routine has "beta" already multiplied in - - Action linear in links h and of form: - - beta S = beta Sum_p (1 - 1/Nc Re Tr Plaq ) - - Writing Sigma = 1/Nc (beta Sigma') where sum over staples is "Sigma' " - - beta S = const - beta/Nc Re Tr h Sigma' - = const - Re Tr h Sigma - - Decompose h and Sigma into (1, sigma_j) ; h_i real, h^2=1, Sigma_i complex - arbitrary. - - Tr h Sigma = h_i Sigma_j Tr (sigma_i sigma_j) = h_i Sigma_j 2 delta_ij - Re Tr h Sigma = 2 h_j Re Sigma_j - - Normalised re Sigma_j = xi u_j - - With u_j a unit vector and U can be in SU(2); - - Re Tr h Sigma = 2 h_j Re Sigma_j = 2 xi (h.u) - - 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] - u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] - - xi = sqrt(Det)/2; - - Write a= u h in SU(2); a has pauli decomp a_j; - - Note: Product b' xi is unvariant because scaling Sigma leaves - normalised vector "u" fixed; Can rescale Sigma so b' = 1. - */ - - //////////////////////////////////////////////////////// - // Real part of Pauli decomposition - // Note a subgroup can project to zero in cold start - //////////////////////////////////////////////////////// - su2Extract(udet, u, V, su2_subgroup); - - ////////////////////////////////////////////////////// - // Normalising this vector if possible; else identity - ////////////////////////////////////////////////////// - LatticeComplex xi(grid); - - LatticeSU2Matrix lident(grid); - - SU2Matrix ident = Complex(1.0); - SU2Matrix pauli1; - SU<2>::generator(0, pauli1); - SU2Matrix pauli2; - SU<2>::generator(1, pauli2); - SU2Matrix pauli3; - SU<2>::generator(2, pauli3); - pauli1 = timesI(pauli1) * 2.0; - pauli2 = timesI(pauli2) * 2.0; - pauli3 = timesI(pauli3) * 2.0; - - LatticeComplex cone(grid); - LatticeReal adet(grid); - adet = abs(toReal(udet)); - lident = Complex(1.0); - cone = Complex(1.0); - Real machine_epsilon = 1.0e-7; - u = where(adet > machine_epsilon, u, lident); - udet = where(adet > machine_epsilon, udet, cone); - - xi = 0.5 * sqrt(udet); // 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] - u = 0.5 * u * - pow(xi, -1.0); // u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] - - // Debug test for sanity - uinv = adj(u); - b = u * uinv - 1.0; - assert(norm2(b) < 1.0e-4); - - /* - Measure: Haar measure dh has d^4a delta(1-|a^2|) - In polars: - da = da0 r^2 sin theta dr dtheta dphi delta( 1 - r^2 -a0^2) - = da0 r^2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r)(sqrt(1-a0^) + - r) ) - = da0 r/2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r) ) - - Action factor Q(h) dh = e^-S[h] dh = e^{ xi Tr uh} dh // beta enters - through xi - = e^{2 xi (h.u)} dh - = e^{2 xi h0u0}.e^{2 xi h1u1}.e^{2 xi - h2u2}.e^{2 xi h3u3} dh - - Therefore for each site, take xi for that site - i) generate |a0|<1 with dist - (1-a0^2)^0.5 e^{2 xi a0 } da0 - - Take alpha = 2 xi = 2 xi [ recall 2 beta/Nc unmod staple norm]; hence 2.0/Nc - factor in Chroma ] - A. Generate two uniformly distributed pseudo-random numbers R and R', R'', - R''' in the unit interval; - B. Set X = -(ln R)/alpha, X' =-(ln R')/alpha; - C. Set C = cos^2(2pi R"), with R" another uniform random number in [0,1] ; - D. Set A = XC; - E. Let d = X'+A; - F. If R'''^2 :> 1 - 0.5 d, go back to A; - G. Set a0 = 1 - d; - - Note that in step D setting B ~ X - A and using B in place of A in step E will - generate a second independent a 0 value. - */ - - ///////////////////////////////////////////////////////// - // count the number of sites by picking "1"'s out of hat - ///////////////////////////////////////////////////////// - Integer hit = 0; - LatticeReal rtmp(grid); - rtmp = where(wheremask, rones, rzeros); - RealD numSites = sum(rtmp); - RealD numAccepted; - LatticeInteger Accepted(grid); - Accepted = Zero(); - LatticeInteger newlyAccepted(grid); - - std::vector xr(4, grid); - std::vector a(4, grid); - LatticeReal d(grid); - d = Zero(); - LatticeReal alpha(grid); - - // std::cout< 1 - 0.5 d, go back to A; - LatticeReal thresh(grid); - thresh = 1.0 - d * 0.5; - xrsq = xr[0] * xr[0]; - LatticeInteger ione(grid); - ione = 1; - LatticeInteger izero(grid); - izero = Zero(); - - newlyAccepted = where(xrsq < thresh, ione, izero); - Accepted = where(newlyAccepted, newlyAccepted, Accepted); - Accepted = where(wheremask, Accepted, izero); - - // FIXME need an iSum for integer to avoid overload on return type?? - rtmp = where(Accepted, rones, rzeros); - numAccepted = sum(rtmp); - - hit++; - - } while ((numAccepted < numSites) && (hit < nheatbath)); - - // G. Set a0 = 1 - d; - a[0] = Zero(); - a[0] = where(wheremask, 1.0 - d, a[0]); - - ////////////////////////////////////////// - // ii) generate a_i uniform on two sphere radius (1-a0^2)^0.5 - ////////////////////////////////////////// - - LatticeReal a123mag(grid); - a123mag = sqrt(abs(1.0 - a[0] * a[0])); - - LatticeReal cos_theta(grid); - LatticeReal sin_theta(grid); - LatticeReal phi(grid); - - random(pRNG, phi); - phi = phi * twopi; // uniform in [0,2pi] - random(pRNG, cos_theta); - cos_theta = (cos_theta * 2.0) - 1.0; // uniform in [-1,1] - sin_theta = sqrt(abs(1.0 - cos_theta * cos_theta)); - - a[1] = a123mag * sin_theta * cos(phi); - a[2] = a123mag * sin_theta * sin(phi); - a[3] = a123mag * cos_theta; - - ua = toComplex(a[0]) * ident + toComplex(a[1]) * pauli1 + - toComplex(a[2]) * pauli2 + toComplex(a[3]) * pauli3; - - b = 1.0; - b = where(wheremask, uinv * ua, b); - su2Insert(b, V, su2_subgroup); - - // mask the assignment back based on Accptance - link = where(Accepted, V * link, link); - - ////////////////////////////// - // Debug Checks - // SU2 check - LatticeSU2Matrix check(grid); // rotated matrix after hb - u = Zero(); - check = ua * adj(ua) - 1.0; - check = where(Accepted, check, u); - assert(norm2(check) < 1.0e-4); - - check = b * adj(b) - 1.0; - check = where(Accepted, check, u); - assert(norm2(check) < 1.0e-4); - - LatticeMatrix Vcheck(grid); - Vcheck = Zero(); - Vcheck = where(Accepted, V * adj(V) - 1.0, Vcheck); - // std::cout< +static void su2Extract(Lattice > &Determinant, + Lattice > &subgroup, + const Lattice > &source, + int su2_index) { + GridBase *grid(source.Grid()); + conformable(subgroup, source); + conformable(subgroup, Determinant); + int i0, i1; + su2SubGroupIndex(i0, i1, su2_index); + + autoView(subgroup_v, subgroup, AcceleratorWrite); + autoView(source_v, source, AcceleratorRead); + autoView(Determinant_v, Determinant, AcceleratorWrite); + accelerator_for(ss, grid->oSites(), 1, { + subgroup_v[ss]()()(0, 0) = source_v[ss]()()(i0, i0); + subgroup_v[ss]()()(0, 1) = source_v[ss]()()(i0, i1); + subgroup_v[ss]()()(1, 0) = source_v[ss]()()(i1, i0); + subgroup_v[ss]()()(1, 1) = source_v[ss]()()(i1, i1); + + iSU2Matrix Sigma = subgroup_v[ss]; + + Sigma = Sigma - adj(Sigma) + trace(adj(Sigma)); + + subgroup_v[ss] = Sigma; + + // this should be purely real + Determinant_v[ss] = + Sigma()()(0, 0) * Sigma()()(1, 1) - Sigma()()(0, 1) * Sigma()()(1, 0); + }); +} + +////////////////////////////////////////////////////////////////////////////////////////// +// Set matrix to one and insert a pauli subgroup +////////////////////////////////////////////////////////////////////////////////////////// +template +static void su2Insert(const Lattice > &subgroup, + Lattice > &dest, int su2_index) { + GridBase *grid(dest.Grid()); + conformable(subgroup, dest); + int i0, i1; + su2SubGroupIndex(i0, i1, su2_index); + + dest = 1.0; // start out with identity + autoView(dest_v, dest, AcceleratorWrite); + autoView(subgroup_v, subgroup, AcceleratorRead); + accelerator_for(ss, grid->oSites(), 1, { + dest_v[ss]()()(i0, i0) = subgroup_v[ss]()()(0, 0); + dest_v[ss]()()(i0, i1) = subgroup_v[ss]()()(0, 1); + dest_v[ss]()()(i1, i0) = subgroup_v[ss]()()(1, 0); + dest_v[ss]()()(i1, i1) = subgroup_v[ss]()()(1, 1); + }); +} + +/////////////////////////////////////////////// +// Generate e^{ Re Tr Staple Link} dlink +// +// *** Note Staple should be appropriate linear compbination between all +// staples. +// *** If already by beta pass coefficient 1.0. +// *** This routine applies the additional 1/Nc factor that comes after trace +// in action. +// +/////////////////////////////////////////////// +template +static void SubGroupHeatBath( + GridSerialRNG &sRNG, GridParallelRNG &pRNG, + RealD beta, // coeff multiplying staple in action (with no 1/Nc) + LatticeMatrix &link, + const LatticeMatrix &barestaple, // multiplied by action coeffs so th + int su2_subgroup, int nheatbath, LatticeInteger &wheremask) { + GridBase *grid = link.Grid(); + + const RealD twopi = 2.0 * M_PI; + + LatticeMatrix staple(grid); + + staple = barestaple * (beta / ncolour); + + LatticeMatrix V(grid); + V = link * staple; + + // Subgroup manipulation in the lie algebra space + LatticeSU2Matrix u( + grid); // Kennedy pendleton "u" real projected normalised Sigma + LatticeSU2Matrix uinv(grid); + LatticeSU2Matrix ua(grid); // a in pauli form + LatticeSU2Matrix b(grid); // rotated matrix after hb + + // Some handy constant fields + LatticeComplex ones(grid); + ones = 1.0; + LatticeComplex zeros(grid); + zeros = Zero(); + LatticeReal rones(grid); + rones = 1.0; + LatticeReal rzeros(grid); + rzeros = Zero(); + LatticeComplex udet(grid); // determinant of real(staple) + LatticeInteger mask_true(grid); + mask_true = 1; + LatticeInteger mask_false(grid); + mask_false = 0; + + /* + PLB 156 P393 (1985) (Kennedy and Pendleton) + + Note: absorb "beta" into the def of sigma compared to KP paper; staple + passed to this routine has "beta" already multiplied in + + Action linear in links h and of form: + + beta S = beta Sum_p (1 - 1/Nc Re Tr Plaq ) + + Writing Sigma = 1/Nc (beta Sigma') where sum over staples is "Sigma' " + + beta S = const - beta/Nc Re Tr h Sigma' + = const - Re Tr h Sigma + + Decompose h and Sigma into (1, sigma_j) ; h_i real, h^2=1, Sigma_i complex + arbitrary. + + Tr h Sigma = h_i Sigma_j Tr (sigma_i sigma_j) = h_i Sigma_j 2 delta_ij + Re Tr h Sigma = 2 h_j Re Sigma_j + + Normalised re Sigma_j = xi u_j + + With u_j a unit vector and U can be in SU(2); + + Re Tr h Sigma = 2 h_j Re Sigma_j = 2 xi (h.u) + + 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] + u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] + + xi = sqrt(Det)/2; + + Write a= u h in SU(2); a has pauli decomp a_j; + + Note: Product b' xi is unvariant because scaling Sigma leaves + normalised vector "u" fixed; Can rescale Sigma so b' = 1. + */ + + //////////////////////////////////////////////////////// + // Real part of Pauli decomposition + // Note a subgroup can project to zero in cold start + //////////////////////////////////////////////////////// + su2Extract(udet, u, V, su2_subgroup); + + ////////////////////////////////////////////////////// + // Normalising this vector if possible; else identity + ////////////////////////////////////////////////////// + LatticeComplex xi(grid); + + LatticeSU2Matrix lident(grid); + + SU2Matrix ident = Complex(1.0); + SU2Matrix pauli1; + GaugeGroup<2, GroupName::SU>::generator(0, pauli1); + SU2Matrix pauli2; + GaugeGroup<2, GroupName::SU>::generator(1, pauli2); + SU2Matrix pauli3; + GaugeGroup<2, GroupName::SU>::generator(2, pauli3); + pauli1 = timesI(pauli1) * 2.0; + pauli2 = timesI(pauli2) * 2.0; + pauli3 = timesI(pauli3) * 2.0; + + LatticeComplex cone(grid); + LatticeReal adet(grid); + adet = abs(toReal(udet)); + lident = Complex(1.0); + cone = Complex(1.0); + Real machine_epsilon = 1.0e-7; + u = where(adet > machine_epsilon, u, lident); + udet = where(adet > machine_epsilon, udet, cone); + + xi = 0.5 * sqrt(udet); // 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] + u = 0.5 * u * pow(xi, -1.0); // u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] + + // Debug test for sanity + uinv = adj(u); + b = u * uinv - 1.0; + assert(norm2(b) < 1.0e-4); + + /* + Measure: Haar measure dh has d^4a delta(1-|a^2|) + In polars: + da = da0 r^2 sin theta dr dtheta dphi delta( 1 - r^2 -a0^2) + = da0 r^2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r)(sqrt(1-a0^) + + r) ) + = da0 r/2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r) ) + + Action factor Q(h) dh = e^-S[h] dh = e^{ xi Tr uh} dh // beta + enters through xi = e^{2 xi (h.u)} dh = e^{2 xi h0u0}.e^{2 xi h1u1}.e^{2 + xi h2u2}.e^{2 xi h3u3} dh + + Therefore for each site, take xi for that site + i) generate |a0|<1 with dist + (1-a0^2)^0.5 e^{2 xi a0 } da0 + + Take alpha = 2 xi = 2 xi [ recall 2 beta/Nc unmod staple norm]; + hence 2.0/Nc factor in Chroma ] A. Generate two uniformly distributed + pseudo-random numbers R and R', R'', R''' in the unit interval; B. Set X = + -(ln R)/alpha, X' =-(ln R')/alpha; C. Set C = cos^2(2pi R"), with R" + another uniform random number in [0,1] ; D. Set A = XC; E. Let d = X'+A; + F. If R'''^2 :> 1 - 0.5 d, go back to A; + G. Set a0 = 1 - d; + + Note that in step D setting B ~ X - A and using B in place of A in step E + will generate a second independent a 0 value. + */ + + ///////////////////////////////////////////////////////// + // count the number of sites by picking "1"'s out of hat + ///////////////////////////////////////////////////////// + Integer hit = 0; + LatticeReal rtmp(grid); + rtmp = where(wheremask, rones, rzeros); + RealD numSites = sum(rtmp); + RealD numAccepted; + LatticeInteger Accepted(grid); + Accepted = Zero(); + LatticeInteger newlyAccepted(grid); + + std::vector xr(4, grid); + std::vector a(4, grid); + LatticeReal d(grid); + d = Zero(); + LatticeReal alpha(grid); + + // std::cout< 1 - 0.5 d, go back to A; + LatticeReal thresh(grid); + thresh = 1.0 - d * 0.5; + xrsq = xr[0] * xr[0]; + LatticeInteger ione(grid); + ione = 1; + LatticeInteger izero(grid); + izero = Zero(); + + newlyAccepted = where(xrsq < thresh, ione, izero); + Accepted = where(newlyAccepted, newlyAccepted, Accepted); + Accepted = where(wheremask, Accepted, izero); + + // FIXME need an iSum for integer to avoid overload on return type?? + rtmp = where(Accepted, rones, rzeros); + numAccepted = sum(rtmp); + + hit++; + + } while ((numAccepted < numSites) && (hit < nheatbath)); + + // G. Set a0 = 1 - d; + a[0] = Zero(); + a[0] = where(wheremask, 1.0 - d, a[0]); + + ////////////////////////////////////////// + // ii) generate a_i uniform on two sphere radius (1-a0^2)^0.5 + ////////////////////////////////////////// + + LatticeReal a123mag(grid); + a123mag = sqrt(abs(1.0 - a[0] * a[0])); + + LatticeReal cos_theta(grid); + LatticeReal sin_theta(grid); + LatticeReal phi(grid); + + random(pRNG, phi); + phi = phi * twopi; // uniform in [0,2pi] + random(pRNG, cos_theta); + cos_theta = (cos_theta * 2.0) - 1.0; // uniform in [-1,1] + sin_theta = sqrt(abs(1.0 - cos_theta * cos_theta)); + + a[1] = a123mag * sin_theta * cos(phi); + a[2] = a123mag * sin_theta * sin(phi); + a[3] = a123mag * cos_theta; + + ua = toComplex(a[0]) * ident + toComplex(a[1]) * pauli1 + + toComplex(a[2]) * pauli2 + toComplex(a[3]) * pauli3; + + b = 1.0; + b = where(wheremask, uinv * ua, b); + su2Insert(b, V, su2_subgroup); + + // mask the assignment back based on Accptance + link = where(Accepted, V * link, link); + + ////////////////////////////// + // Debug Checks + // SU2 check + LatticeSU2Matrix check(grid); // rotated matrix after hb + u = Zero(); + check = ua * adj(ua) - 1.0; + check = where(Accepted, check, u); + assert(norm2(check) < 1.0e-4); + + check = b * adj(b) - 1.0; + check = where(Accepted, check, u); + assert(norm2(check) < 1.0e-4); + + LatticeMatrix Vcheck(grid); + Vcheck = Zero(); + Vcheck = where(Accepted, V * adj(V) - 1.0, Vcheck); + // std::cout< +static void testGenerators(GroupName::SU) { + Matrix ta; + Matrix tb; + std::cout << GridLogMessage + << "Fundamental - Checking trace ta tb is 0.5 delta_ab" + << std::endl; + for (int a = 0; a < AdjointDimension; a++) { + for (int b = 0; b < AdjointDimension; b++) { + generator(a, ta); + generator(b, tb); + Complex tr = TensorRemove(trace(ta * tb)); + std::cout << GridLogMessage << "(" << a << "," << b << ") = " << tr << std::endl; - std::cout << GridLogMessage << ta << std::endl; - } - } - - - - 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 < AdjointDimension; a++) { - for (int b = 0; b < AdjointDimension; 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 < AdjointDimension; 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 < AdjointDimension; a++) { - generator(a, ta); - Complex tr = TensorRemove(trace(ta)); - std::cout << GridLogMessage << a << " " << std::endl; - assert(abs(tr) < 1.0e-6); + 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; } - - // reunitarise?? - template - static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, double scale = 1.0) - { - GridBase *grid = out.Grid(); - - typedef typename LatticeMatrixType::vector_type vector_type; - typedef typename LatticeMatrixType::scalar_type scalar_type; - - typedef iSinglet vTComplexType; - - typedef Lattice LatticeComplexType; - typedef typename GridTypeMapper::scalar_object MatrixType; - - LatticeComplexType ca(grid); - LatticeMatrixType lie(grid); - LatticeMatrixType la(grid); - ComplexD ci(0.0, scale); - // ComplexD cone(1.0, 0.0); - MatrixType ta; - - lie = Zero(); - - for (int a = 0; a < AdjointDimension; a++) { - random(pRNG, ca); - - ca = (ca + conjugate(ca)) * 0.5; - ca = ca - 0.5; - - generator(a, ta); - - la = ci * ca * ta; - - lie = lie + la; // e^{i la ta} - - } - taExp(lie, out); + std::cout << GridLogMessage << "Fundamental - Checking if hermitian" + << std::endl; + for (int a = 0; a < AdjointDimension; a++) { + generator(a, ta); + std::cout << GridLogMessage << a << std::endl; + assert(norm2(ta - adj(ta)) < 1.0e-6); } + std::cout << GridLogMessage << std::endl; - static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG, - LatticeMatrix &out, - Real scale = 1.0) { - GridBase *grid = out.Grid(); - LatticeReal ca(grid); - LatticeMatrix la(grid); - Complex ci(0.0, scale); - Matrix ta; - - out = Zero(); - for (int a = 0; a < AdjointDimension; a++) { - gaussian(pRNG, ca); - generator(a, ta); - la = toComplex(ca) * ta; - out += la; - } - out *= ci; + std::cout << GridLogMessage << "Fundamental - Checking if traceless" + << std::endl; + for (int a = 0; a < AdjointDimension; 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; +} - static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h, - LatticeMatrix &out, - Real scale = 1.0) { - conformable(h, out); - GridBase *grid = out.Grid(); - LatticeMatrix la(grid); - Matrix ta; - - out = Zero(); - for (int a = 0; a < AdjointDimension; a++) { - generator(a, ta); - la = peekColour(h, a) * timesI(ta) * scale; - out += la; - } - } /* * Fundamental rep gauge xform */ - template - static void GaugeTransformFundamental( Fundamental &ferm, GaugeMat &g){ - GridBase *grid = ferm._grid; - conformable(grid,g._grid); - ferm = g*ferm; - } +template +static void GaugeTransformFundamental(Fundamental &ferm, GaugeMat &g) { + GridBase *grid = ferm._grid; + conformable(grid, g._grid); + ferm = g * ferm; +} /* * Adjoint rep gauge xform */ - template - static void GaugeTransform( GaugeField &Umu, GaugeMat &g){ - GridBase *grid = Umu.Grid(); - conformable(grid,g.Grid()); +template +static void GaugeTransform(GaugeField &Umu, GaugeMat &g) { + GridBase *grid = Umu.Grid(); + conformable(grid, g.Grid()); - GaugeMat U(grid); - GaugeMat ag(grid); ag = adj(g); + GaugeMat U(grid); + GaugeMat ag(grid); + ag = adj(g); - for(int mu=0;mu(Umu,mu); - U = g*U*Cshift(ag, mu, 1); - PokeIndex(Umu,U,mu); - } - } - template - static void GaugeTransform( std::vector &U, GaugeMat &g){ - GridBase *grid = g.Grid(); - GaugeMat ag(grid); ag = adj(g); - for(int mu=0;mu - static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, GaugeMat &g){ - LieRandomize(pRNG,g,1.0); - GaugeTransform(Umu,g); - } - - // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 ) - // inverse operation: FundamentalLieAlgebraMatrix - static void projectOnAlgebra(LatticeAlgebraVector &h_out, const LatticeMatrix &in, Real scale = 1.0) { - conformable(h_out, in); - h_out = Zero(); - Matrix Ta; - - for (int a = 0; a < AdjointDimension; a++) { - generator(a, Ta); - pokeColour(h_out, - 2.0 * (trace(timesI(Ta) * in)) * scale, a); - } - } - - template - static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) { - typedef typename GaugeField::vector_type vector_type; - typedef iSUnMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType Umu(out.Grid()); - for (int mu = 0; mu < Nd; mu++) { - LieRandomize(pRNG, Umu, 1.0); - PokeIndex(out, Umu, mu); - } - } - template - static void TepidConfiguration(GridParallelRNG &pRNG,GaugeField &out){ - typedef typename GaugeField::vector_type vector_type; - typedef iSUnMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType Umu(out.Grid()); - for(int mu=0;mu(out,Umu,mu); - } - } - template - static void ColdConfiguration(GaugeField &out){ - typedef typename GaugeField::vector_type vector_type; - typedef iSUnMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType Umu(out.Grid()); - Umu=1.0; - for(int mu=0;mu(out,Umu,mu); - } - } - template - static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){ - ColdConfiguration(out); - } - - template - static void taProj( const LatticeMatrixType &in, LatticeMatrixType &out){ - out = Ta(in); - } - template - static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { - typedef typename LatticeMatrixType::scalar_type ComplexType; - - LatticeMatrixType xn(x.Grid()); - RealD nfac = 1.0; - - xn = x; - ex = xn + ComplexType(1.0); // 1+x - - // Do a 12th order exponentiation - for (int i = 2; i <= 12; ++i) { - nfac = nfac / RealD(i); // 1/2, 1/2.3 ... - xn = xn * x; // x2, x3,x4.... - ex = ex + xn * nfac; // x2/2!, x3/3!.... - } - } -}; - -template -LatticeComplexD Determinant(const Lattice > > > &Umu) -{ - GridBase *grid=Umu.Grid(); - auto lvol = grid->lSites(); - LatticeComplexD ret(grid); - - autoView(Umu_v,Umu,CpuRead); - autoView(ret_v,ret,CpuWrite); - thread_for(site,lvol,{ - Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N,N); - Coordinate lcoor; - grid->LocalIndexToLocalCoor(site, lcoor); - iScalar > > Us; - peekLocalSite(Us, Umu_v, lcoor); - for(int i=0;i -static void ProjectSUn(Lattice > > > &Umu) -{ - Umu = ProjectOnGroup(Umu); - auto det = Determinant(Umu); - - det = conjugate(det); - - for(int i=0;i(Umu,N-1,i); - element = element * det; - PokeIndex(Umu,element,Nc-1,i); + for (int mu = 0; mu < Nd; mu++) { + U = PeekIndex(Umu, mu); + U = g * U * Cshift(ag, mu, 1); + PokeIndex(Umu, U, mu); } } -template -static void ProjectSUn(Lattice >,Nd> > &U) -{ - GridBase *grid=U.Grid(); - // Reunitarise - for(int mu=0;mu(U,mu); - Umu = ProjectOnGroup(Umu); - ProjectSUn(Umu); - PokeIndex(U,Umu,mu); +template +static void GaugeTransform(std::vector &U, GaugeMat &g) { + GridBase *grid = g.Grid(); + GaugeMat ag(grid); + ag = adj(g); + for (int mu = 0; mu < Nd; mu++) { + U[mu] = g * U[mu] * Cshift(ag, mu, 1); } } -// Explicit specialisation for SU(3). -// Explicit specialisation for SU(3). -static void -ProjectSU3 (Lattice > > > &Umu) -{ - GridBase *grid=Umu.Grid(); - const int x=0; - const int y=1; - const int z=2; - // Reunitarise - Umu = ProjectOnGroup(Umu); - autoView(Umu_v,Umu,CpuWrite); - thread_for(ss,grid->oSites(),{ - auto cm = Umu_v[ss]; - cm()()(2,x) = adj(cm()()(0,y)*cm()()(1,z)-cm()()(0,z)*cm()()(1,y)); //x= yz-zy - cm()()(2,y) = adj(cm()()(0,z)*cm()()(1,x)-cm()()(0,x)*cm()()(1,z)); //y= zx-xz - cm()()(2,z) = adj(cm()()(0,x)*cm()()(1,y)-cm()()(0,y)*cm()()(1,x)); //z= xy-yx - Umu_v[ss]=cm; - }); +template +static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, + GaugeMat &g) { + LieRandomize(pRNG, g, 1.0); + GaugeTransform(Umu, g); } -static void ProjectSU3(Lattice >,Nd> > &U) -{ - GridBase *grid=U.Grid(); - // Reunitarise - for(int mu=0;mu(U,mu); - Umu = ProjectOnGroup(Umu); - ProjectSU3(Umu); - PokeIndex(U,Umu,mu); - } -} - -typedef SU<2> SU2; -typedef SU<3> SU3; -typedef SU<4> SU4; -typedef SU<5> SU5; - - -typedef SU FundamentalMatrices; - -NAMESPACE_END(Grid); -#endif diff --git a/Grid/qcd/utils/Sp2n.h b/Grid/qcd/utils/Sp2n.h index 19e30568..6d559c0a 100644 --- a/Grid/qcd/utils/Sp2n.h +++ b/Grid/qcd/utils/Sp2n.h @@ -1,580 +1,334 @@ -#ifndef QCD_UTIL_Sp2n_H -#define QCD_UTIL_Sp2n_H +private: +static int su2subgroups(GroupName::Sp) { return (ncolour/2 * (ncolour/2 - 1)) / 2; } -NAMESPACE_BEGIN(Grid); +// 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 -// Sp(2N) -// ncolour = 2N +template +static void generator(int lieIndex, iGroupMatrix &ta, GroupName::Sp) { + // map lie index into type of generators: diagonal, abcd type, wz type + const int nsp = ncolour/2; + int diagIndex; + int aIndex, bIndex, cIndex, dIndex; + int wIndex, zIndex; // a,b,c,d are N(N-1)/2 and w,z are N + const int mod = nsp * (nsp - 1) * 0.5; + const int offdiag = + 2 * nsp * nsp; // number of generators not in the cartan subalgebra + const int wmod = 4 * mod; + const 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; + } -template -class Sp { -public: - static const int nsp = ncolour/2; - static const int Dimension = ncolour; - static const int AlgebraDimension = nsp*(2*nsp +1); - static int su2subgroups(void) { return (nsp * (nsp - 1)) / 2; } - - - template - using iSp2nMatrix = iScalar > >; - template - using iSU2Matrix = iScalar > >; - template - using iSp2nAlgebraVector = iScalar > >; - - typedef iSp2nMatrix Matrix; - typedef iSp2nMatrix MatrixF; - typedef iSp2nMatrix MatrixD; - - typedef iSp2nMatrix vMatrix; - typedef iSp2nMatrix vMatrixF; - typedef iSp2nMatrix vMatrixD; - - typedef iSp2nAlgebraVector AlgebraVector; - typedef iSp2nAlgebraVector AlgebraVectorF; - typedef iSp2nAlgebraVector AlgebraVectorD; - - typedef iSp2nAlgebraVector vAlgebraVector; - typedef iSp2nAlgebraVector vAlgebraVectorF; - typedef iSp2nAlgebraVector vAlgebraVectorD; - - typedef Lattice LatticeMatrix; - typedef Lattice LatticeMatrixF; - typedef Lattice LatticeMatrixD; - - typedef Lattice LatticeAlgebraVector; - typedef Lattice LatticeAlgebraVectorF; - typedef Lattice LatticeAlgebraVectorD; - - - - // 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 - static void generator(int lieIndex, iSp2nMatrix &ta) { - // map lie index into type of generators: diagonal, abcd type, wz type +} // end of generator - 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 - static void generatorDiagtype(int diagIndex, iSp2nMatrix &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; +template +static void generatorDiagtype(int diagIndex, iGroupMatrix &ta) { + // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra - ta()()(diagIndex,diagIndex) = nrm; - ta()()(diagIndex+nsp,diagIndex+nsp) = -nrm; - } - - template - static void generatorAtype(int aIndex, iSp2nMatrix &ta) { - - // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2) - // with i - static void generatorBtype(int bIndex, iSp2nMatrix &ta) { - - // ta(i,j) = -ta(j,i) = ta(i+N,j+N) = -ta(j+N,i+N) = i / 1/ 2 sqrt(2) - // with i +static void generatorAtype(int aIndex, iGroupMatrix &ta) { + // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2) + // with i - static void generatorCtype(int cIndex, iSp2nMatrix &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; + 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 - static void generatorDtype(int dIndex, iSp2nMatrix &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 = ta * nrm; +} - ta()()(i1,i2+nsp) = i; - ta()()(i2,i1+nsp) = i; - ta()()(i1+nsp,i2) = -i; - ta()()(i2+nsp,i1) = -i; +template +static void generatorBtype(int bIndex, iGroupMatrix &ta) { + // ta(i,j) = -ta(j,i) = ta(i+N,j+N) = -ta(j+N,i+N) = i / 1/ 2 sqrt(2) + // with i - static void generatorWtype(int wIndex, iSp2nMatrix &ta) { - - // ta(i,i+N) = ta(i+N,i) = 1/2 - - ta = Zero(); - RealD nrm = 1.0 / 2; //check + const int nsp=ncolour/2; + int i1, i2; + ta = Zero(); + cplx i(0.0, 1.0); + RealD nrm = 1 / (2 * std::sqrt(2)); + su2SubGroupIndex(i1, i2, bIndex); - ta()()(wIndex,wIndex+nsp) = 1; - ta()()(wIndex+nsp,wIndex) = 1; - - ta = ta * nrm; - } - - template - static void generatorZtype(int zIndex, iSp2nMatrix &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 - //////////////////////////////////////////////////////////////////////// - static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { - assert((su2_index >= 0) && (su2_index < (nsp * (nsp - 1)) / 2)); + ta()()(i1, i2) = i; + ta()()(i2, i1) = -i; + ta()()(i1 + nsp, i2 + nsp) = i; + ta()()(i2 + nsp, i1 + nsp) = -i; - 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 printGenerators(void) { - for (int gen = 0; gen < AlgebraDimension; gen++) { - Matrix ta; - generator(gen, ta); - std::cout << GridLogMessage << "Nc = " << ncolour << std::endl; - std::cout << GridLogMessage << " t_" << gen << std::endl; - std::cout << GridLogMessage << ta << std::endl; - } - } - - - - 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 + ta = ta * nrm; +} + +template +static void generatorCtype(int cIndex, iGroupMatrix &ta) { + // ta(i,j+N) = ta(j,i+N) = ta(i+N,j) = ta(j+N,i) = 1 / 2 sqrt(2) + + const int nsp=ncolour/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 +static void generatorDtype(int dIndex, iGroupMatrix &ta) { + // ta(i,j+N) = ta(j,i+N) = -ta(i+N,j) = -ta(j+N,i) = i / 2 sqrt(2) + + const int nsp=ncolour/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 +static void generatorWtype(int wIndex, iGroupMatrix &ta) { + // ta(i,i+N) = ta(i+N,i) = 1/2 + + const int nsp=ncolour/2; + ta = Zero(); + RealD nrm = 1.0 / 2; // check + + ta()()(wIndex, wIndex + nsp) = 1; + ta()()(wIndex + nsp, wIndex) = 1; + + ta = ta * nrm; +} + +template +static void generatorZtype(int zIndex, iGroupMatrix &ta) { + // ta(i,i+N) = - ta(i+N,i) = i/2 + + const int nsp=ncolour/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 +static void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::Sp) { + const int nsp=ncolour/2; + 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; - 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); - } - + if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6); + if (a != b) assert(abs(tr) < 1.0e-6); } - - - template - static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, double scale = 1.0) - { - GridBase *grid = out.Grid(); + } + 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); + } +} - typedef typename LatticeMatrixType::vector_type vector_type; - typedef typename LatticeMatrixType::scalar_type scalar_type; +public: +template +static void OmegaInvariance(ColourMatrix &in) { + ColourMatrix Omega; + Omega = Zero(); + const int nsp=ncolour/2; - typedef iSinglet vTComplexType; + std::cout << GridLogMessage << "I am a ColourMatrix" << std::endl; - typedef Lattice LatticeComplexType; - typedef typename GridTypeMapper::scalar_object MatrixType; + // 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; + } - LatticeComplexType ca(grid); - LatticeMatrixType lie(grid); - LatticeMatrixType la(grid); - ComplexD ci(0.0, scale); - // ComplexD cone(1.0, 0.0); - MatrixType ta; + 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; + } +} - lie = Zero(); +template +static void OmegaInvariance(GaugeField &in) { + typedef typename GaugeField::vector_type vector_type; + typedef iGroupMatrix vMatrixType; + typedef Lattice LatticeMatrixType; - for (int a = 0; a < AlgebraDimension; a++) { - random(pRNG, ca); + LatticeMatrixType U(in.Grid()); + LatticeMatrixType OmegaLatt(in.Grid()); + LatticeMatrixType identity(in.Grid()); + RealD vol = in.Grid()->gSites(); + ColourMatrix Omega; - ca = (ca + conjugate(ca)) * 0.5; - ca = ca - 0.5; + OmegaLatt = Zero(); + Omega = Zero(); + identity = 1.; - generator(a, ta); + std::cout << GridLogMessage << "I am a GaugeField " << std::endl; - la = ci * ca * ta; + U = PeekIndex(in, 1); - lie = lie + la; // e^{i la ta} + OmegaInvariance(U); +} - } - taExp(lie, out); - } - - - static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG, //same as sun - LatticeMatrix &out, - Real scale = 1.0) { - GridBase *grid = out.Grid(); - LatticeReal ca(grid); - LatticeMatrix la(grid); - Complex ci(0.0, scale); - Matrix ta; +template +static void OmegaInvariance(LatticeColourMatrixD &in) { + const int nsp=ncolour/2; + LatticeColourMatrixD OmegaLatt(in.Grid()); + LatticeColourMatrixD identity(in.Grid()); + RealD vol = in.Grid()->gSites(); + ColourMatrix Omega; - out = Zero(); - for (int a = 0; a < AlgebraDimension; a++) { - gaussian(pRNG, ca); - generator(a, ta); - la = toComplex(ca) * ta; - out += la; - } - out *= ci; - } - - static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h, - LatticeMatrix &out, - Real scale = 1.0) { - conformable(h, out); - GridBase *grid = out.Grid(); - LatticeMatrix la(grid); - Matrix ta; + OmegaLatt = Zero(); + Omega = Zero(); + identity = 1.; - out = Zero(); - for (int a = 0; a < AlgebraDimension; a++) { - generator(a, ta); - la = peekColour(h, a) * timesI(ta) * scale; - out += la; - } - } - - - template - static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { // same as sun - typedef typename LatticeMatrixType::scalar_type ComplexType; + std::cout << GridLogMessage << "I am a LatticeColourMatrix " << std::endl; - LatticeMatrixType xn(x.Grid()); - RealD nfac = 1.0; + for (int i = 0; i < nsp; i++) { + Omega()()(i, nsp + i) = 1.; + Omega()()(nsp + i, i) = -1; + } - xn = x; - ex = xn + ComplexType(1.0); // 1+x + std::cout << GridLogMessage << "Omega = " << Omega()() << std::endl; + OmegaLatt = OmegaLatt + (identity * Omega); - // Do a 12th order exponentiation - for (int i = 2; i <= 12; ++i) { - nfac = nfac / RealD(i); // 1/2, 1/2.3 ... - xn = xn * x; // x2, x3,x4.... - ex = ex + xn * nfac; // x2/2!, x3/3!.... - } - } - - static void projectOnAlgebra(LatticeAlgebraVector &h_out, const LatticeMatrix &in, Real scale = 1.0) { - conformable(h_out, in); - h_out = Zero(); - Matrix Ta; + 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; + } +} - for (int a = 0; a < AlgebraDimension; a++) { - generator(a, Ta); - pokeColour(h_out, - 2.0 * (trace(timesI(Ta) * in)) * scale, a); - } - } - - - template - static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) { - typedef typename GaugeField::vector_type vector_type; - typedef iSp2nMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType Umu(out.Grid()); - for (int mu = 0; mu < Nd; mu++) { - LieRandomize(pRNG, Umu, 1.0); //def - PokeIndex(out, Umu, mu); - } - } - template - static void TepidConfiguration(GridParallelRNG &pRNG,GaugeField &out){ - typedef typename GaugeField::vector_type vector_type; - typedef iSp2nMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType Umu(out.Grid()); - for(int mu=0;mu(out,Umu,mu); - } - } - - template - static void ColdConfiguration(GaugeField &out){ - typedef typename GaugeField::vector_type vector_type; - typedef iSp2nMatrix vMatrixType; - typedef Lattice LatticeMatrixType; - - LatticeMatrixType Umu(out.Grid()); - Umu=1.0; - for(int mu=0;mu(out,Umu,mu); - } - } - template - static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){ - ColdConfiguration(out); - } - - 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 - static void OmegaInvariance(GaugeField &in) - { - typedef typename GaugeField::vector_type vector_type; - typedef iSp2nMatrix vMatrixType; - typedef Lattice 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(in,1); - - OmegaInvariance(U); - - } - - 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; - } - - - } - - -}; // end of class Sp - - - template - static void ProjectSp2n(Lattice > > > &Umu) - { - Umu = ProjectOnSpGroup(Umu); - auto det = Determinant(Umu); // ok ? - - det = conjugate(det); - - for(int i=0;i(Umu,N-1,i); - element = element * det; - PokeIndex(Umu,element,Nc-1,i); - } - } - template - static void ProjectSp2n(Lattice >,Nd> > &U) - { - GridBase *grid=U.Grid(); - for(int mu=0;mu(U,mu); - Umu = ProjectOnSpGroup(Umu); - ProjectSp2n(Umu); - PokeIndex(U,Umu,mu); - } - } - -typedef Sp<2> Sp2; -typedef Sp<4> Sp4; -typedef Sp<6> Sp6; -typedef Sp<8> Sp8; - -NAMESPACE_END(Grid); -#endif diff --git a/Grid/qcd/utils/Sp2nTwoIndex.h b/Grid/qcd/utils/Sp2nTwoIndex.h index 3b685802..98306970 100644 --- a/Grid/qcd/utils/Sp2nTwoIndex.h +++ b/Grid/qcd/utils/Sp2nTwoIndex.h @@ -175,10 +175,10 @@ public: template static void generator(int Index, iSp2nTwoIndexMatrix &i2indTa) { - Vector::template iSp2nMatrix > ta( + Vector > ta( NumGenerators); - Vector::template iSp2nMatrix > eij(Dimension); - typename Sp::template iSp2nMatrix tmp; + Vector > eij(Dimension); + iSp2nMatrix tmp; i2indTa = Zero(); for (int a = 0; a < NumGenerators; a++) @@ -189,7 +189,7 @@ public: for (int a = 0; a < Dimension; a++) { tmp = transpose(ta[Index]) * adj(eij[a]) + adj(eij[a]) * ta[Index]; for (int b = 0; b < Dimension; b++) { - typename Sp::template iSp2nMatrix tmp1 = + iSp2nMatrix tmp1 = tmp * eij[b]; Complex iTr = TensorRemove(timesI(trace(tmp1))); i2indTa()()(a, b) = iTr; diff --git a/Grid/qcd/utils/Utils.h b/Grid/qcd/utils/Utils.h index c654ce18..87196cdc 100644 --- a/Grid/qcd/utils/Utils.h +++ b/Grid/qcd/utils/Utils.h @@ -8,8 +8,7 @@ #include // Include representations -#include -#include +#include #include #include #include diff --git a/configure.ac b/configure.ac index 3a4ed469..04a3cdf4 100644 --- a/configure.ac +++ b/configure.ac @@ -838,7 +838,7 @@ AC_CONFIG_FILES(tests/lanczos/Makefile) AC_CONFIG_FILES(tests/smearing/Makefile) AC_CONFIG_FILES(tests/qdpxx/Makefile) AC_CONFIG_FILES(tests/testu01/Makefile) -AC_CONFIG_FILES(tests/Sp2n/Makefile) +AC_CONFIG_FILES(tests/sp2n/Makefile) AC_CONFIG_FILES(benchmarks/Makefile) AC_CONFIG_FILES(examples/Makefile) AC_OUTPUT diff --git a/scripts/filelist b/scripts/filelist index 4e64cd62..d91da5f1 100755 --- a/scripts/filelist +++ b/scripts/filelist @@ -37,7 +37,6 @@ cd $home/tests dirs=`find . -type d -not -path '*/\.*'` for subdir in $dirs; do cd $home/tests/$subdir - pwd TESTS=`ls T*.cc` TESTLIST=`echo ${TESTS} | sed s/.cc//g ` PREF=`[ $subdir = '.' ] && echo noinst || echo EXTRA` diff --git a/tests/core/Test_lie_generators.cc b/tests/core/Test_lie_generators.cc index e044378c..372e1408 100644 --- a/tests/core/Test_lie_generators.cc +++ b/tests/core/Test_lie_generators.cc @@ -32,7 +32,7 @@ directory #include -#include +#include #include #include diff --git a/tests/core/test_sp2n_lie_gen.cc b/tests/core/test_sp2n_lie_gen.cc deleted file mode 100644 index 6c3766e6..00000000 --- a/tests/core/test_sp2n_lie_gen.cc +++ /dev/null @@ -1,55 +0,0 @@ - - - -#include -#include - -using namespace Grid; - -int main(int argc, char** argv) { - Grid_init(&argc, &argv); - - //std::vector latt({4, 4, 4, 8}); - //GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid( - //latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); - //GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid); - - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "* Generators for Sp(2)" << std::endl; - std::cout << GridLogMessage << "*********************************************" - << std::endl; - - Sp2::printGenerators(); - Sp2::testGenerators(); - - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "* Generators for Sp(4)" << std::endl; - std::cout << GridLogMessage << "*********************************************" - << std::endl; - - Sp4::printGenerators(); - Sp4::testGenerators(); - - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "* Generators for Sp(6)" << std::endl; - std::cout << GridLogMessage << "*********************************************" - << std::endl; - - Sp6::printGenerators(); - Sp6::testGenerators(); - - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "* Generators for Sp(8)" << std::endl; - std::cout << GridLogMessage << "*********************************************" - << std::endl; - - Sp8::printGenerators(); - Sp8::testGenerators(); - - - Grid_finalize(); -} diff --git a/tests/sp2n/Makefile.am b/tests/sp2n/Makefile.am index 0ee7421c..346c8ad0 100644 --- a/tests/sp2n/Makefile.am +++ b/tests/sp2n/Makefile.am @@ -4,4 +4,5 @@ include Make.inc check: tests ./Test_project_on_Sp - ./test_sp2n_lie_gen + ./Test_sp2n_lie_gen + ./Test_Sp_start diff --git a/tests/sp2n/Test_Sp_start.cc b/tests/sp2n/Test_Sp_start.cc index 16a66367..061acbcc 100644 --- a/tests/sp2n/Test_Sp_start.cc +++ b/tests/sp2n/Test_Sp_start.cc @@ -21,7 +21,7 @@ int main (int argc, char **argv) double vol = Umu.Grid()->gSites(); - const int nsp = Sp::nsp; + const int nsp = Nc/2; identity = 1.; Cidentity = 1.; diff --git a/tests/sp2n/Test_project_on_Sp.cc b/tests/sp2n/Test_project_on_Sp.cc index d1cf33bd..36b95cf4 100644 --- a/tests/sp2n/Test_project_on_Sp.cc +++ b/tests/sp2n/Test_project_on_Sp.cc @@ -20,8 +20,7 @@ int main (int argc, char **argv) LatticeColourMatrixD aux(&Grid); LatticeColourMatrixD identity(&Grid); - //const int nsp = Nc / 2; - const int nsp = Sp::nsp; + const int nsp = Nc / 2; identity = 1.0; RealD epsilon = 0.01; diff --git a/tests/sp2n/Test_sp2n_lie_gen.cc b/tests/sp2n/Test_sp2n_lie_gen.cc index 6c3766e6..5c6430a3 100644 --- a/tests/sp2n/Test_sp2n_lie_gen.cc +++ b/tests/sp2n/Test_sp2n_lie_gen.cc @@ -1,8 +1,4 @@ - - - #include -#include using namespace Grid;