Merge pull request #14 from chillenzer/unify_gauge_groups

Unify gauge groups (closes #5)

Grid/qcd/action/gauge/GaugeImplTypes.h

@@ -150,7 +150,7 @@ public:
             P = Ta(P);
             //const int nsp = Nc / 2;
 
-            Sp<Nc>::iSp2nMatrix<Complex> gen;
+            Sp<Nc>::iGroupMatrix<Complex> gen;
             
 
             auto Psum = P;

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 <ayamaguc@staffmail.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: neo <cossu@post.kek.jp>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+
+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<dummy_name, group_name>::value && \
+                      is_su<dummy_name>::value >
+
+#define ONLY_IF_Sp                                                       \
+  typename dummy_name = group_name,                                      \
+           typename = std::enable_if_t <                                 \
+                          std::is_same<dummy_name, group_name>::value && \
+                      is_sp<dummy_name>::value >
+
+NAMESPACE_BEGIN(Grid);
+namespace GroupName {
+class SU {};
+class Sp {};
+}  // namespace GroupName
+
+template <typename group_name>
+struct is_su {
+  static const bool value = false;
+};
+
+template <>
+struct is_su<GroupName::SU> {
+  static const bool value = true;
+};
+
+template <typename group_name>
+struct is_sp {
+  static const bool value = false;
+};
+
+template <>
+struct is_sp<GroupName::Sp> {
+  static const bool value = true;
+};
+
+template <typename group_name>
+constexpr int compute_adjoint_dimension(int ncolour);
+
+template <>
+constexpr int compute_adjoint_dimension<GroupName::SU>(int ncolour) {
+  return ncolour * ncolour - 1;
+}
+
+template <>
+constexpr int compute_adjoint_dimension<GroupName::Sp>(int ncolour) {
+  return ncolour / 2 * (ncolour + 1);
+}
+
+template <int ncolour, class group_name>
+class GaugeGroup {
+ public:
+  static const int Dimension = ncolour;
+  static const int AdjointDimension =
+      compute_adjoint_dimension<group_name>(ncolour);
+  static const int AlgebraDimension =
+      compute_adjoint_dimension<group_name>(ncolour);
+
+  template <typename vtype>
+  using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >;
+  template <typename vtype>
+  using iGroupMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >;
+  template <typename vtype>
+  using iAlgebraVector = iScalar<iScalar<iVector<vtype, AdjointDimension> > >;
+  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<Complex> Matrix;
+  typedef iGroupMatrix<ComplexF> MatrixF;
+  typedef iGroupMatrix<ComplexD> MatrixD;
+
+  typedef iGroupMatrix<vComplex> vMatrix;
+  typedef iGroupMatrix<vComplexF> vMatrixF;
+  typedef iGroupMatrix<vComplexD> vMatrixD;
+
+  // For the projectors to the algebra
+  // these should be real...
+  // keeping complex for consistency with the SIMD vector types
+  typedef iAlgebraVector<Complex> AlgebraVector;
+  typedef iAlgebraVector<ComplexF> AlgebraVectorF;
+  typedef iAlgebraVector<ComplexD> AlgebraVectorD;
+
+  typedef iAlgebraVector<vComplex> vAlgebraVector;
+  typedef iAlgebraVector<vComplexF> vAlgebraVectorF;
+  typedef iAlgebraVector<vComplexD> vAlgebraVectorD;
+
+  typedef Lattice<vMatrix> LatticeMatrix;
+  typedef Lattice<vMatrixF> LatticeMatrixF;
+  typedef Lattice<vMatrixD> LatticeMatrixD;
+
+  typedef Lattice<vAlgebraVector> LatticeAlgebraVector;
+  typedef Lattice<vAlgebraVectorF> LatticeAlgebraVectorF;
+  typedef Lattice<vAlgebraVectorD> LatticeAlgebraVectorD;
+
+  typedef iSU2Matrix<Complex> SU2Matrix;
+  typedef iSU2Matrix<ComplexF> SU2MatrixF;
+  typedef iSU2Matrix<ComplexD> SU2MatrixD;
+
+  typedef iSU2Matrix<vComplex> vSU2Matrix;
+  typedef iSU2Matrix<vComplexF> vSU2MatrixF;
+  typedef iSU2Matrix<vComplexD> vSU2MatrixD;
+
+  typedef Lattice<vSU2Matrix> LatticeSU2Matrix;
+  typedef Lattice<vSU2MatrixF> LatticeSU2MatrixF;
+  typedef Lattice<vSU2MatrixD> LatticeSU2MatrixD;
+
+#include "Grid/qcd/utils/SUn.h"
+#include "Grid/qcd/utils/Sp2n.h"
+
+ public:
+  template <class cplx>
+  static void generator(int lieIndex, iGroupMatrix<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) {
+    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 <typename LatticeMatrixType>
+  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<vector_type> vTComplexType;
+
+    typedef Lattice<vTComplexType> 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 <typename GaugeField>
+  static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
+    typedef typename GaugeField::vector_type vector_type;
+    typedef iGroupMatrix<vector_type> vMatrixType;
+    typedef Lattice<vMatrixType> LatticeMatrixType;
+
+    LatticeMatrixType Umu(out.Grid());
+    for (int mu = 0; mu < Nd; mu++) {
+      LieRandomize(pRNG, Umu, 1.0);
+      PokeIndex<LorentzIndex>(out, Umu, mu);
+    }
+  }
+  template <typename GaugeField>
+  static void TepidConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
+    typedef typename GaugeField::vector_type vector_type;
+    typedef iGroupMatrix<vector_type> vMatrixType;
+    typedef Lattice<vMatrixType> LatticeMatrixType;
+
+    LatticeMatrixType Umu(out.Grid());
+    for (int mu = 0; mu < Nd; mu++) {
+      LieRandomize(pRNG, Umu, 0.01);
+      PokeIndex<LorentzIndex>(out, Umu, mu);
+    }
+  }
+  template <typename GaugeField>
+  static void ColdConfiguration(GaugeField &out) {
+    typedef typename GaugeField::vector_type vector_type;
+    typedef iGroupMatrix<vector_type> vMatrixType;
+    typedef Lattice<vMatrixType> LatticeMatrixType;
+
+    LatticeMatrixType Umu(out.Grid());
+    Umu = 1.0;
+    for (int mu = 0; mu < Nd; mu++) {
+      PokeIndex<LorentzIndex>(out, Umu, mu);
+    }
+  }
+  template <typename GaugeField>
+  static void ColdConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
+    ColdConfiguration(out);
+  }
+
+  template <typename LatticeMatrixType, ONLY_IF_SU>
+  static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) {
+    out = Ta(in);
+  }
+  template <typename LatticeMatrixType>
+  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 <int N>
+LatticeComplexD Determinant(
+    const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &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<iScalar<iMatrix<ComplexD, N> > > 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 <int N>
+static void ProjectSUn(
+    Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
+  Umu = ProjectOnGroup(Umu);
+  auto det = Determinant(Umu);
+
+  det = conjugate(det);
+
+  for (int i = 0; i < N; i++) {
+    auto element = PeekIndex<ColourIndex>(Umu, N - 1, i);
+    element = element * det;
+    PokeIndex<ColourIndex>(Umu, element, Nc - 1, i);
+  }
+}
+template <int N>
+static void ProjectSUn(
+    Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
+  GridBase *grid = U.Grid();
+  // Reunitarise
+  for (int mu = 0; mu < Nd; mu++) {
+    auto Umu = PeekIndex<LorentzIndex>(U, mu);
+    Umu = ProjectOnGroup(Umu);
+    ProjectSUn(Umu);
+    PokeIndex<LorentzIndex>(U, Umu, mu);
+  }
+}
+// Explicit specialisation for SU(3).
+// Explicit specialisation for SU(3).
+static void ProjectSU3(
+    Lattice<iScalar<iScalar<iMatrix<vComplexD, 3> > > > &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<iVector<iScalar<iMatrix<vComplexD, 3> >, Nd> > &U) {
+  GridBase *grid = U.Grid();
+  // Reunitarise
+  for (int mu = 0; mu < Nd; mu++) {
+    auto Umu = PeekIndex<LorentzIndex>(U, mu);
+    Umu = ProjectOnGroup(Umu);
+    ProjectSU3(Umu);
+    PokeIndex<LorentzIndex>(U, Umu, mu);
+  }
+}
+
+template <int ncolour>
+using SU = GaugeGroup<ncolour, GroupName::SU>;
+
+template <int ncolour>
+using Sp = GaugeGroup<ncolour, GroupName::Sp>;
+
+typedef SU<2> SU2;
+typedef SU<3> SU3;
+typedef SU<4> SU4;
+typedef SU<5> SU5;
+
+typedef SU<Nc> FundamentalMatrices;
+
+template <int N>
+static void ProjectSp2n(
+    Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
+  Umu = ProjectOnSpGroup(Umu);
+  auto det = Determinant(Umu);  // ok ?
+
+  det = conjugate(det);
+
+  for (int i = 0; i < N; i++) {
+    auto element = PeekIndex<ColourIndex>(Umu, N - 1, i);
+    element = element * det;
+    PokeIndex<ColourIndex>(Umu, element, Nc - 1, i);
+  }
+}
+template <int N>
+static void ProjectSp2n(
+    Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
+  GridBase *grid = U.Grid();
+  for (int mu = 0; mu < Nd; mu++) {
+    auto Umu = PeekIndex<LorentzIndex>(U, mu);
+    Umu = ProjectOnSpGroup(Umu);
+    ProjectSp2n(Umu);
+    PokeIndex<LorentzIndex>(U, Umu, mu);
+  }
+}
+
+typedef Sp<2> Sp2;
+typedef Sp<4> Sp4;
+typedef Sp<6> Sp6;
+typedef Sp<8> Sp8;
+
+
+NAMESPACE_END(Grid);
+#endif

Grid/qcd/utils/SUn.h

@@ -1,97 +1,13 @@
-/*************************************************************************************
+// This file is #included into the body of the class template definition of
-
+// GaugeGroup. So, image there to be
-Grid physics library, www.github.com/paboyle/Grid
+//
-
+// template <int ncolour, class group_name>
-Source file: ./lib/qcd/utils/SUn.h
+// class GaugeGroup {
-
+//
-Copyright (C) 2015
+// around it.
-
-Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: neo <cossu@post.kek.jp>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-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
-
-NAMESPACE_BEGIN(Grid);
-
-template <int ncolour>
-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 <typename vtype>
-  using iSUnMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >;
-  template <typename vtype>
-  using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >;
-  template <typename vtype>
-  using iSUnAlgebraVector =
-    iScalar<iScalar<iVector<vtype, AdjointDimension> > >;
-
-  //////////////////////////////////////////////////////////////////////////////////////////////////
-  // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix,
-  // SU<2>::LatticeMatrix etc...
-  //////////////////////////////////////////////////////////////////////////////////////////////////
-  typedef iSUnMatrix<Complex> Matrix;
-  typedef iSUnMatrix<ComplexF> MatrixF;
-  typedef iSUnMatrix<ComplexD> MatrixD;
-
-  typedef iSUnMatrix<vComplex> vMatrix;
-  typedef iSUnMatrix<vComplexF> vMatrixF;
-  typedef iSUnMatrix<vComplexD> vMatrixD;
-
-  // For the projectors to the algebra
-  // these should be real...
-  // keeping complex for consistency with the SIMD vector types
-  typedef iSUnAlgebraVector<Complex> AlgebraVector;
-  typedef iSUnAlgebraVector<ComplexF> AlgebraVectorF;
-  typedef iSUnAlgebraVector<ComplexD> AlgebraVectorD;
-
-  typedef iSUnAlgebraVector<vComplex> vAlgebraVector;
-  typedef iSUnAlgebraVector<vComplexF> vAlgebraVectorF;
-  typedef iSUnAlgebraVector<vComplexD> vAlgebraVectorD;
-
-  typedef Lattice<vMatrix> LatticeMatrix;
-  typedef Lattice<vMatrixF> LatticeMatrixF;
-  typedef Lattice<vMatrixD> LatticeMatrixD;
-
-  typedef Lattice<vAlgebraVector> LatticeAlgebraVector;
-  typedef Lattice<vAlgebraVectorF> LatticeAlgebraVectorF;
-  typedef Lattice<vAlgebraVectorD> LatticeAlgebraVectorD;
-
-  typedef iSU2Matrix<Complex> SU2Matrix;
-  typedef iSU2Matrix<ComplexF> SU2MatrixF;
-  typedef iSU2Matrix<ComplexD> SU2MatrixD;
-
-  typedef iSU2Matrix<vComplex> vSU2Matrix;
-  typedef iSU2Matrix<vComplexF> vSU2MatrixF;
-  typedef iSU2Matrix<vComplexD> vSU2MatrixD;
-
-  typedef Lattice<vSU2Matrix> LatticeSU2Matrix;
-  typedef Lattice<vSU2MatrixF> LatticeSU2MatrixF;
-  typedef Lattice<vSU2MatrixD> LatticeSU2MatrixD;
 
+private:
+static int su2subgroups(GroupName::SU) { return (ncolour * (ncolour - 1)) / 2; }
 ////////////////////////////////////////////////////////////////////////
 // There are N^2-1 generators for SU(N).
 //
@@ -141,7 +57,7 @@ public:
 //
 ////////////////////////////////////////////////////////////////////////
 template <class cplx>
-  static void generator(int lieIndex, iSUnMatrix<cplx> &ta) {
+static void generator(int lieIndex, iGroupMatrix<cplx> &ta, GroupName::SU) {
   // map lie index to which type of generator
   int diagIndex;
   int su2Index;
@@ -160,8 +76,8 @@ public:
     generatorSigmaX(su2Index, ta);
 }
 
-  template <class cplx>
+template <class cplx, ONLY_IF_SU>
-  static void generatorSigmaY(int su2Index, iSUnMatrix<cplx> &ta) {
+static void generatorSigmaY(int su2Index, iGroupMatrix<cplx> &ta) {
   ta = Zero();
   int i1, i2;
   su2SubGroupIndex(i1, i2, su2Index);
@@ -170,8 +86,8 @@ public:
   ta = ta * 0.5;
 }
 
-  template <class cplx>
+template <class cplx, ONLY_IF_SU>
-  static void generatorSigmaX(int su2Index, iSUnMatrix<cplx> &ta) {
+static void generatorSigmaX(int su2Index, iGroupMatrix<cplx> &ta) {
   ta = Zero();
   cplx i(0.0, 1.0);
   int i1, i2;
@@ -181,8 +97,8 @@ public:
   ta = ta * 0.5;
 }
 
-  template <class cplx>
+template <class cplx, ONLY_IF_SU>
-  static void generatorDiagonal(int diagIndex, iSUnMatrix<cplx> &ta) {
+static void generatorDiagonal(int diagIndex, iGroupMatrix<cplx> &ta) {
   // diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...)
   ta = Zero();
   int k = diagIndex + 1;                  // diagIndex starts from 0
@@ -194,12 +110,10 @@ public:
   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) {
+static void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::SU) {
   assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2));
 
   int spare = su2_index;
@@ -209,13 +123,14 @@ public:
   i2 = i1 + 1 + spare;
 }
 
+public:
 //////////////////////////////////////////////////////////////////////////////////////////
 // Pull out a subgroup and project on to real coeffs x pauli basis
 //////////////////////////////////////////////////////////////////////////////////////////
-  template <class vcplx>
+template <class vcplx, ONLY_IF_SU>
 static void su2Extract(Lattice<iSinglet<vcplx> > &Determinant,
                        Lattice<iSU2Matrix<vcplx> > &subgroup,
-                         const Lattice<iSUnMatrix<vcplx> > &source,
+                       const Lattice<iGroupMatrix<vcplx> > &source,
                        int su2_index) {
   GridBase *grid(source.Grid());
   conformable(subgroup, source);
@@ -227,7 +142,6 @@ public:
   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);
@@ -248,9 +162,9 @@ public:
 //////////////////////////////////////////////////////////////////////////////////////////
 // Set matrix to one and insert a pauli subgroup
 //////////////////////////////////////////////////////////////////////////////////////////
-  template <class vcplx>
+template <class vcplx, ONLY_IF_SU>
 static void su2Insert(const Lattice<iSU2Matrix<vcplx> > &subgroup,
-                        Lattice<iSUnMatrix<vcplx> > &dest, int su2_index) {
+                      Lattice<iGroupMatrix<vcplx> > &dest, int su2_index) {
   GridBase *grid(dest.Grid());
   conformable(subgroup, dest);
   int i0, i1;
@@ -259,14 +173,12 @@ public:
   dest = 1.0;  // start out with identity
   autoView(dest_v, dest, AcceleratorWrite);
   autoView(subgroup_v, subgroup, AcceleratorRead);
-    accelerator_for(ss, grid->oSites(),1,
+  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);
   });
-
 }
 
 ///////////////////////////////////////////////
@@ -279,12 +191,13 @@ public:
 // in action.
 //
 ///////////////////////////////////////////////
-  static void SubGroupHeatBath(GridSerialRNG &sRNG, GridParallelRNG &pRNG,
+template <ONLY_IF_SU>
+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) 
+    int su2_subgroup, int nheatbath, LatticeInteger &wheremask) {
-  {
   GridBase *grid = link.Grid();
 
   const RealD twopi = 2.0 * M_PI;
@@ -297,7 +210,8 @@ public:
   V = link * staple;
 
   // Subgroup manipulation in the lie algebra space
-    LatticeSU2Matrix u(grid);  // Kennedy pendleton "u" real projected normalised Sigma
+  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
@@ -370,11 +284,11 @@ public:
 
   SU2Matrix ident = Complex(1.0);
   SU2Matrix pauli1;
-    SU<2>::generator(0, pauli1);
+  GaugeGroup<2, GroupName::SU>::generator(0, pauli1);
   SU2Matrix pauli2;
-    SU<2>::generator(1, pauli2);
+  GaugeGroup<2, GroupName::SU>::generator(1, pauli2);
   SU2Matrix pauli3;
-    SU<2>::generator(2, pauli3);
+  GaugeGroup<2, GroupName::SU>::generator(2, pauli3);
   pauli1 = timesI(pauli1) * 2.0;
   pauli2 = timesI(pauli2) * 2.0;
   pauli3 = timesI(pauli3) * 2.0;
@@ -389,8 +303,7 @@ public:
   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 *
+  u = 0.5 * u * pow(xi, -1.0);  //  u   = 1/2xi [ Sig - Sig^dag  + 1 Tr Sigdag]
-      pow(xi, -1.0);  //  u   = 1/2xi [ Sig - Sig^dag  + 1 Tr Sigdag]
 
   // Debug test for sanity
   uinv = adj(u);
@@ -405,29 +318,24 @@ public:
     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
+    Action factor Q(h) dh  = e^-S[h]  dh =  e^{  xi Tr uh} dh    // beta
-      through xi
+    enters through xi =  e^{2 xi (h.u)} dh =  e^{2 xi h0u0}.e^{2 xi h1u1}.e^{2
-      =  e^{2 xi (h.u)} dh
+    xi h2u2}.e^{2 xi h3u3} 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
+    Take alpha = 2 xi  = 2 xi [ recall 2 beta/Nc unmod staple norm];
-      factor in Chroma ]
+    hence 2.0/Nc factor in Chroma ] A. Generate two uniformly distributed
-      A. Generate two uniformly distributed pseudo-random numbers R and R', R'',
+    pseudo-random numbers R and R', R'', R''' in the unit interval; B. Set X =
-      R''' in the unit interval;
+    -(ln R)/alpha, X' =-(ln R')/alpha; C. Set C = cos^2(2pi R"), with R"
-      B. Set X = -(ln R)/alpha, X' =-(ln R')/alpha;
+    another uniform random number in [0,1] ; D. Set A = XC; E. Let d  = X'+A;
-      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
+    Note that in step D setting B ~ X - A and using B in place of A in step E
-      generate a second independent a 0 value.
+    will generate a second independent a 0 value.
   */
 
   /////////////////////////////////////////////////////////
@@ -558,19 +466,8 @@ public:
   /////////////////////////////////
 }
 
