Conjugate residual algorithm; some more unary functions

lib/Make.inc

@@ -1,4 +1,4 @@
 
-HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Comparison.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vRealF.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./Tensors.h ./Threads.h
+HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Comparison.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vRealF.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./Tensors.h ./Threads.h
 
 CCFILES=./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./GridInit.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

lib/Simd.h

@@ -35,6 +35,8 @@ namespace Grid {
   inline RealD conjugate(const RealD  & r){ return r; }
   inline RealD real(const RealD  & r){ return r; }
 
+  inline RealD sqrt(const RealD  & r){ return std::sqrt(r); }
+
   inline ComplexD conjugate(const ComplexD& r){ return(conj(r)); }
   inline ComplexD adj(const ComplexD& r){ return(conjugate(r)); }
   inline ComplexF conjugate(const ComplexF& r ){ return(conj(r)); }
@@ -112,6 +114,7 @@ namespace Grid {
 };
 
 #include <simd/Grid_vector_types.h>
+#include <simd/Grid_vector_unops.h>
 
 namespace Grid {
   // Default precision

lib/Tensors.h

@@ -12,6 +12,7 @@
 #include <tensors/Tensor_peek.h>
 #include <tensors/Tensor_poke.h>
 #include <tensors/Tensor_reality.h>
+#include <tensors/Tensor_unary.h>
 #include <tensors/Tensor_extract_merge.h>
     
 #endif

lib/algorithms/CoarsenedMatrix.h

@@ -0,0 +1,146 @@
+#ifndef  GRID_ALGORITHM_COARSENED_MATRIX_H
+#define  GRID_ALGORITHM_COARSENED_MATRIX_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  class Geometry {
+  public:
+    int npoint;
+    int dimension;
+    std::vector<int> directions   ;
+    std::vector<int> displacements;
+
+    Geometry(int _d) : dimension(_d), npoint(2*_d+1), directions(npoint), displacements(npoint) {
+      for(int d=0;d<dimension;d++){
+	directions[2*d  ] = d;
+	directions[2*d+1] = d;
+	displacements[2*d  ] = +1;
+	displacements[2*d+1] = -1;
+      }
+      directions   [2*dimension]=0;
+      displacements[2*dimension]=0;
+    }
+    
+    std::vector<int> GetDelta(int point) {
+      std::vector<int> delta(dimension,0);
+      delta[directions[point]] = displacements[point];
+      return delta;
+    };
+
+  };
+  
+  // Fine Object == (per site) type of fine field
+  // nbasis      == number of deflation vectors
+  template<class Fobj,class CComplex,int nbasis>
+  class CoarsenedMatrix : public SparseMatrixBase<Lattice<iVector<vComplex,nbasis > > >  {
+  public:
+    
+    typedef iVector<vComplex,nbasis > siteVector;
+    typedef Lattice<iVector<vComplex,nbasis > > CoarseVector;
+    typedef Lattice<iMatrix<vComplex,nbasis > > CoarseMatrix;
+
+    typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
+    typedef Lattice<Fobj >        FineField;
+
+    ////////////////////
+    // Data members
+    ////////////////////
+    Geometry         geom;
+    GridBase *       _grid; 
+    CartesianStencil Stencil; 
+    std::vector<CoarseMatrix> A;
+    std::vector<siteVector,alignedAllocator<siteVector> >   comm_buf;
+      
+    ///////////////////////
+    // Interface
+    ///////////////////////
+    GridBase * Grid(void)         { return _grid; };   // this is all the linalg routines need to know
+
+    RealD M    (const CoarseVector &in, CoarseVector &out){
+      
+      SimpleCompressor<siteVector> compressor;
+      Stencil.HaloExchange(in,comm_buf,compressor);
+
+      //PARALLEL_FOR_LOOP
+      for(int ss=0;ss<Grid()->oSites();ss++){
+        siteVector res = zero;
+	siteVector tmp;
+	siteVector nbr;
+
+	int offset,local,perm;
+	for(int point=0;point<geom.npoint;point++){
+	  offset = Stencil._offsets [point][ss];
+	  local  = Stencil._is_local[point][ss];
+	  perm   = Stencil._permute[point][ss];
+	  
+	  if(local&&perm) { 
+	    permute(nbr,in._odata[offset],perm);
+	  } else if(local) { 
+	    nbr = in._odata[offset];
+	  } else {
+	    nbr = comm_buf[offset];
+	  }
+	  res = res + A[point]._odata[ss]*nbr;
+	}
+	vstream(out._odata[ss],res);
+      }
+      return norm2(out);
+    };
+
+    RealD Mdag (const CoarseVector &in, CoarseVector &out){ 
+      return M(in,out);
+    };
+    // Defer support for further coarsening for now
+    void Mdiag    (const CoarseVector &in,  CoarseVector &out){};
+    void Mdir     (const CoarseVector &in,  CoarseVector &out,int dir, int disp){};
+
+    CoarsenedMatrix(GridCartesian &CoarseGrid) 	: 
+
+      _grid(&CoarseGrid),
+      geom(CoarseGrid._ndimension),
+      Stencil(&CoarseGrid,geom.npoint,Even,geom.directions,geom.displacements),
+      A(geom.npoint,&CoarseGrid)
+    {
+      comm_buf.resize(Stencil._unified_buffer_size);
+    };
+
+    void CoarsenOperator(GridBase *FineGrid,LinearOperatorBase<Lattice<Fobj> > &linop,std::vector<Lattice<Fobj> > & subspace){
+
+      FineField  phi(FineGrid);
+      FineField Mphi(FineGrid);
+      CoarseVector Proj(Grid()); 
+      CoarseScalar InnerProd(Grid()); 
+
+      // Orthogonalise the subblocks over the basis
+      blockOrthogonalise(InnerProd,subspace);
+
+      // Compute the matrix elements of linop between this orthonormal
+      // set of vectors.
+      for(int i=0;i<nbasis;i++){
+	phi=subspace[i];
+	for(int p=0;p<geom.npoint;p++){ 
+	  int dir = geom.directions[p];
+	  int disp= geom.displacements[p];
+
+	  if ( disp==0 )linop.OpDiag(phi,Mphi);
+	  else linop.OpDir(phi,Mphi,dir,disp); 
+
+	  blockProject(Proj,Mphi,subspace);
+
+	  for(int ss=0;ss<Grid()->oSites();ss++){
+	    for(int j=0;j<nbasis;j++){
+	      A[p]._odata[ss](j,i) = Proj._odata[ss](j);
+	    }
+	  }
+
+	}	
+      }
+      std::cout<<"Computed Coarse Operator"<<std::endl;
+    }
+    
+  };
+
+}
+#endif

lib/algorithms/SparseMatrix.h

