Merge branch 'develop' of https://github.com/paboyle/Grid into develop

Grid/communicator/SharedMemoryMPI.cc

@@ -7,6 +7,7 @@
     Copyright (C) 2015
 
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Christoph Lehner <christoph@lhnr.de>
 
     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
@@ -169,6 +170,23 @@ static inline int divides(int a,int b)
 }
 void GlobalSharedMemory::GetShmDims(const Coordinate &WorldDims,Coordinate &ShmDims)
 {
+  ////////////////////////////////////////////////////////////////
+  // Allow user to configure through environment variable
+  ////////////////////////////////////////////////////////////////
+  char* str = getenv(("GRID_SHM_DIMS_" + std::to_string(ShmDims.size())).c_str());
+  if ( str ) {
+    std::vector<int> IntShmDims;
+    GridCmdOptionIntVector(std::string(str),IntShmDims);
+    assert(IntShmDims.size() == WorldDims.size());
+    long ShmSize = 1;
+    for (int dim=0;dim<WorldDims.size();dim++) {
+      ShmSize *= (ShmDims[dim] = IntShmDims[dim]);
+      assert(divides(ShmDims[dim],WorldDims[dim]));
+    }
+    assert(ShmSize == WorldShmSize);
+    return;
+  }
+  
   ////////////////////////////////////////////////////////////////
   // Powers of 2,3,5 only in prime decomposition for now
   ////////////////////////////////////////////////////////////////

Grid/lattice/Lattice_transfer.h

@@ -97,6 +97,20 @@ accelerator_inline void convertType(ComplexF & out, const std::complex<float> &
   out = in;
 }
 
+template<typename T>
+accelerator_inline EnableIf<isGridFundamental<T>> convertType(T & out, const T & in) {
+  out = in;
+}
+
+// This would allow for conversions between GridFundamental types, but is not strictly needed as yet
+/*template<typename T1, typename T2>
+accelerator_inline typename std::enable_if<isGridFundamental<T1>::value && isGridFundamental<T2>::value>::type
+// Or to make this very broad, conversions between anything that's not a GridTensor could be allowed
+//accelerator_inline typename std::enable_if<!isGridTensor<T1>::value && !isGridTensor<T2>::value>::type
+convertType(T1 & out, const T2 & in) {
+  out = in;
+}*/
+
 #ifdef GRID_SIMT
 accelerator_inline void convertType(vComplexF & out, const ComplexF & in) {
   ((ComplexF*)&out)[acceleratorSIMTlane(vComplexF::Nsimd())] = in;
@@ -117,23 +131,18 @@ accelerator_inline void convertType(vComplexD2 & out, const vComplexF & in) {
   Optimization::PrecisionChange::StoD(in.v,out._internal[0].v,out._internal[1].v);
 }
 
-template<typename T1,typename T2,int N>
+template<typename T1,typename T2>
-  accelerator_inline void convertType(iMatrix<T1,N> & out, const iMatrix<T2,N> & in);
+accelerator_inline void convertType(iScalar<T1> & out, const iScalar<T2> & in) {
-template<typename T1,typename T2,int N>
+  convertType(out._internal,in._internal);
-  accelerator_inline void convertType(iVector<T1,N> & out, const iVector<T2,N> & in);
-
-template<typename T1,typename T2, typename std::enable_if<!isGridScalar<T1>::value, T1>::type* = nullptr>
-accelerator_inline void convertType(T1 & out, const iScalar<T2> & in) {
-  convertType(out,in._internal);
 }
 
-template<typename T1, typename std::enable_if<!isGridScalar<T1>::value, T1>::type* = nullptr>
+template<typename T1,typename T2>
-accelerator_inline void convertType(T1 & out, const iScalar<T1> & in) {
+accelerator_inline NotEnableIf<isGridScalar<T1>> convertType(T1 & out, const iScalar<T2> & in) {
   convertType(out,in._internal);
 }
 
 template<typename T1,typename T2>
-accelerator_inline void convertType(iScalar<T1> & out, const T2 & in) {
+accelerator_inline NotEnableIf<isGridScalar<T2>> convertType(iScalar<T1> & out, const T2 & in) {
   convertType(out._internal,in);
 }
 
@@ -150,11 +159,6 @@ accelerator_inline void convertType(iVector<T1,N> & out, const iVector<T2,N> & i
     convertType(out._internal[i],in._internal[i]);
 }
 
-template<typename T, typename std::enable_if<isGridFundamental<T>::value, T>::type* = nullptr>
-accelerator_inline void convertType(T & out, const T & in) {
-  out = in;
-}
-
 template<typename T1,typename T2>
 accelerator_inline void convertType(Lattice<T1> & out, const Lattice<T2> & in) {
   autoView( out_v , out,AcceleratorWrite);

Grid/qcd/action/fermion/WilsonCloverFermion.h

@@ -245,7 +245,7 @@ public:
     return out;
   }
 
-private:
+protected:
   // here fixing the 4 dimensions, make it more general?
 
   RealD csw_r;                                               // Clover coefficient - spatial

Grid/qcd/utils/BaryonUtils.h

@@ -27,7 +27,6 @@
  *************************************************************************************/
 /*  END LEGAL */
 #pragma once
-//#include <Grid/Hadrons/Global.hpp>
 #include <Grid/Eigen/unsupported/CXX11/Tensor>
 
 NAMESPACE_BEGIN(Grid);
@@ -192,14 +191,34 @@ public:
              robj &result);
   template <class mobj, class mobj2, class robj> accelerator_inline
   static void SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,
-             const mobj &Du_tf,
+	     const mobj &Du_tf,
-             const mobj2 &Du_spec,
+	     const mobj2 &Du_spec,
-             const mobj &Dd_tf,
+	     const mobj &Dd_tf,
-             const mobj &Ds_ti,
+	     const mobj &Ds_ti,
-             const Gamma Gamma_H,
+	     const Gamma Gamma_H,
-             const Gamma GammaB_sigma,
+	     const Gamma GammaB_sigma,
-             const Gamma GammaB_nucl,
+	     const Gamma GammaB_nucl,
-             robj &result);
+	     robj &result);
+  template <class mobj, class mobj2, class robj> accelerator_inline
+  static void XiToSigmaQ1EyeSite(const mobj &Dq_loop,
+	     const mobj2 &Dd_spec,
+	     const mobj2 &Ds_spec,
+	     const mobj &Dd_tf,
+	     const mobj &Ds_ti,
+	     const Gamma Gamma_H,
+	     const Gamma GammaB_sigma,
+	     const Gamma GammaB_nucl,
+	     robj &result);
+  template <class mobj, class mobj2, class robj> accelerator_inline
+  static void XiToSigmaQ2EyeSite(const mobj &Dq_loop,
+	     const mobj2 &Dd_spec,
+	     const mobj2 &Ds_spec,
+	     const mobj &Dd_tf,
+	     const mobj &Ds_ti,
+	     const Gamma Gamma_H,
+	     const Gamma GammaB_sigma,
+	     const Gamma GammaB_nucl,
+	     robj &result);
   public:
   template <class mobj>
   static void SigmaToNucleonEye(const PropagatorField &qq_loop,
@@ -213,15 +232,26 @@ public:
          SpinMatrixField &stn_corr);
   template <class mobj>
   static void SigmaToNucleonNonEye(const PropagatorField &qq_ti,
-         const PropagatorField &qq_tf,
+	 const PropagatorField &qq_tf,
-         const mobj &Du_spec,
+	 const mobj &Du_spec,
-         const PropagatorField &qd_tf,
+	 const PropagatorField &qd_tf,
-         const PropagatorField &qs_ti,
+	 const PropagatorField &qs_ti,
-         const Gamma Gamma_H,
+	 const Gamma Gamma_H,
-         const Gamma GammaB_sigma,
+	 const Gamma GammaB_sigma,
-         const Gamma GammaB_nucl,
+	 const Gamma GammaB_nucl,
          const std::string op,
-         SpinMatrixField &stn_corr);
+	 SpinMatrixField &stn_corr);
+  template <class mobj>
+  static void XiToSigmaEye(const PropagatorField &qq_loop,
+	 const mobj &Dd_spec,
+	 const mobj &Ds_spec,
+	 const PropagatorField &qd_tf,
+	 const PropagatorField &qs_ti,
+	 const Gamma Gamma_H,
+	 const Gamma GammaB_sigma,
+	 const Gamma GammaB_nucl,
+         const std::string op,
+	 SpinMatrixField &xts_corr);
 };
 //This computes a baryon contraction on a lattice site, including the spin-trace of the correlation matrix
 template <class FImpl>
