From d0e357ef892fa7296743601f5e5071c57ee1e50e Mon Sep 17 00:00:00 2001
From: paboyle <paboyle@ph.ed.ac.uk>
Date: Mon, 15 Jan 2018 00:23:51 +0000
Subject: [PATCH] CLeanup and no QCD namespace

---
 lib/algorithms/CoarsenedMatrix.h | 770 +++++++++++++++---------------
 lib/algorithms/FFT.h             | 375 +++++++--------
 lib/algorithms/LinearOperator.h  | 776 +++++++++++++++----------------
 lib/algorithms/Preconditioner.h  |  30 +-
 lib/algorithms/SparseMatrix.h    |  72 +--
 5 files changed, 1012 insertions(+), 1011 deletions(-)

diff --git a/lib/algorithms/CoarsenedMatrix.h b/lib/algorithms/CoarsenedMatrix.h
index 8af8d7ac..8b476603 100644
--- a/lib/algorithms/CoarsenedMatrix.h
+++ b/lib/algorithms/CoarsenedMatrix.h
@@ -1,4 +1,4 @@
-    /*************************************************************************************
+/*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -26,455 +26,455 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
     See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
+*************************************************************************************/
+/*  END LEGAL */
 #ifndef  GRID_ALGORITHM_COARSENED_MATRIX_H
 #define  GRID_ALGORITHM_COARSENED_MATRIX_H
 
 