@@ -18,8 +18,9 @@ namespace Grid {
 	Field tmp (in._grid);
 	ni=M(in,tmp);
 	no=Mdag(tmp,out);
-	std::cout << "MdagM "<< ni<<" "<<no<<std::endl;
       }
+      virtual  void Mdiag    (const Field &in, Field &out)=0;
+      virtual  void Mdir     (const Field &in, Field &out,int dir, int disp)=0;
     };
 
   /////////////////////////////////////////////////////////////////////////////////////////////

lib/algorithms/iterative/SchurRedBlack.h

@@ -108,7 +108,7 @@ namespace Grid {
       RealD ns = norm2(in);
       RealD nr = norm2(resid);
 
-      std::cout << "SchurRedBlackDiagMooee solver true unprec resid "<< sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+      std::cout << "SchurRedBlackDiagMooee solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
     }     
   };
 

lib/lattice/Lattice_base.h

@@ -299,6 +299,7 @@ PARALLEL_FOR_LOOP
 #include <lattice/Lattice_reality.h>
 #include <lattice/Lattice_coordinate.h>
 #include <lattice/Lattice_rng.h>
+#include <lattice/Lattice_unary.h>
 #include <lattice/Lattice_transfer.h>
 
 

lib/lattice/Lattice_local.h

@@ -25,8 +25,7 @@ PARALLEL_FOR_LOOP
     
     // localInnerProduct
     template<class vobj>
-    inline auto localInnerProduct (const Lattice<vobj> &lhs,const Lattice<vobj> &rhs)
+    inline auto localInnerProduct (const Lattice<vobj> &lhs,const Lattice<vobj> &rhs) -> Lattice<typename vobj::tensor_reduced>
-      -> Lattice<typename vobj::tensor_reduced>
     {
       Lattice<typename vobj::tensor_reduced> ret(rhs._grid);
 PARALLEL_FOR_LOOP

lib/lattice/Lattice_reality.h

@@ -27,9 +27,9 @@ PARALLEL_FOR_LOOP
         return ret;
     };
 
