From fd057c838f4788f809b49ce137b33acd6b3a4377 Mon Sep 17 00:00:00 2001
From: Alessandro Lupo <Alessandro.Lupo@ed.ac.uk>
Date: Fri, 10 Mar 2023 12:10:46 +0000
Subject: [PATCH] add ProjectOnGaugeGroup and ProjectGn to allow future
 templating in GaugeImplTypes

---
 Grid/qcd/utils/GaugeGroup.h | 27 ++++++++++++++++++++++++
 Grid/qcd/utils/SUn.impl     | 41 +++++++++++++++++++++++++++++++++++++
 Grid/qcd/utils/Sp2n.impl    | 41 +++++++++++++++++++++++++++++++++++++
 3 files changed, 109 insertions(+)

diff --git a/Grid/qcd/utils/GaugeGroup.h b/Grid/qcd/utils/GaugeGroup.h
index f0442e2f..dea8f985 100644
--- a/Grid/qcd/utils/GaugeGroup.h
+++ b/Grid/qcd/utils/GaugeGroup.h
@@ -244,6 +244,29 @@ class GaugeGroup {
       ex = ex + xn * nfac;     // x2/2!, x3/3!....
     }
   }
+    
+  template <int N>        // reunitarise, resimplectify...
+  static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
+      ProjectOnGaugeGroup(U, group_name());
+  }
+       
+  template <int N>        // reunitarise, resimplectify...
+  static void ProjectOnGaugeGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
+      ProjectOnGaugeGroup(Umu, group_name());
+  }
+    
+    
+  template <int N>        // reunitarise, resimplectify... previously ProjectSUn
+  static void ProjectGn(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
+      ProjectGn(Umu, group_name());
+  }
+
+  template <int N>        // reunitarise, resimplectify... previously ProjectSUn
+  static void ProjectGn(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
+      ProjectGn(U, group_name());
+      
+  }
+
 };
 
 template <int N>
@@ -271,6 +294,8 @@ LatticeComplexD Determinant(
   });
   return ret;
 }
+
+
 template <int N>
 static void ProjectSUn(
     Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
@@ -285,6 +310,7 @@ static void ProjectSUn(
     PokeIndex<ColourIndex>(Umu, element, Nc - 1, i);
   }
 }
+
 template <int N>
 static void ProjectSUn(
     Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
@@ -297,6 +323,7 @@ static void ProjectSUn(
     PokeIndex<LorentzIndex>(U, Umu, mu);
   }
 }
+
 // Explicit specialisation for SU(3).
 // Explicit specialisation for SU(3).
 static void ProjectSU3(
diff --git a/Grid/qcd/utils/SUn.impl b/Grid/qcd/utils/SUn.impl
index eef039c7..18be96b1 100644
--- a/Grid/qcd/utils/SUn.impl
+++ b/Grid/qcd/utils/SUn.impl
@@ -621,6 +621,47 @@ static void projectOnAlgebra(LatticeAlgebraVector &h_out,
   }
 }
 
+template <int N>
+static void ProjectOnGaugeGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::SU) {
+  Umu = ProjectOnGroup(Umu);
+ }
+
+template <int N>
+static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U, GroupName::SU) {
+  GridBase *grid = U.Grid();
+  // Reunitarise
+  for (int mu = 0; mu < Nd; mu++) {
+    auto Umu = PeekIndex<LorentzIndex>(U, mu);
+    Umu = ProjectOnGroup(Umu);
+  }
+}
+
+template <int N>
+static void ProjectGn(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::SU) {
+  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 ProjectGn(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U, GroupName::SU) {
+  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);
+  }
+}
+
 /*
  * Fundamental rep gauge xform
  */
diff --git a/Grid/qcd/utils/Sp2n.impl b/Grid/qcd/utils/Sp2n.impl
index 411235ec..44771a14 100644
--- a/Grid/qcd/utils/Sp2n.impl
+++ b/Grid/qcd/utils/Sp2n.impl
@@ -366,6 +366,47 @@ static void TepidConfiguration(GridParallelRNG &pRNG, GaugeField &out, GroupName
   }
 }
 
+template <int N>
+static void ProjectOnGaugeGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::Sp) {
+  Umu = ProjectOnSpGroup(Umu);
+ }
+
+template <int N>
+static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U, GroupName::Sp) {
+  GridBase *grid = U.Grid();
+  // Reunitarise
+  for (int mu = 0; mu < Nd; mu++) {
+    auto Umu = PeekIndex<LorentzIndex>(U, mu);
+    Umu = ProjectOnSpGroup(Umu);
+  }
+}
+
+template <int N>
+static void ProjectGn(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::Sp) {
+  Umu = ProjectOnSpGroup(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 ProjectGn(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U, GroupName::Sp) {
+  GridBase *grid = U.Grid();
+  // Reunitarise
+  for (int mu = 0; mu < Nd; mu++) {
+    auto Umu = PeekIndex<LorentzIndex>(U, mu);
+    Umu = ProjectOnSpGroup(Umu);
+    ProjectSUn(Umu);
+    PokeIndex<LorentzIndex>(U, Umu, mu);
+  }
+}
+
 public:
 template <ONLY_IF_Sp>
 static void OmegaInvariance(ColourMatrix &in) {