projection on Sp2n algebra, to be used instead of Ta

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);

Grid/qcd/utils/Sp2n.impl

@@ -414,13 +414,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.;
-  //     Omega()()(2*ncolour-1-i, i) = -1;
-  // }
   for (int i = 0; i < nsp; i++) {
     Omega()()(i, nsp + i) = 1.;
     Omega()()(nsp + i, i) = -1;
@@ -454,8 +447,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);
@@ -473,8 +464,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;
@@ -496,3 +485,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;
+}

Grid/tensors/Tensor_Ta.h

@@ -66,6 +66,88 @@ 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
+
+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;
+    // following HiRep's suN_utils.c does a sort of Gram-Schmidt
+    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 
@@ -137,8 +219,6 @@ accelerator_inline iMatrix<vtype,N> ProjectOnGroup(const iMatrix<vtype,N> &arg)
 
 // re-do for sp2n
 
-// Ta cannot be defined here for Sp2n because I need the generators from the Sp class
-// It is defined in gauge impl types
 
 template<class vtype> accelerator_inline iScalar<vtype> ProjectOnSpGroup(const iScalar<vtype>&r)
 {

tests/sp2n/Test_project_on_Sp.cc

@@ -20,6 +20,8 @@ int main (int argc, char **argv)
     LatticeColourMatrixD aux(&Grid);
     LatticeColourMatrixD identity(&Grid);
     
+    // Will test resimplectification (from ProjectOnGaugeGroup, ProjectOnSpGroup) and projection on the algebra (from ProjectSp2nAlgebra)
+    
     const int nsp = Nc / 2;
     
     identity = 1.0;
@@ -29,15 +31,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 +54,27 @@ 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 << "Simplectify" << std::endl;
-
-    std::cout <<GridLogMessage << std::endl;
-    std::cout << GridLogMessage << "Projecting on Sp2n " << 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 simplectification " << 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 << "Det after simplectification = " << norm2( Determinant(U) ) / vol << std::endl;
-    std::cout <<GridLogMessage << std::endl;
+    assert( norm2(aux) - 1 < 1e-8);
-    std::cout <<GridLogMessage << std::endl;
     
     // actual sp2n check
-    std::cout << GridLogMessage << "Checking invariance after projection "<< std::endl;
+    std::cout << GridLogMessage << "Checking invariance after simplectification. Will not kill if it fails, but the next check will " << std::endl;
     Sp<Nc>::OmegaInvariance(U);
     
     // 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;
@@ -121,34 +112,110 @@ int main (int argc, char **argv)
     
     std::cout << GridLogMessage << "ok" << std::endl;
     
-    // an explicit check for sp2
+    std::cout << GridLogMessage << "Check the ProjectOnGaugeGroup function to simplectify" << std::endl;
-    /*
+    U = U + 2932.111*identity;
-    if (Nc == 2)
+    std::cout << GridLogMessage << "Resimplectify deformed matrix" << std::endl;
-    {
+    Sp<Nc>::ProjectOnGaugeGroup(U);
-        assert(Nc==2);
+    aux = U*adj(U) - identity;
-        ColourMatrix A;
+    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
-        A = Zero();
+    assert( norm2(aux) < 1e-8);
-
+    Sp<Nc>::OmegaInvariance(U);
-        Complex a(25041994., 12.);
+    
-        Complex b(39., 0.22);
+    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
-        Complex d(10000., -2222.3333);
+    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 << "ok" << std::endl;
+    
+    std::cout << GridLogMessage << "Now check projection on the algebra" << 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 << "Project on sp2n algebra" << 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();