-    template<class vobj> inline auto real(const Lattice<vobj> &z) -> Lattice<decltype(real(z._odata[0]))>
+    template<class vobj> inline auto real(const Lattice<vobj> &z) -> Lattice<vobj>
     {
-      Lattice<decltype(real(z._odata[0]))> ret(z._grid);
+      Lattice<vobj> ret(z._grid);
 PARALLEL_FOR_LOOP
         for(int ss=0;ss<z._grid->oSites();ss++){
             ret._odata[ss] = real(z._odata[ss]);
@@ -37,9 +37,9 @@ PARALLEL_FOR_LOOP
       return ret;
     }
 
-    template<class vobj> inline auto imag(const Lattice<vobj> &z) -> Lattice<decltype(imag(z._odata[0]))>
+    template<class vobj> inline auto imag(const Lattice<vobj> &z) -> Lattice<vobj>
     {
-      Lattice<decltype(imag(z._odata[0]))> ret(z._grid);
+      Lattice<vobj> ret(z._grid);
 PARALLEL_FOR_LOOP
         for(int ss=0;ss<z._grid->oSites();ss++){
             ret._odata[ss] = imag(z._odata[ss]);

lib/lattice/Lattice_transfer.h

@@ -56,10 +56,10 @@ PARALLEL_FOR_LOOP
   }
   
 
-template<class vobj,int nbasis>
+template<class vobj,class CComplex,int nbasis>
-inline void projectBlockBasis(Lattice<iVector<vComplex,nbasis > > &coarseData,
+inline void blockProject(Lattice<iVector<CComplex,nbasis > > &coarseData,
-			      const             Lattice<vobj>   &fineData,
+			 const             Lattice<vobj>   &fineData,
-			      const std::vector<Lattice<vobj> > &Basis)
+			 const std::vector<Lattice<vobj> > &Basis)
 {
   GridBase * fine  = fineData._grid;
   GridBase * coarse= coarseData._grid;
@@ -69,13 +69,14 @@ inline void projectBlockBasis(Lattice<iVector<vComplex,nbasis > > &coarseData,
   assert( nbasis == Basis.size() );
   subdivides(coarse,fine); 
   for(int i=0;i<nbasis;i++){
-    conformable(Basis,fineData);
+    conformable(Basis[i],fineData);
   }
 
   std::vector<int>  block_r      (_ndimension);
   
   for(int d=0 ; d<_ndimension;d++){
     block_r[d] = fine->_rdimensions[d] / coarse->_rdimensions[d];
+    assert(block_r[d]*coarse->_rdimensions[d] == fine->_rdimensions[d]);
   }
 
   coarseData=zero;
@@ -92,20 +93,149 @@ inline void projectBlockBasis(Lattice<iVector<vComplex,nbasis > > &coarseData,
 
     for(int i=0;i<nbasis;i++) {
       
-      coarseData._odata[sc][i]=coarseData._odata[sc][i]
+      coarseData._odata[sc](i)=coarseData._odata[sc](i)
-	+ innerProduct(Basis[i]._odata[sf],fineData._odata[sf]);
+	+ TensorRemove( innerProduct(Basis[i]._odata[sf],fineData._odata[sf]));
 
     }
   }
   return;
+}
+
+template<class vobj,class CComplex>
+inline void blockZAXPY(Lattice<vobj> &fineZ,
+		       const Lattice<CComplex> &coarseA,
+		       const Lattice<vobj> &fineX,
+		       const Lattice<vobj> &fineY)
+{
+  GridBase * fine  = fineZ._grid;
+  GridBase * coarse= coarseA._grid;
+
+  fineZ.checkerboard=fineX.checkerboard;
+  subdivides(coarse,fine); // require they map
+  conformable(fineX,fineY);
+  conformable(fineX,fineZ);
+
+  int _ndimension = coarse->_ndimension;
   
+  std::vector<int>  block_r      (_ndimension);
+
+  // FIXME merge with subdivide checking routine as this is redundant
+  for(int d=0 ; d<_ndimension;d++){
+    block_r[d] = fine->_rdimensions[d] / coarse->_rdimensions[d];
+    assert(block_r[d]*coarse->_rdimensions[d]==fine->_rdimensions[d]);
+  }
+
+PARALLEL_FOR_LOOP
+  for(int sf=0;sf<fine->oSites();sf++){
+    
+    int sc;
+    std::vector<int> coor_c(_ndimension);
+    std::vector<int> coor_f(_ndimension);
+
+    GridBase::CoorFromIndex(coor_f,sf,fine->_rdimensions);
+    for(int d=0;d<_ndimension;d++) coor_c[d]=coor_f[d]/block_r[d];
+    GridBase::IndexFromCoor(coor_c,sc,coarse->_rdimensions);
+
+    // z = A x + y
+    fineZ._odata[sf]=coarseA._odata[sc]*fineX._odata[sf]+fineY._odata[sf];
+
+  }
+
+  return;
+}
+template<class vobj,class CComplex>
+  inline void blockInnerProduct(Lattice<CComplex> &CoarseInner,
+				const Lattice<vobj> &fineX,
+				const Lattice<vobj> &fineY)
+{
+  typedef decltype(innerProduct(fineX._odata[0],fineY._odata[0])) dotp;
+
+  GridBase *coarse(CoarseInner._grid);
+  GridBase *fine  (fineX._grid);
+
+  Lattice<dotp> fine_inner(fine);
+  Lattice<dotp> coarse_inner(coarse);
+
+  fine_inner = localInnerProduct(fineX,fineY);
+  blockSum(coarse_inner,fine_inner);
+  for(int ss=0;ss<coarse->oSites();ss++){
+    CoarseInner._odata[ss] = coarse_inner._odata[ss];
+  }
+}
+template<class vobj,class CComplex>
+inline void blockNormalise(Lattice<CComplex> &ip,Lattice<vobj> &fineX)
+{
+  GridBase *coarse = ip._grid;
+  blockInnerProduct(ip,fineX,fineX);
+  ip = rsqrt(ip);
+  blockZAXPY(fineX,ip,fineX,fineX);
+}
+// useful in multigrid project;
+// Generic name : Coarsen?
+template<class vobj>
+inline void blockSum(Lattice<vobj> &coarseData,const Lattice<vobj> &fineData)
+{
+  GridBase * fine  = fineData._grid;
+  GridBase * coarse= coarseData._grid;
+
+  subdivides(coarse,fine); // require they map
+
+  int _ndimension = coarse->_ndimension;
+  
+  std::vector<int>  block_r      (_ndimension);
+  
+  for(int d=0 ; d<_ndimension;d++){
+    block_r[d] = fine->_rdimensions[d] / coarse->_rdimensions[d];
+  }
+
+  coarseData=zero;
+  for(int sf=0;sf<fine->oSites();sf++){
+    
+    int sc;
+    std::vector<int> coor_c(_ndimension);
+    std::vector<int> coor_f(_ndimension);
+
+    GridBase::CoorFromIndex(coor_f,sf,fine->_rdimensions);
+    for(int d=0;d<_ndimension;d++) coor_c[d]=coor_f[d]/block_r[d];
+    GridBase::IndexFromCoor(coor_c,sc,coarse->_rdimensions);
+
+    coarseData._odata[sc]=coarseData._odata[sc]+fineData._odata[sf];
+
+  }
+  return;
 }
 
 
-template<class vobj,int nbasis>
+template<class vobj,class CComplex>
-inline void promoteBlockBasis(const Lattice<iVector<vComplex,nbasis > > &coarseData,
+inline void blockOrthogonalise(Lattice<CComplex> &ip,std::vector<Lattice<vobj> > &Basis)
-			                        Lattice<vobj>   &fineData,
+{
-			      const std::vector<Lattice<vobj> > &Basis)
+  GridBase *coarse = ip._grid;
+  GridBase *fine   = Basis[0]._grid;
+
+  int       nbasis = Basis.size() ;
+  int  _ndimension = coarse->_ndimension;
+
+  // checks
+  subdivides(coarse,fine); 
+  for(int i=0;i<nbasis;i++){
+    conformable(Basis[i]._grid,fine);
+  }
+
+  for(int v=0;v<nbasis;v++) {
+    for(int u=0;u<v;u++) {
+      //Inner product & remove component 
+      blockInnerProduct(ip,Basis[u],Basis[v]);
+      ip = -ip;
+      blockZAXPY<vobj,CComplex> (Basis[v],ip,Basis[u],Basis[v]);
+    }
+    blockNormalise(ip,Basis[v]);
+  }
+}
+
+template<class vobj,class CComplex,int nbasis>
+inline void blockPromote(const Lattice<iVector<CComplex,nbasis > > &coarseData,
+			 Lattice<vobj>   &fineData,
+			 const std::vector<Lattice<vobj> > &Basis)
 {
   GridBase * fine  = fineData._grid;
   GridBase * coarse= coarseData._grid;
@@ -146,40 +276,6 @@ inline void promoteBlockBasis(const Lattice<iVector<vComplex,nbasis > > &coarseD
   
 }
 
-// useful in multigrid project;
-// Generic name : Coarsen?
-template<class vobj>
-inline void sumBlocks(Lattice<vobj> &coarseData,const Lattice<vobj> &fineData)
-{
-  GridBase * fine  = fineData._grid;
-  GridBase * coarse= coarseData._grid;
-
-  subdivides(coarse,fine); // require they map
-
-  int _ndimension = coarse->_ndimension;
-  
-  std::vector<int>  block_r      (_ndimension);
-  
-  for(int d=0 ; d<_ndimension;d++){
-    block_r[d] = fine->_rdimensions[d] / coarse->_rdimensions[d];
-  }
-
-  coarseData=zero;
-  for(int sf=0;sf<fine->oSites();sf++){
-    
-    int sc;
-    std::vector<int> coor_c(_ndimension);
-    std::vector<int> coor_f(_ndimension);
-
-    GridBase::CoorFromIndex(coor_f,sf,fine->_rdimensions);
-    for(int d=0;d<_ndimension;d++) coor_c[d]=coor_f[d]/block_r[d];
-    GridBase::IndexFromCoor(coor_c,sc,coarse->_rdimensions);
-
-    coarseData._odata[sc]=coarseData._odata[sc]+fineData._odata[sf];
-
-  }
-  return;
-}
 
 }
 #endif

lib/lattice/Lattice_unary.h

@@ -0,0 +1,32 @@
+#ifndef GRID_LATTICE_UNARY_H
+#define GRID_LATTICE_UNARY_H
+
+namespace Grid {
+
+  //////////////////////////////////////////////////////////////////////////////////////////////////////
+  //  avoid copy back routines for mult, mac, sub, add
+  //////////////////////////////////////////////////////////////////////////////////////////////////////
+  template<class obj> Lattice<obj> sqrt(const Lattice<obj> &rhs){
+    Lattice<obj> ret(rhs._grid);
+    ret.checkerboard = rhs.checkerboard;
+    conformable(ret,rhs);
+PARALLEL_FOR_LOOP
+    for(int ss=0;ss<rhs._grid->oSites();ss++){
+      ret._odata[ss]=sqrt(rhs._odata[ss]);
+    }
+    return ret;
+  }
+  template<class obj> Lattice<obj> rsqrt(const Lattice<obj> &rhs){
+    Lattice<obj> ret(rhs._grid);
+    ret.checkerboard = rhs.checkerboard;
+    conformable(ret,rhs);
+PARALLEL_FOR_LOOP
+    for(int ss=0;ss<rhs._grid->oSites();ss++){
+      ret._odata[ss]=rsqrt(rhs._odata[ss]);
+    }
+    return ret;
+  }
+
+
+}
+#endif

lib/qcd/action/fermion/CayleyFermion5D.cc

@@ -165,6 +165,27 @@ namespace QCD {
     }
   }
 
+  void  CayleyFermion5D::Mdir (const LatticeFermion &psi, LatticeFermion &chi,int dir,int disp){
+    LatticeFermion tmp(psi._grid);
+    // Assemble the 5d matrix
+    for(int s=0;s<Ls;s++){
+      if ( s==0 ) {
+	//	tmp = bs psi[s] + cs[s] psi[s+1}
+	//      tmp+= -mass*cs[s] psi[s+1}
+	axpby_ssp_pminus(tmp,beo[s],psi,-ceo[s],psi ,s, s+1);
+	axpby_ssp_pplus(tmp,1.0,tmp,mass*ceo[s],psi,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pminus(tmp,beo[s],psi,mass*ceo[s],psi,s,0);
+	axpby_ssp_pplus(tmp,1.0,tmp,-ceo[s],psi,s,s-1);
+      } else {
+	axpby_ssp_pminus(tmp,beo[s],psi,-ceo[s],psi,s,s+1);
+	axpby_ssp_pplus (tmp,1.0,tmp,-ceo[s],psi,s,s-1);
+      }
+    }
+    // Apply 4d dslash fragment
+    DhopDir(tmp,chi,dir,disp);
+  }
+
   void CayleyFermion5D::MooeeDag    (const LatticeFermion &psi, LatticeFermion &chi)
   {
     for (int s=0;s<Ls;s++){

lib/qcd/action/fermion/CayleyFermion5D.h

@@ -21,6 +21,10 @@ namespace Grid {
       virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
       virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
       virtual void   Instantiatable(void)=0;
+
+      // Efficient support for multigrid coarsening
+      virtual void  Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
+
       //    protected:
       RealD mass;
 

lib/qcd/action/fermion/ContinuedFractionFermion5D.cc

@@ -90,6 +90,18 @@ namespace Grid {
       // Can ignore "dag"
       return M(psi,chi);
     }
+    void  ContinuedFractionFermion5D::Mdir (const LatticeFermion &psi, LatticeFermion &chi,int dir,int disp){
+      DhopDir(psi,chi,dir,disp); // Dslash on diagonal. g5 Dslash is hermitian
+      int sign=1;
+      for(int s=0;s<Ls;s++){
+	if ( s==(Ls-1) ){
+	  ag5xpby_ssp(chi,Beta[s]*ZoloHiInv,chi,0.0,chi,s,s);
+	} else {
+	  ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*ZoloHiInv,chi,0.0,chi,s,s);
+	}
+	sign=-sign; 
+      }
+    }
     void   ContinuedFractionFermion5D::Meooe       (const LatticeFermion &psi, LatticeFermion &chi)
     {
       // Apply 4d dslash

lib/qcd/action/fermion/ContinuedFractionFermion5D.h

@@ -24,6 +24,9 @@ namespace Grid {
       //      virtual void   Instantiatable(void)=0;
       virtual void   Instantiatable(void) =0;
 
+      // Efficient support for multigrid coarsening
+      virtual void  Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
+
       // Constructors
       ContinuedFractionFermion5D(LatticeGaugeField &_Umu,
 				 GridCartesian         &FiveDimGrid,

lib/qcd/action/fermion/FermionOperator.h

@@ -40,6 +40,9 @@ namespace Grid {
       virtual void DhopOE(const FermionField &in, FermionField &out,int dag)=0;
       virtual void DhopEO(const FermionField &in, FermionField &out,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  DhopDir(const FermionField &in, FermionField &out,int dir,int disp)=0; // implemented by WilsonFermion and WilsonFermion5D
 
     };
 

lib/qcd/action/fermion/PartialFractionFermion5D.cc

@@ -2,6 +2,22 @@
 namespace Grid {
   namespace QCD {
 
+    void  PartialFractionFermion5D::Mdir (const LatticeFermion &psi, LatticeFermion &chi,int dir,int disp){
+      // this does both dag and undag but is trivial; make a common helper routing
+
+      int sign = 1;
+
+      DhopDir(psi,chi,dir,disp);
+
+      int nblock=(Ls-1)/2;
+      for(int b=0;b<nblock;b++){
+	int s = 2*b;
+	ag5xpby_ssp(chi,-scale,chi,0.0,chi,s,s); 
+	ag5xpby_ssp(chi, scale,chi,0.0,chi,s+1,s+1); 
+      }
+      ag5xpby_ssp(chi,p[nblock]*scale/amax,chi,0.0,chi,Ls-1,Ls-1);
+
+    }
     void   PartialFractionFermion5D::Meooe_internal(const LatticeFermion &psi, LatticeFermion &chi,int dag)
     {
       // this does both dag and undag but is trivial; make a common helper routing

lib/qcd/action/fermion/PartialFractionFermion5D.h

@@ -30,6 +30,9 @@ namespace Grid {
 
       virtual void   Instantiatable(void) =0; // ensure no make-eee
 
+      // Efficient support for multigrid coarsening
+      virtual void  Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
+
       // Constructors
       PartialFractionFermion5D(LatticeGaugeField &_Umu,
 				    GridCartesian         &FiveDimGrid,

lib/qcd/action/fermion/WilsonFermion.cc

@@ -93,6 +93,26 @@ void WilsonFermion::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out)
   out.checkerboard = in.checkerboard;
   MooeeInv(in,out);
 }
+void WilsonFermion::Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp)
+{
+  DhopDir(in,out,dir,disp);
+}
+void WilsonFermion::DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp){
+  WilsonCompressor compressor(DaggerNo);
+  Stencil.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
+  
+  assert( (disp==1)||(disp==-1) );
+
+  int skip = (disp==1) ? 0 : 1;
+
+  int dirdisp = dir+skip*4;
+
+PARALLEL_FOR_LOOP
+  for(int sss=0;sss<in._grid->oSites();sss++){
+    DiracOpt::DhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp);
+  }
+
+};
 
 void WilsonFermion::DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField & U,
 				const LatticeFermion &in, LatticeFermion &out,int dag)

lib/qcd/action/fermion/WilsonFermion.h

@@ -26,7 +26,7 @@ namespace Grid {
       void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out);
       virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out); // remain virtual so we 
       virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out); // can derive Clover
-      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out); // from Wilson bas
+      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out); // from Wilson base
       virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
 
       // non-hermitian hopping term; half cb or both
@@ -34,6 +34,10 @@ namespace Grid {
       void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag);
       void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag);
 
+      // Multigrid assistance
+      void   Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
+      void DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
+
       ///////////////////////////////////////////////////////////////
       // Extra methods added by derived
       ///////////////////////////////////////////////////////////////

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -82,6 +82,28 @@ void WilsonFermion5D::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGau
     pokeIndex<LorentzIndex>(Uds,U,mu+4);
   }
 }
+void WilsonFermion5D::DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp)
+{
+  assert( (disp==1)||(disp==-1) );
+
+  WilsonCompressor compressor(DaggerNo);
+  Stencil.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
+  
+  int skip = (disp==1) ? 0 : 1;
+
+  int dirdisp = dir+skip*4;
+
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<Umu._grid->oSites();ss++){
+    for(int s=0;s<Ls;s++){
+      int sU=ss;
+      int sF = s+Ls*sU; 
+      DiracOpt::DhopDir(Stencil,Umu,comm_buf,sF,sU,in,out,dirdisp);
+    }
+  }
+
+};
+
 void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
 				   LatticeDoubledGaugeField & U,
 				   const LatticeFermion &in, LatticeFermion &out,int dag)