-namespace Grid {
+NAMESPACE_BEGIN(Grid);
 
-  class Geometry {
-    //    int dimension;
-  public:
-    int npoint;
-    std::vector<int> directions   ;
-    std::vector<int> displacements;
+class Geometry {
+  //    int dimension;
+public:
+  int npoint;
+  std::vector<int> directions   ;
+  std::vector<int> displacements;
 
   Geometry(int _d)  {
   
-      int base = (_d==5) ? 1:0;
+    int base = (_d==5) ? 1:0;
 
-      // make coarse grid stencil for 4d , not 5d
-      if ( _d==5 ) _d=4;
+    // make coarse grid stencil for 4d , not 5d
+    if ( _d==5 ) _d=4;
 
-      npoint = 2*_d+1;
-      directions.resize(npoint);
-      displacements.resize(npoint);
-      for(int d=0;d<_d;d++){
-	directions[2*d  ] = d+base;
-	directions[2*d+1] = d+base;
-	displacements[2*d  ] = +1;
-	displacements[2*d+1] = -1;
-      }
-      directions   [2*_d]=0;
-      displacements[2*_d]=0;
+    npoint = 2*_d+1;
+    directions.resize(npoint);
+    displacements.resize(npoint);
+    for(int d=0;d<_d;d++){
+      directions[2*d  ] = d+base;
+      directions[2*d+1] = d+base;
+      displacements[2*d  ] = +1;
+      displacements[2*d+1] = -1;
+    }
+    directions   [2*_d]=0;
+    displacements[2*_d]=0;
       
-      //// report back
-      std::cout<<GridLogMessage<<"directions    :";
-      for(int d=0;d<npoint;d++) std::cout<< directions[d]<< " ";
-      std::cout <<std::endl;
-      std::cout<<GridLogMessage<<"displacements :";
-      for(int d=0;d<npoint;d++) std::cout<< displacements[d]<< " ";
-      std::cout<<std::endl;
-    }
+    //// report back
+    std::cout<<GridLogMessage<<"directions    :";
+    for(int d=0;d<npoint;d++) std::cout<< directions[d]<< " ";
+    std::cout <<std::endl;
+    std::cout<<GridLogMessage<<"displacements :";
+    for(int d=0;d<npoint;d++) std::cout<< displacements[d]<< " ";
+    std::cout<<std::endl;
+  }
   
-    /*
-      // Original cleaner code
-    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;
-    };
-    */    
-
+  /*
+  // Original cleaner code
+  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;
   };
+  */    
+
+};
   
-  template<class Fobj,class CComplex,int nbasis>
-  class Aggregation   {
-  public:
-    typedef iVector<CComplex,nbasis >             siteVector;
-    typedef Lattice<siteVector>                 CoarseVector;
-    typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
+template<class Fobj,class CComplex,int nbasis>
+class Aggregation   {
+public:
+  typedef iVector<CComplex,nbasis >             siteVector;
+  typedef Lattice<siteVector>                 CoarseVector;
+  typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
 
-    typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
-    typedef Lattice<Fobj >        FineField;
+  typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
+  typedef Lattice<Fobj >        FineField;
 
-    GridBase *CoarseGrid;
-    GridBase *FineGrid;
-    std::vector<Lattice<Fobj> > subspace;
-    int checkerboard;
+  GridBase *CoarseGrid;
+  GridBase *FineGrid;
+  std::vector<Lattice<Fobj> > subspace;
+  int checkerboard;
 
   Aggregation(GridBase *_CoarseGrid,GridBase *_FineGrid,int _checkerboard) : 
     CoarseGrid(_CoarseGrid),
-      FineGrid(_FineGrid),
-      subspace(nbasis,_FineGrid),
-      checkerboard(_checkerboard)
-	{
-	};
+    FineGrid(_FineGrid),
+    subspace(nbasis,_FineGrid),
+    checkerboard(_checkerboard)
+  {
+  };
   
-    void Orthogonalise(void){
-      CoarseScalar InnerProd(CoarseGrid); 
-      std::cout << GridLogMessage <<" Gramm-Schmidt pass 1"<<std::endl;
-      blockOrthogonalise(InnerProd,subspace);
-      std::cout << GridLogMessage <<" Gramm-Schmidt pass 2"<<std::endl;
-      blockOrthogonalise(InnerProd,subspace);
-      //      std::cout << GridLogMessage <<" Gramm-Schmidt checking orthogonality"<<std::endl;
-      //      CheckOrthogonal();
-    } 
-    void CheckOrthogonal(void){
-      CoarseVector iProj(CoarseGrid); 
-      CoarseVector eProj(CoarseGrid); 
-      for(int i=0;i<nbasis;i++){
-	blockProject(iProj,subspace[i],subspace);
-	eProj=zero; 
-	parallel_for(int ss=0;ss<CoarseGrid->oSites();ss++){
-	  eProj._odata[ss](i)=CComplex(1.0);
-	}
-	eProj=eProj - iProj;
-	std::cout<<GridLogMessage<<"Orthog check error "<<i<<" " << norm2(eProj)<<std::endl;
+  void Orthogonalise(void){
+    CoarseScalar InnerProd(CoarseGrid); 
+    std::cout << GridLogMessage <<" Gramm-Schmidt pass 1"<<std::endl;
+    blockOrthogonalise(InnerProd,subspace);
+    std::cout << GridLogMessage <<" Gramm-Schmidt pass 2"<<std::endl;
+    blockOrthogonalise(InnerProd,subspace);
+    //      std::cout << GridLogMessage <<" Gramm-Schmidt checking orthogonality"<<std::endl;
+    //      CheckOrthogonal();
+  } 
+  void CheckOrthogonal(void){
+    CoarseVector iProj(CoarseGrid); 
+    CoarseVector eProj(CoarseGrid); 
+    for(int i=0;i<nbasis;i++){
+      blockProject(iProj,subspace[i],subspace);
+      eProj=zero; 
+      parallel_for(int ss=0;ss<CoarseGrid->oSites();ss++){
+	eProj._odata[ss](i)=CComplex(1.0);
       }
-      std::cout<<GridLogMessage <<"CheckOrthog done"<<std::endl;
+      eProj=eProj - iProj;
+      std::cout<<GridLogMessage<<"Orthog check error "<<i<<" " << norm2(eProj)<<std::endl;
     }
-    void ProjectToSubspace(CoarseVector &CoarseVec,const FineField &FineVec){
-      blockProject(CoarseVec,FineVec,subspace);
-    }
-    void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){
-      FineVec.checkerboard = subspace[0].checkerboard;
-      blockPromote(CoarseVec,FineVec,subspace);
-    }
-    void CreateSubspaceRandom(GridParallelRNG &RNG){
-      for(int i=0;i<nbasis;i++){
-	random(RNG,subspace[i]);
-	std::cout<<GridLogMessage<<" norm subspace["<<i<<"] "<<norm2(subspace[i])<<std::endl;
-      }
-      Orthogonalise();
+    std::cout<<GridLogMessage <<"CheckOrthog done"<<std::endl;
+  }
+  void ProjectToSubspace(CoarseVector &CoarseVec,const FineField &FineVec){
+    blockProject(CoarseVec,FineVec,subspace);
+  }
+  void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){
+    FineVec.checkerboard = subspace[0].checkerboard;
+    blockPromote(CoarseVec,FineVec,subspace);
+  }
+  void CreateSubspaceRandom(GridParallelRNG &RNG){
+    for(int i=0;i<nbasis;i++){
+      random(RNG,subspace[i]);
+      std::cout<<GridLogMessage<<" norm subspace["<<i<<"] "<<norm2(subspace[i])<<std::endl;
     }
+    Orthogonalise();
+  }
 
-    /*
+  /*
     virtual void CreateSubspaceLanczos(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) 
     {
-      // Run a Lanczos with sloppy convergence
-	const int Nstop = nn;
-	const int Nk = nn+20;
-	const int Np = nn+20;
-	const int Nm = Nk+Np;
-	const int MaxIt= 10000;
-	RealD resid = 1.0e-3;
+    // Run a Lanczos with sloppy convergence
+    const int Nstop = nn;
+    const int Nk = nn+20;
+    const int Np = nn+20;
+    const int Nm = Nk+Np;
+    const int MaxIt= 10000;
+    RealD resid = 1.0e-3;
 
-	Chebyshev<FineField> Cheb(0.5,64.0,21);
-	ImplicitlyRestartedLanczos<FineField> IRL(hermop,Cheb,Nstop,Nk,Nm,resid,MaxIt);
-	//	IRL.lock = 1;
+    Chebyshev<FineField> Cheb(0.5,64.0,21);
+    ImplicitlyRestartedLanczos<FineField> IRL(hermop,Cheb,Nstop,Nk,Nm,resid,MaxIt);
+    //	IRL.lock = 1;
 
-	FineField noise(FineGrid); gaussian(RNG,noise);
-	FineField tmp(FineGrid); 
-	std::vector<RealD>     eval(Nm);
-	std::vector<FineField> evec(Nm,FineGrid);
+    FineField noise(FineGrid); gaussian(RNG,noise);
+    FineField tmp(FineGrid); 
+    std::vector<RealD>     eval(Nm);
+    std::vector<FineField> evec(Nm,FineGrid);
 
-	int Nconv;
-	IRL.calc(eval,evec,
-		 noise,
-		 Nconv);
+    int Nconv;
+    IRL.calc(eval,evec,
+    noise,
+    Nconv);
 
-    	// pull back nn vectors
-	for(int b=0;b<nn;b++){
+    // pull back nn vectors
+    for(int b=0;b<nn;b++){
 
-	  subspace[b]   = evec[b];
+    subspace[b]   = evec[b];
 
-	  std::cout << GridLogMessage <<"subspace["<<b<<"] = "<<norm2(subspace[b])<<std::endl;
+    std::cout << GridLogMessage <<"subspace["<<b<<"] = "<<norm2(subspace[b])<<std::endl;
 
-	  hermop.Op(subspace[b],tmp); 
-	  std::cout<<GridLogMessage << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(tmp)<<std::endl;
+    hermop.Op(subspace[b],tmp); 
+    std::cout<<GridLogMessage << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(tmp)<<std::endl;
 
-	  noise = tmp -  sqrt(eval[b])*subspace[b] ;
+    noise = tmp -  sqrt(eval[b])*subspace[b] ;
 
-	  std::cout<<GridLogMessage << " lambda_"<<b<<" = "<< eval[b] <<"  ;  [ M - Lambda ]_"<<b<<" vec_"<<b<<"  = " <<norm2(noise)<<std::endl;
+    std::cout<<GridLogMessage << " lambda_"<<b<<" = "<< eval[b] <<"  ;  [ M - Lambda ]_"<<b<<" vec_"<<b<<"  = " <<norm2(noise)<<std::endl;
 
-	  noise = tmp +  eval[b]*subspace[b] ;
+    noise = tmp +  eval[b]*subspace[b] ;
 
-	  std::cout<<GridLogMessage << " lambda_"<<b<<" = "<< eval[b] <<"  ;  [ M - Lambda ]_"<<b<<" vec_"<<b<<"  = " <<norm2(noise)<<std::endl;
+    std::cout<<GridLogMessage << " lambda_"<<b<<" = "<< eval[b] <<"  ;  [ M - Lambda ]_"<<b<<" vec_"<<b<<"  = " <<norm2(noise)<<std::endl;
 
-	}
-	Orthogonalise();
-	for(int b=0;b<nn;b++){
-	  std::cout << GridLogMessage <<"subspace["<<b<<"] = "<<norm2(subspace[b])<<std::endl;
-	}
     }
-    */
-    virtual void CreateSubspace(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) {
+    Orthogonalise();
+    for(int b=0;b<nn;b++){
+    std::cout << GridLogMessage <<"subspace["<<b<<"] = "<<norm2(subspace[b])<<std::endl;
+    }
+    }
+  */
+  virtual void CreateSubspace(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) {
 
-      RealD scale;
+    RealD scale;
 
-      ConjugateGradient<FineField> CG(1.0e-2,10000);
-      FineField noise(FineGrid);
-      FineField Mn(FineGrid);
+    ConjugateGradient<FineField> CG(1.0e-2,10000);
+    FineField noise(FineGrid);
+    FineField Mn(FineGrid);
 
-      for(int b=0;b<nn;b++){
+    for(int b=0;b<nn;b++){
 	
-	gaussian(RNG,noise);
+      gaussian(RNG,noise);
+      scale = std::pow(norm2(noise),-0.5); 
+      noise=noise*scale;
+
+      hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise   ["<<b<<"] <n|MdagM|n> "<<norm2(Mn)<<std::endl;
+
+      for(int i=0;i<1;i++){
+
+	CG(hermop,noise,subspace[b]);
+
+	noise = subspace[b];
 	scale = std::pow(norm2(noise),-0.5); 
 	noise=noise*scale;
 
-	hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise   ["<<b<<"] <n|MdagM|n> "<<norm2(Mn)<<std::endl;
-
-	for(int i=0;i<1;i++){
-
-	  CG(hermop,noise,subspace[b]);
-
-	  noise = subspace[b];
-	  scale = std::pow(norm2(noise),-0.5); 
-	  noise=noise*scale;
-
-	}
-
-	hermop.Op(noise,Mn); std::cout<<GridLogMessage << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(Mn)<<std::endl;
-	subspace[b]   = noise;
-
       }
 
-      Orthogonalise();
+      hermop.Op(noise,Mn); std::cout<<GridLogMessage << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(Mn)<<std::endl;
+      subspace[b]   = noise;
 
     }
-  };
-  // 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<CComplex,nbasis > > >  {
-  public:
+
+    Orthogonalise();
+
+  }
+};
+// 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<CComplex,nbasis > > >  {
+public:
     
-    typedef iVector<CComplex,nbasis >             siteVector;
-    typedef Lattice<siteVector>                 CoarseVector;
-    typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
+  typedef iVector<CComplex,nbasis >             siteVector;
+  typedef Lattice<siteVector>                 CoarseVector;
+  typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
 
-    typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
-    typedef Lattice<Fobj >        FineField;
+  typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
+  typedef Lattice<Fobj >        FineField;
 
-    ////////////////////
-    // Data members
-    ////////////////////
-    Geometry         geom;
-    GridBase *       _grid; 
-    CartesianStencil<siteVector,siteVector> Stencil; 
+  ////////////////////
+  // Data members
+  ////////////////////
+  Geometry         geom;
+  GridBase *       _grid; 
+  CartesianStencil<siteVector,siteVector> Stencil; 
 
-    std::vector<CoarseMatrix> A;
+  std::vector<CoarseMatrix> A;
 
       
-    ///////////////////////
-    // Interface
-    ///////////////////////
-    GridBase * Grid(void)         { return _grid; };   // this is all the linalg routines need to know
+  ///////////////////////
+  // Interface
+  ///////////////////////
+  GridBase * Grid(void)         { return _grid; };   // this is all the linalg routines need to know
 
-    RealD M (const CoarseVector &in, CoarseVector &out){
+  RealD M (const CoarseVector &in, CoarseVector &out){
 
-      conformable(_grid,in._grid);
-      conformable(in._grid,out._grid);
+    conformable(_grid,in._grid);
+    conformable(in._grid,out._grid);
 
-      SimpleCompressor<siteVector> compressor;
-      Stencil.HaloExchange(in,compressor);
+    SimpleCompressor<siteVector> compressor;
+    Stencil.HaloExchange(in,compressor);
 
-      parallel_for(int ss=0;ss<Grid()->oSites();ss++){
-        siteVector res = zero;
-	siteVector nbr;
-	int ptype;
-	StencilEntry *SE;
-	for(int point=0;point<geom.npoint;point++){
+    parallel_for(int ss=0;ss<Grid()->oSites();ss++){
+      siteVector res = zero;
+      siteVector nbr;
+      int ptype;
+      StencilEntry *SE;
+      for(int point=0;point<geom.npoint;point++){
 
-	  SE=Stencil.GetEntry(ptype,point,ss);
+	SE=Stencil.GetEntry(ptype,point,ss);
 	  
-	  if(SE->_is_local&&SE->_permute) { 
-	    permute(nbr,in._odata[SE->_offset],ptype);
-	  } else if(SE->_is_local) { 
-	    nbr = in._odata[SE->_offset];
-	  } else {
-	    nbr = Stencil.CommBuf()[SE->_offset];
-	  }
-	  res = res + A[point]._odata[ss]*nbr;
+	if(SE->_is_local&&SE->_permute) { 
+	  permute(nbr,in._odata[SE->_offset],ptype);
+	} else if(SE->_is_local) { 
+	  nbr = in._odata[SE->_offset];
+	} else {
+	  nbr = Stencil.CommBuf()[SE->_offset];
 	}
-	vstream(out._odata[ss],res);
+	res = res + A[point]._odata[ss]*nbr;
       }
-      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)
-    {
-    };
-
-    void CoarsenOperator(GridBase *FineGrid,LinearOperatorBase<Lattice<Fobj> > &linop,
-			 Aggregation<Fobj,CComplex,nbasis> & Subspace){
-
-      FineField iblock(FineGrid); // contributions from within this block
-      FineField oblock(FineGrid); // contributions from outwith this block
-
-      FineField     phi(FineGrid);
-      FineField     tmp(FineGrid);
-      FineField     zz(FineGrid); zz=zero;
-      FineField    Mphi(FineGrid);
-
-      Lattice<iScalar<vInteger> > coor(FineGrid);
-
-      CoarseVector iProj(Grid()); 
-      CoarseVector oProj(Grid()); 
-      CoarseScalar InnerProd(Grid()); 
-
-      // Orthogonalise the subblocks over the basis
-      blockOrthogonalise(InnerProd,Subspace.subspace);
-
-      // Compute the matrix elements of linop between this orthonormal
-      // set of vectors.
-      int self_stencil=-1;
-      for(int p=0;p<geom.npoint;p++){ 
-	A[p]=zero;
-	if( geom.displacements[p]==0){
-	  self_stencil=p;
-	}
-      }
-      assert(self_stencil!=-1);
-
-      for(int i=0;i<nbasis;i++){
-	phi=Subspace.subspace[i];
-	
-	std::cout<<GridLogMessage<<"("<<i<<").."<<std::endl;
-
-	for(int p=0;p<geom.npoint;p++){ 
-
-	  int dir   = geom.directions[p];
-	  int disp  = geom.displacements[p];
-
-	  Integer block=(FineGrid->_rdimensions[dir])/(Grid()->_rdimensions[dir]);
-
-	  LatticeCoordinate(coor,dir);
-
-	  if ( disp==0 ){
-	    linop.OpDiag(phi,Mphi);
-	  }
-	  else  {
-	    linop.OpDir(phi,Mphi,dir,disp); 
-	  }
-
-	  ////////////////////////////////////////////////////////////////////////
-	  // Pick out contributions coming from this cell and neighbour cell
-	  ////////////////////////////////////////////////////////////////////////
-	  if ( disp==0 ) {
-	    iblock = Mphi;
-	    oblock = zero;
-	  } else if ( disp==1 ) {
-	    oblock = where(mod(coor,block)==(block-1),Mphi,zz);
-	    iblock = where(mod(coor,block)!=(block-1),Mphi,zz);
-	  } else if ( disp==-1 ) {
-	    oblock = where(mod(coor,block)==(Integer)0,Mphi,zz);
-	    iblock = where(mod(coor,block)!=(Integer)0,Mphi,zz);
-	  } else {
-	    assert(0);
-	  }
-
-	  Subspace.ProjectToSubspace(iProj,iblock);
-	  Subspace.ProjectToSubspace(oProj,oblock);
-	  //	  blockProject(iProj,iblock,Subspace.subspace);
-	  //	  blockProject(oProj,oblock,Subspace.subspace);
-	  parallel_for(int ss=0;ss<Grid()->oSites();ss++){
-	    for(int j=0;j<nbasis;j++){
-	      if( disp!= 0 ) {
-		A[p]._odata[ss](j,i) = oProj._odata[ss](j);
-	      }
-	      A[self_stencil]._odata[ss](j,i) =	A[self_stencil]._odata[ss](j,i) + iProj._odata[ss](j);
-	    }
-	  }
-	}
-      }
-
-#if 0
-      ///////////////////////////
-      // test code worth preserving in if block
-      ///////////////////////////
-      std::cout<<GridLogMessage<< " Computed matrix elements "<< self_stencil <<std::endl;
-      for(int p=0;p<geom.npoint;p++){
-	std::cout<<GridLogMessage<< "A["<<p<<"]" << std::endl;
-	std::cout<<GridLogMessage<< A[p] << std::endl;
-      }
-      std::cout<<GridLogMessage<< " picking by block0 "<< self_stencil <<std::endl;
-
-      phi=Subspace.subspace[0];
-      std::vector<int> bc(FineGrid->_ndimension,0);
-
-      blockPick(Grid(),phi,tmp,bc);      // Pick out a block
-      linop.Op(tmp,Mphi);                // Apply big dop
-      blockProject(iProj,Mphi,Subspace.subspace); // project it and print it
-      std::cout<<GridLogMessage<< " Computed matrix elements from block zero only "<<std::endl;
-      std::cout<<GridLogMessage<< iProj <<std::endl;
-      std::cout<<GridLogMessage<<"Computed Coarse Operator"<<std::endl;
-#endif
-      //      ForceHermitian();
-      AssertHermitian();
-      // ForceDiagonal();
+      vstream(out._odata[ss],res);
     }
-    void ForceDiagonal(void) {
-
-
-      std::cout<<GridLogMessage<<"**************************************************"<<std::endl;
-      std::cout<<GridLogMessage<<"****   Forcing coarse operator to be diagonal ****"<<std::endl;
-      std::cout<<GridLogMessage<<"**************************************************"<<std::endl;
-      for(int p=0;p<8;p++){
-	A[p]=zero;
-      }
-
-      GridParallelRNG  RNG(Grid()); RNG.SeedFixedIntegers(std::vector<int>({55,72,19,17,34}));
-      Lattice<iScalar<CComplex> > val(Grid()); random(RNG,val);
-
-      Complex one(1.0);
-
-      iMatrix<CComplex,nbasis> ident;  ident=one;
-
-      val = val*adj(val);
-      val = val + 1.0;
-
-      A[8] = val*ident;
-
-      //      for(int s=0;s<Grid()->oSites();s++) {
-      //	A[8]._odata[s]=val._odata[s];
-      //      }
-    }
-    void ForceHermitian(void) {
-      for(int d=0;d<4;d++){
-	int dd=d+1;
-	A[2*d] = adj(Cshift(A[2*d+1],dd,1));
-      }
-      //      A[8] = 0.5*(A[8] + adj(A[8]));
-    }
-    void AssertHermitian(void) {
-      CoarseMatrix AA    (Grid());
-      CoarseMatrix AAc   (Grid());
-      CoarseMatrix Diff  (Grid());
-      for(int d=0;d<4;d++){
-	
-	int dd=d+1;
-	AAc = Cshift(A[2*d+1],dd,1);
-	AA  = A[2*d];
-	
-	Diff = AA - adj(AAc);
-
-	std::cout<<GridLogMessage<<"Norm diff dim "<<d<<" "<< norm2(Diff)<<std::endl;
-	std::cout<<GridLogMessage<<"Norm dim "<<d<<" "<< norm2(AA)<<std::endl;
-	  
-      }
-      Diff = A[8] - adj(A[8]);
-      std::cout<<GridLogMessage<<"Norm diff local "<< norm2(Diff)<<std::endl;
-      std::cout<<GridLogMessage<<"Norm local "<< norm2(A[8])<<std::endl;
-    }
-    
+    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)
+  {
+  };
+
+  void CoarsenOperator(GridBase *FineGrid,LinearOperatorBase<Lattice<Fobj> > &linop,
+		       Aggregation<Fobj,CComplex,nbasis> & Subspace){
+
+    FineField iblock(FineGrid); // contributions from within this block
+    FineField oblock(FineGrid); // contributions from outwith this block
+
+    FineField     phi(FineGrid);
+    FineField     tmp(FineGrid);
+    FineField     zz(FineGrid); zz=zero;
+    FineField    Mphi(FineGrid);
+
+    Lattice<iScalar<vInteger> > coor(FineGrid);
+
+    CoarseVector iProj(Grid()); 
+    CoarseVector oProj(Grid()); 
+    CoarseScalar InnerProd(Grid()); 
+
+    // Orthogonalise the subblocks over the basis
+    blockOrthogonalise(InnerProd,Subspace.subspace);
+
+    // Compute the matrix elements of linop between this orthonormal
+    // set of vectors.
+    int self_stencil=-1;
+    for(int p=0;p<geom.npoint;p++){ 
+      A[p]=zero;
+      if( geom.displacements[p]==0){
+	self_stencil=p;
+      }
+    }
+    assert(self_stencil!=-1);
+
+    for(int i=0;i<nbasis;i++){
+      phi=Subspace.subspace[i];
+	
+      std::cout<<GridLogMessage<<"("<<i<<").."<<std::endl;
+
+      for(int p=0;p<geom.npoint;p++){ 
+
+	int dir   = geom.directions[p];
+	int disp  = geom.displacements[p];
+
+	Integer block=(FineGrid->_rdimensions[dir])/(Grid()->_rdimensions[dir]);
+
+	LatticeCoordinate(coor,dir);
+
+	if ( disp==0 ){
+	  linop.OpDiag(phi,Mphi);
+	}
+	else  {
+	  linop.OpDir(phi,Mphi,dir,disp); 
+	}
+
+	////////////////////////////////////////////////////////////////////////
+	// Pick out contributions coming from this cell and neighbour cell
+	////////////////////////////////////////////////////////////////////////
+	if ( disp==0 ) {
+	  iblock = Mphi;
+	  oblock = zero;
+	} else if ( disp==1 ) {
+	  oblock = where(mod(coor,block)==(block-1),Mphi,zz);
+	  iblock = where(mod(coor,block)!=(block-1),Mphi,zz);
+	} else if ( disp==-1 ) {
+	  oblock = where(mod(coor,block)==(Integer)0,Mphi,zz);
+	  iblock = where(mod(coor,block)!=(Integer)0,Mphi,zz);
+	} else {
+	  assert(0);
+	}
+
+	Subspace.ProjectToSubspace(iProj,iblock);
+	Subspace.ProjectToSubspace(oProj,oblock);
+	//	  blockProject(iProj,iblock,Subspace.subspace);
+	//	  blockProject(oProj,oblock,Subspace.subspace);
+	parallel_for(int ss=0;ss<Grid()->oSites();ss++){
+	  for(int j=0;j<nbasis;j++){
+	    if( disp!= 0 ) {
+	      A[p]._odata[ss](j,i) = oProj._odata[ss](j);
+	    }
+	    A[self_stencil]._odata[ss](j,i) =	A[self_stencil]._odata[ss](j,i) + iProj._odata[ss](j);
+	  }
+	}
+      }
+    }
+
+#if 0
+    ///////////////////////////
+    // test code worth preserving in if block
+    ///////////////////////////
+    std::cout<<GridLogMessage<< " Computed matrix elements "<< self_stencil <<std::endl;
+    for(int p=0;p<geom.npoint;p++){
+      std::cout<<GridLogMessage<< "A["<<p<<"]" << std::endl;
+      std::cout<<GridLogMessage<< A[p] << std::endl;
+    }
+    std::cout<<GridLogMessage<< " picking by block0 "<< self_stencil <<std::endl;
+
+    phi=Subspace.subspace[0];
+    std::vector<int> bc(FineGrid->_ndimension,0);
+
+    blockPick(Grid(),phi,tmp,bc);      // Pick out a block
+    linop.Op(tmp,Mphi);                // Apply big dop
+    blockProject(iProj,Mphi,Subspace.subspace); // project it and print it
+    std::cout<<GridLogMessage<< " Computed matrix elements from block zero only "<<std::endl;
+    std::cout<<GridLogMessage<< iProj <<std::endl;
+    std::cout<<GridLogMessage<<"Computed Coarse Operator"<<std::endl;
+#endif
+    //      ForceHermitian();
+    AssertHermitian();
+    // ForceDiagonal();
+  }
+  void ForceDiagonal(void) {
+
+
+    std::cout<<GridLogMessage<<"**************************************************"<<std::endl;
+    std::cout<<GridLogMessage<<"****   Forcing coarse operator to be diagonal ****"<<std::endl;
+    std::cout<<GridLogMessage<<"**************************************************"<<std::endl;
+    for(int p=0;p<8;p++){
+      A[p]=zero;
+    }
+
+    GridParallelRNG  RNG(Grid()); RNG.SeedFixedIntegers(std::vector<int>({55,72,19,17,34}));
+    Lattice<iScalar<CComplex> > val(Grid()); random(RNG,val);
+
+    Complex one(1.0);
+
+    iMatrix<CComplex,nbasis> ident;  ident=one;
+
+    val = val*adj(val);
+    val = val + 1.0;
+
+    A[8] = val*ident;
+
+    //      for(int s=0;s<Grid()->oSites();s++) {
+    //	A[8]._odata[s]=val._odata[s];
+    //      }
+  }
+  void ForceHermitian(void) {
+    for(int d=0;d<4;d++){
+      int dd=d+1;
+      A[2*d] = adj(Cshift(A[2*d+1],dd,1));
+    }
+    //      A[8] = 0.5*(A[8] + adj(A[8]));
+  }
+  void AssertHermitian(void) {
+    CoarseMatrix AA    (Grid());
+    CoarseMatrix AAc   (Grid());
+    CoarseMatrix Diff  (Grid());
+    for(int d=0;d<4;d++){
+	
+      int dd=d+1;
+      AAc = Cshift(A[2*d+1],dd,1);
+      AA  = A[2*d];
+	
+      Diff = AA - adj(AAc);
+
+      std::cout<<GridLogMessage<<"Norm diff dim "<<d<<" "<< norm2(Diff)<<std::endl;
+      std::cout<<GridLogMessage<<"Norm dim "<<d<<" "<< norm2(AA)<<std::endl;
+	  
+    }
+    Diff = A[8] - adj(A[8]);
+    std::cout<<GridLogMessage<<"Norm diff local "<< norm2(Diff)<<std::endl;
+    std::cout<<GridLogMessage<<"Norm local "<< norm2(A[8])<<std::endl;
+  }
+    
+};
+
+NAMESPACE_END(Grid);
 #endif
diff --git a/lib/algorithms/FFT.h b/lib/algorithms/FFT.h
index ec558ad9..e92540fb 100644
--- a/lib/algorithms/FFT.h
+++ b/lib/algorithms/FFT.h
@@ -1,5 +1,5 @@
 
-    /*************************************************************************************
+/*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -24,8 +24,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
     See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
+*************************************************************************************/
+/*  END LEGAL */
 #ifndef _GRID_FFT_H_
 #define _GRID_FFT_H_
 
@@ -38,64 +38,64 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #endif
 
 
-namespace Grid {
+NAMESPACE_BEGIN(Grid);
 
-  template<class scalar> struct FFTW { };
+template<class scalar> struct FFTW { };
 
 #ifdef HAVE_FFTW	
-  template<> struct FFTW<ComplexD> {
-  public:
+template<> struct FFTW<ComplexD> {
+public:
 
-    typedef fftw_complex FFTW_scalar;
-    typedef fftw_plan    FFTW_plan;
+  typedef fftw_complex FFTW_scalar;
+  typedef fftw_plan    FFTW_plan;
 
-    static FFTW_plan fftw_plan_many_dft(int rank, const int *n,int howmany,
-					FFTW_scalar *in, const int *inembed,		
-					int istride, int idist,		
-					FFTW_scalar *out, const int *onembed,		
-					int ostride, int odist,		
-					int sign, unsigned flags) {
-      return ::fftw_plan_many_dft(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,sign,flags);
-    }	  
+  static FFTW_plan fftw_plan_many_dft(int rank, const int *n,int howmany,
+				      FFTW_scalar *in, const int *inembed,		
+				      int istride, int idist,		
+				      FFTW_scalar *out, const int *onembed,		
+				      int ostride, int odist,		
+				      int sign, unsigned flags) {
+    return ::fftw_plan_many_dft(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,sign,flags);
+  }	  
     
-    static void fftw_flops(const FFTW_plan p,double *add, double *mul, double *fmas){
-      ::fftw_flops(p,add,mul,fmas);
-    }
+  static void fftw_flops(const FFTW_plan p,double *add, double *mul, double *fmas){
+    ::fftw_flops(p,add,mul,fmas);
+  }
 
-    inline static void fftw_execute_dft(const FFTW_plan p,FFTW_scalar *in,FFTW_scalar *out) {
-      ::fftw_execute_dft(p,in,out);
-    }
-    inline static void fftw_destroy_plan(const FFTW_plan p) {
-      ::fftw_destroy_plan(p);
-    }
-  };
+  inline static void fftw_execute_dft(const FFTW_plan p,FFTW_scalar *in,FFTW_scalar *out) {
+    ::fftw_execute_dft(p,in,out);
+  }
+  inline static void fftw_destroy_plan(const FFTW_plan p) {
+    ::fftw_destroy_plan(p);
+  }
+};
 
-  template<> struct FFTW<ComplexF> {
-  public:
+template<> struct FFTW<ComplexF> {
+public:
 
-    typedef fftwf_complex FFTW_scalar;
-    typedef fftwf_plan    FFTW_plan;
+  typedef fftwf_complex FFTW_scalar;
+  typedef fftwf_plan    FFTW_plan;
 
-    static FFTW_plan fftw_plan_many_dft(int rank, const int *n,int howmany,
-					FFTW_scalar *in, const int *inembed,		
-					int istride, int idist,		
-					FFTW_scalar *out, const int *onembed,		
-					int ostride, int odist,		
-					int sign, unsigned flags) {
-      return ::fftwf_plan_many_dft(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,sign,flags);
-    }	  
+  static FFTW_plan fftw_plan_many_dft(int rank, const int *n,int howmany,
+				      FFTW_scalar *in, const int *inembed,		
+				      int istride, int idist,		
+				      FFTW_scalar *out, const int *onembed,		
+				      int ostride, int odist,		
+				      int sign, unsigned flags) {
+    return ::fftwf_plan_many_dft(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,sign,flags);
+  }	  
     
-    static void fftw_flops(const FFTW_plan p,double *add, double *mul, double *fmas){
-      ::fftwf_flops(p,add,mul,fmas);
-    }
+  static void fftw_flops(const FFTW_plan p,double *add, double *mul, double *fmas){
+    ::fftwf_flops(p,add,mul,fmas);
+  }
 
-    inline static void fftw_execute_dft(const FFTW_plan p,FFTW_scalar *in,FFTW_scalar *out) {
-      ::fftwf_execute_dft(p,in,out);
-    }
-    inline static void fftw_destroy_plan(const FFTW_plan p) {
-      ::fftwf_destroy_plan(p);
-    }
-  };
+  inline static void fftw_execute_dft(const FFTW_plan p,FFTW_scalar *in,FFTW_scalar *out) {
+    ::fftwf_execute_dft(p,in,out);
+  }
+  inline static void fftw_destroy_plan(const FFTW_plan p) {
+    ::fftwf_destroy_plan(p);
+  }
+};
 
 #endif
 
@@ -104,203 +104,204 @@ namespace Grid {
 #define FFTW_BACKWARD (+1)
 #endif
 
-  class FFT {
-  private:
+class FFT {
+private:
     
-    GridCartesian *vgrid;
-    GridCartesian *sgrid;
+  GridCartesian *vgrid;
+  GridCartesian *sgrid;
     
-    int Nd;
-    double flops;
-    double flops_call;
-    uint64_t usec;
+  int Nd;
+  double flops;
+  double flops_call;
+  uint64_t usec;
     
-    std::vector<int> dimensions;
-    std::vector<int> processors;
-    std::vector<int> processor_coor;
+  std::vector<int> dimensions;
+  std::vector<int> processors;
+  std::vector<int> processor_coor;
     
-  public:
+public:
     
-    static const int forward=FFTW_FORWARD;
-    static const int backward=FFTW_BACKWARD;
+  static const int forward=FFTW_FORWARD;
+  static const int backward=FFTW_BACKWARD;
     
-    double Flops(void) {return flops;}
-    double MFlops(void) {return flops/usec;}
-    double USec(void)   {return (double)usec;}    
+  double Flops(void) {return flops;}
+  double MFlops(void) {return flops/usec;}
+  double USec(void)   {return (double)usec;}    
 
-    FFT ( GridCartesian * grid ) :
+  FFT ( GridCartesian * grid ) :
     vgrid(grid),
     Nd(grid->_ndimension),
     dimensions(grid->_fdimensions),
     processors(grid->_processors),
     processor_coor(grid->_processor_coor)
-    {
-      flops=0;
-      usec =0;
-      std::vector<int> layout(Nd,1);
-      sgrid = new GridCartesian(dimensions,layout,processors);
-    };
+  {
+    flops=0;
+    usec =0;
+    std::vector<int> layout(Nd,1);
+    sgrid = new GridCartesian(dimensions,layout,processors);
+  };
     
-    ~FFT ( void)  {
-      delete sgrid;
-    }
+  ~FFT ( void)  {
+    delete sgrid;
+  }
     
-    template<class vobj>
-    void FFT_dim_mask(Lattice<vobj> &result,const Lattice<vobj> &source,std::vector<int> mask,int sign){
+  template<class vobj>
+  void FFT_dim_mask(Lattice<vobj> &result,const Lattice<vobj> &source,std::vector<int> mask,int sign){
 
-      conformable(result._grid,vgrid);
-      conformable(source._grid,vgrid);
-      Lattice<vobj> tmp(vgrid);
-      tmp = source;
-      for(int d=0;d<Nd;d++){
-	if( mask[d] ) {
-	  FFT_dim(result,tmp,d,sign);
-	  tmp=result;
-	}
+    conformable(result._grid,vgrid);
+    conformable(source._grid,vgrid);
+    Lattice<vobj> tmp(vgrid);
+    tmp = source;
+    for(int d=0;d<Nd;d++){
+      if( mask[d] ) {
+	FFT_dim(result,tmp,d,sign);
+	tmp=result;
       }
     }
+  }
 
-    template<class vobj>
-    void FFT_all_dim(Lattice<vobj> &result,const Lattice<vobj> &source,int sign){
-      std::vector<int> mask(Nd,1);
-      FFT_dim_mask(result,source,mask,sign);
-    }
+  template<class vobj>
+  void FFT_all_dim(Lattice<vobj> &result,const Lattice<vobj> &source,int sign){
+    std::vector<int> mask(Nd,1);
+    FFT_dim_mask(result,source,mask,sign);
+  }
 
 
-    template<class vobj>
-    void FFT_dim(Lattice<vobj> &result,const Lattice<vobj> &source,int dim, int sign){
+  template<class vobj>
+  void FFT_dim(Lattice<vobj> &result,const Lattice<vobj> &source,int dim, int sign){
 #ifndef HAVE_FFTW
-      assert(0);
+    assert(0);
 #else
-      conformable(result._grid,vgrid);
-      conformable(source._grid,vgrid);
+    conformable(result._grid,vgrid);
+    conformable(source._grid,vgrid);
 
-      int L = vgrid->_ldimensions[dim];
-      int G = vgrid->_fdimensions[dim];
+    int L = vgrid->_ldimensions[dim];
+    int G = vgrid->_fdimensions[dim];
       
-      std::vector<int> layout(Nd,1);
-      std::vector<int> pencil_gd(vgrid->_fdimensions);
+    std::vector<int> layout(Nd,1);
+    std::vector<int> pencil_gd(vgrid->_fdimensions);
       
-      pencil_gd[dim] = G*processors[dim];
+    pencil_gd[dim] = G*processors[dim];
       
-      // Pencil global vol LxLxGxLxL per node
-      GridCartesian pencil_g(pencil_gd,layout,processors);
+    // Pencil global vol LxLxGxLxL per node
+    GridCartesian pencil_g(pencil_gd,layout,processors);
       
-      // Construct pencils
-      typedef typename vobj::scalar_object sobj;
-      typedef typename sobj::scalar_type   scalar;
+    // Construct pencils
+    typedef typename vobj::scalar_object sobj;
+    typedef typename sobj::scalar_type   scalar;
       
-      Lattice<sobj> pgbuf(&pencil_g);
+    Lattice<sobj> pgbuf(&pencil_g);
       
 
-      typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
-      typedef typename FFTW<scalar>::FFTW_plan   FFTW_plan;
+    typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
+    typedef typename FFTW<scalar>::FFTW_plan   FFTW_plan;
       
-      int Ncomp = sizeof(sobj)/sizeof(scalar);
-      int Nlow  = 1;
-      for(int d=0;d<dim;d++){
-        Nlow*=vgrid->_ldimensions[d];
-      }
+    int Ncomp = sizeof(sobj)/sizeof(scalar);
+    int Nlow  = 1;
+    for(int d=0;d<dim;d++){
+      Nlow*=vgrid->_ldimensions[d];
+    }
       
-      int rank = 1;  /* 1d transforms */
-      int n[] = {G}; /* 1d transforms of length G */
-      int howmany = Ncomp;
-      int odist,idist,istride,ostride;
-      idist   = odist   = 1;          /* Distance between consecutive FT's */
-      istride = ostride = Ncomp*Nlow; /* distance between two elements in the same FT */
-      int *inembed = n, *onembed = n;
+    int rank = 1;  /* 1d transforms */
+    int n[] = {G}; /* 1d transforms of length G */
+    int howmany = Ncomp;
+    int odist,idist,istride,ostride;
+    idist   = odist   = 1;          /* Distance between consecutive FT's */
+    istride = ostride = Ncomp*Nlow; /* distance between two elements in the same FT */
+    int *inembed = n, *onembed = n;
       
-      scalar div;
-	  if ( sign == backward ) div = 1.0/G;
-	  else if ( sign == forward ) div = 1.0;
-	  else assert(0);
+    scalar div;
+    if ( sign == backward ) div = 1.0/G;
+    else if ( sign == forward ) div = 1.0;
+    else assert(0);
       
-      FFTW_plan p;
-      {
-        FFTW_scalar *in = (FFTW_scalar *)&pgbuf._odata[0];
-        FFTW_scalar *out= (FFTW_scalar *)&pgbuf._odata[0];
-        p = FFTW<scalar>::fftw_plan_many_dft(rank,n,howmany,
-                                             in,inembed,
-                                             istride,idist,
-                                             out,onembed,
-                                             ostride, odist,
-                                             sign,FFTW_ESTIMATE);
-      }
+    FFTW_plan p;
+    {
+      FFTW_scalar *in = (FFTW_scalar *)&pgbuf._odata[0];
+      FFTW_scalar *out= (FFTW_scalar *)&pgbuf._odata[0];
+      p = FFTW<scalar>::fftw_plan_many_dft(rank,n,howmany,
+					   in,inembed,
+					   istride,idist,
+					   out,onembed,
+					   ostride, odist,
+					   sign,FFTW_ESTIMATE);
+    }
       
-      // Barrel shift and collect global pencil
-      std::vector<int> lcoor(Nd), gcoor(Nd);
-      result = source;
-      int pc = processor_coor[dim];
-      for(int p=0;p<processors[dim];p++) {
-        PARALLEL_REGION
+    // Barrel shift and collect global pencil
+    std::vector<int> lcoor(Nd), gcoor(Nd);
+    result = source;
+    int pc = processor_coor[dim];
+    for(int p=0;p<processors[dim];p++) {
+      PARALLEL_REGION
         {
           std::vector<int> cbuf(Nd);
           sobj s;
           
           PARALLEL_FOR_LOOP_INTERN
-          for(int idx=0;idx<sgrid->lSites();idx++) {
-            sgrid->LocalIndexToLocalCoor(idx,cbuf);
-            peekLocalSite(s,result,cbuf);
-	    cbuf[dim]+=((pc+p) % processors[dim])*L;
-	    //            cbuf[dim]+=p*L;
-            pokeLocalSite(s,pgbuf,cbuf);
-          }
+	    for(int idx=0;idx<sgrid->lSites();idx++) {
+	      sgrid->LocalIndexToLocalCoor(idx,cbuf);
+	      peekLocalSite(s,result,cbuf);
+	      cbuf[dim]+=((pc+p) % processors[dim])*L;
+	      //            cbuf[dim]+=p*L;
+	      pokeLocalSite(s,pgbuf,cbuf);
+	    }
         }
-        if (p != processors[dim] - 1)
+      if (p != processors[dim] - 1)
         {
           result = Cshift(result,dim,L);
         }
-      }
+    }
       
-      // Loop over orthog coords
-      int NN=pencil_g.lSites();
-      GridStopWatch timer;
-      timer.Start();
-      PARALLEL_REGION
+    // Loop over orthog coords
+    int NN=pencil_g.lSites();
+    GridStopWatch timer;
+    timer.Start();
+    PARALLEL_REGION
       {
         std::vector<int> cbuf(Nd);
         
         PARALLEL_FOR_LOOP_INTERN
-        for(int idx=0;idx<NN;idx++) {
-          pencil_g.LocalIndexToLocalCoor(idx, cbuf);
-          if ( cbuf[dim] == 0 ) {  // restricts loop to plane at lcoor[dim]==0
-            FFTW_scalar *in = (FFTW_scalar *)&pgbuf._odata[idx];
-            FFTW_scalar *out= (FFTW_scalar *)&pgbuf._odata[idx];
-            FFTW<scalar>::fftw_execute_dft(p,in,out);
-          }
-        }
+	  for(int idx=0;idx<NN;idx++) {
+	    pencil_g.LocalIndexToLocalCoor(idx, cbuf);
+	    if ( cbuf[dim] == 0 ) {  // restricts loop to plane at lcoor[dim]==0
+	      FFTW_scalar *in = (FFTW_scalar *)&pgbuf._odata[idx];
+	      FFTW_scalar *out= (FFTW_scalar *)&pgbuf._odata[idx];
+	      FFTW<scalar>::fftw_execute_dft(p,in,out);
+	    }
+	  }
       }
-      timer.Stop();
+    timer.Stop();
       
-      // performance counting
-      double add,mul,fma;
-      FFTW<scalar>::fftw_flops(p,&add,&mul,&fma);
-      flops_call = add+mul+2.0*fma;
-      usec += timer.useconds();
-      flops+= flops_call*NN;
+    // performance counting
+    double add,mul,fma;
+    FFTW<scalar>::fftw_flops(p,&add,&mul,&fma);
+    flops_call = add+mul+2.0*fma;
+    usec += timer.useconds();
+    flops+= flops_call*NN;
       
-      // writing out result
-      PARALLEL_REGION
+    // writing out result
+    PARALLEL_REGION
       {
         std::vector<int> clbuf(Nd), cgbuf(Nd);
         sobj s;
         
         PARALLEL_FOR_LOOP_INTERN
-        for(int idx=0;idx<sgrid->lSites();idx++) {
-          sgrid->LocalIndexToLocalCoor(idx,clbuf);
-          cgbuf = clbuf;
-          cgbuf[dim] = clbuf[dim]+L*pc;
-          peekLocalSite(s,pgbuf,cgbuf);
-          pokeLocalSite(s,result,clbuf);
-        }
+	  for(int idx=0;idx<sgrid->lSites();idx++) {
+	    sgrid->LocalIndexToLocalCoor(idx,clbuf);
+	    cgbuf = clbuf;
+	    cgbuf[dim] = clbuf[dim]+L*pc;
+	    peekLocalSite(s,pgbuf,cgbuf);
+	    pokeLocalSite(s,result,clbuf);
+	  }
       }
-      result = result*div;
+    result = result*div;
       
-      // destroying plan
-      FFTW<scalar>::fftw_destroy_plan(p);
+    // destroying plan
+    FFTW<scalar>::fftw_destroy_plan(p);
 #endif
-    }
-  };
-}
+  }
+};
+
+NAMESPACE_END(Grid);
 
 #endif
diff --git a/lib/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h
index 26746e6e..4bc169eb 100644
--- a/lib/algorithms/LinearOperator.h
+++ b/lib/algorithms/LinearOperator.h
@@ -1,4 +1,4 @@
-    /*************************************************************************************
+/*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -24,431 +24,431 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
     See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
+*************************************************************************************/
+/*  END LEGAL */
 #ifndef  GRID_ALGORITHM_LINEAR_OP_H
 #define  GRID_ALGORITHM_LINEAR_OP_H
 
-namespace Grid {
+NAMESPACE_BEGIN(Grid);
 
-  /////////////////////////////////////////////////////////////////////////////////////////////
-  // LinearOperators Take a something and return a something.
-  /////////////////////////////////////////////////////////////////////////////////////////////
-  //
-  // Hopefully linearity is satisfied and the AdjOp is indeed the Hermitian conjugateugate (transpose if real):
-  //SBase
-  //   i)  F(a x + b y) = aF(x) + b F(y).
-  //  ii)  <x|Op|y> = <y|AdjOp|x>^\ast
-  //
-  // Would be fun to have a test linearity & Herm Conj function!
-  /////////////////////////////////////////////////////////////////////////////////////////////
-    template<class Field> class LinearOperatorBase {
-    public:
+/////////////////////////////////////////////////////////////////////////////////////////////
+// LinearOperators Take a something and return a something.
+/////////////////////////////////////////////////////////////////////////////////////////////
+//
+// Hopefully linearity is satisfied and the AdjOp is indeed the Hermitian conjugateugate (transpose if real):
+//SBase
+//   i)  F(a x + b y) = aF(x) + b F(y).
+//  ii)  <x|Op|y> = <y|AdjOp|x>^\ast
+//
+// Would be fun to have a test linearity & Herm Conj function!
+/////////////////////////////////////////////////////////////////////////////////////////////
+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
+  // 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
 
-      virtual void Op     (const Field &in, Field &out) = 0; // Abstract base
-      virtual void AdjOp  (const Field &in, Field &out) = 0; // Abstract base
-      virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2)=0;
-      virtual void HermOp(const Field &in, Field &out)=0;
-    };
+  virtual void Op     (const Field &in, Field &out) = 0; // Abstract base
+  virtual void AdjOp  (const Field &in, Field &out) = 0; // Abstract base
+  virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2)=0;
+  virtual void HermOp(const Field &in, Field &out)=0;
+};
 
 
-  /////////////////////////////////////////////////////////////////////////////////////////////
-  // By sharing the class for Sparse Matrix across multiple operator wrappers, we can share code
-  // between RB and non-RB variants. Sparse matrix is like the fermion action def, and then
-  // the wrappers implement the specialisation of "Op" and "AdjOp" to the cases minimising
-  // replication of code.
-  //
-  // I'm not entirely happy with implementation; to share the Schur code between herm and non-herm
-  // while still having a "OpAndNorm" in the abstract base I had to implement it in both cases
-  // with an assert trap in the non-herm. This isn't right; there must be a better C++ way to
-  // do it, but I fear it required multiple inheritance and mixed in abstract base classes
-  /////////////////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////////////////////
+// By sharing the class for Sparse Matrix across multiple operator wrappers, we can share code
+// between RB and non-RB variants. Sparse matrix is like the fermion action def, and then
+// the wrappers implement the specialisation of "Op" and "AdjOp" to the cases minimising
+// replication of code.
+//
+// I'm not entirely happy with implementation; to share the Schur code between herm and non-herm
+// while still having a "OpAndNorm" in the abstract base I had to implement it in both cases
+// with an assert trap in the non-herm. This isn't right; there must be a better C++ way to
+// do it, but I fear it required multiple inheritance and mixed in abstract base classes
+/////////////////////////////////////////////////////////////////////////////////////////////
 
-    ////////////////////////////////////////////////////////////////////
-    // Construct herm op from non-herm matrix
-    ////////////////////////////////////////////////////////////////////
-    template<class Matrix,class Field>
-    class MdagMLinearOperator : public LinearOperatorBase<Field> {
-      Matrix &_Mat;
-    public:
-    MdagMLinearOperator(Matrix &Mat): _Mat(Mat){};
+////////////////////////////////////////////////////////////////////
+// Construct herm op from non-herm matrix
+////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field>
+class MdagMLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat;
+public:
+  MdagMLinearOperator(Matrix &Mat): _Mat(Mat){};
 
-      // Support for coarsening to a multigrid
-      void OpDiag (const Field &in, Field &out) {
-	_Mat.Mdiag(in,out);
-      }
-      void OpDir  (const Field &in, Field &out,int dir,int disp) {
-	_Mat.Mdir(in,out,dir,disp);
-      }
-      void Op     (const Field &in, Field &out){
-	_Mat.M(in,out);
-      }
-      void AdjOp     (const Field &in, Field &out){
-	_Mat.Mdag(in,out);
-      }
-      void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-	_Mat.MdagM(in,out,n1,n2);
-      }
-      void HermOp(const Field &in, Field &out){
-	RealD n1,n2;
-	HermOpAndNorm(in,out,n1,n2);
-      }
-    };
+  // Support for coarsening to a multigrid
+  void OpDiag (const Field &in, Field &out) {
+    _Mat.Mdiag(in,out);
+  }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) {
+    _Mat.Mdir(in,out,dir,disp);
+  }
+  void Op     (const Field &in, Field &out){
+    _Mat.M(in,out);
+  }
+  void AdjOp     (const Field &in, Field &out){
+    _Mat.Mdag(in,out);
+  }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+    _Mat.MdagM(in,out,n1,n2);
+  }
+  void HermOp(const Field &in, Field &out){
+    RealD n1,n2;
+    HermOpAndNorm(in,out,n1,n2);
+  }
+};
 
-    ////////////////////////////////////////////////////////////////////
-    // Construct herm op and shift it for mgrid smoother
-    ////////////////////////////////////////////////////////////////////
-    template<class Matrix,class Field>
-    class ShiftedMdagMLinearOperator : public LinearOperatorBase<Field> {
-      Matrix &_Mat;
-      RealD _shift;
-    public:
-    ShiftedMdagMLinearOperator(Matrix &Mat,RealD shift): _Mat(Mat), _shift(shift){};
-      // Support for coarsening to a multigrid
-      void OpDiag (const Field &in, Field &out) {
-	_Mat.Mdiag(in,out);
-	assert(0);
-      }
-      void OpDir  (const Field &in, Field &out,int dir,int disp) {
-	_Mat.Mdir(in,out,dir,disp);
-	assert(0);
-      }
-      void Op     (const Field &in, Field &out){
-	_Mat.M(in,out);
-	assert(0);
-      }
-      void AdjOp     (const Field &in, Field &out){
-	_Mat.Mdag(in,out);
-	assert(0);
-      }
-      void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-	_Mat.MdagM(in,out,n1,n2);
-	out = out + _shift*in;
+////////////////////////////////////////////////////////////////////
+// Construct herm op and shift it for mgrid smoother
+////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field>
+class ShiftedMdagMLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat;
+  RealD _shift;
+public:
+  ShiftedMdagMLinearOperator(Matrix &Mat,RealD shift): _Mat(Mat), _shift(shift){};
+  // Support for coarsening to a multigrid
+  void OpDiag (const Field &in, Field &out) {
+    _Mat.Mdiag(in,out);
+    assert(0);
+  }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) {
+    _Mat.Mdir(in,out,dir,disp);
+    assert(0);
+  }
+  void Op     (const Field &in, Field &out){
+    _Mat.M(in,out);
+    assert(0);
+  }
+  void AdjOp     (const Field &in, Field &out){
+    _Mat.Mdag(in,out);
+    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);
-	n1=real(dot);
-	n2=norm2(out);
-      }
-      void HermOp(const Field &in, Field &out){
-	RealD n1,n2;
-	HermOpAndNorm(in,out,n1,n2);
-      }
-    };
+    ComplexD dot;	
+    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);
+  }
+};
 
