expand ProjecOnGaugeGroup, change ProjectOnSp2nAlgebra into SpTa, fixing some of its issues

Grid/lattice/Lattice_ET.h

@@ -345,7 +345,7 @@ GridUnopClass(UnaryNot, Not(a));
 GridUnopClass(UnaryTrace, trace(a));
 GridUnopClass(UnaryTranspose, transpose(a));
 GridUnopClass(UnaryTa, Ta(a));
-GridUnopClass(UnaryProjectSp2nAlgebra, ProjectSp2nAlgebra(a));
+GridUnopClass(UnarySpTa, SpTa(a));
 GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a));
 GridUnopClass(UnaryProjectOnSpGroup, ProjectOnSpGroup(a));
 GridUnopClass(UnaryTimesI, timesI(a));
@@ -458,7 +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(SpTa, UnarySpTa);
 GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup);
 GRID_DEF_UNOP(ProjectOnSpGroup, UnaryProjectOnSpGroup);
 GRID_DEF_UNOP(timesI, UnaryTimesI);

Grid/qcd/utils/GaugeGroup.h

@@ -300,9 +300,10 @@ class GaugeGroup {
     ColdConfiguration(out);
   }
 
-  template <typename LatticeMatrixType, ONLY_IF_SU>
+  template <typename LatticeMatrixType>
   static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) {
-    out = Ta(in);
+    taProj(in, out, group_name());
+    //out = Ta(in);
   }
     
   template <typename LatticeMatrixType>
@@ -333,6 +334,32 @@ class GaugeGroup {
       ProjectOnGaugeGroup(Umu, group_name());
   }
     
+    template <class vtype>
+    static void ProjectOnGaugeGroup(iScalar<vtype> &r) {
+      ProjectOnGaugeGroup(r, group_name());
+    }
+
+    template <class vtype, int N>
+    static void ProjectOnGaugeGroup(iVector<vtype,N> &r) {
+      r = ProjectOnGaugeGroup(r, group_name());
+    }
+
+    template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
+    static void ProjectOnGaugeGroup(iMatrix<vtype,N> &arg) {
+      arg = ProjectOnGaugeGroup(arg, group_name());
+    }
+
+    
+/*
+  template <int N>        // reunitarise, resimplectify...
+  static void ProjectOnGaugeGroup(iVector<iScalar<iMatrix<vComplexD, N> >, Nd> &U) {
+      ProjectOnGaugeGroup(U, group_name());
+  }
+         
+  template <int N>        // reunitarise, resimplectify...
+  static void ProjectOnGaugeGroup(iScalar<iScalar<iMatrix<vComplexD, N> > > &Umu) {
+      ProjectOnGaugeGroup(Umu, group_name());
+  }*/
 
   template <int N>       // reunitarise, resimplectify... previously ProjectSUn
   static void ProjectGn(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
@@ -349,7 +376,7 @@ class GaugeGroup {
    }
 
     template <int N>    // reunitarise, resimplectify... previously ProjectSUn
-    static void ProjectGn(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U, GroupName::SU) {
+    static void ProjectGn(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
       // Reunitarise
       for (int mu = 0; mu < Nd; mu++) {
         auto Umu = PeekIndex<LorentzIndex>(U, mu);

Grid/qcd/utils/SUn.impl

@@ -508,6 +508,7 @@ static void testGenerators(GroupName::SU) {
   std::cout << GridLogMessage << std::endl;
 }
 
+
 template <int N>
 static void ProjectOnGaugeGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::SU) {
   Umu = ProjectOnGroup(Umu);
@@ -522,6 +523,44 @@ static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >,
   }
 }
 
+template <class vtype>
+static void ProjectOnGaugeGroup(iScalar<vtype> &r, GroupName::SU) {
+  r = ProjectOnGroup(r);
+}
+
+template <class vtype, int N>
+static void ProjectOnGaugeGroup(iVector<vtype,N> &r, GroupName::SU) {
+  r = ProjectOnGroup(r);
+}
+
+template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
+static void ProjectOnGaugeGroup(iMatrix<vtype,N> &arg, GroupName::SU) {
+  arg = ProjectOnGroup(arg);
+}
+
+
+
+
+/*
+template <int N>    // non-lattice objects
+static void ProjectOnGaugeGroup(iScalar<iScalar<iMatrix<vComplexD, N> > >  &Umu, GroupName::SU) {
+  Umu = ProjectOnGroup(Umu);
+ }
+
+template <int N>
+static void ProjectOnGaugeGroup(iVector<iScalar<iMatrix<vComplexD, N> >, Nd>  &U, GroupName::SU) {
+  // Reunitarise
+  for (int mu = 0; mu < Nd; mu++) {
+    auto Umu = PeekIndex<LorentzIndex>(U, mu);
+    Umu = ProjectOnGroup(Umu);
+  }
+}*/
+
+template <typename LatticeMatrixType>
+static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out, GroupName::SU) {
+  out = Ta(in);
+}
+
 /*
  * Fundamental rep gauge xform
  */

Grid/qcd/utils/Sp2n.impl

@@ -267,6 +267,41 @@ static void ProjectOnGaugeGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >,
   }
 }
 