lib/qcd/action/fermion/WilsonFermion5D.h

@@ -57,6 +57,10 @@ namespace Grid {
       void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag);
       void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag);
 
+      // add a DhopComm
+      // -- suboptimal interface will presently trigger multiple comms.
+      void DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
+
       ///////////////////////////////////////////////////////////////
       // New methods added 
       ///////////////////////////////////////////////////////////////

lib/qcd/action/fermion/WilsonKernels.cc

@@ -13,7 +13,6 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     vHalfSpinColourVector Uchi;
     int offset,local,perm, ptype;
 
-    //#define VERBOSE( A)  if ( ss<10 ) { std::cout << "site " <<ss << " " #A " neigh " << offset << " perm "<< perm <<std::endl;}    
 
     // Xp
     int ss = sF;
@@ -33,12 +32,6 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     mult(&Uchi(),&U._odata[sU](Xp),&chi());
     spReconXp(result,Uchi);
 
-    //    std::cout << "XP_RECON"<<std::endl;
-    //    std::cout << result()(0)(0) <<" "<<result()(0)(1) <<" "<<result()(0)(2) <<std::endl;
-    //    std::cout << result()(1)(0) <<" "<<result()(1)(1) <<" "<<result()(1)(2) <<std::endl;
-    //    std::cout << result()(2)(0) <<" "<<result()(2)(1) <<" "<<result()(2)(2) <<std::endl;
-    //    std::cout << result()(3)(0) <<" "<<result()(3)(1) <<" "<<result()(3)(2) <<std::endl;
-
     // Yp
     offset = st._offsets [Yp][ss];
     local  = st._is_local[Yp][ss];
