Merge pull request #36 from LupoA/sp2n/gpu-bugfix

Sp2n/gpu bugfix [close #30]

Grid/qcd/action/gauge/GaugeImplTypes.h

@@ -145,7 +145,7 @@ public:
     accelerator_for(ss, P.Grid()->oSites(),1,{
       for (int mu = 0; mu < Nd; mu++) {
           U_v[ss](mu) = Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu);
-          Group::ProjectOnGeneralGroup(U_v[ss](mu));
+          U_v[ss](mu) = Group::ProjectOnGeneralGroup(U_v[ss](mu));
       }
     });
    //auto end = std::chrono::high_resolution_clock::now();

Grid/qcd/utils/GaugeGroup.h

@@ -326,33 +326,33 @@ class GaugeGroup {
   static void ProjectOnGeneralGroup(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >, Nd> > &U) {
     for (int mu = 0; mu < Nd; mu++) {
       auto Umu = PeekIndex<LorentzIndex>(U, mu);
-      ProjectOnGeneralGroup(Umu);
+      Umu = ProjectOnGeneralGroup(Umu);
     }
   }
        
   template <int N>
-  static void ProjectOnGeneralGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
+  static Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > ProjectOnGeneralGroup(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
-      ProjectOnGeneralGroup(Umu, group_name());
+      return ProjectOnGeneralGroup(Umu, group_name());
   }
     
   template <class vtype>
-  static void ProjectOnGeneralGroup(iScalar<vtype> &r) {
+  accelerator_inline static iScalar<vtype> ProjectOnGeneralGroup(const iScalar<vtype> &r) {
-      ProjectOnGeneralGroup(r, group_name());
+      return ProjectOnGeneralGroup(r, group_name());
   }
 
   template <class vtype, int N>
-  static void ProjectOnGeneralGroup(iVector<vtype,N> &r) {
+  accelerator_inline static iVector<vtype,N> ProjectOnGeneralGroup(const iVector<vtype,N> &r) {
-      ProjectOnGeneralGroup(r, group_name());
+      return ProjectOnGeneralGroup(r, group_name());
   }
 
   template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
-  static void ProjectOnGeneralGroup(iMatrix<vtype,N> &arg) {
+  accelerator_inline static iMatrix<vtype,N> ProjectOnGeneralGroup(const iMatrix<vtype,N> &arg) {
-      ProjectOnGeneralGroup(arg, group_name());
+      return ProjectOnGeneralGroup(arg, group_name());
   }
 
   template <int N>       // Projects on SU(N), Sp(2N)
   static void ProjectOnSpecialGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu) {
-       ProjectOnGeneralGroup(Umu);
+       Umu = ProjectOnGeneralGroup(Umu);
        auto det = Determinant(Umu);
 
        det = conjugate(det);

Grid/qcd/utils/SUn.impl

@@ -539,23 +539,24 @@ Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > Inverse(const Lattice<iScala
 }
 
 template <int N>
-static void ProjectOnGeneralGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::SU) {
+static Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > >
-  Umu = ProjectOnGroup(Umu);
+ProjectOnGeneralGroup(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::SU) {
+  return ProjectOnGroup(Umu);
 }
 
 template <class vtype>
-static void ProjectOnGeneralGroup(iScalar<vtype> &r, GroupName::SU) {
+accelerator_inline static iScalar<vtype> ProjectOnGeneralGroup(const iScalar<vtype> &r, GroupName::SU) {
-  r = ProjectOnGroup(r);
+  return ProjectOnGroup(r);
 }
 
 template <class vtype, int N>
-static void ProjectOnGeneralGroup(iVector<vtype,N> &r, GroupName::SU) {
+accelerator_inline static iVector<vtype,N> ProjectOnGeneralGroup(const iVector<vtype,N> &r, GroupName::SU) {
-  r = ProjectOnGroup(r);
+  return ProjectOnGroup(r);
 }
 
 template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
-static void ProjectOnGeneralGroup(iMatrix<vtype,N> &arg, GroupName::SU) {
+accelerator_inline static iMatrix<vtype,N> ProjectOnGeneralGroup(const iMatrix<vtype,N> &arg, GroupName::SU) {
-  arg = ProjectOnGroup(arg);
+  return ProjectOnGroup(arg);
 }
 
 template <typename LatticeMatrixType>

Grid/qcd/utils/Sp2n.impl

@@ -254,23 +254,24 @@ static void testGenerators(GroupName::Sp) {
 }
 
 template <int N>
-static void ProjectOnGeneralGroup(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::Sp) {
+static Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > >
-  Umu = ProjectOnSpGroup(Umu);
+ProjectOnGeneralGroup(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::Sp) {
+  return ProjectOnSpGroup(Umu);
 }
 
 template <class vtype>
-static void ProjectOnGeneralGroup(iScalar<vtype> &r, GroupName::Sp) {
+accelerator_inline static iScalar<vtype> ProjectOnGeneralGroup(const iScalar<vtype> &r, GroupName::Sp) {
-  r = ProjectOnSpGroup(r);
+  return ProjectOnSpGroup(r);
 }
 
 template <class vtype, int N>
-static void ProjectOnGeneralGroup(iVector<vtype,N> &r, GroupName::Sp) {
+accelerator_inline static iVector<vtype,N> ProjectOnGeneralGroup(const iVector<vtype,N> &r, GroupName::Sp) {
-  r = ProjectOnSpGroup(r);
+  return ProjectOnSpGroup(r);
 }
 
 template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
-static void ProjectOnGeneralGroup(iMatrix<vtype,N> &arg, GroupName::Sp) {
+accelerator_inline static iMatrix<vtype,N> ProjectOnGeneralGroup(const iMatrix<vtype,N> &arg, GroupName::Sp) {
-  arg = ProjectOnSpGroup(arg);
+  return ProjectOnSpGroup(arg);
 }
 
 template <typename LatticeMatrixType>   

tests/sp2n/Test_2as_base.cc

@@ -50,7 +50,6 @@ static void run_base_checks() {
       }
 
       RealD realA;
-      RealD realB;
 
       std::cout << GridLogMessage << "checking base has symmetry " << S << std::endl;
       for (int a=0; a < this_irrep_dim; a++)
@@ -101,7 +100,6 @@ static void run_base_checks() {
       }
       for (int gen_id = 0; gen_id < this_algebra_dim; gen_id++)
       {
-          Complex iTr;
           sum = 0;
           sum_im = 0;
           std::cout << GridLogMessage <<  "generator number " << gen_id << std::endl;

tests/sp2n/Test_project_on_Sp.cc

@@ -93,7 +93,7 @@ void test_group_projections(T U) {
   std::cout << GridLogMessage << "Apply to deformed matrix" << std::endl;
 
   U = U + Delta * identity;
-  Sp<Nc>::ProjectOnGeneralGroup(U);
+  U = Sp<Nc>::ProjectOnGeneralGroup(U);
   assert(is_element_of_sp2n_group(U));
 
   name = "ProjectOnSpecialGroup";

tests/sp2n/Test_sp2n_lie_gen.cc

@@ -1,53 +1,65 @@
 #include <Grid/Grid.h>
 
+#include <iostream>
+
 using namespace Grid;
 
-template<int ncolour>
+template <int ngroup>
-void run_checks(bool print_generators = 0) {
+std::ostream& operator<<(std::ostream& o, Sp<ngroup> g) {
-    std::cout << GridLogMessage << "*********************************************"
+  return o << "Sp(" << ngroup << ") Fundamental";
-                << std::endl;
-    std::cout << GridLogMessage << "* Generators for Sp(" << ncolour << ")" << "Fundamental" << std::endl;
-    std::cout << GridLogMessage << "*********************************************"
-                << std::endl;
-    
-    if (print_generators)
-    {
-        Sp<ncolour>::printGenerators();
-    }
-    Sp<ncolour>::testGenerators();
-    
-    if (Sp_TwoIndex<ncolour, Symmetric>::Dimension > 1) {
-        std::cout << GridLogMessage << "*********************************************"
-                    << std::endl;
-        std::cout << GridLogMessage << "* Generators for Sp(" << ncolour << ")" << "TwoIndex Symmetric: " << std::endl;
-        std::cout << GridLogMessage << "*********************************************"
-                    << std::endl;
-        if (print_generators) {
-            Sp_TwoIndex<ncolour, Symmetric>::printGenerators();
-        }
-        Sp_TwoIndex<ncolour, Symmetric>::testGenerators();
-    }
-    
-    if (Sp_TwoIndex<ncolour, AntiSymmetric>::Dimension > 1) {
-        std::cout << GridLogMessage << "*********************************************"
-                    << std::endl;
-        std::cout << GridLogMessage << "* Generators for Sp(" << ncolour << ")" << "TwoIndex AntiSymmetric: " << std::endl;
-        std::cout << GridLogMessage << "*********************************************"
-                    << std::endl;
-        if (print_generators) {
-            Sp_TwoIndex<ncolour, AntiSymmetric>::printGenerators();
-        }
-        Sp_TwoIndex<ncolour, AntiSymmetric>::testGenerators();
-    }
 }
-  
+
+template <int ngroup, TwoIndexSymmetry S>
+std::ostream& operator<<(std::ostream& o, Sp_TwoIndex<ngroup, S> g) {
+  return o << "Sp(" << ngroup << ") TwoIndex "
+           << (S == Symmetric ? "Symmetric" : "AntiSymmetric");
+}
+
+template <class Group>
+void run_check_on(bool print_generators = false) {
+  std::cout << GridLogMessage << "*********************************************"
+            << std::endl;
+  std::cout << GridLogMessage << "* Generators for " << Group() << std::endl;
+  std::cout << GridLogMessage << "*********************************************"
+            << std::endl;
+
+  if (print_generators) {
+    Group::printGenerators();
+  }
+  Group::testGenerators();
+}
+
+template <int ngroup>
+void run_checks() {
+  run_check_on<Sp<ngroup>>();
+  run_check_on<Sp_TwoIndex<ngroup, Symmetric>>();
+  run_check_on<Sp_TwoIndex<ngroup, AntiSymmetric>>();
+}
+
+template <>
+void run_checks<2>() {
+  // Print generators because they are small enough to be actually helpful.
+  run_check_on<Sp<2>>(true);
+  run_check_on<Sp_TwoIndex<2, Symmetric>>(true);
+  // The AntiSymmetric representation is 0 dimensional. This makes problems in
+  // device code.
+}
+
+template <>
+void run_checks<4>() {
+  // Print generators because they are small enough to be actually helpful.
+  run_check_on<Sp<4>>(true);
+  run_check_on<Sp_TwoIndex<4, Symmetric>>(true);
+  run_check_on<Sp_TwoIndex<4, AntiSymmetric>>(true);
+}
+
 int main(int argc, char** argv) {
   Grid_init(&argc, &argv);
-    
+
-  run_checks<2>(1); //  check and print Nc=2
+  run_checks<2>();
-  run_checks<4>(1); //  check and print Nc=4
+  run_checks<4>();
-  run_checks<6>();  //  check Nc=6
+  run_checks<6>();
-  run_checks<8>();  //  check Nc=8
+  run_checks<8>();
 
   Grid_finalize();
 }