+template <class vtype>
+static void ProjectOnGaugeGroup(iScalar<vtype> &r, GroupName::Sp) {
+  r = ProjectOnSpGroup(r);
+}
+
+template <class vtype, int N>
+static void ProjectOnGaugeGroup(iVector<vtype,N> &r, GroupName::Sp) {
+  r = ProjectOnSpGroup(r);
+}
+
+template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
+static void ProjectOnGaugeGroup(iMatrix<vtype,N> &arg, GroupName::Sp) {
+  arg = ProjectOnSpGroup(arg);
+}
+
+/*
+template <int N>    // non-lattice objects
+static void ProjectOnGaugeGroup(iScalar<iScalar<iMatrix<vComplexD, N> > >  &Umu, GroupName::Sp) {
+  Umu = ProjectOnSpGroup(Umu);
+ }
+
+template <int N>
+static void ProjectOnGaugeGroup(iVector<iScalar<iMatrix<vComplexD, N> >, Nd>  &U, GroupName::Sp) {
+  // Reunitarise
+  for (int mu = 0; mu < Nd; mu++) {
+    auto Umu = PeekIndex<LorentzIndex>(U, mu);
+    Umu = ProjectOnSpGroup(Umu);
+  }
+}*/
+
+template <typename LatticeMatrixType>   
+static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out, GroupName::Sp) {
+  out = SpTa(in);
+}
+
 public:
 template <ONLY_IF_Sp>
 static void OmegaInvariance(ColourMatrix &in) {

Grid/tensors/Tensor_Ta.h

@@ -66,23 +66,23 @@ 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. We do a Gram-Schmidt
+// for sp2n can't be as simple as Ta. We do a Gram-Schmidt
 
-template<class vtype> accelerator_inline iScalar<vtype> ProjectSp2nAlgebra(const iScalar<vtype>&r)
+template<class vtype> accelerator_inline iScalar<vtype> SpTa(const iScalar<vtype>&r)
 {
   iScalar<vtype> ret;
-  ret._internal = ProjectSp2nAlgebra(r._internal);
+  ret._internal = SpTa(r._internal);
-  return ret;
+  return Ta(ret);
 }
-template<class vtype,int N> accelerator_inline iVector<vtype,N> ProjectSp2nAlgebra(const iVector<vtype,N>&r)
+template<class vtype,int N> accelerator_inline iVector<vtype,N> SpTa(const iVector<vtype,N>&r)
 {
   iVector<vtype,N> ret;
   for(int i=0;i<N;i++){
-    ret._internal[i] = ProjectSp2nAlgebra(r._internal[i]);
+    ret._internal[i] = SpTa(r._internal[i]);
   }
-  return ret;
+  return Ta(ret);
 }
-template<class vtype,int N> accelerator_inline iMatrix<vtype,N> ProjectSp2nAlgebra(const iMatrix<vtype,N> &arg)
+template<class vtype,int N> accelerator_inline iMatrix<vtype,N> SpTa(const iMatrix<vtype,N> &arg)
 {
     iMatrix<vtype,N> ret;
     vtype nrm;
@@ -142,7 +142,7 @@ template<class vtype,int N> accelerator_inline iMatrix<vtype,N> ProjectSp2nAlgeb
     }
     
   
-  return ret;
+  return Ta(ret);
 }
 
 

tests/sp2n/Test_project_on_Sp.cc

@@ -17,6 +17,7 @@ int main (int argc, char **argv)
     
     LatticeGaugeField Umu(&Grid);
     LatticeColourMatrixD U(&Grid);
+    LatticeColourMatrixD Up(&Grid);
     LatticeColourMatrixD aux(&Grid);
     LatticeColourMatrixD identity(&Grid);
     
@@ -59,6 +60,11 @@ 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 << "#   #   #   #" << 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;
@@ -123,7 +129,7 @@ int main (int argc, char **argv)
     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++)
@@ -200,23 +206,27 @@ int main (int argc, char **argv)
             assert(  amizeroo.real() < 10e-6 );
         }
     }
-    
+    std::cout <<GridLogMessage << std::endl;
-    std::cout << GridLogMessage << "Testing ProjectSp2nAlgebra" << 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 + 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 << "Apply SpTa to deformed matrix" << std::endl;
-    std::cout << GridLogMessage << "Apply ProjectSp2nAlgebra to deformed matrix" << std::endl;
+    U = SpTa(U);
-    U = ProjectSp2nAlgebra(U);
+    
-    aux = U*adj(U) - identity;
+    aux = U - adj(U);
-    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
+    std::cout << GridLogMessage << "SpTa ::: T - Tda = " << norm2(aux) << std::endl;
-    assert( norm2(aux) < 1e-8);
+    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;
@@ -258,6 +268,71 @@ int main (int argc, char **argv)
         }
     }
     
+    //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 + 666.*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 );
+        }
+    }
     
     
     Grid_finalize();