@@ -93,8 +86,7 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     perm   = st._permute[Xm][ss];
     ptype  = st._permute_type[Xm];
 
-    if ( local && perm ) 
+    if ( local && perm ) {
-    {
       spProjXm(tmp,in._odata[offset]);
       permute(chi,tmp,ptype);
     } else if ( local ) {
@@ -104,12 +96,6 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     }
     mult(&Uchi(),&U._odata[sU](Xm),&chi());
     accumReconXm(result,Uchi);
-    //  std::cout << "XM_RECON_ACCUM"<<std::endl;
-    //    std::cout << result()(0)(0) <<" "<<result()(0)(1) <<" "<<result()(0)(2) <<std::endl;
-    //    std::cout << result()(1)(0) <<" "<<result()(1)(1) <<" "<<result()(1)(2) <<std::endl;
-    //    std::cout << result()(2)(0) <<" "<<result()(2)(1) <<" "<<result()(2)(2) <<std::endl;
-    //    std::cout << result()(3)(0) <<" "<<result()(3)(1) <<" "<<result()(3)(2) <<std::endl;
-
 
     // Ym
     offset = st._offsets [Ym][ss];
@@ -308,4 +294,136 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
 
     vstream(out._odata[ss],result*(-0.5));
 }
+
+void DiracOpt::DhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
+			std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
+		       int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp)
+{
+    vHalfSpinColourVector  tmp;    
+    vHalfSpinColourVector  chi;    
+    vSpinColourVector result;
+    vHalfSpinColourVector Uchi;
+    int offset,local,perm, ptype;
+    int ss=sF;
+    
+    offset = st._offsets [dirdisp][ss];
+    local  = st._is_local[dirdisp][ss];
+    perm   = st._permute[dirdisp][ss];
+    ptype  = st._permute_type[dirdisp];
+
+    // Xp
+    if(dirdisp==Xp){
+      if ( local && perm ) {
+	spProjXp(tmp,in._odata[offset]);
+	permute(chi,tmp,ptype);
+      } else if ( local ) {
+	spProjXp(chi,in._odata[offset]);
+      } else { 
+	chi=buf[offset];
+      }
+      mult(&Uchi(),&U._odata[sU](Xp),&chi());
+      spReconXp(result,Uchi);
+    }
+
+    // Yp
+    if ( dirdisp==Yp ){
+      if ( local && perm ) {
+	spProjYp(tmp,in._odata[offset]);
+	permute(chi,tmp,ptype);
+      } else if ( local ) {
+	spProjYp(chi,in._odata[offset]);
+      } else { 
+	chi=buf[offset];
+      }
+      mult(&Uchi(),&U._odata[sU](Yp),&chi());
+      spReconYp(result,Uchi);
+    }
+
+    // Zp
+    if ( dirdisp ==Zp ){
+      if ( local && perm ) {
+	spProjZp(tmp,in._odata[offset]);
+	permute(chi,tmp,ptype);
+      } else if ( local ) {
+	spProjZp(chi,in._odata[offset]);
+      } else { 
+	chi=buf[offset];
+      }
+      mult(&Uchi(),&U._odata[sU](Zp),&chi());
+      spReconZp(result,Uchi);
+    }
+
+    // Tp
+    if ( dirdisp ==Tp ){
+      if ( local && perm ) {
+	spProjTp(tmp,in._odata[offset]);
+	permute(chi,tmp,ptype);
+      } else if ( local ) {
+	spProjTp(chi,in._odata[offset]);
+      } else { 
+	chi=buf[offset];
+      }
+      mult(&Uchi(),&U._odata[sU](Tp),&chi());
+      spReconTp(result,Uchi);
+    }
+
+    // Xm
+    if ( dirdisp==Xm ){
+      if ( local && perm ) {
+	spProjXm(tmp,in._odata[offset]);
+	permute(chi,tmp,ptype);
+      } else if ( local ) {
+	spProjXm(chi,in._odata[offset]);
+      } else { 
+	chi=buf[offset];
+      }
+      mult(&Uchi(),&U._odata[sU](Xm),&chi());
+      spReconXm(result,Uchi);
+    }
+
+    // Ym
+    if ( dirdisp == Ym ){
+      if ( local && perm ) {
+	spProjYm(tmp,in._odata[offset]);
+	permute(chi,tmp,ptype);
+      } else if ( local ) {
+	spProjYm(chi,in._odata[offset]);
+      } else { 
+	chi=buf[offset];
+      }
+      mult(&Uchi(),&U._odata[sU](Ym),&chi());
+      spReconYm(result,Uchi);
+    }
+
+    // Zm
+    if ( dirdisp == Zm ){
+      if ( local && perm ) {
+	spProjZm(tmp,in._odata[offset]);
+	permute(chi,tmp,ptype);
+      } else if ( local ) {
+	spProjZm(chi,in._odata[offset]);
+      } else { 
+	chi=buf[offset];
+      }
+      mult(&Uchi(),&U._odata[sU](Zm),&chi());
+      spReconZm(result,Uchi);
+    }
+
+    // Tm
+    if ( dirdisp==Tm ) {
+      if ( local && perm ) {
+	spProjTm(tmp,in._odata[offset]);
+	permute(chi,tmp,ptype);
+      } else if ( local ) {
+	spProjTm(chi,in._odata[offset]);
+      } else { 
+	chi=buf[offset];
+      }
+      mult(&Uchi(),&U._odata[sU](Tm),&chi());
+      spReconTm(result,Uchi);
+    }
+
+    vstream(out._odata[ss],result*(-0.5));
+}
+
 }}

