From 283f14b7c1151538a5558bfb45fa3b18b52fef5b Mon Sep 17 00:00:00 2001
From: Alessandro Lupo <Alessandro.Lupo@ed.ac.uk>
Date: Tue, 12 Oct 2021 12:12:16 +0100
Subject: [PATCH] fix sp2n projection

---
 Grid/tensors/Tensor_Ta.h         | 133 ++++++++++++++-----------------
 tests/sp2n/Test_project_on_Sp.cc | 122 +++++++++++++++++++++++++++-
 2 files changed, 181 insertions(+), 74 deletions(-)

diff --git a/Grid/tensors/Tensor_Ta.h b/Grid/tensors/Tensor_Ta.h
index f53296bc..7b4bccca 100644
--- a/Grid/tensors/Tensor_Ta.h
+++ b/Grid/tensors/Tensor_Ta.h
@@ -161,88 +161,75 @@ accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg
   iMatrix<vtype,N> ret(arg);
   vtype nrm;
   vtype inner;
-  for(int c1=0;c1<N;c1++){
+  for(int c1=0;c1<N/2;c1++){
 
-       // Normalises row c1
-        zeroit(inner);
-        for(int c2=0;c2<N;c2++)
-          inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
+    // Normalises row c1
+    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;
+    nrm = sqrt(inner);
+    nrm = 1.0/nrm;
+    for(int c2=0;c2<N;c2++)
+      ret._internal[c1][c2]*= nrm;
 
-        // Compute row c1+N2/2: c1+N/2 = - \Omega c1*
-        for(int c2=0;c2<N/2;c2++)
-          //tmp = conjugate(ret._internal[c1][c2]);
-          //ret._internal[c1+N/2][c2+N/2] = tmp;
-            ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]);
+    // Compute row c1+N2/2: c1+N/2 = - \Omega c1*
+    for(int c2=0;c2<N/2;c2++)
+      ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]);
 
-        for(int c2=N/2;c2<N;c2++)
-          //tmp = conjugate(ret._internal[c1][c2]);
-         // ret._internal[c1+N/2][c2-N/2] = -tmp;
-            ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);
+    for(int c2=N/2;c2<N;c2++)
+      ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);
 
-        // Remove c1 from rows c1+1...N/2-1
-        for (int b=c1+1; b<N/2; ++b){
-          decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
-          zeroit(pr);
-          for(int c=0; c<N; ++c)
-            pr += conjugate(ret._internal[c1][c])*ret._internal[b][c];
+    // Remove c1 from rows c1+1...N/2-1
+    for (int b=c1+1; b<N/2; ++b){
+      decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
+      zeroit(pr);
+      for(int c=0; c<N; ++c)
+        pr += conjugate(ret._internal[c1][c])*ret._internal[b][c];
 
-          for(int c=0; c<N; ++c){
-            ret._internal[b][c] -= pr * ret._internal[c1][c];
-          }
-        }
-
-        // Remove c1+N/2 from rows c1+1...N/2-1
-        for (int b=c1+1; b<N; ++b){
-          decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
-          zeroit(pr);
-          for(int c=0; c<N; ++c)
-            pr += conjugate(ret._internal[c1+N/2][c])*ret._internal[b][c];
-
-          for(int c=0; c<N; ++c){
-            ret._internal[b][c] -= pr * ret._internal[c1+N/2][c];
-          }
-        }
-      } //end for loop over c1
-
-      // Compute the last row
-      {
-      // Normalise last row
-        int c1 = N/2-1;
-        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;
-        }
-      // Compute row cN-1+N2/2: cN-1+N/2 = - \Omega cN-1*
-        for(int c2=0;c2<N/2;c2++)
-        {
-          //tmp = conjugate(ret._internal[c1][c2]);
-          // ret._internal[c1+N/2][c2+N/2] = tmp;
-             ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]);
-        }
-        for(int c2=N/2;c2<N;c2++)
-        {
-          //tmp = conjugate(ret._internal[c1][c2]);
-          //ret._internal[c1+N/2][c2-N/2] = -tmp;
-            ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);
-        }
+      for(int c=0; c<N; ++c){
+        ret._internal[b][c] -= pr * ret._internal[c1][c];
       }
