Flop count set as in DiRAC-ITT-2020 (mistaken 20% low, but must maintain consistency)

Added untraced baryon contraction code

Grid/lattice/Lattice.h

@@ -36,7 +36,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/lattice/Lattice_local.h>
 #include <Grid/lattice/Lattice_reduction.h>
 #include <Grid/lattice/Lattice_peekpoke.h>
-//#include <Grid/lattice/Lattice_reality.h>
+#include <Grid/lattice/Lattice_reality.h>
 #include <Grid/lattice/Lattice_real_imag.h>
 #include <Grid/lattice/Lattice_comparison_utils.h>
 #include <Grid/lattice/Lattice_comparison.h>

Grid/lattice/Lattice_ET.h

@@ -342,19 +342,14 @@ inline void ExpressionViewClose(LatticeTrinaryExpression<Op, T1, T2, T3> &expr)
 
 GridUnopClass(UnarySub, -a);
 GridUnopClass(UnaryNot, Not(a));
-GridUnopClass(UnaryAdj, adj(a));
-GridUnopClass(UnaryConj, conjugate(a));
 GridUnopClass(UnaryTrace, trace(a));
 GridUnopClass(UnaryTranspose, transpose(a));
 GridUnopClass(UnaryTa, Ta(a));
 GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a));
-GridUnopClass(UnaryToReal, toReal(a));
-GridUnopClass(UnaryToComplex, toComplex(a));
 GridUnopClass(UnaryTimesI, timesI(a));
 GridUnopClass(UnaryTimesMinusI, timesMinusI(a));
 GridUnopClass(UnaryAbs, abs(a));
 GridUnopClass(UnarySqrt, sqrt(a));
-GridUnopClass(UnaryRsqrt, rsqrt(a));
 GridUnopClass(UnarySin, sin(a));
 GridUnopClass(UnaryCos, cos(a));
 GridUnopClass(UnaryAsin, asin(a));
@@ -456,20 +451,17 @@ GridTrinOpClass(TrinaryWhere,
 GRID_DEF_UNOP(operator-, UnarySub);
 GRID_DEF_UNOP(Not, UnaryNot);
 GRID_DEF_UNOP(operator!, UnaryNot);
-GRID_DEF_UNOP(adj, UnaryAdj);
-GRID_DEF_UNOP(conjugate, UnaryConj);
+//GRID_DEF_UNOP(adj, UnaryAdj);
+//GRID_DEF_UNOP(conjugate, UnaryConj);
 GRID_DEF_UNOP(trace, UnaryTrace);
 GRID_DEF_UNOP(transpose, UnaryTranspose);
 GRID_DEF_UNOP(Ta, UnaryTa);
 GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup);
-GRID_DEF_UNOP(toReal, UnaryToReal);
-GRID_DEF_UNOP(toComplex, UnaryToComplex);
 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
                                // abs-fabs-dabs-labs thing
 GRID_DEF_UNOP(sqrt, UnarySqrt);
-GRID_DEF_UNOP(rsqrt, UnaryRsqrt);
 GRID_DEF_UNOP(sin, UnarySin);
 GRID_DEF_UNOP(cos, UnaryCos);
 GRID_DEF_UNOP(asin, UnaryAsin);
@@ -494,27 +486,27 @@ GRID_DEF_TRINOP(where, TrinaryWhere);
 /////////////////////////////////////////////////////////////
 template <class Op, class T1>
 auto closure(const LatticeUnaryExpression<Op, T1> &expr)
-  -> Lattice<decltype(expr.op.func(vecEval(0, expr.arg1)))> 
+  -> Lattice<typename std::remove_const<decltype(expr.op.func(vecEval(0, expr.arg1)))>::type > 
 {
-  Lattice<decltype(expr.op.func(vecEval(0, expr.arg1)))> ret(expr);
+  Lattice<typename std::remove_const<decltype(expr.op.func(vecEval(0, expr.arg1)))>::type > ret(expr);
   return ret;
 }
 template <class Op, class T1, class T2>
 auto closure(const LatticeBinaryExpression<Op, T1, T2> &expr)
-  -> Lattice<decltype(expr.op.func(vecEval(0, expr.arg1),vecEval(0, expr.arg2)))> 
+  -> Lattice<typename std::remove_const<decltype(expr.op.func(vecEval(0, expr.arg1),vecEval(0, expr.arg2)))>::type >
 {
-  Lattice<decltype(expr.op.func(vecEval(0, expr.arg1),vecEval(0, expr.arg2)))> ret(expr);
+  Lattice<typename std::remove_const<decltype(expr.op.func(vecEval(0, expr.arg1),vecEval(0, expr.arg2)))>::type > ret(expr);
   return ret;
 }
 template <class Op, class T1, class T2, class T3>
 auto closure(const LatticeTrinaryExpression<Op, T1, T2, T3> &expr)
-  -> Lattice<decltype(expr.op.func(vecEval(0, expr.arg1),
+  -> Lattice<typename std::remove_const<decltype(expr.op.func(vecEval(0, expr.arg1),
 				   vecEval(0, expr.arg2),
-				   vecEval(0, expr.arg3)))> 
+				   vecEval(0, expr.arg3)))>::type >
 {
-  Lattice<decltype(expr.op.func(vecEval(0, expr.arg1),
+  Lattice<typename std::remove_const<decltype(expr.op.func(vecEval(0, expr.arg1),
 				vecEval(0, expr.arg2),
-			        vecEval(0, expr.arg3)))>  ret(expr);
+			        vecEval(0, expr.arg3)))>::type >  ret(expr);
   return ret;
 }
 #define EXPRESSION_CLOSURE(function)					\

Grid/lattice/Lattice_basis.h

@@ -62,7 +62,7 @@ void basisRotate(VField &basis,Matrix& Qt,int j0, int j1, int k0,int k1,int Nm)
     basis_v.push_back(basis[k].View(AcceleratorWrite));
   }
 
-#if ( (!defined(GRID_SYCL)) && (!defined(GRID_CUDA)) && (!defined(GRID_HIP)) )
+#if ( (!defined(GRID_SYCL)) && (!defined(GRID_CUDA)) )
   int max_threads = thread_max();
   Vector < vobj > Bt(Nm * max_threads);
   thread_region