lib/qcd/action/fermion/WilsonKernels.h

@@ -20,6 +20,9 @@ namespace Grid {
       static void DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
 			      std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
 			      int sF,int sU,const LatticeFermion &in, LatticeFermion &out);
+      static void DhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U,
+			   std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
+			  int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp);
 
     };
 

lib/qcd/action/fermion/g5HermitianLinop.h

@@ -11,11 +11,22 @@ class Gamma5HermitianLinearOperator : public LinearOperatorBase<Field> {
 public:
   Gamma5HermitianLinearOperator(Matrix &Mat): _Mat(Mat){};
   void Op     (const Field &in, Field &out){
-    _Mat.M(in,out);
+    HermOp(in,out);
   }
   void AdjOp     (const Field &in, Field &out){
-    _Mat.M(in,out);
+    HermOp(in,out);
   }
+  void OpDiag (const Field &in, Field &out) {
+    Field tmp(in._grid);
+    _Mat.Mdiag(in,tmp);
+    G5R5(out,tmp);
+  }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) {
+    Field tmp(in._grid);
+    _Mat.Mdir(in,tmp,dir,disp);
+    G5R5(out,tmp);
+  }
+
   void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
 
     HermOp(in,out);

lib/simd/Grid_avx.h

@@ -314,9 +314,9 @@ namespace Optimization {
   template<>
     inline Grid::ComplexF Reduce<Grid::ComplexF, __m256>::operator()(__m256 in){
     __m256 v1,v2;
-    Optimization::permute(v1,in,0); // sse 128; paired complex single
+    Optimization::permute(v1,in,0); // avx 256; quad complex single
     v1 = _mm256_add_ps(v1,in);
-    Optimization::permute(v2,v1,1); // avx 256; quad complex single
+    Optimization::permute(v2,v1,1); 
     v1 = _mm256_add_ps(v1,v2);
     u256f conv; conv.v = v1;
     return Grid::ComplexF(conv.f[0],conv.f[1]);
@@ -367,7 +367,6 @@ namespace Optimization {
     assert(0);
   }
   
-  
 }
 
 //////////////////////////////////////////////////////////////////////////////////////

lib/simd/Grid_vector_types.h

@@ -78,11 +78,18 @@ namespace Grid {
     typedef typename RealPart < Scalar_type >::type Real; 
     typedef Vector_type     vector_type;
     typedef Scalar_type     scalar_type;
+
+
+    typedef union conv_t_union {
+	Vector_type v;
+	Scalar_type s[sizeof(Vector_type)/sizeof(Scalar_type)];
+      conv_t_union(){};
+      } conv_t;
+    
    
     Vector_type v;
     
     static inline int Nsimd(void) { return sizeof(Vector_type)/sizeof(Scalar_type);}
-
     
     Grid_simd& operator=(const Grid_simd&& rhs){v=rhs.v;return *this;};
     Grid_simd& operator=(const Grid_simd& rhs){v=rhs.v;return *this;}; //faster than not declaring it and leaving to the compiler
@@ -192,6 +199,27 @@ namespace Grid {
       return *this;
     }
 
+
+    ///////////////////////////////////////
+    // Not all functions are supported
+    // through SIMD and must breakout to 
+    // scalar type and back again. This
+    // provides support
+    ///////////////////////////////////////
+
+    template<class functor> friend inline Grid_simd SimdApply (const functor &func,const Grid_simd &v) {
+
+      Grid_simd ret;
+      Grid_simd::conv_t conv;
+
+      conv.v = v.v;
+      for(int i=0;i<Nsimd();i++){
+	conv.s[i]=func(conv.s[i]);
+      }
+      ret.v = conv.v;
+      return ret;
+    }
+
     ////////////////////////////////////////////////////////////////////
     // General permute; assumes vector length is same across 
     // all subtypes; may not be a good assumption, but could

lib/simd/Grid_vector_unops.h

@@ -0,0 +1,60 @@
+#ifndef GRID_VECTOR_UNOPS
+#define GRID_VECTOR_UNOPS
+
+namespace Grid { 
+
+  template<class scalar> struct SqrtRealFunctor {
+    scalar operator()(const scalar &a) const {
+      return sqrt(real(a));
+    }
+  };
+
+  template<class scalar> struct RSqrtRealFunctor {
+    scalar operator()(const scalar &a)  const {
+      return scalar(1.0/sqrt(real(a)));
+    }
+  };
+
+  template<class scalar> struct CosRealFunctor {
+    scalar operator()(const scalar &a)  const {
+      return cos(real(a));
+    }
+  };
+
+  template<class scalar> struct SinRealFunctor {
+    scalar operator()(const scalar &a)  const {
+      return sin(real(a));
+    }
+  };
+
+  template<class scalar> struct PowRealFunctor {
+    double y;
+  PowRealFunctor(double _y) : y(_y) {};
+    scalar operator()(const scalar &a)  const {
+      return pow(real(a),y);
+    }
+  };
+
+  template < class S, class V > 
+  inline Grid_simd<S,V> sqrt(const Grid_simd<S,V> &r) {
+    return SimdApply(SqrtRealFunctor<S>(),r);
+  }
+  template < class S, class V > 
+  inline Grid_simd<S,V> rsqrt(const Grid_simd<S,V> &r) {
+    return SimdApply(RSqrtRealFunctor<S>(),r);
+  }
+  template < class S, class V > 
+  inline Grid_simd<S,V> cos(const Grid_simd<S,V> &r) {
+    return SimdApply(CosRealFunctor<S>(),r);
+  }
+  template < class S, class V > 
+  inline Grid_simd<S,V> sin(const Grid_simd<S,V> &r) {
+    return SimdApply(CosRealFunctor<S>(),r);
+  }
+  template < class S, class V > 
+  inline Grid_simd<S,V> pow(const Grid_simd<S,V> &r,double y) {
+    return SimdApply(PowRealFunctor<S>(y),r);
+  }
+
+}
+#endif

lib/tensors/Tensor_class.h

@@ -123,6 +123,13 @@ public:
   typedef iScalar<tensor_reduced_v> tensor_reduced;
   typedef iVector<recurse_scalar_object,N> scalar_object;
 
+  template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > strong_inline auto operator = (T arg) -> iVector<vtype,N>
+    { 
+      zeroit(*this);
+      for(int i=0;i<N;i++)
+	_internal[i] = arg;
+      return *this;
+    }
 
   enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
   iVector(const Zero &z){ *this = zero; };

lib/tensors/Tensor_unary.h

@@ -0,0 +1,37 @@
+#ifndef GRID_TENSOR_UNARY_H
+#define GRID_TENSOR_UNARY_H
+namespace Grid {
+
+#define UNARY_REAL(func)\
+template<class obj> inline auto func(const iScalar<obj> &z) -> iScalar<obj>\
+{\
+    iScalar<obj> ret;\
+    ret._internal = func( (z._internal));\
+    return ret;\
+}\
+template<class obj,int N> inline auto func(const iVector<obj,N> &z) -> iVector<obj,N>\
+{\
+    iVector<obj,N> ret;\
+    for(int c1=0;c1<N;c1++){\
+      ret._internal[c1] = func( (z._internal[c1]));\
+    }\
+    return ret;\
+}\
+template<class obj,int N> inline auto func(const iMatrix<obj,N> &z) -> iMatrix<obj,N>\
+{\
+    iMatrix<obj,N> ret;\
+    for(int c1=0;c1<N;c1++){\
+    for(int c2=0;c2<N;c2++){\
+      ret._internal[c1][c2] = func( (z._internal[c1][c2]));\
+    }}\
+    return ret;\
+}
+
+UNARY_REAL(sqrt);
+UNARY_REAL(rsqrt);
+UNARY_REAL(sin);
+UNARY_REAL(cos);
+
+
+}
+#endif

tests/Make.inc

@@ -1,5 +1,5 @@
 
-bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_even_odd Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_even_odd Test_gamma Test_main Test_nersc_io Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_even_odd
+bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_gamma Test_main Test_nersc_io Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
 
 
 Test_GaugeAction_SOURCES=Test_GaugeAction.cc
@@ -10,10 +10,18 @@ Test_cayley_cg_SOURCES=Test_cayley_cg.cc
 Test_cayley_cg_LDADD=-lGrid
 
 
+Test_cayley_coarsen_support_SOURCES=Test_cayley_coarsen_support.cc
+Test_cayley_coarsen_support_LDADD=-lGrid
+
+
 Test_cayley_even_odd_SOURCES=Test_cayley_even_odd.cc
 Test_cayley_even_odd_LDADD=-lGrid
 
 
+Test_cf_coarsen_support_SOURCES=Test_cf_coarsen_support.cc
+Test_cf_coarsen_support_LDADD=-lGrid
+
+
 Test_cf_cr_unprec_SOURCES=Test_cf_cr_unprec.cc
 Test_cf_cr_unprec_LDADD=-lGrid
 

tests/Test_GaugeAction.cc

@@ -115,7 +115,7 @@ int main (int argc, char ** argv)
   Complex l  = TensorRemove(Tl);
   std::cout << "calculated link trace " <<l*LinkTraceScale<<std::endl;
 
-  sumBlocks(cPlaq,Plaq);
+  blockSum(cPlaq,Plaq);
   TComplex TcP = sum(cPlaq);
   Complex ll= TensorRemove(TcP);
   std::cout << "coarsened plaquettes sum to " <<ll*PlaqScale<<std::endl;

tests/Test_cayley_coarsen_support.cc

@@ -0,0 +1,108 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+
+template<class d>
+struct scal {
+  d internal;
+};
+
+  Gamma::GammaMatrix Gmu [] = {
+    Gamma::GammaX,
+    Gamma::GammaY,
+    Gamma::GammaZ,
+    Gamma::GammaT
+  };
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+
+
+  const int Ls=8;
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  // Construct a coarsened grid
+  std::vector<int> clatt = GridDefaultLatt();
+  for(int d=0;d<clatt.size();d++){
+    clatt[d] = clatt[d]/2;
+  }
+  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());;
+  GridCartesian *Coarse5d =  SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
+
+  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);
+
+  LatticeFermion    src(FGrid); random(RNG5,src);
+  LatticeFermion result(FGrid); result=zero;
+  LatticeFermion    ref(FGrid); ref=zero;
+  LatticeFermion    tmp(FGrid);
+  LatticeFermion    err(FGrid);
+  LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
+
+  std::vector<LatticeColourMatrix> U(4,UGrid);
+
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = peekIndex<LorentzIndex>(Umu,mu);
+  }
+  
+  RealD mass=0.5;
+  RealD M5=1.8;
+
+  DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+  Gamma5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
+
+  HermIndefOp.Op(src,ref);
+  HermIndefOp.OpDiag(src,result);
+  
+  for(int d=0;d<4;d++){
+    HermIndefOp.OpDir(src,tmp,d,+1); result=result+tmp; 
+    std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
+    HermIndefOp.OpDir(src,tmp,d,-1); result=result+tmp;
+    std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
+  }
+  err = result-ref;
+  std::cout<<"Error "<<norm2(err)<<std::endl;
+
+  const int nbasis = 8;
+  std::vector<LatticeFermion> subspace(nbasis,FGrid);
+  
+  for(int b=0;b<nbasis;b++){
+    random(RNG5,subspace[b]);
+  }
+  std::cout << "Computed randoms"<< std::endl;
+
+  CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOp(*Coarse5d);
+
+  LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace);
+
+  typedef   Lattice<iVector<vComplex,nbasis > > coarse_vec;
+  
+  coarse_vec c_src (Coarse5d);  c_src= zero;
+  coarse_vec c_res (Coarse5d);
+  
+  Complex one(1.0);
+  c_src = one;  // 1 in every element for vector 1.
+
+  //  TODO
+  // -- promote from subspace, check we get the vector we wanted
+  // -- apply ldop; check we get the same as inner product of M times big vec
+  // -- pick blocks one by one. Evaluate matrix elements.
+
+  std::cout << "Multiplying by LittleDiracOp "<< std::endl;
+  LittleDiracOp.M(c_src,c_res);
+
+  std::cout << "Done "<< std::endl;
+  Grid_finalize();
+}

