From 4e31e4e094530456653479f82b3990078eaa1e7d Mon Sep 17 00:00:00 2001
From: Alessandro Lupo <Alessandro.Lupo@ed.ac.uk>
Date: Wed, 13 Oct 2021 15:07:23 +0100
Subject: [PATCH] Better  tests

---
 Grid/qcd/utils/Sp2n.h            |  78 ++++++++++++++++++++++++
 configure.ac                     |   1 -
 tests/sp2n/Test_Sp_start.cc      | 100 +++++++++++++++++--------------
 tests/sp2n/Test_project_on_Sp.cc |  55 +++++++++--------
 4 files changed, 161 insertions(+), 73 deletions(-)

diff --git a/Grid/qcd/utils/Sp2n.h b/Grid/qcd/utils/Sp2n.h
index 1dafd771..ad6933e3 100644
--- a/Grid/qcd/utils/Sp2n.h
+++ b/Grid/qcd/utils/Sp2n.h
@@ -433,7 +433,85 @@ public:
       ColdConfiguration(out);
     }
     
+    static void OmegaInvariance(ColourMatrix &in)
+    {
+        
+        ColourMatrix Omega;
+        Omega = Zero();
+        
+        for (int i = 0; i < ncolour; i++)
+        {
+            Omega()()(i, 2*ncolour-1-i) = 1.;
+            Omega()()(2*ncolour-1-i, i) = -1;
+        }
+
+        auto diff = Omega - (in * Omega * transpose(in) );
+        auto sdiff = norm2(diff);
+        if (norm2(sdiff) < 1e-8)
+        {
+            std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl;
+        } else {
+            std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT left invariant by " << norm2(sdiff) << std::endl;
+        }
+        
+    }
     
+    template<typename GaugeField>
+    static void OmegaInvariance(GaugeField &in)
+    {
+        typedef typename GaugeField::vector_type vector_type;
+        typedef iSp2nMatrix<vector_type> vMatrixType;
+        typedef Lattice<vMatrixType> LatticeMatrixType;
+        
+        LatticeMatrixType U(in.Grid());
+        LatticeMatrixType OmegaLatt(in.Grid());
+        LatticeMatrixType identity(in.Grid());
+        RealD vol = in.Grid()->gSites();
+        ColourMatrix Omega;
+        
+        OmegaLatt = Zero();
+        Omega = Zero();
+        identity = 1.;
+        
+        U = PeekIndex<LorentzIndex>(in,1);
+        
+        OmegaInvariance(U);
+        
+    }
+    
+       static void OmegaInvariance(LatticeColourMatrixD &in)
+       {
+                      
+           LatticeColourMatrixD OmegaLatt(in.Grid());
+           LatticeColourMatrixD identity(in.Grid());
+           RealD vol = in.Grid()->gSites();
+           ColourMatrix Omega;
+           
+           OmegaLatt = Zero();
+           Omega = Zero();
+           identity = 1.;
+           
+           for (int i = 0; i < ncolour; i++)
+           {
+               Omega()()(i, 2*ncolour-1-i) = 1.;
+               Omega()()(2*ncolour-1-i, i) = -1;
+           }
+    
+           OmegaLatt = OmegaLatt + (identity*Omega);
+           
+           auto diff = OmegaLatt - (in*Omega*transpose(in));
+           auto sdiff = sum(diff);
+           //assert( norm2(sdiff) < 1e-8 );
+           if (norm2(sdiff) < 1e-8)
+           {
+               std::cout << GridLogMessage << "Symplectic condition satisfied: Omega invariant" << std::endl;
+           } else {
+               std::cout << GridLogMessage << "WARNING!!!!!! Matrix Omega NOT left invariant by " << norm2(sdiff) << std::endl;
+           }
+           
+           
+       }
+        
     
 }; // end of class Sp
     
diff --git a/configure.ac b/configure.ac
index caadca41..c6a85588 100644
--- a/configure.ac
+++ b/configure.ac
@@ -787,7 +787,6 @@ os (target)                 : $target_os
 compiler vendor             : ${ax_cv_cxx_compiler_vendor}
 compiler version            : ${ax_cv_gxx_version}
 ----- BUILD OPTIONS -----------------------------------
-#gauge group                 : `if test "${ac_ENABLE_SP}x" == "yesx"; then echo Sp2n; else echo SUn; fi`
 Nc                          : ${ac_Nc}
 SIMD                        : ${ac_SIMD}${SIMD_GEN_WIDTH_MSG}
 Threading                   : ${ac_openmp}
