diff --git a/Grid/algorithms/CoarsenedMatrix.h b/Grid/algorithms/CoarsenedMatrix.h
index a6b01986..4c26f799 100644
--- a/Grid/algorithms/CoarsenedMatrix.h
+++ b/Grid/algorithms/CoarsenedMatrix.h
@@ -541,17 +541,14 @@ public:
   ///////////////////////
   GridBase * Grid(void)         { return _grid; };   // this is all the linalg routines need to know
 
-  RealD M (const CoarseVector &in, CoarseVector &out){
-
+  void M (const CoarseVector &in, CoarseVector &out)
+  {
     conformable(_grid,in.Grid());
     conformable(in.Grid(),out.Grid());
 
-    //    RealD Nin = norm2(in);
     SimpleCompressor<siteVector> compressor;
 
-    double comms_usec = -usecond();
     Stencil.HaloExchange(in,compressor);
-    comms_usec += usecond();
 
     auto in_v = in.View();
     auto out_v = out.View();
@@ -565,12 +562,7 @@ public:
     typedef decltype(coalescedRead(in_v[0])) calcVector;
     typedef decltype(coalescedRead(in_v[0](0))) calcComplex;
 
-    GridStopWatch ArithmeticTimer;
     int osites=Grid()->oSites();
-    //    double flops = osites*Nsimd*nbasis*nbasis*8.0*geom.npoint;
-    //    double bytes = osites*nbasis*nbasis*geom.npoint*sizeof(CComplex);
-    double usecs =-usecond();
-    // assert(geom.npoint==9);
 
     accelerator_for(sss, Grid()->oSites()*nbasis, Nsimd, {
       int ss = sss/nbasis;
@@ -598,23 +590,9 @@ public:
       }
       coalescedWrite(out_v[ss](b),res,lane);
     });
-    usecs +=usecond();
-
-    double nrm_usec=-usecond();
-    RealD Nout= norm2(out);
-    nrm_usec+=usecond();
-
-    /*
-        std::cout << GridLogMessage << "\tNorm        " << nrm_usec << " us" <<std::endl;
-        std::cout << GridLogMessage << "\tHalo        " << comms_usec << " us" <<std::endl;
-        std::cout << GridLogMessage << "\tMatrix      " << usecs << " us" <<std::endl;
-        std::cout << GridLogMessage << "\t  mflop/s   " << flops/usecs<<std::endl;
-        std::cout << GridLogMessage << "\t  MB/s      " << bytes/usecs<<std::endl;
-    */
-    return Nout;
   };
 
-  RealD Mdag (const CoarseVector &in, CoarseVector &out)
+  void Mdag (const CoarseVector &in, CoarseVector &out)
   {
     if(hermitian) {
       // corresponds to Petrov-Galerkin coarsening
@@ -625,7 +603,6 @@ public:
       G5C(tmp, in); 
       M(tmp, out);
       G5C(out, out);
-      return norm2(out);
     }
   };
   void MdirComms(const CoarseVector &in)
@@ -870,8 +847,6 @@ public:
 	    auto A_self  = A[self_stencil].View();
 
 	    accelerator_for(ss, Grid()->oSites(), Fobj::Nsimd(),{ coalescedWrite(A_p[ss](j,i),oZProj_v(ss)); });
-	    //      if( disp!= 0 ) { accelerator_for(ss, Grid()->oSites(), Fobj::Nsimd(),{ coalescedWrite(A_p[ss](j,i),oZProj_v(ss)); });}
-	    //	    accelerator_for(ss, Grid()->oSites(), Fobj::Nsimd(),{ coalescedWrite(A_self[ss](j,i),A_self(ss)(j,i)+iZProj_v(ss)); });
 
 	  }
 	}
diff --git a/Grid/algorithms/LinearOperator.h b/Grid/algorithms/LinearOperator.h
index c41f8eef..1add212c 100644
--- a/Grid/algorithms/LinearOperator.h
+++ b/Grid/algorithms/LinearOperator.h
@@ -43,7 +43,6 @@ NAMESPACE_BEGIN(Grid);
 /////////////////////////////////////////////////////////////////////////////////////////////
 template<class Field> class LinearOperatorBase {
 public:
-
   // Support for coarsening to a multigrid
   virtual void OpDiag (const Field &in, Field &out) = 0; // Abstract base
   virtual void OpDir  (const Field &in, Field &out,int dir,int disp) = 0; // Abstract base
@@ -94,7 +93,10 @@ public:
     _Mat.Mdag(in,out);
   }
   void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-    _Mat.MdagM(in,out,n1,n2);
+    _Mat.MdagM(in,out);
+    ComplexD dot = innerProduct(in,out);
+    n1=real(dot);
+    n2=norm2(out);
   }
   void HermOp(const Field &in, Field &out){
     _Mat.MdagM(in,out);
@@ -131,17 +133,14 @@ public:
     assert(0);
   }
   void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-    _Mat.MdagM(in,out,n1,n2);
-    out = out + _shift*in;
-
-    ComplexD dot;	
-    dot= innerProduct(in,out);
+    HermOp(in,out);
+    ComplexD dot = innerProduct(in,out);
     n1=real(dot);
     n2=norm2(out);
   }
   void HermOp(const Field &in, Field &out){
-    RealD n1,n2;
-    HermOpAndNorm(in,out,n1,n2);
+    _Mat.MdagM(in,out);
+    out = out + _shift*in;
   }
 };
 
@@ -170,7 +169,7 @@ public:
     _Mat.M(in,out);
   }
   void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-    _Mat.M(in,out);
+    HermOp(in,out);
     ComplexD dot= innerProduct(in,out); n1=real(dot);
     n2=norm2(out);
   }
@@ -208,339 +207,305 @@ public:
   }
 };
 
-    //////////////////////////////////////////////////////////
-    // Even Odd Schur decomp operators; there are several
-    // ways to introduce the even odd checkerboarding
-    //////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////
+// Even Odd Schur decomp operators; there are several
+// ways to introduce the even odd checkerboarding
+//////////////////////////////////////////////////////////
 