tests/Test_cf_coarsen_support.cc

@@ -0,0 +1,87 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+template<class d>
+struct scal {
+  d internal;
+};
+
+  Gamma::GammaMatrix Gmu [] = {
+    Gamma::GammaX,
+    Gamma::GammaY,
+    Gamma::GammaZ,
+    Gamma::GammaT
+  };
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int Ls=9;
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  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);
+
+  LatticeFermion    src(FGrid); random(RNG5,src);
+  LatticeFermion result(FGrid); result=zero;
+  LatticeFermion    ref(FGrid); ref=zero;
+  LatticeFermion    tmp(FGrid);
+  LatticeFermion    err(FGrid);
+  LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
+
+  std::vector<LatticeColourMatrix> U(4,UGrid);
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = peekIndex<LorentzIndex>(Umu,mu);
+  }
+  
+  RealD mass=0.1;
+  RealD M5=1.8;
+
+  {
+    OverlapWilsonContFracTanhFermion Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
+    HermitianLinearOperator<OverlapWilsonContFracTanhFermion,LatticeFermion> HermIndefOp(Dcf);
+
+    HermIndefOp.Op(src,ref);
+    HermIndefOp.OpDiag(src,result);
+    
+    for(int d=0;d<4;d++){
+      HermIndefOp.OpDir(src,tmp,d,+1); result=result+tmp; 
+      std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
+      HermIndefOp.OpDir(src,tmp,d,-1); result=result+tmp;
+      std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
+    }
+    err = result-ref;
+    std::cout<<"Error "<<norm2(err)<<std::endl;
+  }
+
+  {
+    OverlapWilsonPartialFractionTanhFermion Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
+    HermitianLinearOperator<OverlapWilsonPartialFractionTanhFermion,LatticeFermion> HermIndefOp(Dpf);
+    
+    HermIndefOp.Op(src,ref);
+    HermIndefOp.OpDiag(src,result);
+    
+    for(int d=0;d<4;d++){
+      HermIndefOp.OpDir(src,tmp,d,+1); result=result+tmp; 
+      std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
+      HermIndefOp.OpDir(src,tmp,d,-1); result=result+tmp;
+      std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
+    }
+
+    err = result-ref;
+    std::cout<<"Error "<<norm2(err)<<std::endl;
+  }
+
+
+  Grid_finalize();
+}