-    ////////////////////////////////////////////////////////////////////
-    // Wrap an already herm matrix
-    ////////////////////////////////////////////////////////////////////
-    template<class Matrix,class Field>
-    class HermitianLinearOperator : public LinearOperatorBase<Field> {
-      Matrix &_Mat;
-    public:
-    HermitianLinearOperator(Matrix &Mat): _Mat(Mat){};
-      // Support for coarsening to a multigrid
-      void OpDiag (const Field &in, Field &out) {
-	_Mat.Mdiag(in,out);
-      }
-      void OpDir  (const Field &in, Field &out,int dir,int disp) {
-	_Mat.Mdir(in,out,dir,disp);
-      }
-      void Op     (const Field &in, Field &out){
-	_Mat.M(in,out);
-      }
-      void AdjOp     (const Field &in, Field &out){
-	_Mat.M(in,out);
-      }
-      void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-	_Mat.M(in,out);
+////////////////////////////////////////////////////////////////////
+// Wrap an already herm matrix
+////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field>
+class HermitianLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat;
+public:
+  HermitianLinearOperator(Matrix &Mat): _Mat(Mat){};
+  // Support for coarsening to a multigrid
+  void OpDiag (const Field &in, Field &out) {
+    _Mat.Mdiag(in,out);
+  }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) {
+    _Mat.Mdir(in,out,dir,disp);
+  }
+  void Op     (const Field &in, Field &out){
+    _Mat.M(in,out);
+  }
+  void AdjOp     (const Field &in, Field &out){
+    _Mat.M(in,out);
+  }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+    _Mat.M(in,out);
 	