-  static void printGenerators(void) {
+template <ONLY_IF_SU>
-    for (int gen = 0; gen < AdjointDimension; gen++) {
+static void testGenerators(GroupName::SU) {
-      Matrix ta;
-      generator(gen, ta);
-      std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
-                << std::endl;
-      std::cout << GridLogMessage << ta << std::endl;
-    }
-  }
-
-
-
-  static void testGenerators(void) {
   Matrix ta;
   Matrix tb;
   std::cout << GridLogMessage
@@ -608,83 +505,10 @@ public:
   std::cout << GridLogMessage << std::endl;
 }
 
-  // reunitarise??
-  template <typename LatticeMatrixType>
-  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<vector_type> vTComplexType;
-
-    typedef Lattice<vTComplexType> 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;
-    }
-  }
 /*
  * Fundamental rep gauge xform
  */
-  template<typename Fundamental,typename GaugeMat>
+template <typename Fundamental, typename GaugeMat, ONLY_IF_SU>
 static void GaugeTransformFundamental(Fundamental &ferm, GaugeMat &g) {
   GridBase *grid = ferm._grid;
   conformable(grid, g._grid);
@@ -694,13 +518,14 @@ public:
  * Adjoint rep gauge xform
  */
 
-  template<typename GaugeField,typename GaugeMat>
+template <typename GaugeField, typename GaugeMat, ONLY_IF_SU>
 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 ag(grid);
+  ag = adj(g);
 
   for (int mu = 0; mu < Nd; mu++) {
     U = PeekIndex<LorentzIndex>(Umu, mu);
@@ -708,186 +533,18 @@ public:
     PokeIndex<LorentzIndex>(Umu, U, mu);
   }
 }
