From 026e736dfabdb256056a5fc6dc31d443556a5bfc Mon Sep 17 00:00:00 2001
From: Alessandro Lupo <Alessandro.Lupo@ed.ac.uk>
Date: Mon, 3 Apr 2023 16:31:19 +0100
Subject: [PATCH] Projection on algebra can now be templated. Fix #12

---
 Grid/lattice/Lattice_ET.h        |   2 +
 Grid/qcd/utils/Sp2n.impl         |  24 +++-
 Grid/tensors/Tensor_Ta.h         |  80 ++++++++++++
 tests/sp2n/Test_project_on_Sp.cc | 211 +++++++++++++++++++++++--------
 4 files changed, 261 insertions(+), 56 deletions(-)

diff --git a/Grid/lattice/Lattice_ET.h b/Grid/lattice/Lattice_ET.h
index 9042eb8f..661d7441 100644
--- a/Grid/lattice/Lattice_ET.h
+++ b/Grid/lattice/Lattice_ET.h
@@ -345,6 +345,7 @@ GridUnopClass(UnaryNot, Not(a));
 GridUnopClass(UnaryTrace, trace(a));
 GridUnopClass(UnaryTranspose, transpose(a));
 GridUnopClass(UnaryTa, Ta(a));
+GridUnopClass(UnaryProjectSp2nAlgebra, ProjectSp2nAlgebra(a));
 GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a));
 GridUnopClass(UnaryProjectOnSpGroup, ProjectOnSpGroup(a));
 GridUnopClass(UnaryTimesI, timesI(a));
@@ -457,6 +458,7 @@ GRID_DEF_UNOP(operator!, UnaryNot);
 GRID_DEF_UNOP(trace, UnaryTrace);
 GRID_DEF_UNOP(transpose, UnaryTranspose);
 GRID_DEF_UNOP(Ta, UnaryTa);
+GRID_DEF_UNOP(ProjectSp2nAlgebra, UnaryProjectSp2nAlgebra);
 GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup);
 GRID_DEF_UNOP(ProjectOnSpGroup, UnaryProjectOnSpGroup);
 GRID_DEF_UNOP(timesI, UnaryTimesI);