@@ -1200,7 +1230,6 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,
 
   Real ee;
 
-
   for (int ie_n=0; ie_n < 6 ; ie_n++){
     int a_n    = (ie_n < 3 ? ie_n       : (6-ie_n)%3 ); //epsilon[ie_n][0]; //a
     int b_n    = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b
@@ -1354,4 +1383,194 @@ void BaryonUtils<FImpl>::SigmaToNucleonNonEye(const PropagatorField &qq_ti,
   }
 }
 
+
+/***********************************************************************
+ * The following code is for Xi -> Sigma rare hypeon decays            *
+ **********************************************************************/
+
+/* Dq_loop is a quark line from t_H to t_H
+ * Dd_spec is a quark line from t_i to t_f
+ * Ds_spec is a quark line from t_i to t_f
+ * Dd_tf is a quark line from t_f to t_H
+ * Ds_ti is a quark line from t_i to t_H */
+template <class FImpl>
+template <class mobj, class mobj2, class robj> accelerator_inline
+void BaryonUtils<FImpl>::XiToSigmaQ1EyeSite(const mobj &Dq_loop,
+						 const mobj2 &Dd_spec,
+						 const mobj2 &Ds_spec,
+						 const mobj &Dd_tf,
+						 const mobj &Ds_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_xi,
+		                 		 const Gamma GammaB_sigma,
+						 robj &result)
+{
+
+  Gamma g5(Gamma::Algebra::Gamma5); 
+
+  auto DdG = Dd_spec * GammaB_sigma;
+  auto GDs = GammaB_xi * Ds_spec;
+  // Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
+  auto DsGDd = Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5;
+  // DsGDd * GammaB
+  auto DsGDdG = DsGDd * GammaB_sigma;
+  // GammaB * DsGDd
+  auto GDsGDd = GammaB_xi * DsGDd;
+  // GammaB * DsGDd * GammaB
+  auto GDsGDdG = GDsGDd * GammaB_sigma;
+  // \gamma_\mu^L * Dq_loop 
+  auto trGDq = TensorRemove(trace(Gamma_H * Dq_loop)); 
+
+  Real ee;
+
+  for (int ie_s=0; ie_s < 6 ; ie_s++){
+      int a_s = (ie_s < 3 ? ie_s       : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a'
+      int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b'
+      int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c'
+      int eSgn_s = (ie_s < 3 ? 1 : -1);
+    for (int ie_x=0; ie_x < 6 ; ie_x++){
+      int a_x = (ie_x < 3 ? ie_x       : (6-ie_x)%3 ); //epsilon[ie_x][0]; //a'
+      int b_x = (ie_x < 3 ? (ie_x+1)%3 : (8-ie_x)%3 ); //epsilon[ie_x][1]; //b'
+      int c_x = (ie_x < 3 ? (ie_x+2)%3 : (7-ie_x)%3 ); //epsilon[ie_x][2]; //c'
+      int eSgn_x = (ie_x < 3 ? 1 : -1);
+      ee = Real(eSgn_s * eSgn_x);
+      auto ee_GD = ee * trGDq;
+      for (int alpha_x=0; alpha_x<Ns; alpha_x++){
+      for (int beta_s=0; beta_s<Ns; beta_s++){
+        auto GDsGDdG_ab_ba = GDsGDd()(alpha_x,beta_s)(b_x,a_s);
+        auto Ds_ab_ab = Ds_spec()(alpha_x,beta_s)(a_x,b_s);
+        auto DsGDdG_ab_aa = DsGDd()(alpha_x,beta_s)(a_x,a_s);
+        auto GDs_ab_bb = GDs()(alpha_x,beta_s)(b_x,b_s);
+        for (int gamma_x=0; gamma_x<Ns; gamma_x++){
+          auto DdG_cb_ca = DdG()(gamma_x,beta_s)(c_x,a_s);
+          for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+            result()(gamma_x,gamma_s)() -= ee_GD * Dd_spec()(gamma_x,gamma_s)(c_x,c_s) * GDsGDdG_ab_ba * Ds_ab_ab;
+            result()(gamma_x,gamma_s)() += ee_GD * DdG_cb_ca * GDsGDd()(alpha_x,gamma_s)(b_x,c_s) * Ds_ab_ab;
+            result()(gamma_x,gamma_s)() += ee_GD * Dd_spec()(gamma_x,gamma_s)(c_x,c_s) * DsGDdG_ab_aa * GDs_ab_bb;
+            result()(gamma_x,gamma_s)() -= ee_GD * DdG_cb_ca * DsGDd()(alpha_x,gamma_s)(a_x,c_s) * GDs_ab_bb;
+          }
+	}
+      }}
+    }
+  }
+}
+
+
+//Equivalent to "One-trace"
+/* Dq_loop is a quark line from t_H to t_H
+ * Dd_spec is a quark line from t_i to t_f
+ * Ds_spec is a quark line from t_i to t_f
+ * Dd_tf is a quark line from t_f to t_H
+ * Ds_ti is a quark line from t_i to t_H */
+template <class FImpl>
+template <class mobj, class mobj2, class robj> accelerator_inline
+void BaryonUtils<FImpl>::XiToSigmaQ2EyeSite(const mobj &Dq_loop,
+						 const mobj2 &Dd_spec,
+						 const mobj2 &Ds_spec,
+						 const mobj &Dd_tf,
+						 const mobj &Ds_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_xi,
+		                 		 const Gamma GammaB_sigma,
+						 robj &result)
+{
+
+  Gamma g5(Gamma::Algebra::Gamma5); 
+
+  auto DdG = Dd_spec * GammaB_sigma;
+  auto GDs = GammaB_xi * Ds_spec;
+  // Ds * \gamma_\mu^L * Dq_loop * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
+  auto DsGDqGDd = Ds_ti * Gamma_H * Dq_loop * Gamma_H * g5 * adj(Dd_tf) * g5;
+  // DsGDd * GammaB
+  auto DsGDqGDdG = DsGDqGDd * GammaB_sigma;
+  // GammaB * DsGDd
+  auto GDsGDqGDd = GammaB_xi * DsGDqGDd;
+  // GammaB * DsGDd * GammaB
+  auto GDsGDqGDdG = GDsGDqGDd * GammaB_sigma;
+
+  Real ee;
+
+  for (int ie_s=0; ie_s < 6 ; ie_s++){
+      int a_s = (ie_s < 3 ? ie_s       : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a'
+      int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b'
+      int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c'
+      int eSgn_s = (ie_s < 3 ? 1 : -1);
+    for (int ie_x=0; ie_x < 6 ; ie_x++){
+      int a_x = (ie_x < 3 ? ie_x       : (6-ie_x)%3 ); //epsilon[ie_x][0]; //a'
+      int b_x = (ie_x < 3 ? (ie_x+1)%3 : (8-ie_x)%3 ); //epsilon[ie_x][1]; //b'
+      int c_x = (ie_x < 3 ? (ie_x+2)%3 : (7-ie_x)%3 ); //epsilon[ie_x][2]; //c'
+      int eSgn_x = (ie_x < 3 ? 1 : -1);
+      ee = Real(eSgn_s * eSgn_x);
+      for (int alpha_x=0; alpha_x<Ns; alpha_x++){
+      for (int beta_s=0; beta_s<Ns; beta_s++){
+        auto GDsGDqGDdG_ab_ba = GDsGDqGDdG()(alpha_x,beta_s)(b_x,a_s);
+        auto Ds_ab_ab = Ds_spec()(alpha_x,beta_s)(a_x,b_s);
+        auto DsGDqGDdG_ab_aa = GDsGDqGDdG()(alpha_x,beta_s)(a_x,a_s);
+        auto GDs_ab_bb = GDs()(alpha_x,beta_s)(b_x,b_s);
+          for (int gamma_x=0; gamma_x<Ns; gamma_x++){
+            auto DdG_cb_ca = DdG()(gamma_x, beta_s)(c_x,a_s);
+          for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+            result()(gamma_x,gamma_s)() -= ee * Dd_spec()(gamma_x,gamma_s)(c_x,c_s) * GDsGDqGDdG_ab_ba * Ds_ab_ab;
+            result()(gamma_x,gamma_s)() += ee * DdG_cb_ca * GDsGDqGDd()(alpha_x,gamma_s)(b_x,c_s) * Ds_ab_ab;
+            result()(gamma_x,gamma_s)() += ee * Dd_spec()(gamma_x,gamma_s)(c_x,c_s) * DsGDqGDdG_ab_aa * GDs_ab_bb;
+            result()(gamma_x,gamma_s)() -= ee * DdG_cb_ca * DsGDqGDd()(alpha_x,gamma_s)(a_x,c_s) * GDs_ab_bb;
+          }}
+      }}
+    }
+  }
+}
+
+template<class FImpl>
+template <class mobj>
+void BaryonUtils<FImpl>::XiToSigmaEye(const PropagatorField &qq_loop,
+						 const mobj &Dd_spec,
+						 const mobj &Ds_spec,
+						 const PropagatorField &qd_tf,
+						 const PropagatorField &qs_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_xi,
+		                 		 const Gamma GammaB_sigma,
+						 const std::string op,
+						 SpinMatrixField &xts_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 = qs_ti.Grid();
+
+  autoView( vcorr   , xts_corr , AcceleratorWrite);
+  autoView( vq_loop , qq_loop  , AcceleratorRead);
+  autoView( vd_tf   , qd_tf    , AcceleratorRead);
+  autoView( vs_ti   , qs_ti    , AcceleratorRead);
+
+  Vector<mobj> my_Dq_spec{Dd_spec,Ds_spec};
+  mobj * Dq_spec_p = &my_Dq_spec[0];
+
+  if(op == "Q1"){
+    accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
+      auto Dq_loop = vq_loop(ss);
+      auto Dd_tf   = vd_tf(ss);
+      auto Ds_ti   = vs_ti(ss);
+      typedef decltype(coalescedRead(vcorr[0])) spinor;
+      spinor result=Zero();
+      XiToSigmaQ1EyeSite(Dq_loop,Dq_spec_p[0],Dq_spec_p[1],Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result);
+      coalescedWrite(vcorr[ss],result);
+    }  );//end loop over lattice sites
+  } else if(op == "Q2"){
+    accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
+      auto Dq_loop = vq_loop(ss);
+      auto Dd_tf   = vd_tf(ss);
+      auto Ds_ti   = vs_ti(ss);
+      typedef decltype(coalescedRead(vcorr[0])) spinor;
+      spinor result=Zero();
+      XiToSigmaQ2EyeSite(Dq_loop,Dq_spec_p[0],Dq_spec_p[1],Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result);
+      coalescedWrite(vcorr[ss],result);
+    }  );//end loop over lattice sites
+  } else {
+    assert(0 && "Weak Operator not correctly specified");
+  }
+}
+
+
 NAMESPACE_END(Grid);