diff --git a/tests/sp2n/Test_Sp_start.cc b/tests/sp2n/Test_Sp_start.cc
index eaeb59f4..da19779f 100644
--- a/tests/sp2n/Test_Sp_start.cc
+++ b/tests/sp2n/Test_Sp_start.cc
@@ -14,54 +14,79 @@ int main (int argc, char **argv)
     GridRedBlackCartesian     RBGrid(&Grid);
     
     LatticeGaugeField Umu(&Grid);
-    LatticeGaugeField identity(&Grid);
     LatticeColourMatrixD U(&Grid);
-    LatticeColourMatrixD Cidentity(&Grid);
     LatticeColourMatrixD aux(&Grid);
-    
-    identity = 1.0;
-    Cidentity = 1.0;
-    Complex i(0., 1.);
+    LatticeGaugeField identity(&Grid);
+    LatticeColourMatrixD Cidentity(&Grid);
     
     double vol = Umu.Grid()->gSites();
-    
+
     const int nsp = Nc / 2;
+    identity = 1.;
+    Cidentity = 1.;
+    
     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);
     
     std::cout << GridLogMessage << std::endl;
-    
-    
+    std::cout << GridLogMessage << "-------" << std::endl;
+    std::cout << GridLogMessage << "Checking Cold Configuration " << std::endl;
+    std::cout << GridLogMessage << std::endl;
     Sp<nsp>::ColdConfiguration(pRNG,Umu);
+    Sp<nsp>::OmegaInvariance(Umu);
+    std::cout << GridLogMessage << std::endl;
+
     Umu = Umu - identity;
     U = PeekIndex<LorentzIndex>(Umu,1);
     
     std::cout << GridLogMessage << "U - 1 = " << norm2(sum(U)) << std::endl;
     assert ( norm2(sum(U)) < 1e-6);
     std::cout << GridLogMessage << std::endl;
-    
-    std::cout << GridLogMessage << "Hot Start " << std::endl;
+    std::cout << GridLogMessage << "-------" << std::endl;
+    std::cout << GridLogMessage << "Checking Hot Configuration" << std::endl;
+    std::cout << GridLogMessage << std::endl;
     Sp<nsp>::HotConfiguration(pRNG,Umu);
+    Sp<nsp>::OmegaInvariance(Umu);
+    std::cout << GridLogMessage << std::endl;
+    U = PeekIndex<LorentzIndex>(Umu,1);
+    std::cout << GridLogMessage << "Checking unitarity " << std::endl;
+    aux = U*adj(U) - Cidentity;
+    assert( norm2(aux) < 1e-6 );
+    std::cout << GridLogMessage << "ok" << std::endl;
+    std::cout << GridLogMessage << std::endl;
+    std::cout << GridLogMessage << "Checking determinant " << std::endl;
+    std::cout << GridLogMessage << "Det = " <<  norm2( Determinant(U) )/vol << std::endl;
+    std::cout << GridLogMessage << std::endl;
+    std::cout << GridLogMessage << "-------" << std::endl;
+    std::cout << GridLogMessage << "Checking Tepid Configuration" << std::endl;
+    std::cout << GridLogMessage << std::endl;
+    Sp<nsp>::TepidConfiguration(pRNG,Umu);
+    Sp<nsp>::OmegaInvariance(Umu);
+    std::cout << GridLogMessage << std::endl;
     U = PeekIndex<LorentzIndex>(Umu,1);
-    
     std::cout << GridLogMessage << "Checking unitarity " << std::endl;
-    
     aux = U*adj(U) - Cidentity;
-
     assert( norm2(aux) < 1e-6 );
-    
     std::cout << GridLogMessage << "ok" << std::endl;
     std::cout << GridLogMessage << std::endl;
-    
     std::cout << GridLogMessage << "Checking determinant " << std::endl;
     std::cout << GridLogMessage << "Det = " <<  norm2( Determinant(U) )/vol << std::endl;
     std::cout << GridLogMessage << std::endl;
     
-    std::cout << GridLogMessage << "Checking the matrix is in Sp2n " << 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;
+    std::cout << GridLogMessage  << std::endl;
+    
+    Complex i(0., 1.);
+    
+    
+    Sp<nsp>::HotConfiguration(pRNG,Umu);
+    U = PeekIndex<LorentzIndex>(Umu,0);
     for (int c1 = 0; c1 < nsp; c1++) //check on W
     {
         for (int c2 = 0; c2 < nsp; c2++)
@@ -71,17 +96,14 @@ int main (int argc, char **argv)
             auto Ww = conjugate( Wstar );
             auto amizero = sum(W - Ww);
             auto amizeroo = TensorRemove(amizero);
-            std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
+            //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 << "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++)
@@ -91,27 +113,18 @@ int main (int argc, char **argv)
             auto minusXx = conjugate(minusXstar);
             auto amizero = sum (X + minusXx);
             auto amizeroo = TensorRemove(amizero);
-            std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
+            //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 << "ok " << std::endl;
         }
     }
     
