better projection test

Grid/tensors/Tensor_Ta.h

@@ -212,7 +212,6 @@ accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg
   iMatrix<vtype,N> ret(arg);
   vtype nrm;
   vtype inner;
-  vtype tmp;
   
   for(int c1=0;c1<N/2;c1++)
   {
@@ -250,20 +249,15 @@ accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg
       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] = conjugate(ret._internal[c1][c2]);          // down right in the new matrix = (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] = -conjugate(ret._internal[c1][c2]);;     // down left in the new matrix = -(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;
 }

tests/sp2n/Test_project_on_Sp.cc

@@ -2,11 +2,121 @@
 
 using namespace Grid;
 
+int SpGroupQuiz (const LatticeColourMatrixD Uin)
+{
+    double vol = Uin.Grid()->gSites();
+    const int nsp = Nc / 2;
+    LatticeColourMatrixD Omega(Uin.Grid());
+    Sp<Nc>::Omega(Omega);
+    LatticeColourMatrixD aux(Uin.Grid());
+    LatticeColourMatrixD identity(Uin.Grid());
+    Complex i(0., 1.);
+    identity = 1;
+    
+    std::cout << GridLogMessage << "Check matrix is non-zero " << std::endl;
+    assert(norm2(Uin) > 1e-8);
+    aux = Uin*adj(Uin) - identity;
+    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux)/vol << std::endl;
+    assert(norm2(aux) < 1e-8);
+    aux = Omega - (Uin * Omega * transpose(Uin));
+    std::cout << GridLogMessage << "Omega - U Omega transpose(U) = " << norm2(aux)/vol << std::endl;
+    assert(norm2(aux) < 1e-8);
+    //Sp<Nc>::OmegaInvariance(Uin)
+    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
+    std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
+    std::cout << GridLogMessage << "      (  -X^*  W^* )  " << std::endl;
+    for (int c1 = 0; c1 < nsp; c1++) //check on W
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto W = PeekIndex<ColourIndex>(Uin,c1,c2);
+            auto Wstar =  PeekIndex<ColourIndex>(Uin,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>(Uin,c1,c2+nsp);
+            auto minusXstar = PeekIndex<ColourIndex>(Uin,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 << "|Det| = " << norm2( Determinant(Uin) ) / vol << std::endl;
+    assert( norm2( Determinant(Uin) ) / vol - 1 < 1e-8);
+    return 0;
+}
+
+int SpAntiHermitianAlgebraQuiz (const LatticeColourMatrixD Uin)
+{
+    double vol = Uin.Grid()->gSites();
+    const int nsp = Nc / 2;
+    LatticeColourMatrixD Omega(Uin.Grid());
+    Sp<Nc>::Omega(Omega);
+    LatticeColourMatrixD aux(Uin.Grid());
+    LatticeColourMatrixD identity(Uin.Grid());
+    Complex i(0., 1.);
+    identity = 1;
+    
+    std::cout << GridLogMessage << "Check matrix is non-zero " << std::endl;
+    assert(norm2(Uin) > 1e-8);
+    aux = Uin - adj(Uin);
+    std::cout << GridLogMessage << "SpTa ::: T - Tda = " << norm2(aux)/vol << std::endl;
+    aux = Uin + adj(Uin);
+    std::cout << GridLogMessage << "SpTa ::: T + Tda = " << norm2(aux)/vol << std::endl;
+    assert( norm2(aux) - 1 < 1e-8);
+    std::cout << GridLogMessage << "Check that Omega T Omega + conj(T) = 0 " << std::endl;
+    aux = Omega*Uin*Omega + conjugate(Uin);
+    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;
+    for (int c1 = 0; c1 < nsp; c1++) //check on W
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto W = PeekIndex<ColourIndex>(Uin,c1,c2);
+            auto Wstar =  PeekIndex<ColourIndex>(Uin,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>(Uin,c1,c2+nsp);
+            auto minusXstar = PeekIndex<ColourIndex>(Uin,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 0;
+}
+
 int main (int argc, char **argv)
 {
     Grid_init(&argc,&argv);
     
-    
     Coordinate latt_size   = GridDefaultLatt();
     Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
     Coordinate mpi_layout  = GridDefaultMpi();
@@ -14,27 +124,18 @@ int main (int argc, char **argv)
     GridCartesian             Grid(latt_size,simd_layout,mpi_layout);
     GridRedBlackCartesian     RBGrid(&Grid);
     
-    
     LatticeGaugeField Umu(&Grid);
     LatticeColourMatrixD U(&Grid);
-    LatticeColourMatrixD Up(&Grid);
     LatticeColourMatrixD aux(&Grid);
     LatticeColourMatrixD identity(&Grid);
-    
+    LatticeColourMatrixD Omega(&Grid);
-    // Will test resimplectification-related functionalities (from ProjectOnGaugeGroup, ProjectOnSpGroup, ProjectGn) and projection on the algebra (from ProjectSp2nAlgebra)
+    Sp<Nc>::Omega(Omega);
-    // 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;
-    RealD epsilon = 0.01;
+    RealD epsilon = 0.1;
-    RealD Delta = 666.;
     Complex i(0., 1.);
-    RealD u = 0.;
     double vol = Umu.Grid()->gSites();
+    const int nsp = Nc / 2;
     
     std::vector<int> pseeds({1,2,3,4,5});
     GridParallelRNG  pRNG(&Grid); pRNG.SeedFixedIntegers(pseeds);
@@ -42,236 +143,76 @@ 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;
+    aux = U*adj(U) - identity;
     std::cout <<GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
     assert ( norm2(aux) < 1e-8 );
-    std::cout <<GridLogMessage << std::endl;
+    if (Nc != 2)    // Sp2 = SU2
-    if (Nc != 2)
     {
-        std::cout << GridLogMessage << "This matrix should not leave Omega invariant, expect a warning" << std::endl;
+        std::cout << GridLogMessage <<  "Checking matrix does NOT leave Omega invariant, to avoid a trivial test " << std::endl;
+        aux = Omega - (U * Omega * transpose(U));
+        assert ( norm2(aux) > 1e-8 );
     }
-    Sp<Nc>::OmegaInvariance(U);
-    std::cout <<GridLogMessage << std::endl;
-    
-    U = U + epsilon*identity;
-    aux = U*adj(U) - identity;
     
+    U = U + epsilon*identity - i*identity;
     std::cout << GridLogMessage << "Unitary matrix deformed " << std::endl;
+    auto det =  sum( Determinant(U) );
+    std::cout << GridLogMessage << "Re(Det) = " << real(det) / vol << std::endl;
+    det = det*i;
+    std::cout << GridLogMessage << "Im(Det) = " << real(det) / vol << std::endl;
+    aux = U*adj(U) - identity;
     std::cout << GridLogMessage << "now U adjU - 1 = " << norm2(aux) << std::endl;
     
-    std::cout <<GridLogMessage << std::endl;
+    //  ProjectOnSpGroup
-    std::cout << GridLogMessage << "#   #   #   #" << std::endl;
-    std::cout << GridLogMessage << "Group" << std::endl;
-    std::cout << GridLogMessage << "#   #   #   #" << std::endl;
-    std::cout <<GridLogMessage << std::endl;
-    //  Testing ProjectOnSpGroup
     std::cout << GridLogMessage << "Testing ProjectOnSpGroup" << std::endl;
     std::cout << GridLogMessage << "Apply ProjectOnSpGroup to deformed matrix" << std::endl;
     U = ProjectOnSpGroup(U);
+    std::cout << GridLogMessage << "Run checks:" << std::endl;
+    SpGroupQuiz(U);
+    det =  sum( Determinant(U) );
+    std::cout << GridLogMessage << "Re(Det) after ProjectOnSpGroup (nothing to assert) = " << real(det) / vol << std::endl;
+    det = det*i;
+    std::cout << GridLogMessage << "Im(Det) after ProjectOnSpGroup (nothing to assert) = " << real(det) / vol << std::endl;
     
-    aux = U*adj(U) - identity;
+    //  ProjectOnGaugeGroup
-    std::cout << GridLogMessage << "Unitary check after ProjectOnSpGroup " << std::endl;
-    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
-    assert( norm2(aux) < 1e-8);
-    
-    // actual sp2n check
-    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 << "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 );
-        }
-    }
-    
     SU<Nc>::HotConfiguration(pRNG,Umu); //refresh
-    U = PeekIndex<LorentzIndex>(Umu,1);
+    U = PeekIndex<LorentzIndex>(Umu,0);
-    
     std::cout << GridLogMessage << "Testing ProjectOnGaugeGroup" << std::endl;
-    U = U + Delta*identity;
+    U = U + epsilon*identity - i*identity;
     std::cout << GridLogMessage << "Apply ProjectOnGaugeGroup to deformed matrix" << std::endl;
     Sp<Nc>::ProjectOnGaugeGroup(U);
-    aux = U*adj(U) - identity;
+    std::cout << GridLogMessage <<"Run checks:" << std::endl;
-    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
+    SpGroupQuiz(U);
-    assert( norm2(aux) < 1e-8);
+    det =  sum( Determinant(U) );
-    Sp<Nc>::OmegaInvariance(U);
+    std::cout << GridLogMessage << "Re(Det) after ProjectOnGaugeGroup (nothing to assert) = " << real(det) / vol << std::endl;
-    
+    det = det*i;
-    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
+    std::cout << GridLogMessage << "Im(Det) after ProjectOnGaugeGroup (nothing to assert) = " << real(det) / vol << std::endl;
-    std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
-    std::cout << GridLogMessage << "      (  -X^*  W^* )  " << 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 );
-        }
-    }
     