-	ComplexD dot= innerProduct(in,out); n1=real(dot);
-	n2=norm2(out);
-      }
-      void HermOp(const Field &in, Field &out){
-	_Mat.M(in,out);
-      }
-    };
+    ComplexD dot= innerProduct(in,out); n1=real(dot);
+    n2=norm2(out);
+  }
+  void HermOp(const Field &in, Field &out){
+    _Mat.M(in,out);
+  }
+};
 
-    //////////////////////////////////////////////////////////
-    // 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);
-	ni=Mpc(in,tmp);
-	no=MpcDag(tmp,out);
-      }
-      virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-	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);
-      }
-    };
-    template<class Matrix,class Field>
-      class SchurDiagMooeeOperator :  public SchurOperatorBase<Field> {
-    protected:
-      Matrix &_Mat;
-    public:
-      SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){};
-      virtual  RealD Mpc      (const Field &in, Field &out) {
-	Field tmp(in._grid);
-//	std::cout <<"grid pointers: in._grid="<< in._grid << " out._grid=" << out._grid << "  _Mat.Grid=" << _Mat.Grid() << " _Mat.RedBlackGrid=" << _Mat.RedBlackGrid() << std::endl;
+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);
+    ni=Mpc(in,tmp);
+    no=MpcDag(tmp,out);
+  }
+  virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+    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);
+  }
+};
+template<class Matrix,class Field>
+class SchurDiagMooeeOperator :  public SchurOperatorBase<Field> {
+protected:
+  Matrix &_Mat;
+public:
+  SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){};
+  virtual  RealD Mpc      (const Field &in, Field &out) {
+    Field tmp(in._grid);
+    //	std::cout <<"grid pointers: in._grid="<< in._grid << " out._grid=" << out._grid << "  _Mat.Grid=" << _Mat.Grid() << " _Mat.RedBlackGrid=" << _Mat.RedBlackGrid() << std::endl;
 
-	_Mat.Meooe(in,tmp);
-	_Mat.MooeeInv(tmp,out);
-	_Mat.Meooe(out,tmp);
+    _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.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.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){};
+    _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);
+  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);
+    _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);
+    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);
+    _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){};
+    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);
+  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);
+    _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);
+    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);
+    _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;
-    public:
-      SchurStaggeredOperator (Matrix &Mat): _Mat(Mat){};
-      virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-	GridLogIterative.TimingMode(1);
-	std::cout << GridLogIterative << " HermOpAndNorm "<<std::endl;
-	n2 = Mpc(in,out);
-	std::cout << GridLogIterative << " HermOpAndNorm.Mpc "<<std::endl;
-	ComplexD dot= innerProduct(in,out);
-	std::cout << GridLogIterative << " HermOpAndNorm.innerProduct "<<std::endl;
-	n1 = real(dot);
-      }
-      virtual void HermOp(const Field &in, Field &out){
-	std::cout << GridLogIterative << " HermOp "<<std::endl;
-	Mpc(in,out);
-      }
-      virtual  RealD Mpc      (const Field &in, Field &out) {
-	Field tmp(in._grid);
-	Field tmp2(in._grid);
+    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;
+public:
+  SchurStaggeredOperator (Matrix &Mat): _Mat(Mat){};
+  virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+    GridLogIterative.TimingMode(1);
+    std::cout << GridLogIterative << " HermOpAndNorm "<<std::endl;
+    n2 = Mpc(in,out);
+    std::cout << GridLogIterative << " HermOpAndNorm.Mpc "<<std::endl;
+    ComplexD dot= innerProduct(in,out);
+    std::cout << GridLogIterative << " HermOpAndNorm.innerProduct "<<std::endl;
+    n1 = real(dot);
+  }
+  virtual void HermOp(const Field &in, Field &out){
+    std::cout << GridLogIterative << " HermOp "<<std::endl;
+    Mpc(in,out);
+  }
+  virtual  RealD Mpc      (const Field &in, Field &out) {
+    Field tmp(in._grid);
+    Field tmp2(in._grid);
 
-	std::cout << GridLogIterative << " HermOp.Mpc "<<std::endl;
-	_Mat.Mooee(in,out);
-	_Mat.Mooee(out,tmp);
-	std::cout << GridLogIterative << " HermOp.MooeeMooee "<<std::endl;
+    std::cout << GridLogIterative << " HermOp.Mpc "<<std::endl;
+    _Mat.Mooee(in,out);
+    _Mat.Mooee(out,tmp);
+    std::cout << GridLogIterative << " HermOp.MooeeMooee "<<std::endl;
 
-	_Mat.Meooe(in,out);
-	_Mat.Meooe(out,tmp2);
-	std::cout << GridLogIterative << " HermOp.MeooeMeooe "<<std::endl;
+    _Mat.Meooe(in,out);
+    _Mat.Meooe(out,tmp2);
+    std::cout << GridLogIterative << " HermOp.MeooeMeooe "<<std::endl;
 
-	RealD nn=axpy_norm(out,-1.0,tmp2,tmp);
-	std::cout << GridLogIterative << " HermOp.axpy_norm "<<std::endl;
-	return nn;
-      }
-      virtual  RealD MpcDag   (const Field &in, Field &out){
-	return Mpc(in,out);
-      }
-      virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) {
-	assert(0);// Never need with staggered
-      }
-    };
-    template<class Matrix,class Field> using SchurStagOperator = SchurStaggeredOperator<Matrix,Field>;
+    RealD nn=axpy_norm(out,-1.0,tmp2,tmp);
+    std::cout << GridLogIterative << " HermOp.axpy_norm "<<std::endl;
+    return nn;
+  }
+  virtual  RealD MpcDag   (const Field &in, Field &out){
+    return Mpc(in,out);
+  }
+  virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) {
+    assert(0);// Never need with staggered
+  }
+};
+template<class Matrix,class Field> using SchurStagOperator = SchurStaggeredOperator<Matrix,Field>;
 
 
-    /////////////////////////////////////////////////////////////
-    // Base classes for functions of operators
-    /////////////////////////////////////////////////////////////
-    template<class Field> class OperatorFunction {
-    public:
-      virtual void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) = 0;
-    };
+/////////////////////////////////////////////////////////////
+// Base classes for functions of operators
+/////////////////////////////////////////////////////////////
+template<class Field> class OperatorFunction {
+public:
+  virtual void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) = 0;
+};
 