-      // assuming the determinant is ok
-      return ret;
+    }
+
+    // Remove c1+N/2 from rows c1+1...N/2-1
+    for (int b=c1+1; b<N; ++b){
+      decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
+      zeroit(pr);
+      for(int c=0; c<N; ++c)
+        pr += conjugate(ret._internal[c1+N/2][c])*ret._internal[b][c];
+
+      for(int c=0; c<N; ++c){
+        ret._internal[b][c] -= pr * ret._internal[c1+N/2][c];
+      }
+    }
   }
 
+  // Compute the last row
+  {
+  // Normalise last row
+    int c1 = N/2-1;
+    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;
+
+  // Compute row cN-1+N2/2: cN-1+N/2 = - \Omega cN-1*
+    for(int c2=0;c2<N/2;c2++)
+      ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]);
+
+    for(int c2=N/2;c2<N;c2++)
+      ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);
+  }
+  // assuming the determinant is ok
+  return ret;
+}
+
+
 NAMESPACE_END(Grid);
 
 #endif
diff --git a/tests/sp2n/Test_project_on_Sp.cc b/tests/sp2n/Test_project_on_Sp.cc
index 61ec6d0c..d494d828 100644
--- a/tests/sp2n/Test_project_on_Sp.cc
+++ b/tests/sp2n/Test_project_on_Sp.cc
@@ -17,6 +17,14 @@ int main (int argc, char **argv)
     
     LatticeGaugeField Umu(&Grid);
     LatticeColourMatrixD U(&Grid);
+    LatticeColourMatrixD aux(&Grid);
+    LatticeColourMatrixD identity(&Grid);
+    
+    identity = 1.0;
+    RealD epsilon = 0.01;
+    Complex i(0., 1.);
+    RealD u = 0.;
+    double vol = Umu.Grid()->gSites();
     
     std::vector<int> pseeds({1,2,3,4,5});
     std::vector<int> sseeds({6,7,8,9,10});
@@ -24,10 +32,122 @@ int main (int argc, char **argv)
     GridSerialRNG    sRNG;       sRNG.SeedFixedIntegers(sseeds);
     
     SU<Nc>::HotConfiguration(pRNG,Umu);
-    U = PeekIndex<LorentzIndex>(Umu,0);
+    U = PeekIndex<LorentzIndex>(Umu,2);
     
+    // it is unitary
+    aux = U*adj(U) - identity;
+    std::cout <<GridLogMessage << std::endl;
+    std::cout << GridLogMessage << "Unitary check of starting random SUn matrix " << std::endl;
+    std::cout <<GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
+    std::cout <<GridLogMessage << std::endl;
+    std::cout <<GridLogMessage << std::endl;
     
+    U = U + epsilon*identity;
+    aux = U*adj(U) - identity;
+    
+    std::cout << GridLogMessage << "Unitary matrix modified " << 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;
+
     U = ProjectOnSpGroup(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 << "U adjU - 1 = " << norm2(aux) << std::endl;
+    assert( norm2(aux) < 1e-6);
+    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;
+    // checks on elements
+    // an Sp2n matrix is a block matrix of the type
+    //
+    //
+    //     U = (    W       X   )
+    //         (    -X^*    W^* )
+    //
+    // will check this
+    //
+    std::cout << GridLogMessage << "Checking that U[ i , j ] - conj U[ i+N/2, j+N/2 ] = 0 for i = 0, ... , N/2 " << std::endl;
+    std::cout <<GridLogMessage << std::endl;
+    const int nsp = Nc / 2;
+    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);
+            std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+            std::cout << GridLogMessage << "ok " << std::endl;
+        }
+        
+    }
+    std::cout <<GridLogMessage << std::endl;
+    std::cout << GridLogMessage << "Checking that U [ i , j+N/2 ] + conj U [ i+N/2, j+N/2 ] = 0 for i = 0, ... , N/2 " << std::endl;
+    std::cout <<GridLogMessage << std::endl;
+    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);
+            std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+            std::cout << GridLogMessage << "ok " << std::endl;
+        }
+    }
+    
+    
+    // an explicit check for sp2
+    if (Nc == 2)
+    {
+        assert(Nc==2);
+        ColourMatrix A;
+        A = Zero();
+
+        Complex a(25041994., 12.);
+        Complex b(39., 0.22);
+        Complex d(10000., -2222.3333);
+    
+        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();