Grid/simd/Grid_vector_types.h

@@ -208,8 +208,8 @@ struct RealPart<complex<T> > {
 //////////////////////////////////////
 // type alias used to simplify the syntax of std::enable_if
 template <typename T> using Invoke = typename T::type;
-template <typename Condition, typename ReturnType> using EnableIf    = Invoke<std::enable_if<Condition::value, ReturnType> >;
+template <typename Condition, typename ReturnType = void> using EnableIf    = Invoke<std::enable_if<Condition::value, ReturnType> >;
-template <typename Condition, typename ReturnType> using NotEnableIf = Invoke<std::enable_if<!Condition::value, ReturnType> >;
+template <typename Condition, typename ReturnType = void> using NotEnableIf = Invoke<std::enable_if<!Condition::value, ReturnType> >;
 
 ////////////////////////////////////////////////////////
 // Check for complexity with type traits

Grid/tensors/Tensor_exp.h

@@ -28,7 +28,7 @@ Author: neo <cossu@post.kek.jp>
 #ifndef GRID_MATH_EXP_H
 #define GRID_MATH_EXP_H
 
-#define DEFAULT_MAT_EXP 12
+#define DEFAULT_MAT_EXP 20
 
 NAMESPACE_BEGIN(Grid);
 

Grid/util/Init.cc

@@ -140,7 +140,7 @@ void GridCmdOptionCSL(std::string str,std::vector<std::string> & vec)
 }
 
 template<class VectorInt>