-    template<class Field> class LinearFunction {
-    public:
-      virtual void operator() (const Field &in, Field &out) = 0;
-    };
+template<class Field> class LinearFunction {
+public:
+  virtual void operator() (const Field &in, Field &out) = 0;
+};
 
-    template<class Field> class IdentityLinearFunction : public LinearFunction<Field> {
-    public:
-      void operator() (const Field &in, Field &out){
-	out = in;
-      };
-    };
-
-
-    /////////////////////////////////////////////////////////////
-    // Base classes for Multishift solvers for operators
-    /////////////////////////////////////////////////////////////
-    template<class Field> class OperatorMultiFunction {
-    public:
-      virtual void operator() (LinearOperatorBase<Field> &Linop, const Field &in, std::vector<Field> &out) = 0;
-    };
-
-    // FIXME : To think about
-
-    // Chroma functionality list defining LinearOperator
-    /*
-     virtual void operator() (T& chi, const T& psi, enum PlusMinus isign) const = 0;
-     virtual void operator() (T& chi, const T& psi, enum PlusMinus isign, Real epsilon) const
-     virtual const Subset& subset() const = 0;
-     virtual unsigned long nFlops() const { return 0; }
-     virtual void deriv(P& ds_u, const T& chi, const T& psi, enum PlusMinus isign) const
-     class UnprecLinearOperator : public DiffLinearOperator<T,P,Q>
-       const Subset& subset() const {return all;}
-     };
-    */
-
-  ////////////////////////////////////////////////////////////////////////////////////////////
-  // Hermitian operator Linear function and operator function
-  ////////////////////////////////////////////////////////////////////////////////////////////
-    template<class Field>
-      class HermOpOperatorFunction : public OperatorFunction<Field> {
-      void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
-	Linop.HermOp(in,out);
-      };
-    };
-
-    template<typename Field>
-      class PlainHermOp : public LinearFunction<Field> {
-    public:
-      LinearOperatorBase<Field> &_Linop;
-      
-      PlainHermOp(LinearOperatorBase<Field>& linop) : _Linop(linop) 
-      {}
-      
-      void operator()(const Field& in, Field& out) {
-	_Linop.HermOp(in,out);
-      }
-    };
-
-    template<typename Field>
-    class FunctionHermOp : public LinearFunction<Field> {
-    public:
-      OperatorFunction<Field>   & _poly;
-      LinearOperatorBase<Field> &_Linop;
-      
-      FunctionHermOp(OperatorFunction<Field> & poly,LinearOperatorBase<Field>& linop) 
-	: _poly(poly), _Linop(linop) {};
-      
-      void operator()(const Field& in, Field& out) {
-	_poly(_Linop,in,out);
-      }
-    };
-
-  template<class Field>
-  class Polynomial : public OperatorFunction<Field> {
-  private:
-    std::vector<RealD> Coeffs;
-  public:
-    Polynomial(std::vector<RealD> &_Coeffs) : Coeffs(_Coeffs) { };
-
-    // Implement the required interface
-    void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
-
-      Field AtoN(in._grid);
-      Field Mtmp(in._grid);
-      AtoN = in;
-      out = AtoN*Coeffs[0];
-      for(int n=1;n<Coeffs.size();n++){
-	Mtmp = AtoN;
-	Linop.HermOp(Mtmp,AtoN);
-	out=out+AtoN*Coeffs[n];
-      }
-    };
+template<class Field> class IdentityLinearFunction : public LinearFunction<Field> {
+public:
+  void operator() (const Field &in, Field &out){
+    out = in;
   };
+};
 