-    template<class Field>
-    class SchurOperatorBase :  public LinearOperatorBase<Field> {
-    public:
-      virtual  RealD Mpc      (const Field &in, Field &out) =0;
-      virtual  RealD MpcDag   (const Field &in, Field &out) =0;
-      virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) {
-      Field tmp(in.Grid());
-      tmp.Checkerboard() = in.Checkerboard();
-	ni=Mpc(in,tmp);
-	no=MpcDag(tmp,out);
-      }
-      virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-      out.Checkerboard() = in.Checkerboard();
-	MpcDagMpc(in,out,n1,n2);
-      }
-      virtual void HermOp(const Field &in, Field &out){
-	RealD n1,n2;
-	HermOpAndNorm(in,out,n1,n2);
-      }
-      void Op     (const Field &in, Field &out){
-	Mpc(in,out);
-      }
-      void AdjOp     (const Field &in, Field &out){ 
-	MpcDag(in,out);
-      }
-      // Support for coarsening to a multigrid
-      void OpDiag (const Field &in, Field &out) {
-	assert(0); // must coarsen the unpreconditioned system
-      }
-      void OpDir  (const Field &in, Field &out,int dir,int disp) {
-	assert(0);
-      }
-      void OpDirAll  (const Field &in, std::vector<Field> &out){
-	assert(0);
-      };
-    };
-    template<class Matrix,class Field>
-    class SchurDiagMooeeOperator :  public SchurOperatorBase<Field> {
-    public:
-      Matrix &_Mat;
-      SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){};
-      virtual  RealD Mpc      (const Field &in, Field &out) {
-      Field tmp(in.Grid());
-      tmp.Checkerboard() = !in.Checkerboard();
-
-	_Mat.Meooe(in,tmp);
-	_Mat.MooeeInv(tmp,out);
-	_Mat.Meooe(out,tmp);
-
-	_Mat.Mooee(in,out);
-	return axpy_norm(out,-1.0,tmp,out);
-      }
-      virtual  RealD MpcDag   (const Field &in, Field &out){
-	Field tmp(in.Grid());
-
-	_Mat.MeooeDag(in,tmp);
-        _Mat.MooeeInvDag(tmp,out);
-	_Mat.MeooeDag(out,tmp);
-
-	_Mat.MooeeDag(in,out);
-	return axpy_norm(out,-1.0,tmp,out);
-      }
-    };
-    template<class Matrix,class Field>
-      class SchurDiagOneOperator :  public SchurOperatorBase<Field> {
-    protected:
-      Matrix &_Mat;
-    public:
-      SchurDiagOneOperator (Matrix &Mat): _Mat(Mat){};
-
-      virtual  RealD Mpc      (const Field &in, Field &out) {
-	Field tmp(in.Grid());
-
-	_Mat.Meooe(in,out);
-	_Mat.MooeeInv(out,tmp);
-	_Mat.Meooe(tmp,out);
-	_Mat.MooeeInv(out,tmp);
-
-	return axpy_norm(out,-1.0,tmp,in);
-      }
-      virtual  RealD MpcDag   (const Field &in, Field &out){
-	Field tmp(in.Grid());
-
-	_Mat.MooeeInvDag(in,out);
-	_Mat.MeooeDag(out,tmp);
-	_Mat.MooeeInvDag(tmp,out);
-	_Mat.MeooeDag(out,tmp);
-
-	return axpy_norm(out,-1.0,tmp,in);
-      }
-    };
-    template<class Matrix,class Field>
-      class SchurDiagTwoOperator :  public SchurOperatorBase<Field> {
-    protected:
-      Matrix &_Mat;
-    public:
-      SchurDiagTwoOperator (Matrix &Mat): _Mat(Mat){};
-
-      virtual  RealD Mpc      (const Field &in, Field &out) {
-	Field tmp(in.Grid());
-
-	_Mat.MooeeInv(in,out);
-	_Mat.Meooe(out,tmp);
-	_Mat.MooeeInv(tmp,out);
-	_Mat.Meooe(out,tmp);
-
-	return axpy_norm(out,-1.0,tmp,in);
-      }
-      virtual  RealD MpcDag   (const Field &in, Field &out){
-	Field tmp(in.Grid());
-
-	_Mat.MeooeDag(in,out);
-	_Mat.MooeeInvDag(out,tmp);
-	_Mat.MeooeDag(tmp,out);
-	_Mat.MooeeInvDag(out,tmp);
-
-	return axpy_norm(out,-1.0,tmp,in);
-      }
-    };
-
-    template<class Field>
-    class NonHermitianSchurOperatorBase :  public LinearOperatorBase<Field> 
-    {
-      public:
-        virtual RealD Mpc      (const Field& in, Field& out) = 0;
-        virtual RealD MpcDag   (const Field& in, Field& out) = 0;
-        virtual void  MpcDagMpc(const Field& in, Field& out, RealD& ni, RealD& no) {
-          Field tmp(in.Grid());
-          tmp.Checkerboard() = in.Checkerboard();
-	        ni = Mpc(in,tmp);
-	        no = MpcDag(tmp,out);
-        }
-        virtual void HermOpAndNorm(const Field& in, Field& out, RealD& n1, RealD& n2) {
-          assert(0);
-        }
-        virtual void HermOp(const Field& in, Field& out) {
-          assert(0);
-        }
-        void Op(const Field& in, Field& out) {
-          Mpc(in, out);
-        }
-        void AdjOp(const Field& in, Field& out) { 
-          MpcDag(in, out);
-        }
-        // Support for coarsening to a multigrid
-        void OpDiag(const Field& in, Field& out) {
-          assert(0); // must coarsen the unpreconditioned system
-        }
-        void OpDir(const Field& in, Field& out, int dir, int disp) {
-          assert(0);
-        }
-        void OpDirAll(const Field& in, std::vector<Field>& out){
-          assert(0);
-        };
-    };
-
-    template<class Matrix, class Field>
-    class NonHermitianSchurDiagMooeeOperator :  public NonHermitianSchurOperatorBase<Field> 
-    {
-      public:
-        Matrix& _Mat;
-        NonHermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){};
-        virtual RealD Mpc(const Field& in, Field& out) {
-          Field tmp(in.Grid());
-          tmp.Checkerboard() = !in.Checkerboard();
-
-  	      _Mat.Meooe(in, tmp);
-	        _Mat.MooeeInv(tmp, out);
-	        _Mat.Meooe(out, tmp);
-
-	        _Mat.Mooee(in, out);
-	
-          return axpy_norm(out, -1.0, tmp, out);
-        }
-        virtual RealD MpcDag(const Field& in, Field& out) {
-	        Field tmp(in.Grid());
-
-	        _Mat.MeooeDag(in, tmp);
-          _Mat.MooeeInvDag(tmp, out);
-	        _Mat.MeooeDag(out, tmp);
-
-	        _Mat.MooeeDag(in, out);
-	
-          return axpy_norm(out, -1.0, tmp, out);
-      }
-    };
-    
-    template<class Matrix,class Field>
-    class NonHermitianSchurDiagOneOperator : public NonHermitianSchurOperatorBase<Field> 
-    {
-      protected:
-        Matrix &_Mat;
-    
-      public:
-        NonHermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){};
-        virtual RealD Mpc(const Field& in, Field& out) {
-	        Field tmp(in.Grid());
-
-	        _Mat.Meooe(in, out);
-	        _Mat.MooeeInv(out, tmp);
-	        _Mat.Meooe(tmp, out);
-	        _Mat.MooeeInv(out, tmp);
-
-	        return axpy_norm(out, -1.0, tmp, in);
-        }
-        virtual RealD MpcDag(const Field& in, Field& out) {
-	        Field tmp(in.Grid());
-
-	        _Mat.MooeeInvDag(in, out);
-	        _Mat.MeooeDag(out, tmp);
-	        _Mat.MooeeInvDag(tmp, out);
-	        _Mat.MeooeDag(out, tmp);
-
-	        return axpy_norm(out, -1.0, tmp, in);
-        }
-    };
-
-    template<class Matrix, class Field>
-    class NonHermitianSchurDiagTwoOperator : public NonHermitianSchurOperatorBase<Field> 
-    {
-      protected:
-        Matrix& _Mat;
-    
-      public:
-        NonHermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){};
-
-        virtual RealD Mpc(const Field& in, Field& out) {
-          Field tmp(in.Grid());
-
-	        _Mat.MooeeInv(in, out);
-	        _Mat.Meooe(out, tmp);
-	        _Mat.MooeeInv(tmp, out);
-	        _Mat.Meooe(out, tmp);
-
-	        return axpy_norm(out, -1.0, tmp, in);
-        }
-        virtual RealD MpcDag(const Field& in, Field& out) {
-	        Field tmp(in.Grid());
-
-          _Mat.MeooeDag(in, out);
-          _Mat.MooeeInvDag(out, tmp);
-          _Mat.MeooeDag(tmp, out);
-          _Mat.MooeeInvDag(out, tmp);
-
-          return axpy_norm(out, -1.0, tmp, in);
-        }
-    };
-
-    ///////////////////////////////////////////////////////////////////////////////////////////////////
-    // Left  handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) psi = eta  -->  ( 1 - Moo^-1 Moe Mee^-1 Meo ) psi = Moo^-1 eta
-    // Right handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) Moo^-1 Moo psi = eta  -->  ( 1 - Moe Mee^-1 Meo Moo^-1) phi=eta ; psi = Moo^-1 phi
-    ///////////////////////////////////////////////////////////////////////////////////////////////////
-    template<class Matrix,class Field> using SchurDiagOneRH = SchurDiagTwoOperator<Matrix,Field> ;
-    template<class Matrix,class Field> using SchurDiagOneLH = SchurDiagOneOperator<Matrix,Field> ;
-    ///////////////////////////////////////////////////////////////////////////////////////////////////
-    //  Staggered use
-    ///////////////////////////////////////////////////////////////////////////////////////////////////
-    template<class Matrix,class Field>
-      class SchurStaggeredOperator :  public SchurOperatorBase<Field> {
-    protected:
-      Matrix &_Mat;
-      Field tmp;
-      RealD mass;
-      double tMpc;
-      double tIP;
-      double tMeo;
-      double taxpby_norm;
-      uint64_t ncall;
-public:
-      void Report(void)
-      {
-	std::cout << GridLogMessage << " HermOpAndNorm.Mpc "<< tMpc/ncall<<" usec "<<std::endl;
-	std::cout << GridLogMessage << " HermOpAndNorm.IP  "<< tIP /ncall<<" usec "<<std::endl;
-	std::cout << GridLogMessage << " Mpc.MeoMoe        "<< tMeo/ncall<<" usec "<<std::endl;
-	std::cout << GridLogMessage << " Mpc.axpby_norm    "<< taxpby_norm/ncall<<" usec "<<std::endl;
-      }
-      SchurStaggeredOperator (Matrix &Mat): _Mat(Mat), tmp(_Mat.RedBlackGrid()) 
-      { 
-	assert( _Mat.isTrivialEE() );
-	mass = _Mat.Mass();
-	tMpc=0;
-	tIP =0;
-        tMeo=0;
-        taxpby_norm=0;
-	ncall=0;
-      }
+template<class Field>
+class SchurOperatorBase :  public LinearOperatorBase<Field> {
+ public:
+  virtual  void Mpc      (const Field &in, Field &out) =0;
+  virtual  void MpcDag   (const Field &in, Field &out) =0;
+  virtual  void MpcDagMpc(const Field &in, Field &out) {
+    Field tmp(in.Grid());
+    tmp.Checkerboard() = in.Checkerboard();
+    Mpc(in,tmp);
+    MpcDag(tmp,out);
+  }
   virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-	ncall++;
-	tMpc-=usecond();
-    n2 = Mpc(in,out);
-	tMpc+=usecond();
-	tIP-=usecond();
-    ComplexD dot= innerProduct(in,out);
-	tIP+=usecond();
-    n1 = real(dot);
+    out.Checkerboard() = in.Checkerboard();
+    MpcDagMpc(in,out);
+    ComplexD dot= innerProduct(in,out); 
+    n1=real(dot);
+    n2=norm2(out);
   }
   virtual void HermOp(const Field &in, Field &out){
-	ncall++;
-	tMpc-=usecond();
-	_Mat.Meooe(in,out);
-	_Mat.Meooe(out,tmp);
-	tMpc+=usecond();
-	taxpby_norm-=usecond();
-	axpby(out,-1.0,mass*mass,tmp,in);
-	taxpby_norm+=usecond();
+    out.Checkerboard() = in.Checkerboard();
+    MpcDagMpc(in,out);
   }
-  virtual  RealD Mpc      (const Field &in, Field &out) 
-  {
+  void Op     (const Field &in, Field &out){
+    Mpc(in,out);
+  }
+  void AdjOp     (const Field &in, Field &out){ 
+    MpcDag(in,out);
+  }
+  // Support for coarsening to a multigrid
+  void OpDiag (const Field &in, Field &out) {
+    assert(0); // must coarsen the unpreconditioned system
+  }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) {
+    assert(0);
+  }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){
+    assert(0);
+  };
+};
+template<class Matrix,class Field>
+  class SchurDiagMooeeOperator :  public SchurOperatorBase<Field> {
+ public:
+    Matrix &_Mat;
+    SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){};
+    virtual  void Mpc      (const Field &in, Field &out) {
+      Field tmp(in.Grid());
+      tmp.Checkerboard() = !in.Checkerboard();
+      
+      _Mat.Meooe(in,tmp);
+      _Mat.MooeeInv(tmp,out);
+      _Mat.Meooe(out,tmp);
+      _Mat.Mooee(in,out);
+      axpy(out,-1.0,tmp,out);
+    }
+    virtual void MpcDag   (const Field &in, Field &out){
+      Field tmp(in.Grid());
+	
+      _Mat.MeooeDag(in,tmp);
+      _Mat.MooeeInvDag(tmp,out);
+      _Mat.MeooeDag(out,tmp);
+      _Mat.MooeeDag(in,out);
+      axpy(out,-1.0,tmp,out);
+    }
+};
+template<class Matrix,class Field>
+  class SchurDiagOneOperator :  public SchurOperatorBase<Field> {
+ protected:
+    Matrix &_Mat;
+ public:
+    SchurDiagOneOperator (Matrix &Mat): _Mat(Mat){};
+    
+    virtual void Mpc      (const Field &in, Field &out) {
+      Field tmp(in.Grid());
 
+      _Mat.Meooe(in,out);
+      _Mat.MooeeInv(out,tmp);
+      _Mat.Meooe(tmp,out);
+      _Mat.MooeeInv(out,tmp);
+      axpy(out,-1.0,tmp,in);
+    }
+    virtual void MpcDag   (const Field &in, Field &out){
+      Field tmp(in.Grid());
+      
+      _Mat.MooeeInvDag(in,out);
+      _Mat.MeooeDag(out,tmp);
+      _Mat.MooeeInvDag(tmp,out);
+      _Mat.MeooeDag(out,tmp);
+      axpy(out,-1.0,tmp,in);
+    }
+};
+template<class Matrix,class Field>
+  class SchurDiagTwoOperator :  public SchurOperatorBase<Field> {
+ protected:
+    Matrix &_Mat;
+ public:
+    SchurDiagTwoOperator (Matrix &Mat): _Mat(Mat){};
+    
+    virtual void Mpc      (const Field &in, Field &out) {
+      Field tmp(in.Grid());
+      
+      _Mat.MooeeInv(in,out);
+      _Mat.Meooe(out,tmp);
+      _Mat.MooeeInv(tmp,out);
+      _Mat.Meooe(out,tmp);
+      
+      axpy(out,-1.0,tmp,in);
+    }
+    virtual  void MpcDag   (const Field &in, Field &out){
+      Field tmp(in.Grid());
+
+      _Mat.MeooeDag(in,out);
+      _Mat.MooeeInvDag(out,tmp);
+      _Mat.MeooeDag(tmp,out);
+      _Mat.MooeeInvDag(out,tmp);
+
+      axpy(out,-1.0,tmp,in);
+    }
+};
+
+template<class Field>
+class NonHermitianSchurOperatorBase :  public LinearOperatorBase<Field> 
+{
+ public:
+  virtual void  Mpc      (const Field& in, Field& out) = 0;
+  virtual void  MpcDag   (const Field& in, Field& out) = 0;
+  virtual void  MpcDagMpc(const Field& in, Field& out) {
+    Field tmp(in.Grid());
+    tmp.Checkerboard() = in.Checkerboard();
+    Mpc(in,tmp);
+    MpcDag(tmp,out);
+  }
+  virtual void HermOpAndNorm(const Field& in, Field& out, RealD& n1, RealD& n2) {
+    assert(0);
+  }
+  virtual void HermOp(const Field& in, Field& out) {
+    assert(0);
+  }
+  void Op(const Field& in, Field& out) {
+    Mpc(in, out);
+  }
+  void AdjOp(const Field& in, Field& out) { 
+    MpcDag(in, out);
+  }
+  // Support for coarsening to a multigrid
+  void OpDiag(const Field& in, Field& out) {
+    assert(0); // must coarsen the unpreconditioned system
+  }
+  void OpDir(const Field& in, Field& out, int dir, int disp) {
+    assert(0);
+  }
+  void OpDirAll(const Field& in, std::vector<Field>& out){
+    assert(0);
+  };
+};
+
+template<class Matrix, class Field>
+class NonHermitianSchurDiagMooeeOperator :  public NonHermitianSchurOperatorBase<Field> 
+{
+ public:
+  Matrix& _Mat;
+ NonHermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){};
+  virtual void Mpc(const Field& in, Field& out) {
+    Field tmp(in.Grid());
+    tmp.Checkerboard() = !in.Checkerboard();
+    
+    _Mat.Meooe(in, tmp);
+    _Mat.MooeeInv(tmp, out);
+    _Mat.Meooe(out, tmp);
+    
+    _Mat.Mooee(in, out);
+    
+    axpy(out, -1.0, tmp, out);
+  }
+  virtual void MpcDag(const Field& in, Field& out) {
+    Field tmp(in.Grid());
+    
+    _Mat.MeooeDag(in, tmp);
+    _Mat.MooeeInvDag(tmp, out);
+    _Mat.MeooeDag(out, tmp);
+	  
+    _Mat.MooeeDag(in, out);
+    
+    axpy(out, -1.0, tmp, out);
+  }
+};
+    
+template<class Matrix,class Field>
+class NonHermitianSchurDiagOneOperator : public NonHermitianSchurOperatorBase<Field> 
+{
+ protected:
+  Matrix &_Mat;
+  
+ public:
+  NonHermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){};
+  virtual void Mpc(const Field& in, Field& out) {
+    Field tmp(in.Grid());
+	  
+    _Mat.Meooe(in, out);
+    _Mat.MooeeInv(out, tmp);
+    _Mat.Meooe(tmp, out);
+    _Mat.MooeeInv(out, tmp);
+
+    axpy(out, -1.0, tmp, in);
+  }
+  virtual void MpcDag(const Field& in, Field& out) {
+    Field tmp(in.Grid());
+    
+    _Mat.MooeeInvDag(in, out);
+    _Mat.MeooeDag(out, tmp);
+    _Mat.MooeeInvDag(tmp, out);
+    _Mat.MeooeDag(out, tmp);
+    
+    axpy(out, -1.0, tmp, in);
+  }
+};
+
+template<class Matrix, class Field>
+class NonHermitianSchurDiagTwoOperator : public NonHermitianSchurOperatorBase<Field> 
+{
+ protected:
+  Matrix& _Mat;
+  
+ public:
+ NonHermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){};
+
+  virtual void Mpc(const Field& in, Field& out) {
+    Field tmp(in.Grid());
+    
+    _Mat.MooeeInv(in, out);
+    _Mat.Meooe(out, tmp);
+    _Mat.MooeeInv(tmp, out);
+    _Mat.Meooe(out, tmp);
+
+    axpy(out, -1.0, tmp, in);
+  }
+  virtual void MpcDag(const Field& in, Field& out) {
+    Field tmp(in.Grid());
+    
+    _Mat.MeooeDag(in, out);
+    _Mat.MooeeInvDag(out, tmp);
+    _Mat.MeooeDag(tmp, out);
+    _Mat.MooeeInvDag(out, tmp);
+
+    axpy(out, -1.0, tmp, in);
+  }
+};
+
+///////////////////////////////////////////////////////////////////////////////////////////////////
+// Left  handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) psi = eta  -->  ( 1 - Moo^-1 Moe Mee^-1 Meo ) psi = Moo^-1 eta
+// Right handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) Moo^-1 Moo psi = eta  -->  ( 1 - Moe Mee^-1 Meo Moo^-1) phi=eta ; psi = Moo^-1 phi
+///////////////////////////////////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field> using SchurDiagOneRH = SchurDiagTwoOperator<Matrix,Field> ;
+template<class Matrix,class Field> using SchurDiagOneLH = SchurDiagOneOperator<Matrix,Field> ;
+///////////////////////////////////////////////////////////////////////////////////////////////////
+//  Staggered use
+///////////////////////////////////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field>
+class SchurStaggeredOperator :  public SchurOperatorBase<Field> {
+ protected:
+  Matrix &_Mat;
+  Field tmp;
+  RealD mass;
+ public:
+  SchurStaggeredOperator (Matrix &Mat): _Mat(Mat), tmp(_Mat.RedBlackGrid()) 
+  { 
+    assert( _Mat.isTrivialEE() );
+    mass = _Mat.Mass();
+  }
+  virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+    Mpc(in,out);
+    ComplexD dot= innerProduct(in,out);
+    n1 = real(dot);
+    n2 =0.0;
+  }
+  virtual void HermOp(const Field &in, Field &out){
+    Mpc(in,out);
+    //    _Mat.Meooe(in,out);
+    //    _Mat.Meooe(out,tmp);
+    //    axpby(out,-1.0,mass*mass,tmp,in);
+  }
+  virtual  void Mpc      (const Field &in, Field &out) 
+  {
     Field tmp(in.Grid());
     Field tmp2(in.Grid());
+	
+    //    _Mat.Mooee(in,out);
+    //    _Mat.Mooee(out,tmp);
 
-    //    std::cout << GridLogIterative << " HermOp.Mpc "<<std::endl;
-    _Mat.Mooee(in,out);
-    _Mat.Mooee(out,tmp);
-    //    std::cout << GridLogIterative << " HermOp.MooeeMooee "<<std::endl;
-
-    tMeo-=usecond();
     _Mat.Meooe(in,out);
     _Mat.Meooe(out,tmp);
-    tMeo+=usecond();
-    taxpby_norm-=usecond();
-    RealD nn=axpby_norm(out,-1.0,mass*mass,tmp,in);
-    taxpby_norm+=usecond();
-    return nn;
+    axpby(out,-1.0,mass*mass,tmp,in);
   }