@@ -161,11 +161,12 @@ void basisRotateJ(Field &result,std::vector<Field> &basis,Eigen::MatrixXd& Qt,in
   double * Qt_j = & Qt_jv[0];
   for(int k=0;k<Nm;++k) Qt_j[k]=Qt(j,k);
 
+  auto basis_vp=& basis_v[0];
   autoView(result_v,result,AcceleratorWrite);
   accelerator_for(ss, grid->oSites(),vobj::Nsimd(),{
     auto B=coalescedRead(zz);
     for(int k=k0; k<k1; ++k){
-      B +=Qt_j[k] * coalescedRead(basis_v[k][ss]);
+      B +=Qt_j[k] * coalescedRead(basis_vp[k][ss]);
     }
     coalescedWrite(result_v[ss], B);
   });

Grid/lattice/Lattice_reality.h

@@ -45,8 +45,8 @@ template<class vobj> inline Lattice<vobj> adj(const Lattice<vobj> &lhs){
   autoView( ret_v, ret, AcceleratorWrite);
 
   ret.Checkerboard()=lhs.Checkerboard();
-  accelerator_for( ss, lhs_v.size(), vobj::Nsimd(), {
-    coalescedWrite(ret_v[ss], adj(lhs_v(ss)));
+  accelerator_for( ss, lhs_v.size(), 1, {
+     ret_v[ss] = adj(lhs_v[ss]);
   });
   return ret;
 };
@@ -64,6 +64,53 @@ template<class vobj> inline Lattice<vobj> conjugate(const Lattice<vobj> &lhs){
   return ret;
 };
 
+template<class vobj> inline Lattice<typename vobj::Complexified> toComplex(const Lattice<vobj> &lhs){
+  Lattice<typename vobj::Complexified> ret(lhs.Grid());
+
+  autoView( lhs_v, lhs, AcceleratorRead);
+  autoView( ret_v, ret, AcceleratorWrite);
+
+  ret.Checkerboard() = lhs.Checkerboard();
+  accelerator_for( ss, lhs_v.size(), 1, {
+    ret_v[ss] = toComplex(lhs_v[ss]);
+  });
+  return ret;
+};
+template<class vobj> inline Lattice<typename vobj::Realified> toReal(const Lattice<vobj> &lhs){
+  Lattice<typename vobj::Realified> ret(lhs.Grid());
+
+  autoView( lhs_v, lhs, AcceleratorRead);
+  autoView( ret_v, ret, AcceleratorWrite);
+
+  ret.Checkerboard() = lhs.Checkerboard();
+  accelerator_for( ss, lhs_v.size(), 1, {
+    ret_v[ss] = toReal(lhs_v[ss]);
+  });
+  return ret;
+};
+
+
+template<class Expression,typename std::enable_if<is_lattice_expr<Expression>::value,void>::type * = nullptr> 
+auto toComplex(const Expression &expr)  -> decltype(closure(expr)) 
+{
+  return toComplex(closure(expr));
+}
+template<class Expression,typename std::enable_if<is_lattice_expr<Expression>::value,void>::type * = nullptr> 
+auto toReal(const Expression &expr)  -> decltype(closure(expr)) 
+{
+  return toReal(closure(expr));
+}
+template<class Expression,typename std::enable_if<is_lattice_expr<Expression>::value,void>::type * = nullptr> 
+auto adj(const Expression &expr)  -> decltype(closure(expr)) 
+{
+  return adj(closure(expr));
+}
+template<class Expression,typename std::enable_if<is_lattice_expr<Expression>::value,void>::type * = nullptr> 
+auto conjugate(const Expression &expr)  -> decltype(closure(expr)) 
+{
+  return conjugate(closure(expr));
+}
+
 NAMESPACE_END(Grid);
 
 #endif

Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h

@@ -133,14 +133,14 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
   pickCheckerboard(Even, CloverTermEven, CloverTerm);
   pickCheckerboard(Odd, CloverTermOdd, CloverTerm);
 
-  pickCheckerboard(Even, CloverTermDagEven, closure(adj(CloverTerm)));
-  pickCheckerboard(Odd, CloverTermDagOdd, closure(adj(CloverTerm)));
+  pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm));
+  pickCheckerboard(Odd, CloverTermDagOdd, adj(CloverTerm));
 
   pickCheckerboard(Even, CloverTermInvEven, CloverTermInv);
   pickCheckerboard(Odd, CloverTermInvOdd, CloverTermInv);
 
-  pickCheckerboard(Even, CloverTermInvDagEven, closure(adj(CloverTermInv)));
-  pickCheckerboard(Odd, CloverTermInvDagOdd, closure(adj(CloverTermInv)));
+  pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv));
+  pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv));
 }
 
 template <class Impl>

Grid/qcd/utils/BaryonUtils.h

@@ -51,7 +51,7 @@ public:
 
   private: 
   template <class mobj, class robj>