-}
+
+/////////////////////////////////////////////////////////////
+// Base classes for Multishift solvers for operators
+/////////////////////////////////////////////////////////////
+template<class Field> class OperatorMultiFunction {
+public:
+  virtual void operator() (LinearOperatorBase<Field> &Linop, const Field &in, std::vector<Field> &out) = 0;
+};
+
+// FIXME : To think about
+
+// Chroma functionality list defining LinearOperator
+/*
+  virtual void operator() (T& chi, const T& psi, enum PlusMinus isign) const = 0;
+  virtual void operator() (T& chi, const T& psi, enum PlusMinus isign, Real epsilon) const
+  virtual const Subset& subset() const = 0;
+  virtual unsigned long nFlops() const { return 0; }
+  virtual void deriv(P& ds_u, const T& chi, const T& psi, enum PlusMinus isign) const
+  class UnprecLinearOperator : public DiffLinearOperator<T,P,Q>
+  const Subset& subset() const {return all;}
+  };
+*/
+
+////////////////////////////////////////////////////////////////////////////////////////////
+// Hermitian operator Linear function and operator function
+////////////////////////////////////////////////////////////////////////////////////////////
+template<class Field>
+class HermOpOperatorFunction : public OperatorFunction<Field> {
+  void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
+    Linop.HermOp(in,out);
+  };
+};
+
+template<typename Field>
+class PlainHermOp : public LinearFunction<Field> {
+public:
+  LinearOperatorBase<Field> &_Linop;
+      
+  PlainHermOp(LinearOperatorBase<Field>& linop) : _Linop(linop) 
+  {}
+      
+  void operator()(const Field& in, Field& out) {
+    _Linop.HermOp(in,out);
+  }
+};
+
+template<typename Field>
+class FunctionHermOp : public LinearFunction<Field> {
+public:
+  OperatorFunction<Field>   & _poly;
+  LinearOperatorBase<Field> &_Linop;
+      
+  FunctionHermOp(OperatorFunction<Field> & poly,LinearOperatorBase<Field>& linop) 
+    : _poly(poly), _Linop(linop) {};
+      
+  void operator()(const Field& in, Field& out) {
+    _poly(_Linop,in,out);
+  }
+};
+
+template<class Field>
+class Polynomial : public OperatorFunction<Field> {
+private:
+  std::vector<RealD> Coeffs;
+public:
+  Polynomial(std::vector<RealD> &_Coeffs) : Coeffs(_Coeffs) { };
+
+  // Implement the required interface
+  void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
+
+    Field AtoN(in._grid);
+    Field Mtmp(in._grid);
+    AtoN = in;
+    out = AtoN*Coeffs[0];
+    for(int n=1;n<Coeffs.size();n++){
+      Mtmp = AtoN;
+      Linop.HermOp(Mtmp,AtoN);
+      out=out+AtoN*Coeffs[n];
+    }
+  };
+};
+
+NAMESPACE_END(Grid);
 
 #endif