+    //  ProjectGn
+    SU<Nc>::HotConfiguration(pRNG,Umu); //refresh
+    U = PeekIndex<LorentzIndex>(Umu,0);
     std::cout << GridLogMessage << "Testing ProjectGn" << std::endl;
-    U = U + Delta*identity;
+    U = U + epsilon*identity - i*identity;
     std::cout << GridLogMessage << "Apply ProjectGn to deformed matrix" << std::endl;
     Sp<Nc>::ProjectGn(U);
-    aux = U*adj(U) - identity;
+    std::cout << GridLogMessage << "Run checks:" << std::endl;
-    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
+    SpGroupQuiz(U);
-    assert( norm2(aux) < 1e-8);
+    det =  sum( Determinant(U) );
-    std::cout << GridLogMessage << "Det after ProjectGn = " << norm2( Determinant(U) ) / vol << std::endl;
+    std::cout << GridLogMessage << "Re(Det) after ProjectGn (constrained) = " << real(det) / vol << std::endl;
-    assert( norm2(aux) - 1 < 1e-8);
+    assert ( (real(det) / vol) - 1 < 1e-8 );
-    Sp<Nc>::OmegaInvariance(U);
+    det = det*i;
-    
+    std::cout << GridLogMessage << "Im(Det) after ProjectGn (constrained) = " << real(det) / vol << std::endl;
-    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
+    assert ( real(det) / vol < 1e-8 );
-    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 << 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;
     
+    //  SpTa
     SU<Nc>::HotConfiguration(pRNG,Umu);//refresh
-    U = PeekIndex<LorentzIndex>(Umu,2);
+    U = PeekIndex<LorentzIndex>(Umu,0);
-
+    std::cout << GridLogMessage << "Testing SpTa" << std::endl;
-    U = U + Delta*identity;
+    U = U + epsilon*identity - i*identity;
-    std::cout << GridLogMessage << "Matrix deformed " << std::endl;
     std::cout << GridLogMessage << "Apply SpTa to deformed matrix" << std::endl;
     U = SpTa(U);
-    
+    std::cout << GridLogMessage << "Run checks:" << std::endl;
-    aux = U - adj(U);
+    SpAntiHermitianAlgebraQuiz(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;
-    assert( norm2(aux) - 1 < 1e-8);
-    
-    std::cout << GridLogMessage << "Check that Omega T Omega + conj(T) = 0 " << std::endl;
-
-    LatticeColourMatrixD Omega(&Grid);
-    Sp<Nc>::Omega(Omega);
-    aux = Omega*U*Omega + conjugate(U);
-    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 );
-        }
-    }
         
     Grid_finalize();