fix projection

Grid/qcd/utils/Sp2n.h

@@ -493,13 +493,14 @@ public:
            
            for (int i = 0; i < ncolour; i++)
            {
-               Omega()()(i, 2*ncolour-1-i) = 1.;
-               Omega()()(2*ncolour-1-i, i) = -1;
+               Omega()()(i, ncolour+i) = 1.;
+               Omega()()(ncolour+i, i) = -1;
            }
            
+           std::cout << GridLogMessage << "Omega = " << Omega()() << std::endl;
            OmegaLatt = OmegaLatt + (identity*Omega);
            
-           auto diff = OmegaLatt - (in*Omega*transpose(in));
+           auto diff = OmegaLatt - (in*OmegaLatt*transpose(in));
            auto sdiff = sum(diff);
            //assert( norm2(sdiff) < 1e-8 );
            if (norm2(sdiff) < 1e-8)

Grid/tensors/Tensor_Ta.h

@@ -156,6 +156,75 @@ template<class vtype,int N> accelerator_inline iVector<vtype,N> ProjectOnSpGroup
 // int N is 2n in Sp(2n)
 template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
 accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg)
+{
+  // need a check for the group type?
+  iMatrix<vtype,N> ret(arg);
+  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;
+}
+
+
+
+
+
+
+/*
+template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
+accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg)
 {
   // need a check for the group type?
   iMatrix<vtype,N> ret(arg);
@@ -175,7 +244,7 @@ accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg
 
     // 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]);
+      ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]); // lupo: c1,c2? dont we want c1 c2+N/2
 
     for(int c2=N/2;c2<N;c2++)
       ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);
@@ -203,7 +272,7 @@ accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg
         ret._internal[b][c] -= pr * ret._internal[c1+N/2][c];
       }
     }
-  }
+  } //end c1 0, ... N/2-1 loop
 
   // Compute the last row
   {
@@ -227,7 +296,7 @@ accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg
   }
   // assuming the determinant is ok
   return ret;
-}
+}*/
 
 
 NAMESPACE_END(Grid);

tests/sp2n/Test_project_on_Sp.cc

@@ -45,7 +45,7 @@ int main (int argc, char **argv)
     std::cout <<GridLogMessage << std::endl;
     if (Nc != 2)
     {
-        std::cout << "This matrix should not leave Omega invariant, expect a warning" << std::endl;
+        std::cout << GridLogMessage << "This matrix should not leave Omega invariant, expect a warning" << std::endl;
     }
     Sp<nsp>::OmegaInvariance(U);
     std::cout <<GridLogMessage << std::endl;
@@ -56,18 +56,18 @@ 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 << "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;
 
     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 << "U adjU - 1 = " << norm2(aux) << std::endl;
-    assert( norm2(aux) < 1e-6);
+    assert( norm2(aux) < 1e-8);
     std::cout <<GridLogMessage << std::endl;
     std::cout <<GridLogMessage << std::endl;
     
@@ -82,6 +82,8 @@ int main (int argc, char **argv)
     
     // 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;