diff --git a/Grid/qcd/utils/Sp2n.impl b/Grid/qcd/utils/Sp2n.impl
index 77da9ad0..347eb3ce 100644
--- a/Grid/qcd/utils/Sp2n.impl
+++ b/Grid/qcd/utils/Sp2n.impl
@@ -274,8 +274,6 @@ static void OmegaInvariance(ColourMatrix &in) {
   Omega = Zero();
   const int nsp=ncolour/2;
 
-  std::cout << GridLogMessage << "I am a ColourMatrix" << std::endl;
-
   // for (int i = 0; i < ncolour; i++)         wrong?!
   //{
   //     Omega()()(i, 2*ncolour-1-i) = 1.;
@@ -314,8 +312,6 @@ static void OmegaInvariance(GaugeField &in) {
   Omega = Zero();
   identity = 1.;
 
-  std::cout << GridLogMessage << "I am a GaugeField " << std::endl;
-
   U = PeekIndex<LorentzIndex>(in, 1);
 
   OmegaInvariance(U);
@@ -333,8 +329,6 @@ static void OmegaInvariance(LatticeColourMatrixD &in) {
   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;
@@ -356,3 +350,21 @@ static void OmegaInvariance(LatticeColourMatrixD &in) {
   }
 }
 
+template <ONLY_IF_Sp>
+static void Omega(LatticeColourMatrixD &in) {
+  const int nsp=ncolour/2;
+  LatticeColourMatrixD OmegaLatt(in.Grid());
+  LatticeColourMatrixD identity(in.Grid());
+  ColourMatrix Omega;
+
+  OmegaLatt = Zero();
+  Omega = Zero();
+  identity = 1.;
+
+  for (int i = 0; i < nsp; i++) {
+    Omega()()(i, nsp + i) = 1.;
+    Omega()()(nsp + i, i) = -1;
+  }
+  OmegaLatt = OmegaLatt + (identity * Omega);
+  in = OmegaLatt;
+}
diff --git a/Grid/tensors/Tensor_Ta.h b/Grid/tensors/Tensor_Ta.h
index 2ee51fb4..87169ffc 100644
--- a/Grid/tensors/Tensor_Ta.h
+++ b/Grid/tensors/Tensor_Ta.h
@@ -66,6 +66,86 @@ template<class vtype,int N> accelerator_inline iMatrix<vtype,N> Ta(const iMatrix
   return ret;
 }
 
+// for sp2n can't do something as simple as Ta. We do a Gram-Schmidt
+
+template<class vtype> accelerator_inline iScalar<vtype> ProjectSp2nAlgebra(const iScalar<vtype>&r)
+{
+  iScalar<vtype> ret;
+  ret._internal = ProjectSp2nAlgebra(r._internal);
+  return ret;
+}
+template<class vtype,int N> accelerator_inline iVector<vtype,N> ProjectSp2nAlgebra(const iVector<vtype,N>&r)
+{
+  iVector<vtype,N> ret;
+  for(int i=0;i<N;i++){
+    ret._internal[i] = ProjectSp2nAlgebra(r._internal[i]);
+  }
+  return ret;
+}
+template<class vtype,int N> accelerator_inline iMatrix<vtype,N> ProjectSp2nAlgebra(const iMatrix<vtype,N> &arg)
+{
+    iMatrix<vtype,N> ret;
+    vtype nrm;
+    vtype inner;
+    vtype tmp;
+    
+    for(int c1=0;c1<N/2;c1++)
+    {
+        
+        for (int b=0; b<c1; b++)                  // remove the b-rows from U_c1
+        {
+            decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
+            decltype(ret._internal[b][b]*ret._internal[b][b]) prn;
+            zeroit(pr);
+            zeroit(prn);
+            
+            for(int c=0; c<N; c++)
+            {
+                pr += conjugate(ret._internal[c1][c])*ret._internal[b][c];        // <U_c1 | U_b >
+                prn += conjugate(ret._internal[c1][c])*ret._internal[b+N/2][c];   // <U_c1 | U_{b+N} >
+            }
+         
+
+            for(int c=0; c<N; c++)
+            {
+                ret._internal[c1][c] -= (conjugate(pr) * ret._internal[b][c] + conjugate(prn) * ret._internal[b+N/2][c] );    //  U_c1 -= (  <U_c1 | U_b > U_b + <U_c1 | U_{b+N} > U_{b+N}  )
+            }
+        }
+      
+        zeroit(inner);
+        for(int c2=0;c2<N;c2++)
+        {
+            inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
+        }
+        
+        nrm = sqrt(inner);
+        nrm = 1.0/nrm;
+        for(int c2=0;c2<N;c2++)
+        {
+            ret._internal[c1][c2]*= nrm;
+        }
+        
+
+        for(int c2=0;c2<N/2;c2++)
+        {
+            tmp = conjugate(ret._internal[c1][c2]);       // (up-left)* of the old matrix
+            ret._internal[c1+N/2][c2+N/2] = -tmp;          // down right in the new matrix = -(up-left)* of the old matrix
+        }
+           
+        for(int c2=N/2;c2<N;c2++)
+        {
+            tmp = conjugate(ret._internal[c1][c2]);   // (up-right)* of the old
+            ret._internal[c1+N/2][c2-N/2] = tmp;     // down left in the new matrix = (up-right)* of the old
+        }
+      
+        
+    }
+    
+  
+  return ret;
+}
+
+
 
 /////////////////////////////////////////////// 
 // ProjectOnGroup function for scalar, vector, matrix 
diff --git a/tests/sp2n/Test_project_on_Sp.cc b/tests/sp2n/Test_project_on_Sp.cc
index 36b95cf4..eefee4f8 100644
--- a/tests/sp2n/Test_project_on_Sp.cc
+++ b/tests/sp2n/Test_project_on_Sp.cc
@@ -20,6 +20,12 @@ int main (int argc, char **argv)
     LatticeColourMatrixD aux(&Grid);
     LatticeColourMatrixD identity(&Grid);
     
+    // Will test resimplectification-related functionalities (from ProjectOnGaugeGroup, ProjectOnSpGroup, ProjectGn) and projection on the algebra (from ProjectSp2nAlgebra)
+    // we work with matrices with positive determinant so detU = 1 even if in principle ProjectOnGaugeGroup and ProjectOnSpGroup allow for detU=-1
+    // so the checks will be the same for the three functions
+    // NB only ProjectGn is the proper simplectification function
+    
+    
     const int nsp = Nc / 2;
     
     identity = 1.0;
@@ -29,15 +35,12 @@ int main (int argc, char **argv)
     double vol = Umu.Grid()->gSites();
     
     std::vector<int> pseeds({1,2,3,4,5});
-    std::vector<int> sseeds({6,7,8,9,10});
     GridParallelRNG  pRNG(&Grid); pRNG.SeedFixedIntegers(pseeds);
-    GridSerialRNG    sRNG;       sRNG.SeedFixedIntegers(sseeds);
     
     SU<Nc>::HotConfiguration(pRNG,Umu);
-    U = PeekIndex<LorentzIndex>(Umu,2);
+    U = PeekIndex<LorentzIndex>(Umu,0);
     
     aux = U*adj(U) - identity;
-    std::cout <<GridLogMessage << std::endl;
     std::cout << GridLogMessage << "Starting with random SUn matrix " << std::endl;
     std::cout << GridLogMessage << "Unitary check " << std::endl;
     std::cout <<GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
@@ -55,35 +58,23 @@ int main (int argc, char **argv)
     
     std::cout << GridLogMessage << "Unitary matrix deformed " << std::endl;
     std::cout << GridLogMessage << "now U adjU - 1 = " << norm2(aux) << std::endl;
-    std::cout <<GridLogMessage << std::endl;
-
-    std::cout <<GridLogMessage << std::endl;
-    std::cout << GridLogMessage << "Projecting on Sp2n " << std::endl;
-
+    
+    //  Testing ProjectOnSpGroup
+    std::cout << GridLogMessage << "Testing ProjectOnSpGroup" << std::endl;
+    std::cout << GridLogMessage << "Apply ProjectOnSpGroup to deformed matrix" << std::endl;
     U = ProjectOnSpGroup(U);
-    //U = ProjectOnGroup(U);
+
     aux = U*adj(U) - identity;
-    std::cout <<GridLogMessage << std::endl;
-    std::cout <<GridLogMessage << std::endl;
-    std::cout << GridLogMessage << "Unitary check after Sp(2n) projection " << std::endl;
+    std::cout << GridLogMessage << "Unitary check after ProjectOnSpGroup " << std::endl;
     std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
     assert( norm2(aux) < 1e-8);
-    std::cout <<GridLogMessage << std::endl;
-    std::cout <<GridLogMessage << std::endl;
-    
-    // checks on determinant
-    std::cout << GridLogMessage << "Det after Projection on Sp2n = " << norm2( Determinant(U) ) / vol << std::endl;
-    std::cout <<GridLogMessage << std::endl;
-    std::cout <<GridLogMessage << std::endl;
     
     // actual sp2n check
-    std::cout << GridLogMessage << "Checking invariance after projection "<< std::endl;
-    Sp<Nc>::OmegaInvariance(U);
+    std::cout << GridLogMessage << "Checking Omega invariance after ProjectOnSpGroup" << std::endl;
+    Sp<Nc>::OmegaInvariance(U);     // no assertion here, but the next check will kill us if we are not simplectic
     
     // checks on elements
     
-    std::cout <<GridLogMessage << std::endl;
-    std::cout <<GridLogMessage << std::endl;
     std::cout << GridLogMessage << "Checking the structure is " << std::endl;
     std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
     std::cout << GridLogMessage << "      (  -X^*  W^* )  " << std::endl;
@@ -119,38 +110,158 @@ int main (int argc, char **argv)
         }
     }
     