diff --git a/lib/algorithms/Preconditioner.h b/lib/algorithms/Preconditioner.h
index a26016be..65ac3cf9 100644
--- a/lib/algorithms/Preconditioner.h
+++ b/lib/algorithms/Preconditioner.h
@@ -1,4 +1,4 @@
-    /*************************************************************************************
+/*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -23,24 +23,24 @@ Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
     See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
+*************************************************************************************/
+/*  END LEGAL */
 #ifndef GRID_PRECONDITIONER_H
 #define GRID_PRECONDITIONER_H
 
-namespace Grid {
+NAMESPACE_BEGIN(Grid);
 
-  template<class Field> class Preconditioner :  public LinearFunction<Field> { 
-    virtual void operator()(const Field &src, Field & psi)=0;
-  };
+template<class Field> class Preconditioner :  public LinearFunction<Field> { 
+  virtual void operator()(const Field &src, Field & psi)=0;
+};
 
-  template<class Field> class TrivialPrecon :  public Preconditioner<Field> { 
-  public:
-    void operator()(const Field &src, Field & psi){
-      psi = src;
-    }
-    TrivialPrecon(void){};
-  };
+template<class Field> class TrivialPrecon :  public Preconditioner<Field> { 
+public:
+  void operator()(const Field &src, Field & psi){
+    psi = src;
+  }
+  TrivialPrecon(void){};
+};
 
-}
+NAMESPACE_END(Grid);
 #endif
