From 7e3b298d3d712b5c8309454939cbf1cee2e82a33 Mon Sep 17 00:00:00 2001
From: Alessandro Lupo <Alessandro.Lupo@ed.ac.uk>
Date: Mon, 11 Oct 2021 16:28:15 +0100
Subject: [PATCH] project on sp2n

---
 Grid/lattice/Lattice_ET.h |   2 +
 Grid/tensors/Tensor_Ta.h  | 108 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 110 insertions(+)

diff --git a/Grid/lattice/Lattice_ET.h b/Grid/lattice/Lattice_ET.h
index 4a8a7423..9042eb8f 100644
--- a/Grid/lattice/Lattice_ET.h
+++ b/Grid/lattice/Lattice_ET.h
@@ -346,6 +346,7 @@ GridUnopClass(UnaryTrace, trace(a));
 GridUnopClass(UnaryTranspose, transpose(a));
 GridUnopClass(UnaryTa, Ta(a));
 GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a));
+GridUnopClass(UnaryProjectOnSpGroup, ProjectOnSpGroup(a));
 GridUnopClass(UnaryTimesI, timesI(a));
 GridUnopClass(UnaryTimesMinusI, timesMinusI(a));
 GridUnopClass(UnaryAbs, abs(a));
@@ -457,6 +458,7 @@ GRID_DEF_UNOP(trace, UnaryTrace);
 GRID_DEF_UNOP(transpose, UnaryTranspose);
 GRID_DEF_UNOP(Ta, UnaryTa);
 GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup);
+GRID_DEF_UNOP(ProjectOnSpGroup, UnaryProjectOnSpGroup);
 GRID_DEF_UNOP(timesI, UnaryTimesI);
 GRID_DEF_UNOP(timesMinusI, UnaryTimesMinusI);
 GRID_DEF_UNOP(abs, UnaryAbs);  // abs overloaded in cmath C++98; DON'T do the
diff --git a/Grid/tensors/Tensor_Ta.h b/Grid/tensors/Tensor_Ta.h
index 90e57b2b..f53296bc 100644
--- a/Grid/tensors/Tensor_Ta.h
+++ b/Grid/tensors/Tensor_Ta.h
@@ -135,6 +135,114 @@ accelerator_inline iMatrix<vtype,N> ProjectOnGroup(const iMatrix<vtype,N> &arg)
   return ret;
 }
 
+// re-do for sp2n
+
+template<class vtype> accelerator_inline iScalar<vtype> ProjectOnSpGroup(const iScalar<vtype>&r)
+{
+  iScalar<vtype> ret;
+  ret._internal = ProjectOnSpGroup(r._internal);
+  return ret;
+}
+template<class vtype,int N> accelerator_inline iVector<vtype,N> ProjectOnSpGroup(const iVector<vtype,N>&r)
+{
+  iVector<vtype,N> ret;
+  for(int i=0;i<N;i++){
+    ret._internal[i] = ProjectOnSpGroup(r._internal[i]);
+  }
+  return ret;
+}
+
+
+// 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;
+  for(int c1=0;c1<N;c1++){
+
+       // Normalises row c1
+        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;
+
+        // Compute row c1+N2/2: c1+N/2 = - \Omega c1*
+        for(int c2=0;c2<N/2;c2++)
+          //tmp = conjugate(ret._internal[c1][c2]);
+          //ret._internal[c1+N/2][c2+N/2] = tmp;
+            ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]);
+
+        for(int c2=N/2;c2<N;c2++)
+          //tmp = conjugate(ret._internal[c1][c2]);
+         // ret._internal[c1+N/2][c2-N/2] = -tmp;
+            ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);
+
+        // Remove c1 from rows c1+1...N/2-1
+        for (int b=c1+1; b<N/2; ++b){
+          decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
+          zeroit(pr);
+          for(int c=0; c<N; ++c)
+            pr += conjugate(ret._internal[c1][c])*ret._internal[b][c];
+
+          for(int c=0; c<N; ++c){
+            ret._internal[b][c] -= pr * ret._internal[c1][c];
+          }
+        }
+
+        // Remove c1+N/2 from rows c1+1...N/2-1
+        for (int b=c1+1; b<N; ++b){
+          decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
+          zeroit(pr);
+          for(int c=0; c<N; ++c)
+            pr += conjugate(ret._internal[c1+N/2][c])*ret._internal[b][c];
+
+          for(int c=0; c<N; ++c){
+            ret._internal[b][c] -= pr * ret._internal[c1+N/2][c];
+          }
+        }
+      } //end for loop over c1
+
+      // Compute the last row
+      {
+      // Normalise last row
+        int c1 = N/2-1;
+        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;
+        }
+      // Compute row cN-1+N2/2: cN-1+N/2 = - \Omega cN-1*
+        for(int c2=0;c2<N/2;c2++)
+        {
+          //tmp = conjugate(ret._internal[c1][c2]);
+          // ret._internal[c1+N/2][c2+N/2] = tmp;
+             ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]);
+        }
+        for(int c2=N/2;c2<N;c2++)
+        {
+          //tmp = conjugate(ret._internal[c1][c2]);
+          //ret._internal[c1+N/2][c2-N/2] = -tmp;
+            ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);
+        }
+      }
+      // assuming the determinant is ok
+      return ret;
+  }
+
 NAMESPACE_END(Grid);
 
 #endif