-  template<typename GaugeMat>
+template <typename GaugeMat, ONLY_IF_SU>
 static void GaugeTransform(std::vector<GaugeMat> &U, GaugeMat &g) {
   GridBase *grid = g.Grid();
-    GaugeMat ag(grid); ag = adj(g);
+  GaugeMat ag(grid);
+  ag = adj(g);
   for (int mu = 0; mu < Nd; mu++) {
     U[mu] = g * U[mu] * Cshift(ag, mu, 1);
   }
 }
-  template<typename GaugeField,typename GaugeMat>
+template <typename GaugeField, typename GaugeMat, ONLY_IF_SU>
-  static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, GaugeMat &g){
+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 <typename GaugeField>
-  static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
-    typedef typename GaugeField::vector_type vector_type;
-    typedef iSUnMatrix<vector_type> vMatrixType;
-    typedef Lattice<vMatrixType> LatticeMatrixType;
-
-    LatticeMatrixType Umu(out.Grid());
-    for (int mu = 0; mu < Nd; mu++) {
-      LieRandomize(pRNG, Umu, 1.0);
-      PokeIndex<LorentzIndex>(out, Umu, mu);
-    }
-  }
-  template<typename GaugeField>
-  static void TepidConfiguration(GridParallelRNG &pRNG,GaugeField &out){
-    typedef typename GaugeField::vector_type vector_type;
-    typedef iSUnMatrix<vector_type> vMatrixType;
-    typedef Lattice<vMatrixType> LatticeMatrixType;
-
-    LatticeMatrixType Umu(out.Grid());
-    for(int mu=0;mu<Nd;mu++){
-      LieRandomize(pRNG,Umu,0.01);
-      PokeIndex<LorentzIndex>(out,Umu,mu);
-    }
-  }
-  template<typename GaugeField>
-  static void ColdConfiguration(GaugeField &out){
-    typedef typename GaugeField::vector_type vector_type;
-    typedef iSUnMatrix<vector_type> vMatrixType;
-    typedef Lattice<vMatrixType> LatticeMatrixType;
-
-    LatticeMatrixType Umu(out.Grid());
-    Umu=1.0;
-    for(int mu=0;mu<Nd;mu++){
-      PokeIndex<LorentzIndex>(out,Umu,mu);
-    }
-  }
-  template<typename GaugeField>
-  static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){
-    ColdConfiguration(out);
-  }
-
-  template<typename LatticeMatrixType>
-  static void taProj( const LatticeMatrixType &in,  LatticeMatrixType &out){
-    out = Ta(in);
-  }
-  template <typename LatticeMatrixType>
-  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<int N>
-LatticeComplexD Determinant(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &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<iScalar<iMatrix<ComplexD, N> > > 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<int N>
-static void ProjectSUn(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
-{
-  Umu      = ProjectOnGroup(Umu);
-  auto det = Determinant(Umu);
-
-  det = conjugate(det);
-
-  for(int i=0;i<N;i++){
-    auto element = PeekIndex<ColourIndex>(Umu,N-1,i);
-    element = element * det;
-    PokeIndex<ColourIndex>(Umu,element,Nc-1,i);
-  }
-}
-template<int N>
-static void ProjectSUn(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >,Nd> > &U)
-{
-  GridBase *grid=U.Grid();
-  // Reunitarise
-  for(int mu=0;mu<Nd;mu++){
-    auto Umu = PeekIndex<LorentzIndex>(U,mu);
-    Umu      = ProjectOnGroup(Umu);
-    ProjectSUn(Umu);
-    PokeIndex<LorentzIndex>(U,Umu,mu);
-  }
-}
-// Explicit specialisation for SU(3).
-// Explicit specialisation for SU(3).
-static void
-ProjectSU3 (Lattice<iScalar<iScalar<iMatrix<vComplexD, 3> > > > &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<iVector<iScalar<iMatrix<vComplexD, 3> >,Nd> > &U)
-{
-  GridBase *grid=U.Grid();
-  // Reunitarise
-  for(int mu=0;mu<Nd;mu++){
-    auto Umu = PeekIndex<LorentzIndex>(U,mu);
-    Umu      = ProjectOnGroup(Umu);
-    ProjectSU3(Umu);
-    PokeIndex<LorentzIndex>(U,Umu,mu);
-  }
-}
-
-typedef SU<2> SU2;
-typedef SU<3> SU3;
-typedef SU<4> SU4;
-typedef SU<5> SU5;
-
-
-typedef SU<Nc> FundamentalMatrices;
-
-NAMESPACE_END(Grid);
-#endif

Grid/qcd/utils/Sp2n.h

@@ -1,54 +1,6 @@
 
-#ifndef QCD_UTIL_Sp2n_H
+private:
-#define QCD_UTIL_Sp2n_H
+static int su2subgroups(GroupName::Sp) { return (ncolour/2 * (ncolour/2 - 1)) / 2; }
-
-NAMESPACE_BEGIN(Grid);
-
-// Sp(2N)
-// ncolour = 2N
-
-
-template <int ncolour>
-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 <typename vtype>
-    using iSp2nMatrix = iScalar<iScalar<iMatrix<vtype, Dimension> > >;
-    template <typename vtype>
-    using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >;
-    template <typename vtype>
-    using iSp2nAlgebraVector = iScalar<iScalar<iVector<vtype, AlgebraDimension> > >;
-    
-    typedef iSp2nMatrix<Complex> Matrix;
-    typedef iSp2nMatrix<ComplexF> MatrixF;
-    typedef iSp2nMatrix<ComplexD> MatrixD;
-    
-    typedef iSp2nMatrix<vComplex> vMatrix;
-    typedef iSp2nMatrix<vComplexF> vMatrixF;
-    typedef iSp2nMatrix<vComplexD> vMatrixD;
-    
-    typedef iSp2nAlgebraVector<Complex> AlgebraVector;
-    typedef iSp2nAlgebraVector<ComplexF> AlgebraVectorF;
-    typedef iSp2nAlgebraVector<ComplexD> AlgebraVectorD;
-    
-    typedef iSp2nAlgebraVector<vComplex> vAlgebraVector;
-    typedef iSp2nAlgebraVector<vComplexF> vAlgebraVectorF;
-    typedef iSp2nAlgebraVector<vComplexD> vAlgebraVectorD;
-    
-    typedef Lattice<vMatrix> LatticeMatrix;
-    typedef Lattice<vMatrixF> LatticeMatrixF;
-    typedef Lattice<vMatrixD> LatticeMatrixD;
-    
-    typedef Lattice<vAlgebraVector> LatticeAlgebraVector;
-    typedef Lattice<vAlgebraVectorF> LatticeAlgebraVectorF;
-    typedef Lattice<vAlgebraVectorD> LatticeAlgebraVectorD;
-    
-    
 
 // Sp(2N) has N(2N+1) = 2N^2+N generators
 //
@@ -60,17 +12,19 @@ public:
 //     there are 6 types named a,b,c,d and w,z
 //     abcd are N(N-1)/2 each while wz are N each
 
- template <class cplx>
+template <class cplx, ONLY_IF_Sp>
-    static void generator(int lieIndex, iSp2nMatrix<cplx> &ta) {
+static void generator(int lieIndex, iGroupMatrix<cplx> &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
-        int mod = nsp * (nsp-1) * 0.5;
+  const int mod = nsp * (nsp - 1) * 0.5;
-        int offdiag = 2*nsp*nsp; // number of generators not in the cartan subalgebra
+  const int offdiag =
-        int wmod = 4*mod;
+      2 * nsp * nsp;  // number of generators not in the cartan subalgebra
-        int zmod = wmod+nsp;
+  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;
@@ -104,13 +58,15 @@ public:
     generatorBtype(bIndex, ta);
     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;
     cIndex = lieIndex - 2 * mod;
     generatorCtype(cIndex, ta);
     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;
     dIndex = lieIndex - 3 * mod;
     generatorDtype(dIndex, ta);
@@ -119,11 +75,11 @@ public:
 
 }  // end of generator
 
-    template <class cplx>
+template <class cplx, ONLY_IF_Sp>
-    static void generatorDiagtype(int diagIndex, iSp2nMatrix<cplx> &ta) {
+static void generatorDiagtype(int diagIndex, iGroupMatrix<cplx> &ta) {
-        
   // ta(i,i) = - ta(i+N,i+N) = 1/2 for each i index of the cartan subalgebra
 
+  const int nsp=ncolour/2;
   ta = Zero();
   RealD nrm = 1.0 / 2;
 
@@ -131,13 +87,13 @@ public:
   ta()()(diagIndex + nsp, diagIndex + nsp) = -nrm;
 }
 
-    template <class cplx>
+template <class cplx, ONLY_IF_Sp>
-    static void generatorAtype(int aIndex, iSp2nMatrix<cplx> &ta) {
+static void generatorAtype(int aIndex, iGroupMatrix<cplx> &ta) {
-        
   // ta(i,j) = ta(j,i) = -ta(i+N,j+N) = -ta(j+N,i+N) = 1 / 2 sqrt(2)
   // with i<j and i=0,...,N-2
   // follows that j=i+1, ... , N
   int i1, i2;
+  const int nsp=ncolour/2;
   ta = Zero();
   RealD nrm = 1 / (2 * std::sqrt(2));
 
@@ -150,20 +106,19 @@ public:
   ta = ta * nrm;
 }
 
-    template <class cplx>
+template <class cplx, ONLY_IF_Sp>
-    static void generatorBtype(int bIndex, iSp2nMatrix<cplx> &ta) {
+static void generatorBtype(int bIndex, iGroupMatrix<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
 
+  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()()(i1, i2) = i;
   ta()()(i2, i1) = -i;
   ta()()(i1 + nsp, i2 + nsp) = i;
@@ -172,12 +127,11 @@ public:
   ta = ta * nrm;
 }
 
-    template <class cplx>
+template <class cplx, ONLY_IF_Sp>
-    static void generatorCtype(int cIndex, iSp2nMatrix<cplx> &ta) {
+static void generatorCtype(int cIndex, iGroupMatrix<cplx> &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));
@@ -191,11 +145,11 @@ public:
   ta = ta * nrm;
 }
 
-    template <class cplx>
+template <class cplx, ONLY_IF_Sp>
-    static void generatorDtype(int dIndex, iSp2nMatrix<cplx> &ta) {
+static void generatorDtype(int dIndex, iGroupMatrix<cplx> &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);
@@ -210,11 +164,11 @@ public:
   ta = ta * nrm;
 }
 
-    template <class cplx>
+template <class cplx, ONLY_IF_Sp>
-     static void generatorWtype(int wIndex, iSp2nMatrix<cplx> &ta) {
+static void generatorWtype(int wIndex, iGroupMatrix<cplx> &ta) {
-         
   // ta(i,i+N) =  ta(i+N,i) = 1/2
 
+  const int nsp=ncolour/2;
   ta = Zero();
   RealD nrm = 1.0 / 2;  // check
 
@@ -224,11 +178,11 @@ public:
   ta = ta * nrm;
 }
 
-    template <class cplx>
+template <class cplx, ONLY_IF_Sp>
-     static void generatorZtype(int zIndex, iSp2nMatrix<cplx> &ta) {
+static void generatorZtype(int zIndex, iGroupMatrix<cplx> &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);
@@ -238,11 +192,12 @@ public:
   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) {
+template <ONLY_IF_Sp>
+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;
@@ -252,26 +207,12 @@ public:
   i2 = i1 + 1 + spare;
 }
 
-    
+static void testGenerators(GroupName::Sp) {
-    
-    
-    
-    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;
+  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);
@@ -281,180 +222,33 @@ public:
                 << 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;
+  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;
+  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);
   }
-        
 }
 
-    
+public:
-    template <typename LatticeMatrixType>
+template <ONLY_IF_Sp>
-    static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, double scale = 1.0)
+static void OmegaInvariance(ColourMatrix &in) {
-    {
-      GridBase *grid = out.Grid();
-
-      typedef typename LatticeMatrixType::vector_type vector_type;
-      typedef typename LatticeMatrixType::scalar_type scalar_type;
-
-      typedef iSinglet<vector_type> vTComplexType;
-
-      typedef Lattice<vTComplexType> 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 < AlgebraDimension; 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, //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;
-
-      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;
-
-      out = Zero();
-      for (int a = 0; a < AlgebraDimension; a++) {
-        generator(a, ta);
-        la = peekColour(h, a) * timesI(ta) * scale;
-        out += la;
-      }
-    }
-    
-    
-    template <typename LatticeMatrixType>
-    static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { // same as sun
-      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!....
-      }
-    }
-    
- 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 < AlgebraDimension; a++) {
-     generator(a, Ta);
-     pokeColour(h_out, - 2.0 * (trace(timesI(Ta) * in)) * scale, a);
-   }
- }
-    
-    
-    template <typename GaugeField>
-    static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
-      typedef typename GaugeField::vector_type vector_type;
-      typedef iSp2nMatrix<vector_type> vMatrixType;
-      typedef Lattice<vMatrixType> LatticeMatrixType;
-
-      LatticeMatrixType Umu(out.Grid());
-      for (int mu = 0; mu < Nd; mu++) {
-        LieRandomize(pRNG, Umu, 1.0);   //def
-        PokeIndex<LorentzIndex>(out, Umu, mu);
-      }
-    }
-    template<typename GaugeField>
-    static void TepidConfiguration(GridParallelRNG &pRNG,GaugeField &out){
-      typedef typename GaugeField::vector_type vector_type;
-      typedef iSp2nMatrix<vector_type> vMatrixType;
-      typedef Lattice<vMatrixType> LatticeMatrixType;
-
-      LatticeMatrixType Umu(out.Grid());
-      for(int mu=0;mu<Nd;mu++){
-        LieRandomize(pRNG,Umu,0.01);    //def
-        PokeIndex<LorentzIndex>(out,Umu,mu);
-      }
-    }
-    
-    template<typename GaugeField>
-    static void ColdConfiguration(GaugeField &out){
-      typedef typename GaugeField::vector_type vector_type;
-      typedef iSp2nMatrix<vector_type> vMatrixType;
-      typedef Lattice<vMatrixType> LatticeMatrixType;
-
-      LatticeMatrixType Umu(out.Grid());
-      Umu=1.0;
-      for(int mu=0;mu<Nd;mu++){
-        PokeIndex<LorentzIndex>(out,Umu,mu);
-      }
-    }
-    template<typename GaugeField>
-    static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){
-      ColdConfiguration(out);
-    }
-    
-    static void OmegaInvariance(ColourMatrix &in)
-    {
-        
   ColourMatrix Omega;
   Omega = Zero();
+  const int nsp=ncolour/2;
 
   std::cout << GridLogMessage << "I am a ColourMatrix" << std::endl;
 
@@ -463,28 +257,27 @@ public:
   //     Omega()()(i, 2*ncolour-1-i) = 1.;
   //     Omega()()(2*ncolour-1-i, i) = -1;
   // }
-        for (int i = 0; i < nsp; i++)
+  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)
+  if (norm2(sdiff) < 1e-8) {
-        {
+    std::cout << GridLogMessage
-            std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl;
+              << "Symplectic condition satisfied: Omega invariant" << std::endl;
   } 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, ONLY_IF_Sp>
-    
+static void OmegaInvariance(GaugeField &in) {
-    template<typename GaugeField>
-    static void OmegaInvariance(GaugeField &in)
-    {
   typedef typename GaugeField::vector_type vector_type;
-        typedef iSp2nMatrix<vector_type> vMatrixType;
+  typedef iGroupMatrix<vector_type> vMatrixType;
   typedef Lattice<vMatrixType> LatticeMatrixType;
 
   LatticeMatrixType U(in.Grid());
@@ -502,12 +295,11 @@ public:
   U = PeekIndex<LorentzIndex>(in, 1);
 
   OmegaInvariance(U);
-        
 }
 
-       static void OmegaInvariance(LatticeColourMatrixD &in)
+template <ONLY_IF_Sp>
-       {
+static void OmegaInvariance(LatticeColourMatrixD &in) {
-                      
+  const int nsp=ncolour/2;
   LatticeColourMatrixD OmegaLatt(in.Grid());
   LatticeColourMatrixD identity(in.Grid());
   RealD vol = in.Grid()->gSites();
@@ -519,8 +311,7 @@ public:
 
   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()()(nsp + i, i) = -1;
   }
@@ -531,50 +322,13 @@ public:
   auto diff = OmegaLatt - (in * OmegaLatt * transpose(in));
   auto sdiff = sum(diff);
   // assert( norm2(sdiff) < 1e-8 );
-           if (norm2(sdiff) < 1e-8)
+  if (norm2(sdiff) < 1e-8) {
-           {
+    std::cout << GridLogMessage
-               std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl;
+              << "Symplectic condition satisfied: Omega invariant" << std::endl;
   } 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
-    
- 
- template<int N>
- static void ProjectSp2n(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
- {
-   Umu      = ProjectOnSpGroup(Umu);
-   auto det = Determinant(Umu); // ok ?
-
-   det = conjugate(det);
-
-   for(int i=0;i<N;i++){
-     auto element = PeekIndex<ColourIndex>(Umu,N-1,i);
-     element = element * det;
-     PokeIndex<ColourIndex>(Umu,element,Nc-1,i);
-   }
- }
- template<int N>
- static void ProjectSp2n(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >,Nd> > &U)
- {
-   GridBase *grid=U.Grid();
-   for(int mu=0;mu<Nd;mu++){
-     auto Umu = PeekIndex<LorentzIndex>(U,mu);
-     Umu      = ProjectOnSpGroup(Umu);
-     ProjectSp2n(Umu);
-     PokeIndex<LorentzIndex>(U,Umu,mu);
   }
 }
 
-typedef Sp<2> Sp2;
-typedef Sp<4> Sp4;
-typedef Sp<6> Sp6;
-typedef Sp<8> Sp8;
-
-NAMESPACE_END(Grid);
-#endif

Grid/qcd/utils/Sp2nTwoIndex.h

@@ -175,10 +175,10 @@ public:
 
   template <class cplx>
   static void generator(int Index, iSp2nTwoIndexMatrix<cplx> &i2indTa) {
-    Vector<typename Sp<ncolour>::template iSp2nMatrix<cplx> > ta(
+    Vector<iSp2nMatrix<cplx> > ta(
 								NumGenerators);
-    Vector<typename Sp<ncolour>::template iSp2nMatrix<cplx> > eij(Dimension);
+    Vector<iSp2nMatrix<cplx> > eij(Dimension);
-    typename Sp<ncolour>::template iSp2nMatrix<cplx> tmp;
+    iSp2nMatrix<cplx> 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<ncolour>::template iSp2nMatrix<cplx> tmp1 =
+        iSp2nMatrix<cplx> tmp1 =
 	  tmp * eij[b]; 
         Complex iTr = TensorRemove(timesI(trace(tmp1)));
         i2indTa()()(a, b) = iTr;

Grid/qcd/utils/Utils.h

@@ -8,8 +8,7 @@
 #include <Grid/qcd/utils/ScalarObjs.h>
 
 // Include representations
-#include <Grid/qcd/utils/SUn.h>
+#include <Grid/qcd/utils/GaugeGroup.h>
-#include <Grid/qcd/utils/Sp2n.h>
 #include <Grid/qcd/utils/SUnAdjoint.h>
 #include <Grid/qcd/utils/SUnTwoIndex.h>
 #include <Grid/qcd/utils/Sp2nTwoIndex.h>

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

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`

tests/core/Test_lie_generators.cc

@@ -32,7 +32,7 @@ directory
 
 #include <Grid/qcd/utils/CovariantCshift.h>
 
-#include <Grid/qcd/utils/SUn.h>
+#include <Grid/qcd/utils/GaugeGroup.h>
 #include <Grid/qcd/utils/SUnAdjoint.h>
 #include <Grid/qcd/utils/SUnTwoIndex.h>
 

tests/core/test_sp2n_lie_gen.cc

@@ -1,55 +0,0 @@
-
-
-
-#include <Grid/Grid.h>
-#include <Grid/qcd/utils/Sp2n.h>
-
-using namespace Grid;
-
-int main(int argc, char** argv) {
-  Grid_init(&argc, &argv);
-    
-  //std::vector<int> 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();
-}

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

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<Nc>::nsp;
+    const int nsp = Nc/2;
     identity = 1.;
     Cidentity = 1.;
     

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 = Nc / 2;
-    const int nsp = Sp<Nc>::nsp;
     
     identity = 1.0;
     RealD epsilon = 0.01;

tests/sp2n/Test_sp2n_lie_gen.cc

@@ -1,8 +1,4 @@
-
-
-
 #include <Grid/Grid.h>
-#include <Grid/qcd/utils/Sp2n.h>
 
 using namespace Grid;