-    std::cout << GridLogMessage << "ok" << std::endl;
     
-    // an explicit check for sp2
-    /*
-    if (Nc == 2)
+    std::cout << GridLogMessage << "Testing ProjectOnGaugeGroup" << std::endl;
+    U = U + 2932.111*identity;
+    std::cout << GridLogMessage << "Apply ProjectOnGaugeGroup to deformed matrix" << std::endl;
+    Sp<Nc>::ProjectOnGaugeGroup(U);
+    aux = U*adj(U) - identity;
+    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
+    assert( norm2(aux) < 1e-8);
+    Sp<Nc>::OmegaInvariance(U);
+    
+    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
+    std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
+    std::cout << GridLogMessage << "      (  -X^*  W^* )  " << std::endl;
+    std::cout <<GridLogMessage << std::endl;
+    for (int c1 = 0; c1 < nsp; c1++) //check on W
     {
-        assert(Nc==2);
-        ColourMatrix A;
-        A = Zero();
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto W = PeekIndex<ColourIndex>(U,c1,c2);
+            auto Wstar =  PeekIndex<ColourIndex>(U,c1+nsp,c2+nsp);
+            auto Ww = conjugate( Wstar );
+            auto amizero = sum(W - Ww);
+            auto amizeroo = TensorRemove(amizero);
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+        }
+        
+    }
 
-        Complex a(25041994., 12.);
-        Complex b(39., 0.22);
-        Complex d(10000., -2222.3333);
+    for (int c1 = 0; c1 < nsp ; c1++)
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto X = PeekIndex<ColourIndex>(U,c1,c2+nsp);
+            auto minusXstar = PeekIndex<ColourIndex>(U,c1+nsp,c2);
+            auto minusXx = conjugate(minusXstar);
+            auto amizero = sum (X + minusXx);
+            auto amizeroo = TensorRemove(amizero);
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+        }
+    }
+    
+    std::cout << GridLogMessage << "Testing ProjectGn" << std::endl;
+    U = U + 2932.111*identity;
+    std::cout << GridLogMessage << "Apply ProjectGn to deformed matrix" << std::endl;
+    Sp<Nc>::ProjectGn(U);
+    aux = U*adj(U) - identity;
+    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
+    assert( norm2(aux) < 1e-8);
+    std::cout << GridLogMessage << "Det after ProjectGn = " << norm2( Determinant(U) ) / vol << std::endl;
+    assert( norm2(aux) - 1 < 1e-8);
+    Sp<Nc>::OmegaInvariance(U);
+    
+    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
+    std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
+    std::cout << GridLogMessage << "      (  -X^*  W^* )  " << std::endl;
+    std::cout <<GridLogMessage << std::endl;
+    for (int c1 = 0; c1 < nsp; c1++) //check on W
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto W = PeekIndex<ColourIndex>(U,c1,c2);
+            auto Wstar =  PeekIndex<ColourIndex>(U,c1+nsp,c2+nsp);
+            auto Ww = conjugate( Wstar );
+            auto amizero = sum(W - Ww);
+            auto amizeroo = TensorRemove(amizero);
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+        }
+        
+    }
+
+    for (int c1 = 0; c1 < nsp ; c1++)
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto X = PeekIndex<ColourIndex>(U,c1,c2+nsp);
+            auto minusXstar = PeekIndex<ColourIndex>(U,c1+nsp,c2);
+            auto minusXx = conjugate(minusXstar);
+            auto amizero = sum (X + minusXx);
+            auto amizeroo = TensorRemove(amizero);
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+        }
+    }
+    
+    std::cout << GridLogMessage << "Testing ProjectSp2nAlgebra" << std::endl;
+    
+    U = PeekIndex<LorentzIndex>(Umu,1);
+    U = U + 666.*identity;
+    aux = U*adj(U) - identity;
+    std::cout << GridLogMessage << "Matrix deformed " << std::endl;
+    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
+    std::cout << GridLogMessage << "Apply ProjectSp2nAlgebra to deformed matrix" << std::endl;
+    U = ProjectSp2nAlgebra(U);
+    aux = U*adj(U) - identity;
+    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
+    assert( norm2(aux) < 1e-8);
+    std::cout << GridLogMessage << "Check that Omega U Omega = conj(U)" << std::endl;
+    
+    LatticeColourMatrixD Omega(&Grid);
+    
+    Sp<Nc>::Omega(Omega);
+    aux = Omega*U*Omega - conjugate(U);
+    std::cout << GridLogMessage << "Omega U Omega - conj(U) = " << norm2(aux) << std::endl;
+    assert( norm2(aux) < 1e-8);
+    
+    
+    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
+    std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
+    std::cout << GridLogMessage << "      (  X^*  -W^* )  " << std::endl;
+    std::cout <<GridLogMessage << std::endl;
+    for (int c1 = 0; c1 < nsp; c1++) //check on W
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto W = PeekIndex<ColourIndex>(U,c1,c2);
+            auto Wstar =  PeekIndex<ColourIndex>(U,c1+nsp,c2+nsp);
+            auto Ww = conjugate( Wstar );
+            auto amizero = sum(W + Ww);
+            auto amizeroo = TensorRemove(amizero);
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+        }
+        
+    }
+
+    for (int c1 = 0; c1 < nsp ; c1++)
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto X = PeekIndex<ColourIndex>(U,c1,c2+nsp);
+            auto minusXstar = PeekIndex<ColourIndex>(U,c1+nsp,c2);
+            auto minusXx = conjugate(minusXstar);
+            auto amizero = sum (X - minusXx);
+            auto amizeroo = TensorRemove(amizero);
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+        }
+    }
+    
     
-        A()()(0,0) = a;
-        A()()(0,1) = b;
-        A()()(1,0) = i;
-        A()()(1,1) = d;
-        std::cout <<GridLogMessage << std::endl;
-        std::cout <<GridLogMessage << std::endl;
-        std::cout << GridLogMessage << "An explicit check for Sp2" << std::endl;
-        std::cout <<GridLogMessage << std::endl;
-        std::cout << GridLogMessage << "Building a non unitary matrix by hand with funny entries " << std::endl;
-        std::cout << GridLogMessage << "A = " << A << std::endl;
-        std::cout << GridLogMessage << "Projecting on Sp2 " << std::endl;
-        A = ProjectOnSpGroup(A);
-        std::cout << GridLogMessage << "now A = " << A << std::endl;
-        std::cout << GridLogMessage << "A(0,0) - conjA(1,1) = " << A()()(0,0) - adj ( A()()(1,1) )<< std::endl;
-        std::cout << GridLogMessage << "A(0,1) + conjA(1,0) = " << A()()(0,1) + adj ( A()()(1,0) )<< std::endl;
-    }*/
     
     Grid_finalize();
 
 
 }
+