tests/Test_cf_cr_unprec.cc

@@ -0,0 +1,58 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+template<class d>
+struct scal {
+  d internal;
+};
+
+  Gamma::GammaMatrix Gmu [] = {
+    Gamma::GammaX,
+    Gamma::GammaY,
+    Gamma::GammaZ,
+    Gamma::GammaT
+  };
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int Ls=9;
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  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);
+
+  LatticeFermion    src(FGrid); random(RNG5,src);
+  LatticeFermion result(FGrid); result=zero;
+  LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
+
+  std::vector<LatticeColourMatrix> U(4,UGrid);
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = peekIndex<LorentzIndex>(Umu,mu);
+  }
+  
+  RealD mass=0.1;
+  RealD M5=1.8;
+
+  OverlapWilsonContFracTanhFermion Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
+
+  ConjugateResidual<LatticeFermion> MCR(1.0e-8,10000);
+
+  MdagMLinearOperator<OverlapWilsonContFracTanhFermion,LatticeFermion> HermPosDefOp(Dcf);
+  MCR(HermPosDefOp,src,result);
+
+  HermitianLinearOperator<OverlapWilsonContFracTanhFermion,LatticeFermion> HermIndefOp(Dcf);
+  MCR(HermIndefOp,src,result);
+
+  Grid_finalize();
+}

tests/Test_dwf_cr_unprec.cc

@@ -0,0 +1,63 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+
+template<class d>
+struct scal {
+  d internal;
+};
+
+  Gamma::GammaMatrix Gmu [] = {
+    Gamma::GammaX,
+    Gamma::GammaY,
+    Gamma::GammaZ,
+    Gamma::GammaT
+  };
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+
+
+  const int Ls=8;
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+
+  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);
+
+  LatticeFermion    src(FGrid); random(RNG5,src);
+  LatticeFermion result(FGrid); result=zero;
+  LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
+
+  std::vector<LatticeColourMatrix> U(4,UGrid);
+
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = peekIndex<LorentzIndex>(Umu,mu);
+  }
+
+  ConjugateResidual<LatticeFermion> MCR(1.0e-8,10000);
+  
+  RealD mass=0.5;
+  RealD M5=1.8;
+  DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+
+  MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermOp(Ddwf);
+  MCR(HermOp,src,result);
+
+  Gamma5HermitianLinearOperator<DomainWallFermion,LatticeFermion> g5HermOp(Ddwf);
+  MCR(g5HermOp,src,result);
+
+
+  Grid_finalize();
+}

tests/Test_nersc_io.cc

@@ -83,7 +83,7 @@ int main (int argc, char ** argv)
   Complex l  = TensorRemove(Tl);
   std::cout << "calculated link trace " <<l*LinkTraceScale<<std::endl;
 
-  sumBlocks(cPlaq,Plaq);
+  blockSum(cPlaq,Plaq);
   TComplex TcP = sum(cPlaq);
   Complex ll= TensorRemove(TcP);
   std::cout << "coarsened plaquettes sum to " <<ll*PlaqScale<<std::endl;

tests/Test_wilson_cr_unprec.cc

@@ -0,0 +1,59 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+template<class d>
+struct scal {
+  d internal;
+};
+
+  Gamma::GammaMatrix Gmu [] = {
+    Gamma::GammaX,
+    Gamma::GammaY,
+    Gamma::GammaZ,
+    Gamma::GammaT
+  };
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplexF::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+  GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+
+  std::vector<int> seeds({1,2,3,4});
+  GridParallelRNG          pRNG(&Grid);  pRNG.SeedFixedIntegers(seeds);
+
+  LatticeFermion src(&Grid); random(pRNG,src);
+  RealD nrm = norm2(src);
+  LatticeFermion result(&Grid); result=zero;
+  LatticeGaugeField Umu(&Grid); random(pRNG,Umu);
+
+  std::vector<LatticeColourMatrix> U(4,&Grid);
+
+  double volume=1;
+  for(int mu=0;mu<Nd;mu++){
+    volume=volume*latt_size[mu];
+  }  
+
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = peekIndex<LorentzIndex>(Umu,mu);
+  }
+  
+  RealD mass=0.5;
+  WilsonFermion Dw(Umu,Grid,RBGrid,mass);
+
+  MdagMLinearOperator<WilsonFermion,LatticeFermion> HermOp(Dw);
+
+  ConjugateResidual<LatticeFermion> MCR(1.0e-8,10000);
+
+  MCR(HermOp,src,result);
+
+  Grid_finalize();
+}