-  virtual  RealD MpcDag   (const Field &in, Field &out){
-    return Mpc(in,out);
+  virtual  void MpcDag   (const Field &in, Field &out){
+    Mpc(in,out);
   }
   virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) {
     assert(0);// Never need with staggered
@@ -548,7 +513,6 @@ public:
 };
 template<class Matrix,class Field> using SchurStagOperator = SchurStaggeredOperator<Matrix,Field>;
 
-
 /////////////////////////////////////////////////////////////
 // Base classes for functions of operators
 /////////////////////////////////////////////////////////////
diff --git a/Grid/algorithms/SparseMatrix.h b/Grid/algorithms/SparseMatrix.h
index b959f53c..8a265b3f 100644
--- a/Grid/algorithms/SparseMatrix.h
+++ b/Grid/algorithms/SparseMatrix.h
@@ -38,16 +38,12 @@ template<class Field> class SparseMatrixBase {
 public:
   virtual GridBase *Grid(void) =0;
   // Full checkerboar operations
-  virtual RealD M    (const Field &in, Field &out)=0;
-  virtual RealD Mdag (const Field &in, Field &out)=0;
-  virtual void  MdagM(const Field &in, Field &out,RealD &ni,RealD &no) {
-    Field tmp (in.Grid());
-    ni=M(in,tmp);
-    no=Mdag(tmp,out);
-  }
+  virtual void  M    (const Field &in, Field &out)=0;
+  virtual void  Mdag (const Field &in, Field &out)=0;
   virtual void  MdagM(const Field &in, Field &out) {
-    RealD ni, no;
-    MdagM(in,out,ni,no);
+    Field tmp (in.Grid());
+    M(in,tmp);
+    Mdag(tmp,out);
   }
   virtual  void Mdiag    (const Field &in, Field &out)=0;
   virtual  void Mdir     (const Field &in, Field &out,int dir, int disp)=0;
diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.h b/Grid/qcd/action/fermion/CayleyFermion5D.h
index f27f4c23..c7d68d73 100644
--- a/Grid/qcd/action/fermion/CayleyFermion5D.h
+++ b/Grid/qcd/action/fermion/CayleyFermion5D.h
@@ -40,8 +40,8 @@ public:
 public:
 
   // override multiply
-  virtual RealD  M    (const FermionField &in, FermionField &out);
-  virtual RealD  Mdag (const FermionField &in, FermionField &out);
+  virtual void   M    (const FermionField &in, FermionField &out);
+  virtual void   Mdag (const FermionField &in, FermionField &out);
 
   // half checkerboard operations
   virtual void   Meooe       (const FermionField &in, FermionField &out);
diff --git a/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h b/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h
index 5aa7bfbd..2300afd3 100644
--- a/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h
+++ b/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h
@@ -41,8 +41,8 @@ public:
 public:
 
   // override multiply
-  virtual RealD  M    (const FermionField &in, FermionField &out);
-  virtual RealD  Mdag (const FermionField &in, FermionField &out);
+  virtual void   M    (const FermionField &in, FermionField &out);
+  virtual void   Mdag (const FermionField &in, FermionField &out);
 
   // half checkerboard operaions
   virtual void   Meooe       (const FermionField &in, FermionField &out);
diff --git a/Grid/qcd/action/fermion/DomainWallEOFAFermion.h b/Grid/qcd/action/fermion/DomainWallEOFAFermion.h
index a2d0e733..bcc97176 100644
--- a/Grid/qcd/action/fermion/DomainWallEOFAFermion.h
+++ b/Grid/qcd/action/fermion/DomainWallEOFAFermion.h
@@ -53,8 +53,8 @@ public:
   virtual void  DtildeInv  (const FermionField& in, FermionField& out);
 
   // override multiply
-  virtual RealD M          (const FermionField& in, FermionField& out);
-  virtual RealD Mdag       (const FermionField& in, FermionField& out);
+  virtual void  M          (const FermionField& in, FermionField& out);
+  virtual void  Mdag       (const FermionField& in, FermionField& out);
 
   // half checkerboard operations
   virtual void  Mooee      (const FermionField& in, FermionField& out);
diff --git a/Grid/qcd/action/fermion/FermionOperator.h b/Grid/qcd/action/fermion/FermionOperator.h
index f0c2a039..570e350d 100644
--- a/Grid/qcd/action/fermion/FermionOperator.h
+++ b/Grid/qcd/action/fermion/FermionOperator.h
@@ -58,8 +58,8 @@ public:
   virtual GridBase *GaugeRedBlackGrid(void)   =0;
 
   // override multiply
-  virtual RealD  M    (const FermionField &in, FermionField &out)=0;
-  virtual RealD  Mdag (const FermionField &in, FermionField &out)=0;
+  virtual void  M    (const FermionField &in, FermionField &out)=0;
+  virtual void  Mdag (const FermionField &in, FermionField &out)=0;
 
   // half checkerboard operaions
   virtual void   Meooe       (const FermionField &in, FermionField &out)=0;
@@ -86,15 +86,14 @@ public:
   virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;
   virtual void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;
 
-
   virtual void  Mdiag  (const FermionField &in, FermionField &out) { Mooee(in,out);};   // Same as Mooee applied to both CB's
   virtual void  Mdir   (const FermionField &in, FermionField &out,int dir,int disp)=0;   // case by case Wilson, Clover, Cayley, ContFrac, PartFrac
   virtual void  MdirAll(const FermionField &in, std::vector<FermionField> &out)=0;   // case by case Wilson, Clover, Cayley, ContFrac, PartFrac
 
 
-      virtual void  MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m,std::vector<double> twist) { assert(0);};
+  virtual void  MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m,std::vector<double> twist) { assert(0);};
 