diff --git a/lib/algorithms/SparseMatrix.h b/lib/algorithms/SparseMatrix.h
index 1611a6f4..88e9e7bd 100644
--- a/lib/algorithms/SparseMatrix.h
+++ b/lib/algorithms/SparseMatrix.h
@@ -1,4 +1,4 @@
-    /*************************************************************************************
+/*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -23,49 +23,49 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
     See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
+*************************************************************************************/
+/*  END LEGAL */
 #ifndef  GRID_ALGORITHM_SPARSE_MATRIX_H
 #define  GRID_ALGORITHM_SPARSE_MATRIX_H
 
 
-namespace Grid {
+NAMESPACE_BEGIN(Grid);
 
-  /////////////////////////////////////////////////////////////////////////////////////////////
-  // Interface defining what I expect of a general sparse matrix, such as a Fermion action
-  /////////////////////////////////////////////////////////////////////////////////////////////
-    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 Mdiag    (const Field &in, Field &out)=0;
-      virtual  void Mdir     (const Field &in, Field &out,int dir, int disp)=0;
-    };
+/////////////////////////////////////////////////////////////////////////////////////////////
+// Interface defining what I expect of a general sparse matrix, such as a Fermion action
+/////////////////////////////////////////////////////////////////////////////////////////////
+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 Mdiag    (const Field &in, Field &out)=0;
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp)=0;
+};
 
-  /////////////////////////////////////////////////////////////////////////////////////////////
-  // Interface augmented by a red black sparse matrix, such as a Fermion action
-  /////////////////////////////////////////////////////////////////////////////////////////////
-    template<class Field> class CheckerBoardedSparseMatrixBase : public SparseMatrixBase<Field> {
-    public:
-      virtual GridBase *RedBlackGrid(void)=0;
-      // half checkerboard operaions
-      virtual  void Meooe    (const Field &in, Field &out)=0;
-      virtual  void Mooee    (const Field &in, Field &out)=0;
-      virtual  void MooeeInv (const Field &in, Field &out)=0;
+/////////////////////////////////////////////////////////////////////////////////////////////
+// Interface augmented by a red black sparse matrix, such as a Fermion action
+/////////////////////////////////////////////////////////////////////////////////////////////
+template<class Field> class CheckerBoardedSparseMatrixBase : public SparseMatrixBase<Field> {
+public:
+  virtual GridBase *RedBlackGrid(void)=0;
+  // half checkerboard operaions
+  virtual  void Meooe    (const Field &in, Field &out)=0;
+  virtual  void Mooee    (const Field &in, Field &out)=0;
+  virtual  void MooeeInv (const Field &in, Field &out)=0;
 
-      virtual  void MeooeDag    (const Field &in, Field &out)=0;
-      virtual  void MooeeDag    (const Field &in, Field &out)=0;
-      virtual  void MooeeInvDag (const Field &in, Field &out)=0;
+  virtual  void MeooeDag    (const Field &in, Field &out)=0;
+  virtual  void MooeeDag    (const Field &in, Field &out)=0;
+  virtual  void MooeeInvDag (const Field &in, Field &out)=0;
 
-    };
+};
 
-}
+NAMESPACE_END(Grid);
 
 #endif