+    std::cout << GridLogMessage << "Hot start ok " << std::endl;
     
-    std::vector<int> ppseeds({11,12,13,14,15});
-    std::vector<int> ssseeds({16,17,18,19,20});
-    GridParallelRNG  ppRNG(&Grid); pRNG.SeedFixedIntegers(ppseeds);
-    GridSerialRNG    ssRNG;       sRNG.SeedFixedIntegers(ssseeds);
-    
-    std::cout << GridLogMessage << std::endl;
-    std::cout << GridLogMessage << "Testing gauge implementation " << std::endl;
-    
-    SymplPeriodicGimplR::HotConfiguration(ppRNG, Umu);
-    
-    U = PeekIndex<LorentzIndex>(Umu,1);
-    
+    Sp<nsp>::HotConfiguration(pRNG,Umu);
+    U = PeekIndex<LorentzIndex>(Umu,0);
     for (int c1 = 0; c1 < nsp; c1++) //check on W
     {
         for (int c2 = 0; c2 < nsp; c2++)
@@ -121,17 +134,14 @@ int main (int argc, char **argv)
             auto Ww = conjugate( Wstar );
             auto amizero = sum(W - Ww);
             auto amizeroo = TensorRemove(amizero);
-            std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
+            //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 << "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++)
@@ -141,15 +151,17 @@ int main (int argc, char **argv)
             auto minusXx = conjugate(minusXstar);
             auto amizero = sum (X + minusXx);
             auto amizeroo = TensorRemove(amizero);
-            std::cout << GridLogMessage << "diff(real,im) = " << amizeroo << std::endl;
+            //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 << "ok " << std::endl;
         }
     }
     
+    std::cout << GridLogMessage << "Tepid start ok " << std::endl;
     
+    std::cout << GridLogMessage << std::endl;
     
     Grid_finalize();
 
diff --git a/tests/sp2n/Test_project_on_Sp.cc b/tests/sp2n/Test_project_on_Sp.cc
index d494d828..4e14938e 100644
--- a/tests/sp2n/Test_project_on_Sp.cc
+++ b/tests/sp2n/Test_project_on_Sp.cc
@@ -4,7 +4,7 @@ using namespace Grid;
 
 int main (int argc, char **argv)
 {
-     Grid_init(&argc,&argv);
+    Grid_init(&argc,&argv);
     
     
     Coordinate latt_size   = GridDefaultLatt();
@@ -20,6 +20,8 @@ int main (int argc, char **argv)
     LatticeColourMatrixD aux(&Grid);
     LatticeColourMatrixD identity(&Grid);
     
+    const int nsp = Nc / 2;
+    
     identity = 1.0;
     RealD epsilon = 0.01;
     Complex i(0., 1.);
@@ -34,20 +36,28 @@ int main (int argc, char **argv)
     SU<Nc>::HotConfiguration(pRNG,Umu);
     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 << "Starting with random SUn matrix " << std::endl;
+    std::cout << GridLogMessage << "Unitary check " << std::endl;
     std::cout <<GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
+    assert ( norm2(aux) < 1e-8 );
     std::cout <<GridLogMessage << std::endl;
+    if (Nc != 2)
+    {
+        std::cout << "This matrix should not leave Omega invariant, expect a warning" << std::endl;
+    }
+    Sp<nsp>::OmegaInvariance(U);
     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 << "Unitary matrix deformed " << std::endl;
     std::cout << GridLogMessage << "now U adjU - 1 = " << norm2(aux) << std::endl;
     std::cout <<GridLogMessage << std::endl;
+    std::cout << "This matrix should not leave Omega invariant, expect a warning" << std::endl;
+    Sp<nsp>::OmegaInvariance(U);
     std::cout <<GridLogMessage << std::endl;
     std::cout << GridLogMessage << "Projecting on Sp2n " << std::endl;
 
@@ -65,18 +75,17 @@ int main (int argc, char **argv)
     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<nsp>::OmegaInvariance(U);
+    
     // 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 << "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;
-    const int nsp = Nc / 2;
     for (int c1 = 0; c1 < nsp; c1++) //check on W
     {
         for (int c2 = 0; c2 < nsp; c2++)
@@ -86,17 +95,13 @@ int main (int argc, char **argv)
             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++)
@@ -106,16 +111,16 @@ int main (int argc, char **argv)
             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;
         }
     }
     
+    std::cout << GridLogMessage << "ok" << std::endl;
     
     // an explicit check for sp2
+    /*
     if (Nc == 2)
     {
         assert(Nc==2);
@@ -141,13 +146,7 @@ int main (int argc, char **argv)
         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();