-      virtual void  FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary,std::vector<double> twist) 
+  virtual void  FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary,std::vector<double> twist) 
       {
 	FFT theFFT((GridCartesian *) in.Grid());
 
diff --git a/Grid/qcd/action/fermion/ImprovedStaggeredFermion.h b/Grid/qcd/action/fermion/ImprovedStaggeredFermion.h
index 0cfae7b6..ecf44ed7 100644
--- a/Grid/qcd/action/fermion/ImprovedStaggeredFermion.h
+++ b/Grid/qcd/action/fermion/ImprovedStaggeredFermion.h
@@ -71,8 +71,8 @@ public:
   // override multiply; cut number routines if pass dagger argument
   // and also make interface more uniformly consistent
   //////////////////////////////////////////////////////////////////
-  RealD M(const FermionField &in, FermionField &out);
-  RealD Mdag(const FermionField &in, FermionField &out);
+  void M(const FermionField &in, FermionField &out);
+  void Mdag(const FermionField &in, FermionField &out);
 
   /////////////////////////////////////////////////////////
   // half checkerboard operations
diff --git a/Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h b/Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h
index 0ce1c701..d1bb0e9c 100644
--- a/Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h
+++ b/Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h
@@ -1,4 +1,3 @@
-
 /*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
@@ -74,8 +73,8 @@ public:
   GridBase *FermionRedBlackGrid(void)    { return _FiveDimRedBlackGrid;}
 
   // full checkerboard operations; leave unimplemented as abstract for now
-  RealD  M    (const FermionField &in, FermionField &out);
-  RealD  Mdag (const FermionField &in, FermionField &out);
+  void  M    (const FermionField &in, FermionField &out);
+  void  Mdag (const FermionField &in, FermionField &out);
 
   // half checkerboard operations
   void   Meooe       (const FermionField &in, FermionField &out);
diff --git a/Grid/qcd/action/fermion/MobiusEOFAFermion.h b/Grid/qcd/action/fermion/MobiusEOFAFermion.h
index 6b214233..6e4f79eb 100644
--- a/Grid/qcd/action/fermion/MobiusEOFAFermion.h
+++ b/Grid/qcd/action/fermion/MobiusEOFAFermion.h
@@ -56,8 +56,8 @@ public:
   virtual void  DtildeInv        (const FermionField& in, FermionField& out);
 
   // override multiply
-  virtual RealD M                (const FermionField& in, FermionField& out);
-  virtual RealD Mdag             (const FermionField& in, FermionField& out);
+  virtual void  M                (const FermionField& in, FermionField& out);
+  virtual void  Mdag             (const FermionField& in, FermionField& out);
 
   // half checkerboard operations
   virtual void  Mooee            (const FermionField& in, FermionField& out);
diff --git a/Grid/qcd/action/fermion/PartialFractionFermion5D.h b/Grid/qcd/action/fermion/PartialFractionFermion5D.h
index 928abd3f..54f8547f 100644
--- a/Grid/qcd/action/fermion/PartialFractionFermion5D.h
+++ b/Grid/qcd/action/fermion/PartialFractionFermion5D.h
@@ -47,8 +47,8 @@ public:
   void   M_internal(const FermionField &in, FermionField &out,int dag);
 
   // override multiply
-  virtual RealD  M    (const FermionField &in, FermionField &out);
-  virtual RealD  Mdag (const FermionField &in, FermionField &out);
+  virtual void   M    (const FermionField &in, FermionField &out);
+  virtual void   Mdag (const FermionField &in, FermionField &out);
 
   // half checkerboard operaions
   virtual void   Meooe       (const FermionField &in, FermionField &out);
diff --git a/Grid/qcd/action/fermion/WilsonCloverFermion.h b/Grid/qcd/action/fermion/WilsonCloverFermion.h
index 3847b0d9..4b25d00e 100644
--- a/Grid/qcd/action/fermion/WilsonCloverFermion.h
+++ b/Grid/qcd/action/fermion/WilsonCloverFermion.h
@@ -109,9 +109,8 @@ public:
     ImportGauge(_Umu);
   }
 
-  virtual RealD M(const FermionField &in, FermionField &out);
-  virtual RealD Mdag(const FermionField &in, FermionField &out);
-
+  virtual void M(const FermionField &in, FermionField &out);
+  virtual void Mdag(const FermionField &in, FermionField &out);
   virtual void Mooee(const FermionField &in, FermionField &out);
   virtual void MooeeDag(const FermionField &in, FermionField &out);
   virtual void MooeeInv(const FermionField &in, FermionField &out);
diff --git a/Grid/qcd/action/fermion/WilsonFermion.h b/Grid/qcd/action/fermion/WilsonFermion.h
index 2e0bc9bf..1c4dd3cf 100644
--- a/Grid/qcd/action/fermion/WilsonFermion.h
+++ b/Grid/qcd/action/fermion/WilsonFermion.h
@@ -78,8 +78,8 @@ public:
   // override multiply; cut number routines if pass dagger argument
   // and also make interface more uniformly consistent
   //////////////////////////////////////////////////////////////////
-  virtual RealD M(const FermionField &in, FermionField &out);
-  virtual RealD Mdag(const FermionField &in, FermionField &out);
+  virtual void  M(const FermionField &in, FermionField &out);
+  virtual void  Mdag(const FermionField &in, FermionField &out);
 
   /////////////////////////////////////////////////////////
   // half checkerboard operations
diff --git a/Grid/qcd/action/fermion/WilsonFermion5D.h b/Grid/qcd/action/fermion/WilsonFermion5D.h
index ea71376c..804b1d10 100644
--- a/Grid/qcd/action/fermion/WilsonFermion5D.h
+++ b/Grid/qcd/action/fermion/WilsonFermion5D.h
@@ -1,4 +1,3 @@
-
 /*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
@@ -99,8 +98,8 @@ public:
   GridBase *FermionRedBlackGrid(void)    { return _FiveDimRedBlackGrid;}
 
   // full checkerboard operations; leave unimplemented as abstract for now
-  virtual RealD  M    (const FermionField &in, FermionField &out){assert(0); return 0.0;};
-  virtual RealD  Mdag (const FermionField &in, FermionField &out){assert(0); return 0.0;};
+  virtual void   M    (const FermionField &in, FermionField &out){assert(0);};
+  virtual void   Mdag (const FermionField &in, FermionField &out){assert(0);};
 
   // half checkerboard operations; leave unimplemented as abstract for now
   virtual void   Meooe       (const FermionField &in, FermionField &out){assert(0);};
diff --git a/Grid/qcd/action/fermion/WilsonTMFermion5D.h b/Grid/qcd/action/fermion/WilsonTMFermion5D.h
index 71acf763..982e722a 100644
--- a/Grid/qcd/action/fermion/WilsonTMFermion5D.h
+++ b/Grid/qcd/action/fermion/WilsonTMFermion5D.h
@@ -120,7 +120,8 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
     }
   }
   
-  virtual RealD M(const FermionField &in, FermionField &out) {
+  virtual void M(const FermionField &in, FermionField &out) 
+  {
     out.Checkerboard() = in.Checkerboard();
     this->Dhop(in, out, DaggerNo);
     FermionField tmp(out.Grid());
@@ -129,11 +130,12 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
       ComplexD b(0.0,this->mu[s]);
       axpbg5y_ssp(tmp,a,in,b,in,s,s);
     }
-    return axpy_norm(out, 1.0, tmp, out);
+    axpy(out, 1.0, tmp, out);
   }
   
   // needed for fast PV
-  void update(const std::vector<RealD>& _mass, const std::vector<RealD>& _mu) {
+  void update(const std::vector<RealD>& _mass, const std::vector<RealD>& _mu) 
+  {
     assert(_mass.size() == _mu.size());
     assert(_mass.size() == this->FermionGrid()->_fdimensions[0]);
     this->mass = _mass;
diff --git a/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h
index e379026c..e9675b36 100644
--- a/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h
+++ b/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h
@@ -323,7 +323,7 @@ void CayleyFermion5D<Impl>::MeooeDag5D    (const FermionField &psi, FermionField
 }
 
 template<class Impl>
-RealD CayleyFermion5D<Impl>::M    (const FermionField &psi, FermionField &chi)
+void CayleyFermion5D<Impl>::M    (const FermionField &psi, FermionField &chi)
 {
   FermionField Din(psi.Grid());
   
@@ -335,11 +335,10 @@ RealD CayleyFermion5D<Impl>::M    (const FermionField &psi, FermionField &chi)
   axpby(chi,1.0,1.0,chi,psi); 
   
   M5D(psi,chi);
-  return(norm2(chi));
 }
 
 template<class Impl>
-RealD CayleyFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi)
+void CayleyFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi)
 {
   // Under adjoint
   //D1+        D1- P-    ->   D1+^dag   P+ D2-^dag
@@ -354,7 +353,6 @@ RealD CayleyFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi)
   M5Ddag(psi,chi);
   // ((b D_W + D_w hop terms +1) on s-diag
   axpby (chi,1.0,1.0,chi,psi); 
-  return norm2(chi);
 }
 
 // half checkerboard operations
diff --git a/Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h
index beeb3e00..6687800e 100644
--- a/Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h
+++ b/Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h
@@ -94,7 +94,7 @@ void ContinuedFractionFermion5D<Impl>::SetCoefficientsZolotarev(RealD zolo_hi,Ap
 
 
 template<class Impl>
-RealD  ContinuedFractionFermion5D<Impl>::M           (const FermionField &psi, FermionField &chi)
+void ContinuedFractionFermion5D<Impl>::M           (const FermionField &psi, FermionField &chi)
 {
   int Ls = this->Ls;
 
@@ -116,15 +116,14 @@ RealD  ContinuedFractionFermion5D<Impl>::M           (const FermionField &psi, F
     }
     sign=-sign; 
   }
-  return norm2(chi);
 }
 template<class Impl>
-RealD  ContinuedFractionFermion5D<Impl>::Mdag        (const FermionField &psi, FermionField &chi)
+void ContinuedFractionFermion5D<Impl>::Mdag        (const FermionField &psi, FermionField &chi)
 {
   // This matrix is already hermitian. (g5 Dw) = Dw dag g5 = (g5 Dw)dag
   // The rest of matrix is symmetric.
   // Can ignore "dag"
-  return M(psi,chi);
+  M(psi,chi);
 }
 template<class Impl>
 void  ContinuedFractionFermion5D<Impl>::Mdir (const FermionField &psi, FermionField &chi,int dir,int disp){
diff --git a/Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h b/Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h
index 3684fd6c..64ee4033 100644
--- a/Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h
+++ b/Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h
@@ -89,7 +89,7 @@ void DomainWallEOFAFermion<Impl>::DtildeInv(const FermionField& psi, FermionFiel
 /*****************************************************************************************************/
 
 template<class Impl>
-RealD DomainWallEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
+void DomainWallEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
 {
   FermionField Din(psi.Grid());
 
@@ -97,11 +97,10 @@ RealD DomainWallEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
   this->DW(Din, chi, DaggerNo);
   axpby(chi, 1.0, 1.0, chi, psi);
   this->M5D(psi, chi);
-  return(norm2(chi));
 }
 
 template<class Impl>
-RealD DomainWallEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
+void DomainWallEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
 {
   FermionField Din(psi.Grid());
 
@@ -109,7 +108,6 @@ RealD DomainWallEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& c
   this->MeooeDag5D(Din, chi);
   this->M5Ddag(psi, chi);
   axpby(chi, 1.0, 1.0, chi, psi);
-  return(norm2(chi));
 }
 
 /********************************************************************
diff --git a/Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermion5DImplementation.h
index 23692d49..44a201c1 100644
--- a/Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermion5DImplementation.h
+++ b/Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermion5DImplementation.h
@@ -548,21 +548,24 @@ void ImprovedStaggeredFermion5D<Impl>::MdirAll(const FermionField &in, std::vect
   assert(0);
 }
 template <class Impl>
-RealD ImprovedStaggeredFermion5D<Impl>::M(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion5D<Impl>::M(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   Dhop(in, out, DaggerNo);
-  return axpy_norm(out, mass, in, out);
+  axpy(out, mass, in, out);
 }
 
 template <class Impl>
-RealD ImprovedStaggeredFermion5D<Impl>::Mdag(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion5D<Impl>::Mdag(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   Dhop(in, out, DaggerYes);
-  return axpy_norm(out, mass, in, out);
+  axpy(out, mass, in, out);
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion5D<Impl>::Meooe(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion5D<Impl>::Meooe(const FermionField &in, FermionField &out) 
+{
   if (in.Checkerboard() == Odd) {
     DhopEO(in, out, DaggerNo);
   } else {
@@ -570,7 +573,8 @@ void ImprovedStaggeredFermion5D<Impl>::Meooe(const FermionField &in, FermionFiel
   }
 }
 template <class Impl>
-void ImprovedStaggeredFermion5D<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion5D<Impl>::MeooeDag(const FermionField &in, FermionField &out) 
+{
   if (in.Checkerboard() == Odd) {
     DhopEO(in, out, DaggerYes);
   } else {
@@ -579,27 +583,30 @@ void ImprovedStaggeredFermion5D<Impl>::MeooeDag(const FermionField &in, FermionF
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion5D<Impl>::Mooee(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion5D<Impl>::Mooee(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   typename FermionField::scalar_type scal(mass);
   out = scal * in;
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion5D<Impl>::MooeeDag(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion5D<Impl>::MooeeDag(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   Mooee(in, out);
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion5D<Impl>::MooeeInv(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion5D<Impl>::MooeeInv(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   out = (1.0 / (mass)) * in;
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion5D<Impl>::MooeeInvDag(const FermionField &in,
-						   FermionField &out) {
+void ImprovedStaggeredFermion5D<Impl>::MooeeInvDag(const FermionField &in,FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   MooeeInv(in, out);
 }
diff --git a/Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermionImplementation.h b/Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermionImplementation.h
index 37675da0..57f4cb89 100644
--- a/Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermionImplementation.h
+++ b/Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermionImplementation.h
@@ -171,21 +171,24 @@ void ImprovedStaggeredFermion<Impl>::ImportGauge(const GaugeField &_Uthin,const
 /////////////////////////////
 
 template <class Impl>
-RealD ImprovedStaggeredFermion<Impl>::M(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion<Impl>::M(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   Dhop(in, out, DaggerNo);
-  return axpy_norm(out, mass, in, out);
+  axpy(out, mass, in, out);
 }
 
 template <class Impl>
-RealD ImprovedStaggeredFermion<Impl>::Mdag(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion<Impl>::Mdag(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   Dhop(in, out, DaggerYes);
-  return axpy_norm(out, mass, in, out);
+  axpy(out, mass, in, out);
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::Meooe(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion<Impl>::Meooe(const FermionField &in, FermionField &out) 
+{
   if (in.Checkerboard() == Odd) {
     DhopEO(in, out, DaggerNo);
   } else {
@@ -193,7 +196,8 @@ void ImprovedStaggeredFermion<Impl>::Meooe(const FermionField &in, FermionField
   }
 }
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) 
+{
   if (in.Checkerboard() == Odd) {
     DhopEO(in, out, DaggerYes);
   } else {
@@ -202,27 +206,30 @@ void ImprovedStaggeredFermion<Impl>::MeooeDag(const FermionField &in, FermionFie
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::Mooee(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion<Impl>::Mooee(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   typename FermionField::scalar_type scal(mass);
   out = scal * in;
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   Mooee(in, out);
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) {
+void ImprovedStaggeredFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   out = (1.0 / (mass)) * in;
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::MooeeInvDag(const FermionField &in,
-						 FermionField &out) {
+void ImprovedStaggeredFermion<Impl>::MooeeInvDag(const FermionField &in,FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   MooeeInv(in, out);
 }
@@ -234,7 +241,8 @@ void ImprovedStaggeredFermion<Impl>::MooeeInvDag(const FermionField &in,
 template <class Impl>
 void ImprovedStaggeredFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U, DoubledGaugeField &UUU, 
 						   GaugeField & mat,
-						   const FermionField &A, const FermionField &B, int dag) {
+						   const FermionField &A, const FermionField &B, int dag) 
+{
   assert((dag == DaggerNo) || (dag == DaggerYes));
 
   Compressor compressor;
@@ -284,8 +292,8 @@ void ImprovedStaggeredFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGauge
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
-
+void ImprovedStaggeredFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) 
+{
   conformable(U.Grid(), _grid);
   conformable(U.Grid(), V.Grid());
   conformable(U.Grid(), mat.Grid());
@@ -296,8 +304,8 @@ void ImprovedStaggeredFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionFie
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
-
+void ImprovedStaggeredFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) 
+{
   conformable(U.Grid(), _cbgrid);
   conformable(U.Grid(), V.Grid());
   conformable(U.Grid(), mat.Grid());
@@ -310,8 +318,8 @@ void ImprovedStaggeredFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionF
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
-
+void ImprovedStaggeredFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) 
+{
   conformable(U.Grid(), _cbgrid);
   conformable(U.Grid(), V.Grid());
   conformable(U.Grid(), mat.Grid());
diff --git a/Grid/qcd/action/fermion/implementation/MobiusEOFAFermionImplementation.h b/Grid/qcd/action/fermion/implementation/MobiusEOFAFermionImplementation.h
index 256423e6..9b9db178 100644
--- a/Grid/qcd/action/fermion/implementation/MobiusEOFAFermionImplementation.h
+++ b/Grid/qcd/action/fermion/implementation/MobiusEOFAFermionImplementation.h
@@ -166,7 +166,7 @@ void MobiusEOFAFermion<Impl>::DtildeInv(const FermionField& psi, FermionField& c
 /*****************************************************************************************************/
 
 template<class Impl>
-RealD MobiusEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
+void MobiusEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
 {
   FermionField Din(psi.Grid());
 
@@ -174,11 +174,10 @@ RealD MobiusEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
   this->DW(Din, chi, DaggerNo);
   axpby(chi, 1.0, 1.0, chi, psi);
   this->M5D(psi, chi);
-  return(norm2(chi));
 }
 
 template<class Impl>
-RealD MobiusEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
+void MobiusEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
 {
   FermionField Din(psi.Grid());
 
@@ -186,7 +185,6 @@ RealD MobiusEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
   this->MeooeDag5D(Din, chi);
   this->M5Ddag(psi, chi);
   axpby(chi, 1.0, 1.0, chi, psi);
-  return(norm2(chi));
 }
 
 /********************************************************************
diff --git a/Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h
index edc674cc..0206828b 100644
--- a/Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h
+++ b/Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h
@@ -269,16 +269,14 @@ void   PartialFractionFermion5D<Impl>::M_internal(const FermionField &psi, Fermi
 }
 
 template<class Impl>
-RealD  PartialFractionFermion5D<Impl>::M    (const FermionField &in, FermionField &out)
+void PartialFractionFermion5D<Impl>::M    (const FermionField &in, FermionField &out)
 {
   M_internal(in,out,DaggerNo);
-  return norm2(out);
 }
 template<class Impl>
-RealD  PartialFractionFermion5D<Impl>::Mdag (const FermionField &in, FermionField &out)
+void PartialFractionFermion5D<Impl>::Mdag (const FermionField &in, FermionField &out)
 {
   M_internal(in,out,DaggerYes);
-  return norm2(out);
 }
 
 template<class Impl>
diff --git a/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h
index 9d99d9e7..36447153 100644
--- a/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h
+++ b/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h
@@ -35,7 +35,7 @@ NAMESPACE_BEGIN(Grid);
 
 // *NOT* EO
 template <class Impl>
-RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
 {
   FermionField temp(out.Grid());
 
@@ -47,11 +47,10 @@ RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
   Mooee(in, temp);
 
   out += temp;
-  return norm2(out);
 }
 
 template <class Impl>
-RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
 {
   FermionField temp(out.Grid());
 
@@ -63,7 +62,6 @@ RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
   MooeeDag(in, temp);
 
   out += temp;
-  return norm2(out);
 }
 
 template <class Impl>
diff --git a/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h
index be05fcf8..5267e0c1 100644
--- a/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h
+++ b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h
@@ -102,21 +102,24 @@ void WilsonFermion<Impl>::ImportGauge(const GaugeField &_Umu)
 /////////////////////////////
 
 template <class Impl>
-RealD WilsonFermion<Impl>::M(const FermionField &in, FermionField &out) {
+void WilsonFermion<Impl>::M(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   Dhop(in, out, DaggerNo);
-  return axpy_norm(out, diag_mass, in, out);
+  axpy(out, diag_mass, in, out);
 }
 
 template <class Impl>
-RealD WilsonFermion<Impl>::Mdag(const FermionField &in, FermionField &out) {
+void WilsonFermion<Impl>::Mdag(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   Dhop(in, out, DaggerYes);
-  return axpy_norm(out, diag_mass, in, out);
+  axpy(out, diag_mass, in, out);
 }
 
 template <class Impl>
-void WilsonFermion<Impl>::Meooe(const FermionField &in, FermionField &out) {
+void WilsonFermion<Impl>::Meooe(const FermionField &in, FermionField &out) 
+{
   if (in.Checkerboard() == Odd) {
     DhopEO(in, out, DaggerNo);
   } else {
@@ -125,7 +128,8 @@ void WilsonFermion<Impl>::Meooe(const FermionField &in, FermionField &out) {
 }
 
 template <class Impl>
-void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
+void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) 
+{
   if (in.Checkerboard() == Odd) {
     DhopEO(in, out, DaggerYes);
   } else {
@@ -134,26 +138,30 @@ void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
 }
   
 template <class Impl>
-void WilsonFermion<Impl>::Mooee(const FermionField &in, FermionField &out) {
+void WilsonFermion<Impl>::Mooee(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   typename FermionField::scalar_type scal(diag_mass);
   out = scal * in;
 }
 
 template <class Impl>
-void WilsonFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) {
+void WilsonFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   Mooee(in, out);
 }
 
 template<class Impl>
-void WilsonFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) {
+void WilsonFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   out = (1.0/(diag_mass))*in;
 }
   
 template<class Impl>
-void WilsonFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out) {
+void WilsonFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out) 
+{
   out.Checkerboard() = in.Checkerboard();
   MooeeInv(in,out);
 }
@@ -249,7 +257,8 @@ void WilsonFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U,
 }
 
 template <class Impl>
-void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
+void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) 
+{
   conformable(U.Grid(), _grid);
   conformable(U.Grid(), V.Grid());
   conformable(U.Grid(), mat.Grid());
@@ -260,7 +269,8 @@ void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, cons
 }
 
 template <class Impl>
-void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
+void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) 
+{
   conformable(U.Grid(), _cbgrid);
   conformable(U.Grid(), V.Grid());
   //conformable(U.Grid(), mat.Grid()); not general, leaving as a comment (Guido)
@@ -274,7 +284,8 @@ void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, co
 }
 
 template <class Impl>
-void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
+void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) 
+{
   conformable(U.Grid(), _cbgrid);
   conformable(U.Grid(), V.Grid());
   //conformable(U.Grid(), mat.Grid());
@@ -287,7 +298,8 @@ void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, co
 }
 
 template <class Impl>
-void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out, int dag) {
+void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out, int dag) 
+{
   conformable(in.Grid(), _grid);  // verifies full grid
   conformable(in.Grid(), out.Grid());
 
@@ -297,7 +309,8 @@ void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out, int da
 }
 
 template <class Impl>
-void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out, int dag) {
+void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out, int dag) 
+{
   conformable(in.Grid(), _cbgrid);    // verifies half grid
   conformable(in.Grid(), out.Grid());  // drops the cb check
 
@@ -308,7 +321,8 @@ void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out, int
 }
 
 template <class Impl>
-void WilsonFermion<Impl>::DhopEO(const FermionField &in, FermionField &out,int dag) {
+void WilsonFermion<Impl>::DhopEO(const FermionField &in, FermionField &out,int dag) 
+{
   conformable(in.Grid(), _cbgrid);    // verifies half grid
   conformable(in.Grid(), out.Grid());  // drops the cb check
 
@@ -386,7 +400,8 @@ template <class Impl>
 void WilsonFermion<Impl>::DhopInternalOverlappedComms(StencilImpl &st, LebesgueOrder &lo,
 						      DoubledGaugeField &U,
 						      const FermionField &in,
-						      FermionField &out, int dag) {
+						      FermionField &out, int dag) 
+{
   assert((dag == DaggerNo) || (dag == DaggerYes));
 
   Compressor compressor(dag);
@@ -436,7 +451,8 @@ template <class Impl>
 void WilsonFermion<Impl>::DhopInternalSerial(StencilImpl &st, LebesgueOrder &lo,
                                        DoubledGaugeField &U,
                                        const FermionField &in,
-                                       FermionField &out, int dag) {
+                                       FermionField &out, int dag) 
+{
   assert((dag == DaggerNo) || (dag == DaggerYes));
   Compressor compressor(dag);
   st.HaloExchange(in, compressor);
diff --git a/Grid/qcd/modules/Registration.h b/Grid/qcd/modules/Registration.h
index ec28f020..459e1d0b 100644
--- a/Grid/qcd/modules/Registration.h
+++ b/Grid/qcd/modules/Registration.h
@@ -80,10 +80,10 @@ static Registrar<OneFlavourRatioEOFModule<FermionImplementationPolicy>,
 
 static Registrar< ConjugateGradientModule<WilsonFermionR::FermionField>,   
                   HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CGWFmodXMLInit("ConjugateGradient"); 
-static Registrar< BiCGSTABModule<WilsonFermionR::FermionField>,   
-                  HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CGWFmodXMLInit("BiCGSTAB"); 
-static Registrar< ConjugateResidualModule<WilsonFermionR::FermionField>,   
-                  HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CRWFmodXMLInit("ConjugateResidual"); 
+//static Registrar< BiCGSTABModule<WilsonFermionR::FermionField>,   
+//                  HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CGWFmodXMLInit("BiCGSTAB"); 
+//static Registrar< ConjugateResidualModule<WilsonFermionR::FermionField>,   
+//                  HMC_SolverModuleFactory<solver_string, WilsonFermionR::FermionField, Serialiser> > __CRWFmodXMLInit("ConjugateResidual"); 
 
 // add the staggered, scalar versions here
 
diff --git a/Grid/qcd/smearing/GaugeConfiguration.h b/Grid/qcd/smearing/GaugeConfiguration.h
index f4d00c72..6f2ff2e7 100644
--- a/Grid/qcd/smearing/GaugeConfiguration.h
+++ b/Grid/qcd/smearing/GaugeConfiguration.h
@@ -49,7 +49,7 @@ public:
 
 private:
   const unsigned int smearingLevels;
-  Smear_Stout<Gimpl> StoutSmearing;
+  Smear_Stout<Gimpl> &StoutSmearing;
   std::vector<GaugeField> SmearedSet;
 
   // Member functions
diff --git a/Grid/qcd/utils/CovariantCshift.h b/Grid/qcd/utils/CovariantCshift.h
index ed96f3bf..6ac69150 100644
--- a/Grid/qcd/utils/CovariantCshift.h
+++ b/Grid/qcd/utils/CovariantCshift.h
@@ -52,6 +52,26 @@ namespace PeriodicBC {
     tmp = adj(Link)*field;
     return Cshift(tmp,mu,-1);// moves towards positive mu
   }
+
+  template<class gauge,typename Op, typename T1> auto
+    CovShiftForward(const Lattice<gauge> &Link, 
+		    int mu,
+		    const LatticeUnaryExpression<Op,T1> &expr)
+    -> Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> 
+  {
+    Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> arg(expr);
+    return CovShiftForward(Link,mu,arg);
+  }
+  template<class gauge,typename Op, typename T1> auto
+    CovShiftBackward(const Lattice<gauge> &Link, 
+		     int mu,
+		     const LatticeUnaryExpression<Op,T1> &expr)
+    -> Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> 
+  {
+    Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> arg(expr);
+    return CovShiftForward(Link,mu,arg);
+  }
+
 }
 
 
@@ -122,6 +142,26 @@ namespace ConjugateBC {
     return Cshift(tmp,mu,-1);// moves towards positive mu
   }
 
+  template<class gauge,typename Op, typename T1> auto
+    CovShiftForward(const Lattice<gauge> &Link, 
+		    int mu,
+		    const LatticeUnaryExpression<Op,T1> &expr)
+    -> Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> 
+  {
+    Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> arg(expr);
+    return CovShiftForward(Link,mu,arg);
+  }
+  template<class gauge,typename Op, typename T1> auto
+    CovShiftBackward(const Lattice<gauge> &Link, 
+		     int mu,
+		     const LatticeUnaryExpression<Op,T1> &expr)
+    -> Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> 
+  {
+    Lattice<decltype(expr.op.func(eval(0, expr.arg1)))> arg(expr);
+    return CovShiftForward(Link,mu,arg);
+  }
+
+
 }