-void GridCmdOptionIntVector(std::string &str,VectorInt & vec)
+void GridCmdOptionIntVector(const std::string &str,VectorInt & vec)
 {
   vec.resize(0);
   std::stringstream ss(str);

Grid/util/Init.h

@@ -55,7 +55,7 @@ template<class VectorInt>
 std::string GridCmdVectorIntToString(const VectorInt & vec);
 void GridCmdOptionCSL(std::string str,std::vector<std::string> & vec);
 template<class VectorInt>
-void GridCmdOptionIntVector(std::string &str,VectorInt & vec);
+void GridCmdOptionIntVector(const std::string &str,VectorInt & vec);
 void GridCmdOptionInt(std::string &str,int & val);
 
 

README.md

@@ -149,7 +149,6 @@ If you want to build all the tests at once just use `make tests`.
 - `--enable-numa`: enable NUMA first touch optimisation
 - `--enable-simd=<code>`: setup Grid for the SIMD target `<code>` (default: `GEN`). A list of possible SIMD targets is detailed in a section below.
 - `--enable-gen-simd-width=<size>`: select the size (in bytes) of the generic SIMD vector type (default: 32 bytes).
-- `--enable-precision={single|double}`: set the default precision (default: `double`). **Deprecated option**
 - `--enable-comms=<comm>`: Use `<comm>` for message passing (default: `none`). A list of possible SIMD targets is detailed in a section below.
 - `--enable-rng={sitmo|ranlux48|mt19937}`: choose the RNG (default: `sitmo `).
 - `--disable-timers`: disable system dependent high-resolution timers.

tests/debug/Test_cayley_mres.cc

@@ -33,13 +33,14 @@ using namespace Grid;
 
 
 template<class What> 
-void  TestConserved(What & Ddwf, What & Ddwfrev, 
+void  TestConserved(What & Ddwf,
 		    LatticeGaugeField &Umu,
 		    GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
 		    GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
 		    RealD mass, RealD M5,
 		    GridParallelRNG *RNG4,
-		    GridParallelRNG *RNG5);
+		    GridParallelRNG *RNG5,
+                    What *Ddwfrev=nullptr);
 
   Gamma::Algebra Gmu [] = {
     Gamma::Algebra::GammaX,
@@ -102,10 +103,11 @@ int main (int argc, char ** argv)
   GridRedBlackCartesian * FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGridF);
 
 
-  std::vector<int> seeds4({1,2,3,4});
   std::vector<int> seeds5({5,6,7,8});
   GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
-  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+  GridParallelRNG          RNG4(UGrid);
+  std::vector<int> seeds4({1,2,3,4}); RNG4.SeedFixedIntegers(seeds4);
+  //const std::string seeds4{ "test-gauge-3000" }; RNG4.SeedUniqueString( seeds4 );
 
   LatticeGaugeField Umu(UGrid);
   if( argc > 1 && argv[1][0] != '-' )
@@ -116,8 +118,8 @@ int main (int argc, char ** argv)
   }
   else
   {
-    std::cout<<GridLogMessage <<"Using cold configuration"<<std::endl;
+    std::cout<<GridLogMessage <<"Using hot configuration"<<std::endl;
-    //SU<Nc>::ColdConfiguration(Umu);
+    // SU<Nc>::ColdConfiguration(Umu);
     SU<Nc>::HotConfiguration(RNG4,Umu);
   }
 