-  static void baryon_site(const mobj &D1,
+  static void BaryonSite(const mobj &D1,
 				 const mobj &D2,
 				 const mobj &D3,
 				 const Gamma GammaA_left,
@@ -61,8 +61,18 @@ public:
 				 const int parity,
 				 const bool * wick_contractions,
   				 robj &result);
+  template <class mobj, class robj>
+  static void BaryonSiteMatrix(const mobj &D1,
+         const mobj &D2,
+         const mobj &D3,
+         const Gamma GammaA_left,
+         const Gamma GammaB_left,
+         const Gamma GammaA_right,
+         const Gamma GammaB_right,
+         const bool * wick_contractions,
+           robj &result);
   public:
-  static void Wick_Contractions(std::string qi, 
+  static void WickContractions(std::string qi, 
                  std::string qf, 
                  bool* wick_contractions);
   static void ContractBaryons(const PropagatorField &q1_left,
@@ -75,8 +85,17 @@ public:
 				 const bool* wick_contractions,
 				 const int parity,
 				 ComplexField &baryon_corr);
+  static void ContractBaryonsMatrix(const PropagatorField &q1_left,
+         const PropagatorField &q2_left,
+         const PropagatorField &q3_left,
+         const Gamma GammaA_left,
+         const Gamma GammaB_left,
+         const Gamma GammaA_right,
+         const Gamma GammaB_right,
+         const bool* wick_contractions,
+         SpinMatrixField &baryon_corr);
   template <class mobj, class robj>
-  static void ContractBaryons_Sliced(const mobj &D1,
+  static void ContractBaryonsSliced(const mobj &D1,
 				 const mobj &D2,
 				 const mobj &D3,
 				 const Gamma GammaA_left,
@@ -87,9 +106,20 @@ public:
 				 const int parity,
 				 const int nt,
 				 robj &result);
+  template <class mobj, class robj>
+  static void ContractBaryonsSlicedMatrix(const mobj &D1,
+         const mobj &D2,
+         const mobj &D3,
+         const Gamma GammaA_left,
+         const Gamma GammaB_left,
+         const Gamma GammaA_right,
+         const Gamma GammaB_right,
+         const bool* wick_contractions,
+         const int nt,
+         robj &result);
   private:
   template <class mobj, class mobj2, class robj>
-  static void Baryon_Gamma_3pt_Group1_Site(
+  static void BaryonGamma3ptGroup1Site(
            const mobj &Dq1_ti,
            const mobj2 &Dq2_spec,
            const mobj2 &Dq3_spec,
@@ -101,7 +131,7 @@ public:
            robj &result);
 
   template <class mobj, class mobj2, class robj>
-  static void Baryon_Gamma_3pt_Group2_Site(
+  static void BaryonGamma3ptGroup2Site(
            const mobj2 &Dq1_spec,
            const mobj &Dq2_ti,
            const mobj2 &Dq3_spec,
@@ -113,7 +143,7 @@ public:
            robj &result);
 
   template <class mobj, class mobj2, class robj>
-  static void Baryon_Gamma_3pt_Group3_Site(
+  static void BaryonGamma3ptGroup3Site(
            const mobj2 &Dq1_spec,
            const mobj2 &Dq2_spec,
            const mobj &Dq3_ti,
@@ -125,7 +155,7 @@ public:
            robj &result);
   public:
   template <class mobj>
-  static void Baryon_Gamma_3pt(
+  static void BaryonGamma3pt(
            const PropagatorField &q_ti,
            const mobj &Dq_spec1,
            const mobj &Dq_spec2,
@@ -138,7 +168,7 @@ public:
            SpinMatrixField &stn_corr);
   private: 
   template <class mobj, class mobj2, class robj>
-  static void Sigma_to_Nucleon_Q1_Eye_site(const mobj &Dq_loop,
+  static void SigmaToNucleonQ1EyeSite(const mobj &Dq_loop,
 						 const mobj2 &Du_spec,
 						 const mobj &Dd_tf,
 						 const mobj &Ds_ti,
@@ -147,7 +177,7 @@ public:
 		                 		 const Gamma GammaB_nucl,
 						 robj &result);
   template <class mobj, class mobj2, class robj>
-  static void Sigma_to_Nucleon_Q1_NonEye_site(const mobj &Du_ti,
+  static void SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti,
 						 const mobj &Du_tf,
 						 const mobj2 &Du_spec,
 						 const mobj &Dd_tf,
@@ -159,7 +189,7 @@ public:
 
 
   template <class mobj, class mobj2, class robj>
-  static void Sigma_to_Nucleon_Q2_Eye_site(const mobj &Dq_loop,
+  static void SigmaToNucleonQ2EyeSite(const mobj &Dq_loop,
 						 const mobj2 &Du_spec,
 						 const mobj &Dd_tf,
 						 const mobj &Ds_ti,
@@ -168,7 +198,7 @@ public:
 		                 		 const Gamma GammaB_nucl,
 						 robj &result);
   template <class mobj, class mobj2, class robj>
-  static void Sigma_to_Nucleon_Q2_NonEye_site(const mobj &Du_ti,
+  static void SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,
 						 const mobj &Du_tf,
 						 const mobj2 &Du_spec,
 						 const mobj &Dd_tf,
@@ -179,7 +209,7 @@ public:
 						 robj &result);
   public:
   template <class mobj>
-  static void Sigma_to_Nucleon_Eye(const PropagatorField &qq_loop,
+  static void SigmaToNucleonEye(const PropagatorField &qq_loop,
 				 const mobj &Du_spec,
 				 const PropagatorField &qd_tf,
 				 const PropagatorField &qs_ti,
@@ -189,7 +219,7 @@ public:
 		                 const std::string op,
 				 SpinMatrixField &stn_corr);
   template <class mobj>
-  static void Sigma_to_Nucleon_NonEye(const PropagatorField &qq_ti,
+  static void SigmaToNucleonNonEye(const PropagatorField &qq_ti,
 				 const PropagatorField &qq_tf,
 				 const mobj &Du_spec,
 				 const PropagatorField &qd_tf,
@@ -217,7 +247,7 @@ const Real BaryonUtils<FImpl>::epsilon_sgn[6] = {1.,1.,1.,-1.,-1.,-1.};
 //This is the old version
 template <class FImpl>
 template <class mobj, class robj>
-void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
+void BaryonUtils<FImpl>::BaryonSite(const mobj &D1,
                 const mobj &D2,
                 const mobj &D3,
                          const Gamma GammaA_i,
@@ -329,12 +359,132 @@ void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
     }}
 }
 
+//New version without parity projection or trace
+template <class FImpl>
+template <class mobj, class robj>
+void BaryonUtils<FImpl>::BaryonSiteMatrix(const mobj &D1,
+                const mobj &D2,
+                const mobj &D3,
+                         const Gamma GammaA_i,
+                         const Gamma GammaB_i,
+                         const Gamma GammaA_f,
+                         const Gamma GammaB_f,
+                const bool * wick_contraction,
+                robj &result)
+{
+
+    auto D1_GAi =  D1 * GammaA_i;
+    auto GAf_D1_GAi = GammaA_f * D1_GAi;
+    auto GBf_D1_GAi = GammaB_f * D1_GAi;
+
+    auto D2_GBi = D2 * GammaB_i;
+    auto GBf_D2_GBi = GammaB_f * D2_GBi;
+    auto GAf_D2_GBi = GammaA_f * D2_GBi;
+
+    auto GBf_D3 = GammaB_f * D3;
+    auto GAf_D3 = GammaA_f * D3;
+
+    for (int ie_f=0; ie_f < 6 ; ie_f++){
+        int a_f = epsilon[ie_f][0]; //a
+        int b_f = epsilon[ie_f][1]; //b
+        int c_f = epsilon[ie_f][2]; //c
+    for (int ie_i=0; ie_i < 6 ; ie_i++){
+        int a_i = epsilon[ie_i][0]; //a'
+        int b_i = epsilon[ie_i][1]; //b'
+        int c_i = epsilon[ie_i][2]; //c'
+
+        Real ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i];
+        //This is the \delta_{456}^{123} part
+        if (wick_contraction[0]){
+            for (int rho_i=0; rho_i<Ns; rho_i++){
+            for (int rho_f=0; rho_f<Ns; rho_f++){
+                auto GAf_D1_GAi_rr_cc = GAf_D1_GAi()(rho_f,rho_i)(c_f,c_i);
+                for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+                for (int beta_i=0; beta_i<Ns; beta_i++){
+                    result()(rho_f,rho_i)() += ee  * GAf_D1_GAi_rr_cc
+                                        * D2_GBi    ()(alpha_f,beta_i)(a_f,a_i)
+                                        * GBf_D3    ()(alpha_f,beta_i)(b_f,b_i);
+                }}
+            }}
+        }   
+        //This is the \delta_{456}^{231} part
+        if (wick_contraction[1]){
+            for (int rho_i=0; rho_i<Ns; rho_i++){
+            for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+                auto D1_GAi_ar_ac = D1_GAi()(alpha_f,rho_i)(a_f,c_i);
+                for (int beta_i=0; beta_i<Ns; beta_i++){
+                  auto GBf_D2_GBi_ab_ba = GBf_D2_GBi ()(alpha_f,beta_i)(b_f,a_i);
+                for (int rho_f=0; rho_f<Ns; rho_f++){
+                    result()(rho_f,rho_i)() += ee  * D1_GAi_ar_ac
+                                        * GBf_D2_GBi_ab_ba
+                                        * GAf_D3        ()(rho_f,beta_i)(c_f,b_i);
+                }}
+            }}
+        }   
+        //This is the \delta_{456}^{312} part
+        if (wick_contraction[2]){
+            for (int rho_i=0; rho_i<Ns; rho_i++){
+            for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+                auto GBf_D1_GAi_ar_bc = GBf_D1_GAi()(alpha_f,rho_i)(b_f,c_i);
+                for (int beta_i=0; beta_i<Ns; beta_i++){
+                  auto D3_ab_ab = D3 ()(alpha_f,beta_i)(a_f,b_i);
+                for (int rho_f=0; rho_f<Ns; rho_f++){
+                    result()(rho_f,rho_i)() += ee  * GBf_D1_GAi_ar_bc
+                                        * GAf_D2_GBi    ()(rho_f,beta_i)(c_f,a_i)
+                                        * D3_ab_ab;
+                }}
+            }}
+        }   
+        //This is the \delta_{456}^{132} part
+        if (wick_contraction[3]){
+            for (int rho_i=0; rho_i<Ns; rho_i++){
+            for (int rho_f=0; rho_f<Ns; rho_f++){
+                auto GAf_D1_GAi_rr_cc = GAf_D1_GAi()(rho_f,rho_i)(c_f,c_i);
+                for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+                for (int beta_i=0; beta_i<Ns; beta_i++){
+                    result()(rho_f,rho_i)() -= ee  * GAf_D1_GAi_rr_cc
+                                        * GBf_D2_GBi    ()(alpha_f,beta_i)(b_f,a_i)
+                                        * D3            ()(alpha_f,beta_i)(a_f,b_i);
+                }}
+            }}
+        }   
+        //This is the \delta_{456}^{321} part
+        if (wick_contraction[4]){
+            for (int rho_i=0; rho_i<Ns; rho_i++){
+            for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+                auto GBf_D1_GAi_ar_bc = GBf_D1_GAi()(alpha_f,rho_i)(b_f,c_i);
+                for (int beta_i=0; beta_i<Ns; beta_i++){
+                  auto D2_GBi_ab_aa = D2_GBi()(alpha_f,beta_i)(a_f,a_i);
+                for (int rho_f=0; rho_f<Ns; rho_f++){
+                    result()(rho_f,rho_i)() -= ee  * GBf_D1_GAi_ar_bc
+                                        * D2_GBi_ab_aa
+                                        * GAf_D3    ()(rho_f,beta_i)(c_f,b_i);
+                }}
+            }}
+        }   
+        //This is the \delta_{456}^{213} part
+        if (wick_contraction[5]){
+            for (int rho_i=0; rho_i<Ns; rho_i++){
+            for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+                auto D1_GAi_ar_ac = D1_GAi()(alpha_f,rho_i)(a_f,c_i);
+                for (int beta_i=0; beta_i<Ns; beta_i++){
+                  auto GBf_D3_ab_bb = GBf_D3()(alpha_f,beta_i)(b_f,b_i);
+                for (int rho_f=0; rho_f<Ns; rho_f++){
+                    result()(rho_f,rho_i)() -= ee  * D1_GAi_ar_ac
+                                        * GAf_D2_GBi    ()(rho_f,beta_i)(c_f,a_i)
+                                        * GBf_D3_ab_bb;
+                }}
+            }}
+        }
+    }}
+}
+
 /* Computes which wick contractions should be performed for a    *
  * baryon 2pt function given the initial and finals state quark  *
  * flavours.                                                     *
  * The array wick_contractions must be of length 6               */
 template<class FImpl>
-void BaryonUtils<FImpl>::Wick_Contractions(std::string qi, std::string qf, bool* wick_contractions) {
+void BaryonUtils<FImpl>::WickContractions(std::string qi, std::string qf, bool* wick_contractions) {
     const int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}};
     for (int ie=0; ie < 6 ; ie++) {
         wick_contractions[ie] = (qi.size() == 3 && qf.size() == 3
@@ -365,11 +515,6 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_left,
   assert(Ns==4 && "Baryon code only implemented for N_spin = 4");
   assert(Nc==3 && "Baryon code only implemented for N_colour = 3");
  
-  std::cout << "GammaA (left) " << (GammaA_left.g) <<  std::endl;
-  std::cout << "GammaB (left) " << (GammaB_left.g) <<  std::endl;
-  std::cout << "GammaA (right) " << (GammaA_right.g) <<  std::endl;
-  std::cout << "GammaB (right) " << (GammaB_right.g) <<  std::endl;
- 
   assert(parity==1 || parity == -1 && "Parity must be +1 or -1");
 
   GridBase *grid = q1_left.Grid();
@@ -397,13 +542,62 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_left,
     auto D2 = v2[ss];
     auto D3 = v3[ss];
     vobj result=Zero();
-    baryon_site(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contractions,result);
+    BaryonSite(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contractions,result);
     vbaryon_corr[ss] = result; 
   }  );//end loop over lattice sites
 
   t += usecond();
 
-  std::cout << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s " << std::endl;
+  std::cout << GridLogDebug << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s " << std::endl;
+}
+
+template<class FImpl>
+void BaryonUtils<FImpl>::ContractBaryonsMatrix(const PropagatorField &q1_left,
+             const PropagatorField &q2_left,
+             const PropagatorField &q3_left,
+                         const Gamma GammaA_left,
+                         const Gamma GammaB_left,
+                         const Gamma GammaA_right,
+                         const Gamma GammaB_right,
+             const bool* wick_contractions,
+             SpinMatrixField &baryon_corr)
+{
+
+  assert(Ns==4 && "Baryon code only implemented for N_spin = 4");
+  assert(Nc==3 && "Baryon code only implemented for N_colour = 3");
+ 
+  GridBase *grid = q1_left.Grid();
+  
+  autoView(vbaryon_corr, baryon_corr,CpuWrite);
+  autoView( v1 , q1_left, CpuRead);
+  autoView( v2 , q2_left, CpuRead);
+  autoView( v3 , q3_left, CpuRead);
+
+  // Real bytes =0.;
+  // bytes += grid->oSites() * (432.*sizeof(vComplex) + 126.*sizeof(int) + 36.*sizeof(Real));
+  // for (int ie=0; ie < 6 ; ie++){
+  //   if(ie==0 or ie==3){
+  //      bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contractions[ie];
+  //   }
+  //   else{
+  //      bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contractions[ie];
+  //   }
+  // }
+  // Real t=0.;
+  // t =-usecond();
+
+  accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
+    auto D1 = v1[ss];
+    auto D2 = v2[ss];
+    auto D3 = v3[ss];
+    sobj result=Zero();
+    BaryonSiteMatrix(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,wick_contractions,result);
+    vbaryon_corr[ss] = result; 
+  }  );//end loop over lattice sites
+
+  // t += usecond();
+
+  // std::cout << GridLogDebug << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s " << std::endl;
 
 }
 
@@ -414,7 +608,7 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_left,
  * Wick_Contractions function above                               */
 template <class FImpl>
 template <class mobj, class robj>
-void BaryonUtils<FImpl>::ContractBaryons_Sliced(const mobj &D1,
+void BaryonUtils<FImpl>::ContractBaryonsSliced(const mobj &D1,
 						 const mobj &D2,
 						 const mobj &D3,
 				                 const Gamma GammaA_left,
@@ -430,15 +624,32 @@ void BaryonUtils<FImpl>::ContractBaryons_Sliced(const mobj &D1,
   assert(Ns==4 && "Baryon code only implemented for N_spin = 4");
   assert(Nc==3 && "Baryon code only implemented for N_colour = 3");
  
-  std::cout << "GammaA (left) " << (GammaA_left.g) <<  std::endl;
-  std::cout << "GammaB (left) " << (GammaB_left.g) <<  std::endl;
-  std::cout << "GammaA (right) " << (GammaA_right.g) <<  std::endl;
-  std::cout << "GammaB (right) " << (GammaB_right.g) <<  std::endl;
- 
   assert(parity==1 || parity == -1 && "Parity must be +1 or -1");
 
   for (int t=0; t<nt; t++) {
-    baryon_site(D1[t],D2[t],D3[t],GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contractions,result[t]);
+    BaryonSite(D1[t],D2[t],D3[t],GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contractions,result[t]);
+  }
+}
+
+template <class FImpl>
+template <class mobj, class robj>
+void BaryonUtils<FImpl>::ContractBaryonsSlicedMatrix(const mobj &D1,
+             const mobj &D2,
+             const mobj &D3,
+                         const Gamma GammaA_left,
+                         const Gamma GammaB_left,
+                         const Gamma GammaA_right,
+                         const Gamma GammaB_right,
+             const bool* wick_contractions,
+             const int nt,
+             robj &result)
+{
+
+  assert(Ns==4 && "Baryon code only implemented for N_spin = 4");
+  assert(Nc==3 && "Baryon code only implemented for N_colour = 3");
+
+  for (int t=0; t<nt; t++) {
+    BaryonSiteMatrix(D1[t],D2[t],D3[t],GammaA_left,GammaB_left,GammaA_right,GammaB_right,wick_contractions,result[t]);
   }
 }
 
@@ -454,7 +665,7 @@ void BaryonUtils<FImpl>::ContractBaryons_Sliced(const mobj &D1,
  * Dq4_tf is a quark line from t_f to t_J */
 template<class FImpl>
 template <class mobj, class mobj2, class robj>
-void BaryonUtils<FImpl>::Baryon_Gamma_3pt_Group1_Site(
+void BaryonUtils<FImpl>::BaryonGamma3ptGroup1Site(
                         const mobj &Dq1_ti,
                         const mobj2 &Dq2_spec,
                         const mobj2 &Dq3_spec,
@@ -546,7 +757,7 @@ void BaryonUtils<FImpl>::Baryon_Gamma_3pt_Group1_Site(
  * Dq4_tf is a quark line from t_f to t_J */
 template<class FImpl>
 template <class mobj, class mobj2, class robj>
-void BaryonUtils<FImpl>::Baryon_Gamma_3pt_Group2_Site(
+void BaryonUtils<FImpl>::BaryonGamma3ptGroup2Site(
                         const mobj2 &Dq1_spec,
                         const mobj &Dq2_ti,
                         const mobj2 &Dq3_spec,
@@ -636,7 +847,7 @@ void BaryonUtils<FImpl>::Baryon_Gamma_3pt_Group2_Site(
  * Dq4_tf is a quark line from t_f to t_J */
 template<class FImpl>
 template <class mobj, class mobj2, class robj>
-void BaryonUtils<FImpl>::Baryon_Gamma_3pt_Group3_Site(
+void BaryonUtils<FImpl>::BaryonGamma3ptGroup3Site(
                         const mobj2 &Dq1_spec,
                         const mobj2 &Dq2_spec,
                         const mobj &Dq3_ti,
@@ -728,7 +939,7 @@ void BaryonUtils<FImpl>::Baryon_Gamma_3pt_Group3_Site(
  * https://aportelli.github.io/Hadrons-doc/#/mcontraction        */
 template<class FImpl>
 template <class mobj>
-void BaryonUtils<FImpl>::Baryon_Gamma_3pt(
+void BaryonUtils<FImpl>::BaryonGamma3pt(
                         const PropagatorField &q_ti,
                         const mobj &Dq_spec1,
                         const mobj &Dq_spec2,
@@ -751,7 +962,7 @@ void BaryonUtils<FImpl>::Baryon_Gamma_3pt(
             auto Dq_ti = vq_ti[ss];
             auto Dq_tf = vq_tf[ss];
             sobj result=Zero();
-            Baryon_Gamma_3pt_Group1_Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
+            BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
             vcorr[ss] += result; 
         });//end loop over lattice sites
     } else if (group == 2) {
@@ -759,7 +970,7 @@ void BaryonUtils<FImpl>::Baryon_Gamma_3pt(
             auto Dq_ti = vq_ti[ss];
             auto Dq_tf = vq_tf[ss];
             sobj result=Zero();
-            Baryon_Gamma_3pt_Group2_Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
+            BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
             vcorr[ss] += result; 
         });//end loop over lattice sites
     } else if (group == 3) {
@@ -767,7 +978,7 @@ void BaryonUtils<FImpl>::Baryon_Gamma_3pt(
             auto Dq_ti = vq_ti[ss];
             auto Dq_tf = vq_tf[ss];
             sobj result=Zero();
-            Baryon_Gamma_3pt_Group3_Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
+            BaryonGamma3ptGroup3Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
 
             vcorr[ss] += result; 
         });//end loop over lattice sites
@@ -787,7 +998,7 @@ void BaryonUtils<FImpl>::Baryon_Gamma_3pt(
  * Ds_ti is a quark line from t_i to t_H */
 template <class FImpl>
 template <class mobj, class mobj2, class robj>
-void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q1_Eye_site(const mobj &Dq_loop,
+void BaryonUtils<FImpl>::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop,
 						 const mobj2 &Du_spec,
 						 const mobj &Dd_tf,
 						 const mobj &Ds_ti,
@@ -838,7 +1049,7 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q1_Eye_site(const mobj &Dq_loop,
  * Ds_ti is a quark line from t_i to t_H */
 template <class FImpl>
 template <class mobj, class mobj2, class robj>
-void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q1_NonEye_site(const mobj &Du_ti,
+void BaryonUtils<FImpl>::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti,
 						 const mobj &Du_tf,
 						 const mobj2 &Du_spec,
 						 const mobj &Dd_tf,
@@ -897,7 +1108,7 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q1_NonEye_site(const mobj &Du_ti,
  * Ds_ti is a quark line from t_i to t_H */
 template <class FImpl>
 template <class mobj, class mobj2, class robj>
-void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q2_Eye_site(const mobj &Dq_loop,
+void BaryonUtils<FImpl>::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop,
 						 const mobj2 &Du_spec,
 						 const mobj &Dd_tf,
 						 const mobj &Ds_ti,
@@ -948,7 +1159,7 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q2_Eye_site(const mobj &Dq_loop,
  * Ds_ti is a quark line from t_i to t_H */
 template <class FImpl>
 template <class mobj, class mobj2, class robj>
-void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q2_NonEye_site(const mobj &Du_ti,
+void BaryonUtils<FImpl>::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,
 						 const mobj &Du_tf,
 						 const mobj2 &Du_spec,
 						 const mobj &Dd_tf,
@@ -1002,7 +1213,7 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_Q2_NonEye_site(const mobj &Du_ti,
 
 template<class FImpl>
 template <class mobj>
-void BaryonUtils<FImpl>::Sigma_to_Nucleon_Eye(const PropagatorField &qq_loop,
+void BaryonUtils<FImpl>::SigmaToNucleonEye(const PropagatorField &qq_loop,
 						 const mobj &Du_spec,
 						 const PropagatorField &qd_tf,
 						 const PropagatorField &qs_ti,
@@ -1029,9 +1240,9 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_Eye(const PropagatorField &qq_loop,
     auto Ds_ti = vs_ti[ss];
     sobj result=Zero();
     if(op == "Q1"){
-      Sigma_to_Nucleon_Q1_Eye_site(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
+      SigmaToNucleonQ1EyeSite(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
     } else if(op == "Q2"){
-      Sigma_to_Nucleon_Q2_Eye_site(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
+      SigmaToNucleonQ2EyeSite(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
     } else {
       assert(0 && "Weak Operator not correctly specified");
     }
@@ -1041,7 +1252,7 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_Eye(const PropagatorField &qq_loop,
 
 template<class FImpl>
 template <class mobj>
-void BaryonUtils<FImpl>::Sigma_to_Nucleon_NonEye(const PropagatorField &qq_ti,
+void BaryonUtils<FImpl>::SigmaToNucleonNonEye(const PropagatorField &qq_ti,
 						 const PropagatorField &qq_tf,
 						 const mobj &Du_spec,
 						 const PropagatorField &qd_tf,
@@ -1071,9 +1282,9 @@ void BaryonUtils<FImpl>::Sigma_to_Nucleon_NonEye(const PropagatorField &qq_ti,
     auto Ds_ti = vs_ti[ss];
     sobj result=Zero();
     if(op == "Q1"){
-      Sigma_to_Nucleon_Q1_NonEye_site(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
+      SigmaToNucleonQ1NonEyeSite(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
     } else if(op == "Q2"){
-      Sigma_to_Nucleon_Q2_NonEye_site(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
+      SigmaToNucleonQ2NonEyeSite(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result);
     } else {
       assert(0 && "Weak Operator not correctly specified");
     }

Grid/qcd/utils/SUn.h

@@ -449,7 +449,8 @@ public:
     LatticeReal alpha(grid);
 
     //    std::cout<<GridLogMessage<<"xi "<<xi <<std::endl;
-    alpha = toReal(2.0 * xi);
+    xi = 2.0 *xi;
+    alpha = toReal(xi);
 
     do {
       // A. Generate two uniformly distributed pseudo-random numbers R and R',

Grid/serialisation/JSON_IO.cc

@@ -26,7 +26,7 @@
     *************************************************************************************/
     /*  END LEGAL */
 #include <Grid/Grid.h>
-#ifndef __NVCC__
+#if (!defined(GRID_CUDA)) && (!defined(GRID_HIP))
 
 NAMESPACE_BEGIN(Grid);
 

Grid/simd/Grid_vector_unops.h

@@ -125,14 +125,6 @@ accelerator_inline Grid_simd<S, V> sqrt(const Grid_simd<S, V> &r) {
   return SimdApply(SqrtRealFunctor<S>(), r);
 }
 template <class S, class V>
-accelerator_inline Grid_simd<S, V> rsqrt(const Grid_simd<S, V> &r) {
-  return SimdApply(RSqrtRealFunctor<S>(), r);
-}
-template <class Scalar>
-accelerator_inline Scalar rsqrt(const Scalar &r) {
-  return (RSqrtRealFunctor<Scalar>(), r);
-}
-template <class S, class V>
 accelerator_inline Grid_simd<S, V> cos(const Grid_simd<S, V> &r) {
   return SimdApply(CosRealFunctor<S>(), r);
 }

Grid/tensors/Tensor_Ta.h

@@ -92,17 +92,22 @@ accelerator_inline iMatrix<vtype,N> ProjectOnGroup(const iMatrix<vtype,N> &arg)
 {
   // need a check for the group type?
   iMatrix<vtype,N> ret(arg);
+  vtype rnrm;
   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 = rsqrt(inner);
+    nrm = sqrt(inner);
+    nrm = 1.0/nrm;
     for(int c2=0;c2<N;c2++)
       ret._internal[c1][c2]*= nrm;
       
+    // Remove c1 from rows c1+1...N-1
     for (int b=c1+1; b<N; ++b){
       decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
       zeroit(pr);

Grid/tensors/Tensor_unary.h

@@ -84,7 +84,6 @@ NAMESPACE_BEGIN(Grid);
   }
 
 UNARY(sqrt);
-UNARY(rsqrt);
 UNARY(sin);
 UNARY(cos);
 UNARY(asin);

Grid/threads/Accelerator.h

@@ -333,7 +333,7 @@ inline void *acceleratorAllocDevice(size_t bytes)
   return ptr;
 };
 
-inline void acceleratorFreeShared(void *ptr){ free(ptr);};
+inline void acceleratorFreeShared(void *ptr){ hipFree(ptr);};
 inline void acceleratorFreeDevice(void *ptr){ hipFree(ptr);};
 inline void acceleratorCopyToDevice(void *from,void *to,size_t bytes)  { hipMemcpy(to,from,bytes, hipMemcpyHostToDevice);}
 inline void acceleratorCopyFromDevice(void *from,void *to,size_t bytes){ hipMemcpy(to,from,bytes, hipMemcpyDeviceToHost);}

benchmarks/Benchmark_ITT.cc

@@ -62,7 +62,7 @@ struct time_statistics{
 
 void comms_header(){
   std::cout <<GridLogMessage << " L  "<<"\t"<<" Ls  "<<"\t"
-            <<std::setw(11)<<"bytes"<<"MB/s uni (err/min/max)"<<"\t\t"<<"MB/s bidi (err/min/max)"<<std::endl;
+            <<"bytes\t MB/s uni (err/min/max) \t\t MB/s bidi (err/min/max)"<<std::endl;
 };
 
 Gamma::Algebra Gmu [] = {
@@ -189,11 +189,11 @@ public:
 	//	double rbytes    = dbytes*0.5;
 	double bidibytes = dbytes;
 
-	std::cout<<GridLogMessage << std::setw(4) << lat<<"\t"<<Ls<<"\t"
-		 <<std::setw(11) << bytes<< std::fixed << std::setprecision(1) << std::setw(7)
-		 <<std::right<< xbytes/timestat.mean<<"  "<< xbytes*timestat.err/(timestat.mean*timestat.mean)<< " "
+	std::cout<<GridLogMessage << lat<<"\t"<<Ls<<"\t "
+		 << bytes << " \t "
+		 <<xbytes/timestat.mean<<" \t "<< xbytes*timestat.err/(timestat.mean*timestat.mean)<< " \t "
 		 <<xbytes/timestat.max <<" "<< xbytes/timestat.min  
-		 << "\t\t"<<std::setw(7)<< bidibytes/timestat.mean<< "  " << bidibytes*timestat.err/(timestat.mean*timestat.mean) << " "
+		 << "\t\t"<< bidibytes/timestat.mean<< "  " << bidibytes*timestat.err/(timestat.mean*timestat.mean) << " "
 		 << bidibytes/timestat.max << " " << bidibytes/timestat.min << std::endl;
 	
 	    }
@@ -512,7 +512,6 @@ public:
     NN_global=NN;
     uint64_t SHM=NP/NN;
 
-    Coordinate latt4({local[0]*mpi[0],local[1]*mpi[1],local[2]*mpi[2],local[3]*mpi[3]});
 
     ///////// Welcome message ////////////
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
@@ -701,7 +700,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
   std::cout<<GridLogMessage << " Summary table Ls="<<Ls <<std::endl;
   std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  std::cout<<GridLogMessage << "L \t\t Wilson \t\t DWF4 \t\tt Staggered" <<std::endl;
+  std::cout<<GridLogMessage << "L \t\t Wilson \t\t DWF4 \t\t Staggered" <<std::endl;
   for(int l=0;l<L_list.size();l++){
     std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< wilson[l]<<" \t\t "<<dwf4[l] << " \t\t "<< staggered[l]<<std::endl;
   }
@@ -732,9 +731,9 @@ int main (int argc, char ** argv)
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
     std::cout<<GridLogMessage << " Per Node Summary table Ls="<<Ls <<std::endl;
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << " L \t\t Wilson\t\t DWF4  " <<std::endl;
+    std::cout<<GridLogMessage << " L \t\t Wilson\t\t DWF4\t\t Staggered " <<std::endl;
     for(int l=0;l<L_list.size();l++){
-      std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< wilson[l]/NN<<" \t "<<dwf4[l]/NN<<std::endl;
+      std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< wilson[l]/NN<<" \t "<<dwf4[l]/NN<< " \t "<<staggered[l]/NN<<std::endl;
     }
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
 

tests/core/Test_lie_generators.cc

@@ -69,11 +69,11 @@ int main(int argc, char** argv) {
   std::cout << GridLogMessage << "* Generators for SU(Nc" << std::endl;
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
-  SU<Nc>::printGenerators();
-  std::cout << "Dimension of adjoint representation: "<< SU<Nc>Adjoint::Dimension << std::endl;
-  SU<Nc>Adjoint::printGenerators();
-  SU<Nc>::testGenerators();
-  SU<Nc>Adjoint::testGenerators();
+  SU3::printGenerators();
+  std::cout << "Dimension of adjoint representation: "<< SU3Adjoint::Dimension << std::endl;
+  SU3Adjoint::printGenerators();
+  SU3::testGenerators();
+  SU3Adjoint::testGenerators();
 
   std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
   std::cout<<GridLogMessage<<"* Generators for SU(4)"<<std::endl;
@@ -87,22 +87,22 @@ int main(int argc, char** argv) {
   // Projectors 
   GridParallelRNG gridRNG(grid);
   gridRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
-  SU<Nc>Adjoint::LatticeAdjMatrix Gauss(grid);
-  SU<Nc>::LatticeAlgebraVector ha(grid);
-  SU<Nc>::LatticeAlgebraVector hb(grid);
+  SU3Adjoint::LatticeAdjMatrix Gauss(grid);
+  SU3::LatticeAlgebraVector ha(grid);
+  SU3::LatticeAlgebraVector hb(grid);
   random(gridRNG,Gauss);
 
   std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl;
-  SU<Nc>Adjoint::projectOnAlgebra(ha, Gauss);
+  SU3Adjoint::projectOnAlgebra(ha, Gauss);
   std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl;
   std::cout << GridLogMessage << "Start projector" << std::endl;
-  SU<Nc>Adjoint::projector(hb, Gauss);
+  SU3Adjoint::projector(hb, Gauss);
   std::cout << GridLogMessage << "end projector" << std::endl;
 
   std::cout << GridLogMessage << "ReStart projector" << std::endl;
-  SU<Nc>Adjoint::projector(hb, Gauss);
+  SU3Adjoint::projector(hb, Gauss);
   std::cout << GridLogMessage << "end projector" << std::endl;
-  SU<Nc>::LatticeAlgebraVector diff = ha -hb;
+  SU3::LatticeAlgebraVector diff = ha -hb;
   std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl;
 
 
@@ -114,8 +114,8 @@ int main(int argc, char** argv) {
 
   
   LatticeGaugeField U(grid), V(grid);
-  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U);
-  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V);
+  SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U);
+  SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V);
 
   // Adjoint representation
   // Test group structure
@@ -123,8 +123,8 @@ int main(int argc, char** argv) {
   LatticeGaugeField UV(grid);
   UV = Zero();
   for (int mu = 0; mu < Nd; mu++) {
-    SU<Nc>::LatticeMatrix Umu = peekLorentz(U,mu);
-    SU<Nc>::LatticeMatrix Vmu = peekLorentz(V,mu);
+    SU3::LatticeMatrix Umu = peekLorentz(U,mu);
+    SU3::LatticeMatrix Vmu = peekLorentz(V,mu);
     pokeLorentz(UV,Umu*Vmu, mu);
   }
 
@@ -151,16 +151,16 @@ int main(int argc, char** argv) {
 
   // Check correspondence of algebra and group transformations
   // Create a random vector
-  SU<Nc>::LatticeAlgebraVector h_adj(grid);
+  SU3::LatticeAlgebraVector h_adj(grid);
   typename AdjointRep<Nc>::LatticeMatrix Ar(grid);
   random(gridRNG,h_adj);
   h_adj = real(h_adj);
   SU_Adjoint<Nc>::AdjointLieAlgebraMatrix(h_adj,Ar);
 
   // Re-extract h_adj
-  SU<Nc>::LatticeAlgebraVector h_adj2(grid);
+  SU3::LatticeAlgebraVector h_adj2(grid);
   SU_Adjoint<Nc>::projectOnAlgebra(h_adj2, Ar);
-  SU<Nc>::LatticeAlgebraVector h_diff = h_adj - h_adj2;
+  SU3::LatticeAlgebraVector h_diff = h_adj - h_adj2;
   std::cout << GridLogMessage << "Projections structure check vector difference (Adjoint representation) : " << norm2(h_diff) << std::endl;
 
   // Exponentiate
@@ -183,14 +183,14 @@ int main(int argc, char** argv) {
       
 
   // Construct the fundamental matrix in the group
-  SU<Nc>::LatticeMatrix Af(grid);
-  SU<Nc>::FundamentalLieAlgebraMatrix(h_adj,Af);
-  SU<Nc>::LatticeMatrix Ufund(grid);
+  SU3::LatticeMatrix Af(grid);
+  SU3::FundamentalLieAlgebraMatrix(h_adj,Af);
+  SU3::LatticeMatrix Ufund(grid);
   Ufund  = expMat(Af, 1.0, 16);
   // Check unitarity
-  SU<Nc>::LatticeMatrix uno_f(grid);
+  SU3::LatticeMatrix uno_f(grid);
   uno_f = 1.0;
-  SU<Nc>::LatticeMatrix UnitCheck(grid);
+  SU3::LatticeMatrix UnitCheck(grid);
   UnitCheck = Ufund * adj(Ufund) - uno_f;
   std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck)
             << std::endl;
@@ -260,20 +260,20 @@ int main(int argc, char** argv) {
   std::cout << GridLogMessage << "Test for the Two Index Symmetric projectors"
       << std::endl;
   // Projectors 
-  SU<Nc>TwoIndexSymm::LatticeTwoIndexMatrix Gauss2(grid);
+  SU3TwoIndexSymm::LatticeTwoIndexMatrix Gauss2(grid);
   random(gridRNG,Gauss2);
   
   std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl;
-  SU<Nc>TwoIndexSymm::projectOnAlgebra(ha, Gauss2);
+  SU3TwoIndexSymm::projectOnAlgebra(ha, Gauss2);
   std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl;
   std::cout << GridLogMessage << "Start projector" << std::endl;
-  SU<Nc>TwoIndexSymm::projector(hb, Gauss2);
+  SU3TwoIndexSymm::projector(hb, Gauss2);
   std::cout << GridLogMessage << "end projector" << std::endl;
   
   std::cout << GridLogMessage << "ReStart projector" << std::endl;
-  SU<Nc>TwoIndexSymm::projector(hb, Gauss2);
+  SU3TwoIndexSymm::projector(hb, Gauss2);
   std::cout << GridLogMessage << "end projector" << std::endl;
-  SU<Nc>::LatticeAlgebraVector diff2 = ha - hb;
+  SU3::LatticeAlgebraVector diff2 = ha - hb;
   std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl;
   std::cout << GridLogMessage << "*********************************************"
       << std::endl;
@@ -284,20 +284,20 @@ int main(int argc, char** argv) {
   std::cout << GridLogMessage << "Test for the Two index anti-Symmetric projectors"
       << std::endl;
   // Projectors
-  SU<Nc>TwoIndexAntiSymm::LatticeTwoIndexMatrix Gauss2a(grid);
+  SU3TwoIndexAntiSymm::LatticeTwoIndexMatrix Gauss2a(grid);
   random(gridRNG,Gauss2a);
   
   std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl;
-  SU<Nc>TwoIndexAntiSymm::projectOnAlgebra(ha, Gauss2a);
+  SU3TwoIndexAntiSymm::projectOnAlgebra(ha, Gauss2a);
   std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl;
   std::cout << GridLogMessage << "Start projector" << std::endl;
-  SU<Nc>TwoIndexAntiSymm::projector(hb, Gauss2a);
+  SU3TwoIndexAntiSymm::projector(hb, Gauss2a);
   std::cout << GridLogMessage << "end projector" << std::endl;
   
   std::cout << GridLogMessage << "ReStart projector" << std::endl;
-  SU<Nc>TwoIndexAntiSymm::projector(hb, Gauss2a);
+  SU3TwoIndexAntiSymm::projector(hb, Gauss2a);
   std::cout << GridLogMessage << "end projector" << std::endl;
-  SU<Nc>::LatticeAlgebraVector diff2a = ha - hb;
+  SU3::LatticeAlgebraVector diff2a = ha - hb;
   std::cout << GridLogMessage << "Difference: " << norm2(diff2a) << std::endl;
   std::cout << GridLogMessage << "*********************************************"
       << std::endl;
@@ -311,14 +311,14 @@ int main(int argc, char** argv) {
   // Test group structure
   // (U_f * V_f)_r = U_r * V_r
   LatticeGaugeField U2(grid), V2(grid);
-  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U2);
-  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V2);
+  SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2);
+  SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2);
   
   LatticeGaugeField UV2(grid);
   UV2 = Zero();
   for (int mu = 0; mu < Nd; mu++) {
-    SU<Nc>::LatticeMatrix Umu2 = peekLorentz(U2,mu);
-    SU<Nc>::LatticeMatrix Vmu2 = peekLorentz(V2,mu);
+    SU3::LatticeMatrix Umu2 = peekLorentz(U2,mu);
+    SU3::LatticeMatrix Vmu2 = peekLorentz(V2,mu);
     pokeLorentz(UV2,Umu2*Vmu2, mu);
   }
   
@@ -345,16 +345,16 @@ int main(int argc, char** argv) {
   
   // Check correspondence of algebra and group transformations
   // Create a random vector
-  SU<Nc>::LatticeAlgebraVector h_sym(grid);
+  SU3::LatticeAlgebraVector h_sym(grid);
   typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ar_sym(grid);
   random(gridRNG,h_sym);
   h_sym = real(h_sym);
   SU_TwoIndex<Nc,Symmetric>::TwoIndexLieAlgebraMatrix(h_sym,Ar_sym);
   
   // Re-extract h_sym
-  SU<Nc>::LatticeAlgebraVector h_sym2(grid);
+  SU3::LatticeAlgebraVector h_sym2(grid);
   SU_TwoIndex< Nc, Symmetric>::projectOnAlgebra(h_sym2, Ar_sym);
-  SU<Nc>::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2;
+  SU3::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2;
   std::cout << GridLogMessage << "Projections structure check vector difference (Two Index Symmetric): " << norm2(h_diff_sym) << std::endl;
 
   
@@ -379,11 +379,11 @@ int main(int argc, char** argv) {
   
   
   // Construct the fundamental matrix in the group
-  SU<Nc>::LatticeMatrix Af_sym(grid);
-  SU<Nc>::FundamentalLieAlgebraMatrix(h_sym,Af_sym);
-  SU<Nc>::LatticeMatrix Ufund2(grid);
+  SU3::LatticeMatrix Af_sym(grid);
+  SU3::FundamentalLieAlgebraMatrix(h_sym,Af_sym);
+  SU3::LatticeMatrix Ufund2(grid);
   Ufund2  = expMat(Af_sym, 1.0, 16);
-  SU<Nc>::LatticeMatrix UnitCheck2(grid);
+  SU3::LatticeMatrix UnitCheck2(grid);
   UnitCheck2 = Ufund2 * adj(Ufund2) - uno_f;
   std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2)
       << std::endl;
@@ -421,14 +421,14 @@ int main(int argc, char** argv) {
   // Test group structure
   // (U_f * V_f)_r = U_r * V_r
   LatticeGaugeField U2A(grid), V2A(grid);
-  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U2A);
-  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V2A);
+  SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2A);
+  SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2A);
   
   LatticeGaugeField UV2A(grid);
   UV2A = Zero();
   for (int mu = 0; mu < Nd; mu++) {
-    SU<Nc>::LatticeMatrix Umu2A = peekLorentz(U2,mu);
-    SU<Nc>::LatticeMatrix Vmu2A = peekLorentz(V2,mu);
+    SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu);
+    SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu);
     pokeLorentz(UV2A,Umu2A*Vmu2A, mu);
   }
   
@@ -455,16 +455,16 @@ int main(int argc, char** argv) {
   
   // Check correspondence of algebra and group transformations
   // Create a random vector
-  SU<Nc>::LatticeAlgebraVector h_Asym(grid);
+  SU3::LatticeAlgebraVector h_Asym(grid);
   typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ar_Asym(grid);
   random(gridRNG,h_Asym);
   h_Asym = real(h_Asym);
   SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym);
   
   // Re-extract h_sym
-  SU<Nc>::LatticeAlgebraVector h_Asym2(grid);
+  SU3::LatticeAlgebraVector h_Asym2(grid);
   SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym);
-  SU<Nc>::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2;
+  SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2;
   std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl;
 
   
@@ -489,11 +489,11 @@ int main(int argc, char** argv) {
   
   
   // Construct the fundamental matrix in the group
-  SU<Nc>::LatticeMatrix Af_Asym(grid);
-  SU<Nc>::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym);
-  SU<Nc>::LatticeMatrix Ufund2A(grid);
+  SU3::LatticeMatrix Af_Asym(grid);
+  SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym);
+  SU3::LatticeMatrix Ufund2A(grid);
   Ufund2A  = expMat(Af_Asym, 1.0, 16);
-  SU<Nc>::LatticeMatrix UnitCheck2A(grid);
+  SU3::LatticeMatrix UnitCheck2A(grid);
   UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f;
   std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A)
       << std::endl;

tests/core/Test_unary.cc

@@ -0,0 +1,106 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_quenched_update.cc
+
+    Copyright (C) 2015
+
+Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+ ;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt({8,8,8,8});
+  GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt, 
+							GridDefaultSimd(Nd,vComplexD::Nsimd()),
+							GridDefaultMpi());
+
+  GridCartesian * gridF = SpaceTimeGrid::makeFourDimGrid(latt, 
+							GridDefaultSimd(Nd,vComplexF::Nsimd()),
+							GridDefaultMpi());
+  
+
+  ///////////////////////////////
+  // Configuration of known size
+  ///////////////////////////////
+  LatticeColourMatrixD ident(grid);
+  LatticeColourMatrixD U(grid);
+  LatticeColourMatrixD tmp(grid);
+  LatticeColourMatrixD org(grid);
+  LatticeColourMatrixF UF(gridF);
+
+  LatticeGaugeField Umu(grid);
+
+  ident =1.0;
+
+  // RNG set up for test
+  std::vector<int> pseeds({1,2,3,4,5}); // once I caught a fish alive
+  std::vector<int> sseeds({6,7,8,9,10});// then i let it go again
+  GridParallelRNG  pRNG(grid); pRNG.SeedFixedIntegers(pseeds);
+  GridSerialRNG    sRNG;       sRNG.SeedFixedIntegers(sseeds);
+
+  SU<Nc>::HotConfiguration(pRNG,Umu);
+
+  U = PeekIndex<LorentzIndex>(Umu,0);
+  org=U;
+
+
+  tmp=  U*adj(U) - ident ;
+  RealD Def1 = norm2( tmp );
+  std::cout << " Defect1 "<<Def1<<std::endl;
+
+  tmp = U - org;
+  std::cout << "Diff1 "<<norm2(tmp)<<std::endl;
+  precisionChange(UF,U);
+  precisionChange(U,UF);
+
+  tmp=  U*adj(U) - ident ;
+  RealD Def2 = norm2(  tmp );
+  std::cout << " Defect2 "<<Def2<<std::endl;
+
+  tmp = U - org;
+  std::cout << "Diff2 "<<norm2(tmp)<<std::endl;
+
+  U = ProjectOnGroup(U);
+
+  tmp=  U*adj(U) - ident ;
+  RealD Def3 = norm2(  tmp);
+  std::cout << " Defect3 "<<Def3<<std::endl;
+
+
+  tmp = U - org;
+  std::cout << "Diff3 "<<norm2(tmp)<<std::endl;
+
+
+  Grid_finalize();
+}
+
+
+
+