From ac0a74be0dc1acf7d16de49f05e6ab6874cfceec Mon Sep 17 00:00:00 2001
From: Julian Lenz <j.j.lenz@swansea.ac.uk>
Date: Wed, 3 May 2023 02:32:42 +0100
Subject: [PATCH] Taken care of algebra tests

---
 tests/sp2n/Test_project_on_Sp.cc | 245 ++++++++++++-------------------
 1 file changed, 97 insertions(+), 148 deletions(-)

diff --git a/tests/sp2n/Test_project_on_Sp.cc b/tests/sp2n/Test_project_on_Sp.cc
index 4c2e7b76..046db4bb 100644
--- a/tests/sp2n/Test_project_on_Sp.cc
+++ b/tests/sp2n/Test_project_on_Sp.cc
@@ -68,7 +68,6 @@ void test_group_projections(T U) {
     LatticeColourMatrixD identity(U.Grid());
     identity = 1.0;
 
-    std::cout <<GridLogMessage << std::endl;
     std::cout << GridLogMessage << "#   #   #   #" << std::endl;
     std::cout << GridLogMessage << "Group" << std::endl;
     std::cout << GridLogMessage << "#   #   #   #" << std::endl;
@@ -99,6 +98,102 @@ void test_group_projections(T U) {
     assert(is_element_of_sp2n_group(U));
 }
 
+template<typename T>
+bool has_correct_algebra_block_structure(const T& U) {
+    const int nsp = Nc / 2;
+    Complex i(0., 1.);
+  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 );
+        }
+    }
+    return true;
+}
+
+template<typename T>
+bool is_element_of_sp2n_algebra(T U) {
+ 
+    LatticeColourMatrixD aux(U.Grid());
+    LatticeColourMatrixD identity(U.Grid());
+    identity = 1.0;
+
+    aux = U - adj(U);
+    std::cout << GridLogMessage << "T - Tda = " << norm2(aux) << std::endl;
+    assert( norm2(aux) < 1e-8);
+    aux = U + adj(U);
+    std::cout << GridLogMessage << "T + Tda = " << norm2(aux) << std::endl;
+    assert( norm2(aux) < 1e-8);
+    
+    std::cout << GridLogMessage << "Check that Omega U Omega = conj(U)" << std::endl;
+
+    LatticeColourMatrixD Omega(U.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);
+
+   return has_correct_algebra_block_structure(U); 
+}
+
+
+template< typename T>
+void test_algebra_projections(T U) {
+    RealD Delta = 666.;
+    LatticeColourMatrixD tmp(U.Grid());
+    LatticeColourMatrixD identity(U.Grid());
+    identity = 1.0;
+
+    std::cout << GridLogMessage << "#   #   #   #" << std::endl;
+    std::cout << GridLogMessage << "Algebra" << std::endl;
+    std::cout << GridLogMessage << "#   #   #   #" << std::endl;
+    std::cout <<GridLogMessage << std::endl;
+
+    std::string name="SpTa";
+    std::cout << GridLogMessage << "Testing "<< name << std::endl;
+    std::cout << GridLogMessage << "Apply to deformed matrix" << std::endl;
+
+     U = U + Delta*identity;
+    U = SpTa(U);
+    assert(is_element_of_sp2n_algebra(U));
+    
+name="TaProj";
+    std::cout << GridLogMessage << "Testing "<< name << std::endl;
+    std::cout << GridLogMessage << "Apply to deformed matrix" << std::endl;
+
+     U = U + Delta*identity;
+    Sp<Nc>::taProj(U, tmp);
+    U = tmp;
+    assert(is_element_of_sp2n_algebra(U));
+}    
 
 
 int main (int argc, char **argv)
@@ -139,155 +234,9 @@ int main (int argc, char **argv)
     
     SU<Nc>::HotConfiguration(pRNG,Umu);
     U = PeekIndex<LorentzIndex>(Umu,0);
-    
-    aux = U*adj(U) - identity;
-    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 << GridLogMessage << "This matrix should not leave Omega invariant, expect a warning" << std::endl;
-    }
-    Sp<Nc>::OmegaInvariance(U);
-    std::cout <<GridLogMessage << std::endl;
-    
-    U = U + epsilon*identity;
-    aux = U*adj(U) - identity;
-    
-    std::cout << GridLogMessage << "Unitary matrix deformed " << std::endl;
-    std::cout << GridLogMessage << "now U adjU - 1 = " << norm2(aux) << std::endl;
-    
     test_group_projections(U);
-
-    std::cout <<GridLogMessage << std::endl;
-    std::cout << GridLogMessage << "#   #   #   #" << std::endl;
-    std::cout << GridLogMessage << "Algebra" << std::endl;
-    std::cout << GridLogMessage << "#   #   #   #" << std::endl;
-    std::cout <<GridLogMessage << std::endl;
-    std::cout << GridLogMessage << "Testing SpTa" << std::endl;
-    
     U = PeekIndex<LorentzIndex>(Umu,1);
-    U = U + Delta*identity;
-    std::cout << GridLogMessage << "Matrix deformed " << std::endl;
-    std::cout << GridLogMessage << "Apply SpTa to deformed matrix" << std::endl;
-    U = SpTa(U);
-    
-    aux = U - adj(U);
-    std::cout << GridLogMessage << "SpTa ::: T - Tda = " << norm2(aux) << std::endl;
-    aux = U + adj(U);
-    std::cout << GridLogMessage << "SpTa ::: T + Tda = " << norm2(aux) << std::endl;
-    
-    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 );
-        }
-    }
-    
-    //test Ta
-    /*
-    U = U + 666.*identity;
-    Up = Ta(U);
-    aux = Up - adj(Up);
-    std::cout << GridLogMessage << "TA !!! T - Tda = " << norm2(aux) << std::endl;
-    aux = Up + adj(Up);
-    std::cout << GridLogMessage << "TA !!! T + Tda = " << norm2(aux) << std::endl;*/
-    
-    // test taProj
-    std::cout << GridLogMessage << "Testing taProj" << std::endl;
-    U = U + Delta*identity;
-    std::cout << GridLogMessage << "Matrix deformed " << std::endl;
-    std::cout << GridLogMessage << "Apply taProj to deformed matrix" << std::endl;
-    Sp<Nc>::taProj(U, Up);
-    aux = Up - adj(Up);
-    std::cout << GridLogMessage << "taProj ::: T - Tda = " << norm2(aux) << std::endl;
-    aux = Up + adj(Up);
-    std::cout << GridLogMessage << "taProj ::: T + Tda = " << norm2(aux) << std::endl;
-    
-    std::cout << GridLogMessage << "Check that Omega U Omega = conj(U)" << std::endl;
-    Sp<Nc>::Omega(Omega);
-    aux = Omega*Up*Omega - conjugate(Up);
-    std::cout << GridLogMessage << "Omega U Omega - conj(U) = " << norm2(aux) << std::endl;
-    assert( norm2(aux) < 1e-8);
-    // before it was
-    aux = Omega*U*Omega - conjugate(U);
-    std::cout << GridLogMessage << " before taProj Omega U Omega - conj(U) = " << norm2(aux) << std::endl;
-
-    U = Up;
-    
-    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 );
-        }
-    }
+    test_algebra_projections(U);
         
     Grid_finalize();