@@ -127,7 +129,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage <<"DomainWallFermion test"<<std::endl;
   std::cout<<GridLogMessage <<"======================"<<std::endl;
   DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  TestConserved<DomainWallFermionR>(Ddwf,Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestConserved<DomainWallFermionR>(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   RealD b=1.5;// Scale factor b+c=2, b-c=1
   RealD c=0.5;
@@ -137,13 +139,13 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage <<"MobiusFermion test"<<std::endl;
   std::cout<<GridLogMessage <<"======================"<<std::endl;
   MobiusFermionR Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
-  TestConserved<MobiusFermionR>(Dmob,Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestConserved<MobiusFermionR>(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"ScaledShamirFermion test"<<std::endl;
   std::cout<<GridLogMessage <<"======================"<<std::endl;
   ScaledShamirFermionR Dsham(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,2.0);
-  TestConserved<ScaledShamirFermionR>(Dsham,Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestConserved<ScaledShamirFermionR>(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"ZMobiusFermion test"<<std::endl;
@@ -152,8 +154,7 @@ int main (int argc, char ** argv)
   //  for(int s=0;s<Ls;s++) omegasrev[s]=omegas[Ls-1-s];
   ZMobiusFermionR ZDmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,omegas,b,c);
   ZMobiusFermionR ZDmobrev(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,omegasrev,b,c);
-
+  TestConserved<ZMobiusFermionR>(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5,&ZDmobrev);
-  TestConserved<ZMobiusFermionR>(ZDmob,ZDmobrev,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   Grid_finalize();
 }
@@ -162,20 +163,15 @@ int main (int argc, char ** argv)
 
 template<class Action> 
 void  TestConserved(Action & Ddwf,
-		    Action & Ddwfrev, 
 		    LatticeGaugeField &Umu,
 		    GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
 		    GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
 		    RealD mass, RealD M5,
 		    GridParallelRNG *RNG4,
-		    GridParallelRNG *RNG5)
+		    GridParallelRNG *RNG5,
+                    Action * Ddwfrev)
 {
-  int Ls=Ddwf.Ls;
-
   LatticePropagator phys_src(UGrid);
-
-  std::vector<LatticeColourMatrix> U(4,UGrid);
-  
   LatticePropagator seqsrc(FGrid);
   LatticePropagator prop5(FGrid); 
   LatticePropagator prop5rev(FGrid); 
@@ -194,9 +190,9 @@ void  TestConserved(Action & Ddwf,
   phys_src=Zero();
   pokeSite(kronecker,phys_src,coor);
   
-  MdagMLinearOperator<Action,LatticeFermion> HermOp(Ddwf);
-  MdagMLinearOperator<Action,LatticeFermion> HermOprev(Ddwfrev);
   ConjugateGradient<LatticeFermion> CG(1.0e-16,100000);
+  SchurRedBlackDiagTwoSolve<LatticeFermion> schur(CG);
+  ZeroGuesser<LatticeFermion> zpg;
   for(int s=0;s<Nd;s++){
     for(int c=0;c<Nc;c++){
       LatticeFermion src4  (UGrid); 
@@ -206,20 +202,20 @@ void  TestConserved(Action & Ddwf,
       Ddwf.ImportPhysicalFermionSource(src4,src5);
 
       LatticeFermion result5(FGrid); result5=Zero();
-
+      schur(Ddwf,src5,result5,zpg);
-      // CGNE
+      std::cout<<GridLogMessage<<"spin "<<s<<" color "<<c<<" norm2(sourc5d) "<<norm2(src5)
-      LatticeFermion Mdagsrc5  (FGrid); 
+               <<" norm2(result5d) "<<norm2(result5)<<std::endl;
-      Ddwf.Mdag(src5,Mdagsrc5);
-      CG(HermOp,Mdagsrc5,result5);
       FermToProp<Action>(prop5,result5,s,c);
 
       LatticeFermion result4(UGrid);
       Ddwf.ExportPhysicalFermionSolution(result5,result4);
       FermToProp<Action>(prop4,result4,s,c);
 
-      Ddwfrev.ImportPhysicalFermionSource(src4,src5);
+      if( Ddwfrev ) {
-      Ddwfrev.Mdag(src5,Mdagsrc5);
+        Ddwfrev->ImportPhysicalFermionSource(src4,src5);
-      CG(HermOprev,Mdagsrc5,result5);
+        result5 = Zero();
+        schur(*Ddwfrev,src5,result5,zpg);
+      }
       FermToProp<Action>(prop5rev,result5,s,c);
     }
   }
@@ -251,11 +247,7 @@ void  TestConserved(Action & Ddwf,
       PropToFerm<Action>(src5,seqsrc,s,c);
 
       LatticeFermion result5(FGrid); result5=Zero();
-
+      schur(Ddwf,src5,result5,zpg);
-      // CGNE
-      LatticeFermion Mdagsrc5  (FGrid); 
-      Ddwf.Mdag(src5,Mdagsrc5);
-      CG(HermOp,Mdagsrc5,result5);
 
       LatticeFermion result4(UGrid);
       Ddwf.ExportPhysicalFermionSolution(result5,result4);
@@ -276,10 +268,10 @@ void  TestConserved(Action & Ddwf,
   Ddwf.ContractConservedCurrent(prop5rev,prop5,Vector_mu,phys_src,Current::Vector,Tdir);
   Ddwf.ContractJ5q(prop5,PJ5q);
   
-  PA       = trace(g5*Axial_mu);
+  PA       = trace(g5*Axial_mu);      // Pseudoscalar-Axial conserved current
-  SV       = trace(Vector_mu);
+  SV       = trace(Vector_mu);        // Scalar-Vector conserved current
-  VV       = trace(gT*Vector_mu);
+  VV       = trace(gT*Vector_mu);     // (local) Vector-Vector conserved current
-  PP       = trace(adj(prop4)*prop4);
+  PP       = trace(adj(prop4)*prop4); // Pseudoscalar density
   
   // Spatial sum
   sliceSum(PA,sumPA,Tdir);
@@ -288,15 +280,17 @@ void  TestConserved(Action & Ddwf,
   sliceSum(PP,sumPP,Tdir);
   sliceSum(PJ5q,sumPJ5q,Tdir);
 
-  int Nt=sumPA.size();
+  const int Nt{static_cast<int>(sumPA.size())};
+  std::cout<<GridLogMessage<<"Vector Ward identity by timeslice (~ 0)"<<std::endl;
   for(int t=0;t<Nt;t++){
-    std::cout <<" SV "<<real(TensorRemove(sumSV[t]));
+    std::cout<<GridLogMessage <<" t "<<t<<" SV "<<real(TensorRemove(sumSV[t]))<<" VV "<<real(TensorRemove(sumVV[t]))<<std::endl;
-    std::cout <<" VV "<<real(TensorRemove(sumVV[t]))<<std::endl;
   }
+  std::cout<<GridLogMessage<<"Axial Ward identity by timeslice (defect ~ 0)"<<std::endl;
   for(int t=0;t<Nt;t++){
-    std::cout <<" PAc "<<real(TensorRemove(sumPA[t]));
+    const RealD DmuPAmu{real(TensorRemove(sumPA[t]-sumPA[(t-1+Nt)%Nt]))};
-    std::cout <<" PJ5q "<<real(TensorRemove(sumPJ5q[t]));
+    std::cout<<GridLogMessage<<" t "<<t<<" DmuPAmu "<<DmuPAmu
-    std::cout <<" Ward Identity defect " <<real(TensorRemove(sumPA[t]-sumPA[(t-1+Nt)%Nt] - 2.0*(Ddwf.mass*sumPP[t] + sumPJ5q[t]) ))<<"\n";
+             <<" PP "<<real(TensorRemove(sumPP[t]))<<" PJ5q "<<real(TensorRemove(sumPJ5q[t]))
+             <<" Ward Identity defect " <<(DmuPAmu - 2.*real(TensorRemove(Ddwf.mass*sumPP[t] + sumPJ5q[t])))<<std::endl;
   }
   
   ///////////////////////////////