diff --git a/Grid/algorithms/multigrid/Aggregates.h b/Grid/algorithms/multigrid/Aggregates.h
new file mode 100644
index 00000000..9112b868
--- /dev/null
+++ b/Grid/algorithms/multigrid/Aggregates.h
@@ -0,0 +1,286 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/Aggregates.h
+
+    Copyright (C) 2015
+
+Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+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;
+
+  GridBase *CoarseGrid;
+  GridBase *FineGrid;
+  std::vector<Lattice<Fobj> > subspace;
+  int checkerboard;
+  int Checkerboard(void){return checkerboard;}
+  Aggregation(GridBase *_CoarseGrid,GridBase *_FineGrid,int _checkerboard) : 
+    CoarseGrid(_CoarseGrid),
+    FineGrid(_FineGrid),
+    subspace(nbasis,_FineGrid),
+    checkerboard(_checkerboard)
+  {
+  };
+  
+  void SaveState(std::string file)
+  {
+#ifdef HAVE_LIME
+    emptyUserRecord record;
+    ScidacWriter WR(Grid()->IsBoss());
+    WR.open(file);
+    for(int b=0;b<nbasis;b++){
+      WR.writeScidacFieldRecord(subspace[b],record);
+    }
+    WR.close();
+#endif
+  }
+  void LoadState(std::string file)
+  {
+#ifdef HAVE_LIME
+    emptyUserRecord record;
+    Grid::ScidacReader RD ;
+    RD.open(file);
+    for(int b=0;b<nbasis;b++){
+      RD.readScidacFieldRecord(subspace[b],record);
+    }    
+    RD.close();
+#endif
+  }
+  
+  void Orthogonalise(void){
+    CoarseScalar InnerProd(CoarseGrid); 
+    //    std::cout << GridLogMessage <<" Block Gramm-Schmidt pass 1"<<std::endl;
+    blockOrthogonalise(InnerProd,subspace);
+  } 
+  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);
+  }
+
+  virtual void CreateSubspaceRandom(GridParallelRNG  &RNG) {
+    int nn=nbasis;
+    RealD scale;
+    FineField noise(FineGrid);
+    for(int b=0;b<nn;b++){
+      subspace[b] = Zero();
+      gaussian(RNG,noise);
+      scale = std::pow(norm2(noise),-0.5); 
+      noise=noise*scale;
+      subspace[b] = noise;
+    }
+  }
+  virtual void CreateSubspace(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis)
+  {
+
+    RealD scale;
+
+    ConjugateGradient<FineField> CG(1.0e-2,100,false);
+    FineField noise(FineGrid);
+    FineField Mn(FineGrid);
+
+    for(int b=0;b<nn;b++){
+      
+      subspace[b] = Zero();
+      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 << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(Mn)<<std::endl;
+      subspace[b]   = noise;
+
+    }
+  }
+
+  ////////////////////////////////////////////////////////////////////////////////////////////////
+  // World of possibilities here. But have tried quite a lot of experiments (250+ jobs run on Summit)
+  // and this is the best I found
+  ////////////////////////////////////////////////////////////////////////////////////////////////
+
+  virtual void CreateSubspaceChebyshev(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,
+				       int nn,
+				       double hi,
+				       double lo,
+				       int orderfilter,
+				       int ordermin,
+				       int orderstep,
+				       double filterlo
+				       ) {
+
+    RealD scale;
+
+    FineField noise(FineGrid);
+    FineField Mn(FineGrid);
+    FineField tmp(FineGrid);
+
+    // New normalised noise
+    gaussian(RNG,noise);
+    scale = std::pow(norm2(noise),-0.5); 
+    noise=noise*scale;
+
+    std::cout << GridLogMessage<<" Chebyshev subspace pass-1 : ord "<<orderfilter<<" ["<<lo<<","<<hi<<"]"<<std::endl;
+    std::cout << GridLogMessage<<" Chebyshev subspace pass-2 : nbasis"<<nn<<" min "
+	      <<ordermin<<" step "<<orderstep
+	      <<" lo"<<filterlo<<std::endl;
+
+    // Initial matrix element
+    hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
+
+    int b =0;
+    {
+      // Filter
+      Chebyshev<FineField> Cheb(lo,hi,orderfilter);
+      Cheb(hermop,noise,Mn);
+      // normalise
+      scale = std::pow(norm2(Mn),-0.5); 	Mn=Mn*scale;
+      subspace[b]   = Mn;
+      hermop.Op(Mn,tmp); 
+      std::cout<<GridLogMessage << "filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
+      b++;
+    }
+
+    // Generate a full sequence of Chebyshevs
+    {
+      lo=filterlo;
+      noise=Mn;
+
+      FineField T0(FineGrid); T0 = noise;  
+      FineField T1(FineGrid); 
+      FineField T2(FineGrid);
+      FineField y(FineGrid);
+      
+      FineField *Tnm = &T0;
+      FineField *Tn  = &T1;
+      FineField *Tnp = &T2;
+
+      // Tn=T1 = (xscale M + mscale)in
+      RealD xscale = 2.0/(hi-lo);
+      RealD mscale = -(hi+lo)/(hi-lo);
+      hermop.HermOp(T0,y);
+      T1=y*xscale+noise*mscale;
+
+      for(int n=2;n<=ordermin+orderstep*(nn-2);n++){
+	
+	hermop.HermOp(*Tn,y);
+
+	autoView( y_v , y, AcceleratorWrite);
+	autoView( Tn_v , (*Tn), AcceleratorWrite);
+	autoView( Tnp_v , (*Tnp), AcceleratorWrite);
+	autoView( Tnm_v , (*Tnm), AcceleratorWrite);
+	const int Nsimd = CComplex::Nsimd();
+	accelerator_for(ss, FineGrid->oSites(), Nsimd, {
+	  coalescedWrite(y_v[ss],xscale*y_v(ss)+mscale*Tn_v(ss));
+	  coalescedWrite(Tnp_v[ss],2.0*y_v(ss)-Tnm_v(ss));
+        });
+
+	// Possible more fine grained control is needed than a linear sweep,
+	// but huge productivity gain if this is simple algorithm and not a tunable
+	int m =1;
+	if ( n>=ordermin ) m=n-ordermin;
+	if ( (m%orderstep)==0 ) { 
+	  Mn=*Tnp;
+	  scale = std::pow(norm2(Mn),-0.5);         Mn=Mn*scale;
+	  subspace[b] = Mn;
+	  hermop.Op(Mn,tmp); 
+	  std::cout<<GridLogMessage << n<<" filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
+	  b++;
+	}
+
+	// Cycle pointers to avoid copies
+	FineField *swizzle = Tnm;
+	Tnm    =Tn;
+	Tn     =Tnp;
+	Tnp    =swizzle;
+	  
+      }
+    }
+    assert(b==nn);
+  }
+  virtual void CreateSubspaceChebyshev(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,
+				       int nn,
+				       double hi,
+				       double lo,
+				       int orderfilter
+				       ) {
+
+    RealD scale;
+
+    FineField noise(FineGrid);
+    FineField Mn(FineGrid);
+    FineField tmp(FineGrid);
+
+    // New normalised noise
+    std::cout << GridLogMessage<<" Chebyshev subspace pure noise : ord "<<orderfilter<<" ["<<lo<<","<<hi<<"]"<<std::endl;
+    std::cout << GridLogMessage<<" Chebyshev subspace pure noise  : nbasis "<<nn<<std::endl;
+
+
+    for(int b =0;b<nbasis;b++)
+    {
+      gaussian(RNG,noise);
+      scale = std::pow(norm2(noise),-0.5); 
+      noise=noise*scale;
+
+      // Initial matrix element
+      hermop.Op(noise,Mn);
+      if(b==0) std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
+      // Filter
+      Chebyshev<FineField> Cheb(lo,hi,orderfilter);
+      Cheb(hermop,noise,Mn);
+      // normalise
+      scale = std::pow(norm2(Mn),-0.5); 	Mn=Mn*scale;
+      subspace[b]   = Mn;
+      hermop.Op(Mn,tmp); 
+      std::cout<<GridLogMessage << "filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
+    }
+
+  }
+
+};
+NAMESPACE_END(Grid);
diff --git a/Grid/algorithms/multigrid/CoarsenedMatrix.h b/Grid/algorithms/multigrid/CoarsenedMatrix.h
new file mode 100644
index 00000000..42634004
--- /dev/null
+++ b/Grid/algorithms/multigrid/CoarsenedMatrix.h
@@ -0,0 +1,814 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/CoarsenedMatrix.h
+
+    Copyright (C) 2015
+
+Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef  GRID_ALGORITHM_COARSENED_MATRIX_H
+#define  GRID_ALGORITHM_COARSENED_MATRIX_H
+
+#include <Grid/qcd/QCD.h> // needed for Dagger(Yes|No), Inverse(Yes|No)
+
+NAMESPACE_BEGIN(Grid);
+
+template<class vobj,class CComplex>
+inline void blockMaskedInnerProduct(Lattice<CComplex> &CoarseInner,
+				    const Lattice<decltype(innerProduct(vobj(),vobj()))> &FineMask,
+				    const Lattice<vobj> &fineX,
+				    const Lattice<vobj> &fineY)
+{
+  typedef decltype(innerProduct(vobj(),vobj())) dotp;
+
+  GridBase *coarse(CoarseInner.Grid());
+  GridBase *fine  (fineX.Grid());
+
+  Lattice<dotp> fine_inner(fine); fine_inner.Checkerboard() = fineX.Checkerboard();
+  Lattice<dotp> fine_inner_msk(fine);
+
+  // Multiply could be fused with innerProduct
+  // Single block sum kernel could do both masks.
+  fine_inner = localInnerProduct(fineX,fineY);
+  mult(fine_inner_msk, fine_inner,FineMask);
+  blockSum(CoarseInner,fine_inner_msk);
+}
+
+// Fine Object == (per site) type of fine field
+// nbasis      == number of deflation vectors
+template<class Fobj,class CComplex,int nbasis>
+class CoarsenedMatrix : public CheckerBoardedSparseMatrixBase<Lattice<iVector<CComplex,nbasis > > >  {
+public:
+    
+  typedef iVector<CComplex,nbasis >           siteVector;
+  typedef Lattice<CComplex >                  CoarseComplexField;
+  typedef Lattice<siteVector>                 CoarseVector;
+  typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
+  typedef iMatrix<CComplex,nbasis >  Cobj;
+  typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
+  typedef Lattice<Fobj >        FineField;
+  typedef CoarseVector FermionField;
+
+  // enrich interface, use default implementation as in FermionOperator ///////
+  void Dminus(CoarseVector const& in, CoarseVector& out) { out = in; }
+  void DminusDag(CoarseVector const& in, CoarseVector& out) { out = in; }
+  void ImportPhysicalFermionSource(CoarseVector const& input, CoarseVector& imported) { imported = input; }
+  void ImportUnphysicalFermion(CoarseVector const& input, CoarseVector& imported) { imported = input; }
+  void ExportPhysicalFermionSolution(CoarseVector const& solution, CoarseVector& exported) { exported = solution; };
+  void ExportPhysicalFermionSource(CoarseVector const& solution, CoarseVector& exported) { exported = solution; };
+
+  ////////////////////
+  // Data members
+  ////////////////////
+  Geometry         geom;
+  GridBase *       _grid; 
+  GridBase*        _cbgrid;
+  int hermitian;
+
+  CartesianStencil<siteVector,siteVector,DefaultImplParams> Stencil; 
+  CartesianStencil<siteVector,siteVector,DefaultImplParams> StencilEven;
+  CartesianStencil<siteVector,siteVector,DefaultImplParams> StencilOdd;
+
+  std::vector<CoarseMatrix> A;
+  std::vector<CoarseMatrix> Aeven;
+  std::vector<CoarseMatrix> Aodd;
+
+  CoarseMatrix AselfInv;
+  CoarseMatrix AselfInvEven;
+  CoarseMatrix AselfInvOdd;
+
+  Vector<RealD> dag_factor;
+
+  ///////////////////////
+  // Interface
+  ///////////////////////
+  GridBase * Grid(void)         { return _grid; };   // this is all the linalg routines need to know
+  GridBase * RedBlackGrid()     { return _cbgrid; };
+
+  int ConstEE() { return 0; }
+
+  void M (const CoarseVector &in, CoarseVector &out)
+  {
+    conformable(_grid,in.Grid());
+    conformable(in.Grid(),out.Grid());
+    out.Checkerboard() = in.Checkerboard();
+
+    SimpleCompressor<siteVector> compressor;
+
+    Stencil.HaloExchange(in,compressor);
+    autoView( in_v , in, AcceleratorRead);
+    autoView( out_v , out, AcceleratorWrite);
+    autoView( Stencil_v  , Stencil, AcceleratorRead);
+    int npoint = geom.npoint;
+    typedef LatticeView<Cobj> Aview;
+      
+    Vector<Aview> AcceleratorViewContainer;
+  
+    for(int p=0;p<geom.npoint;p++) AcceleratorViewContainer.push_back(A[p].View(AcceleratorRead));
+    Aview *Aview_p = & AcceleratorViewContainer[0];
+
+    const int Nsimd = CComplex::Nsimd();
+    typedef decltype(coalescedRead(in_v[0])) calcVector;
+    typedef decltype(coalescedRead(in_v[0](0))) calcComplex;
+
+    int osites=Grid()->oSites();
+
+    accelerator_for(sss, Grid()->oSites()*nbasis, Nsimd, {
+      int ss = sss/nbasis;
+      int b  = sss%nbasis;
+      calcComplex res = Zero();
+      calcVector nbr;
+      int ptype;
+      StencilEntry *SE;
+
+      for(int point=0;point<npoint;point++){
+
+	SE=Stencil_v.GetEntry(ptype,point,ss);
+	  
+	if(SE->_is_local) { 
+	  nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
+	} else {
+	  nbr = coalescedRead(Stencil_v.CommBuf()[SE->_offset]);
+	}
+	acceleratorSynchronise();
+
+	for(int bb=0;bb<nbasis;bb++) {
+	  res = res + coalescedRead(Aview_p[point][ss](b,bb))*nbr(bb);
+	}
+      }
+      coalescedWrite(out_v[ss](b),res);
+      });
+
+    for(int p=0;p<geom.npoint;p++) AcceleratorViewContainer[p].ViewClose();
+  };
+
+  void Mdag (const CoarseVector &in, CoarseVector &out)
+  {
+    if(hermitian) {
+      // corresponds to Petrov-Galerkin coarsening
+      return M(in,out);
+    } else {
+      // corresponds to Galerkin coarsening
+      return MdagNonHermitian(in, out);
+    }
+  };
+
+  void MdagNonHermitian(const CoarseVector &in, CoarseVector &out)
+  {
+    conformable(_grid,in.Grid());
+    conformable(in.Grid(),out.Grid());
+    out.Checkerboard() = in.Checkerboard();
+
+    SimpleCompressor<siteVector> compressor;
+
+    Stencil.HaloExchange(in,compressor);
+    autoView( in_v , in, AcceleratorRead);
+    autoView( out_v , out, AcceleratorWrite);
+    autoView( Stencil_v  , Stencil, AcceleratorRead);
+    int npoint = geom.npoint;
+    typedef LatticeView<Cobj> Aview;
+
+    Vector<Aview> AcceleratorViewContainer;
+
+    for(int p=0;p<geom.npoint;p++) AcceleratorViewContainer.push_back(A[p].View(AcceleratorRead));
+    Aview *Aview_p = & AcceleratorViewContainer[0];
+
+    const int Nsimd = CComplex::Nsimd();
+    typedef decltype(coalescedRead(in_v[0])) calcVector;
+    typedef decltype(coalescedRead(in_v[0](0))) calcComplex;
+
+    int osites=Grid()->oSites();
+
+    Vector<int> points(geom.npoint, 0);
+    for(int p=0; p<geom.npoint; p++)
+      points[p] = geom.points_dagger[p];
+
+    auto points_p = &points[0];
+
+    RealD* dag_factor_p = &dag_factor[0];
+
+    accelerator_for(sss, Grid()->oSites()*nbasis, Nsimd, {
+      int ss = sss/nbasis;
+      int b  = sss%nbasis;
+      calcComplex res = Zero();
+      calcVector nbr;
+      int ptype;
+      StencilEntry *SE;
+
+      for(int p=0;p<npoint;p++){
+        int point = points_p[p];
+
+	SE=Stencil_v.GetEntry(ptype,point,ss);
+
+	if(SE->_is_local) {
+	  nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
+	} else {
+	  nbr = coalescedRead(Stencil_v.CommBuf()[SE->_offset]);
+	}
+	acceleratorSynchronise();
+
+	for(int bb=0;bb<nbasis;bb++) {
+	  res = res + dag_factor_p[b*nbasis+bb]*coalescedRead(Aview_p[point][ss](b,bb))*nbr(bb);
+	}
+      }
+      coalescedWrite(out_v[ss](b),res);
+      });
+
+    for(int p=0;p<geom.npoint;p++) AcceleratorViewContainer[p].ViewClose();
+  }
+
+  void MdirComms(const CoarseVector &in)
+  {
+    SimpleCompressor<siteVector> compressor;
+    Stencil.HaloExchange(in,compressor);
+  }
+  void MdirCalc(const CoarseVector &in, CoarseVector &out, int point)
+  {
+    conformable(_grid,in.Grid());
+    conformable(_grid,out.Grid());
+    out.Checkerboard() = in.Checkerboard();
+
+    typedef LatticeView<Cobj> Aview;
+    Vector<Aview> AcceleratorViewContainer;
+    for(int p=0;p<geom.npoint;p++) AcceleratorViewContainer.push_back(A[p].View(AcceleratorRead));
+    Aview *Aview_p = & AcceleratorViewContainer[0];
+
+    autoView( out_v , out, AcceleratorWrite);
+    autoView( in_v  , in, AcceleratorRead);
+    autoView( Stencil_v  , Stencil, AcceleratorRead);
+
+    const int Nsimd = CComplex::Nsimd();
+    typedef decltype(coalescedRead(in_v[0])) calcVector;
+    typedef decltype(coalescedRead(in_v[0](0))) calcComplex;
+
+    accelerator_for(sss, Grid()->oSites()*nbasis, Nsimd, {
+      int ss = sss/nbasis;
+      int b  = sss%nbasis;
+      calcComplex res = Zero();
+      calcVector nbr;
+      int ptype;
+      StencilEntry *SE;
+
+      SE=Stencil_v.GetEntry(ptype,point,ss);
+	  
+      if(SE->_is_local) { 
+	nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
+      } else {
+	nbr = coalescedRead(Stencil_v.CommBuf()[SE->_offset]);
+      }
+      acceleratorSynchronise();
+
+      for(int bb=0;bb<nbasis;bb++) {
+	res = res + coalescedRead(Aview_p[point][ss](b,bb))*nbr(bb);
+      }
+      coalescedWrite(out_v[ss](b),res);
+    });
+    for(int p=0;p<geom.npoint;p++) AcceleratorViewContainer[p].ViewClose();
+  }
+  void MdirAll(const CoarseVector &in,std::vector<CoarseVector> &out)
+  {
+    this->MdirComms(in);
+    int ndir=geom.npoint-1;
+    if ((out.size()!=ndir)&&(out.size()!=ndir+1)) { 
+      std::cout <<"MdirAll out size "<< out.size()<<std::endl;
+      std::cout <<"MdirAll ndir "<< ndir<<std::endl;
+      assert(0);
+    }
+    for(int p=0;p<ndir;p++){
+      MdirCalc(in,out[p],p);
+    }
+  };
+  void Mdir(const CoarseVector &in, CoarseVector &out, int dir, int disp){
+
+    this->MdirComms(in);
+
+    MdirCalc(in,out,geom.point(dir,disp));
+  };
+
+  void Mdiag(const CoarseVector &in, CoarseVector &out)
+  {
+    int point=geom.npoint-1;
+    MdirCalc(in, out, point); // No comms
+  };
+
+  void Mooee(const CoarseVector &in, CoarseVector &out) {
+    MooeeInternal(in, out, DaggerNo, InverseNo);
+  }
+
+  void MooeeInv(const CoarseVector &in, CoarseVector &out) {
+    MooeeInternal(in, out, DaggerNo, InverseYes);
+  }
+
+  void MooeeDag(const CoarseVector &in, CoarseVector &out) {
+    MooeeInternal(in, out, DaggerYes, InverseNo);
+  }
+
+  void MooeeInvDag(const CoarseVector &in, CoarseVector &out) {
+    MooeeInternal(in, out, DaggerYes, InverseYes);
+  }
+
+  void Meooe(const CoarseVector &in, CoarseVector &out) {
+    if(in.Checkerboard() == Odd) {
+      DhopEO(in, out, DaggerNo);
+    } else {
+      DhopOE(in, out, DaggerNo);
+    }
+  }
+
+  void MeooeDag(const CoarseVector &in, CoarseVector &out) {
+    if(in.Checkerboard() == Odd) {
+      DhopEO(in, out, DaggerYes);
+    } else {
+      DhopOE(in, out, DaggerYes);
+    }
+  }
+
+  void Dhop(const CoarseVector &in, CoarseVector &out, int dag) {
+    conformable(in.Grid(), _grid); // verifies full grid
+    conformable(in.Grid(), out.Grid());
+
+    out.Checkerboard() = in.Checkerboard();
+
+    DhopInternal(Stencil, A, in, out, dag);
+  }
+
+  void DhopOE(const CoarseVector &in, CoarseVector &out, int dag) {
+    conformable(in.Grid(), _cbgrid);    // verifies half grid
+    conformable(in.Grid(), out.Grid()); // drops the cb check
+
+    assert(in.Checkerboard() == Even);
+    out.Checkerboard() = Odd;
+
+    DhopInternal(StencilEven, Aodd, in, out, dag);
+  }
+
+  void DhopEO(const CoarseVector &in, CoarseVector &out, int dag) {
+    conformable(in.Grid(), _cbgrid);    // verifies half grid
+    conformable(in.Grid(), out.Grid()); // drops the cb check
+
+    assert(in.Checkerboard() == Odd);
+    out.Checkerboard() = Even;
+
+    DhopInternal(StencilOdd, Aeven, in, out, dag);
+  }
+
+  void MooeeInternal(const CoarseVector &in, CoarseVector &out, int dag, int inv) {
+    out.Checkerboard() = in.Checkerboard();
+    assert(in.Checkerboard() == Odd || in.Checkerboard() == Even);
+
+    CoarseMatrix *Aself = nullptr;
+    if(in.Grid()->_isCheckerBoarded) {
+      if(in.Checkerboard() == Odd) {
+        Aself = (inv) ? &AselfInvOdd : &Aodd[geom.npoint-1];
+        DselfInternal(StencilOdd, *Aself, in, out, dag);
+      } else {
+        Aself = (inv) ? &AselfInvEven : &Aeven[geom.npoint-1];
+        DselfInternal(StencilEven, *Aself, in, out, dag);
+      }
+    } else {
+      Aself = (inv) ? &AselfInv : &A[geom.npoint-1];
+      DselfInternal(Stencil, *Aself, in, out, dag);
+    }
+    assert(Aself != nullptr);
+  }
+
+  void DselfInternal(CartesianStencil<siteVector,siteVector,DefaultImplParams> &st, CoarseMatrix &a,
+                       const CoarseVector &in, CoarseVector &out, int dag) {
+    int point = geom.npoint-1;
+    autoView( out_v, out, AcceleratorWrite);
+    autoView( in_v,  in,  AcceleratorRead);
+    autoView( st_v,  st,  AcceleratorRead);
+    autoView( a_v,   a,   AcceleratorRead);
+
+    const int Nsimd = CComplex::Nsimd();
+    typedef decltype(coalescedRead(in_v[0])) calcVector;
+    typedef decltype(coalescedRead(in_v[0](0))) calcComplex;
+
+    RealD* dag_factor_p = &dag_factor[0];
+
+    if(dag) {
+      accelerator_for(sss, in.Grid()->oSites()*nbasis, Nsimd, {
+        int ss = sss/nbasis;
+        int b  = sss%nbasis;
+        calcComplex res = Zero();
+        calcVector nbr;
+        int ptype;
+        StencilEntry *SE;
+
+        SE=st_v.GetEntry(ptype,point,ss);
+
+        if(SE->_is_local) {
+          nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
+        } else {
+          nbr = coalescedRead(st_v.CommBuf()[SE->_offset]);
+        }
+        acceleratorSynchronise();
+
+        for(int bb=0;bb<nbasis;bb++) {
+          res = res + dag_factor_p[b*nbasis+bb]*coalescedRead(a_v[ss](b,bb))*nbr(bb);
+        }
+        coalescedWrite(out_v[ss](b),res);
+      });
+    } else {
+      accelerator_for(sss, in.Grid()->oSites()*nbasis, Nsimd, {
+        int ss = sss/nbasis;
+        int b  = sss%nbasis;
+        calcComplex res = Zero();
+        calcVector nbr;
+        int ptype;
+        StencilEntry *SE;
+
+        SE=st_v.GetEntry(ptype,point,ss);
+
+        if(SE->_is_local) {
+          nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
+        } else {
+          nbr = coalescedRead(st_v.CommBuf()[SE->_offset]);
+        }
+        acceleratorSynchronise();
+
+        for(int bb=0;bb<nbasis;bb++) {
+          res = res + coalescedRead(a_v[ss](b,bb))*nbr(bb);
+        }
+        coalescedWrite(out_v[ss](b),res);
+      });
+    }
+  }
+
+  void DhopInternal(CartesianStencil<siteVector,siteVector,DefaultImplParams> &st, std::vector<CoarseMatrix> &a,
+                    const CoarseVector &in, CoarseVector &out, int dag) {
+    SimpleCompressor<siteVector> compressor;
+
+    st.HaloExchange(in,compressor);
+    autoView( in_v,  in,  AcceleratorRead);
+    autoView( out_v, out, AcceleratorWrite);
+    autoView( st_v , st,  AcceleratorRead);
+    typedef LatticeView<Cobj> Aview;
+
+    // determine in what order we need the points
+    int npoint = geom.npoint-1;
+    Vector<int> points(npoint, 0);
+    for(int p=0; p<npoint; p++)
+      points[p] = (dag && !hermitian) ? geom.points_dagger[p] : p;
+
+    auto points_p = &points[0];
+
+    Vector<Aview> AcceleratorViewContainer;
+    for(int p=0;p<npoint;p++) AcceleratorViewContainer.push_back(a[p].View(AcceleratorRead));
+    Aview *Aview_p = & AcceleratorViewContainer[0];
+
+    const int Nsimd = CComplex::Nsimd();
+    typedef decltype(coalescedRead(in_v[0])) calcVector;
+    typedef decltype(coalescedRead(in_v[0](0))) calcComplex;
+
+    RealD* dag_factor_p = &dag_factor[0];
+
+    if(dag) {
+      accelerator_for(sss, in.Grid()->oSites()*nbasis, Nsimd, {
+        int ss = sss/nbasis;
+        int b  = sss%nbasis;
+        calcComplex res = Zero();
+        calcVector nbr;
+        int ptype;
+        StencilEntry *SE;
+
+        for(int p=0;p<npoint;p++){
+          int point = points_p[p];
+          SE=st_v.GetEntry(ptype,point,ss);
+
+          if(SE->_is_local) {
+            nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
+          } else {
+            nbr = coalescedRead(st_v.CommBuf()[SE->_offset]);
+          }
+          acceleratorSynchronise();
+
+          for(int bb=0;bb<nbasis;bb++) {
+            res = res + dag_factor_p[b*nbasis+bb]*coalescedRead(Aview_p[point][ss](b,bb))*nbr(bb);
+          }
+        }
+        coalescedWrite(out_v[ss](b),res);
+      });
+    } else {
+      accelerator_for(sss, in.Grid()->oSites()*nbasis, Nsimd, {
+        int ss = sss/nbasis;
+        int b  = sss%nbasis;
+        calcComplex res = Zero();
+        calcVector nbr;
+        int ptype;
+        StencilEntry *SE;
+
+        for(int p=0;p<npoint;p++){
+          int point = points_p[p];
+          SE=st_v.GetEntry(ptype,point,ss);
+
+          if(SE->_is_local) {
+            nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
+          } else {
+            nbr = coalescedRead(st_v.CommBuf()[SE->_offset]);
+          }
+          acceleratorSynchronise();
+
+          for(int bb=0;bb<nbasis;bb++) {
+            res = res + coalescedRead(Aview_p[point][ss](b,bb))*nbr(bb);
+          }
+        }
+        coalescedWrite(out_v[ss](b),res);
+      });
+    }
+
+    for(int p=0;p<npoint;p++) AcceleratorViewContainer[p].ViewClose();
+  }
+  
+  CoarsenedMatrix(GridCartesian &CoarseGrid, int hermitian_=0) 	:
+    _grid(&CoarseGrid),
+    _cbgrid(new GridRedBlackCartesian(&CoarseGrid)),
+    geom(CoarseGrid._ndimension),
+    hermitian(hermitian_),
+    Stencil(&CoarseGrid,geom.npoint,Even,geom.directions,geom.displacements),
+    StencilEven(_cbgrid,geom.npoint,Even,geom.directions,geom.displacements),
+    StencilOdd(_cbgrid,geom.npoint,Odd,geom.directions,geom.displacements),
+    A(geom.npoint,&CoarseGrid),
+    Aeven(geom.npoint,_cbgrid),
+    Aodd(geom.npoint,_cbgrid),
+    AselfInv(&CoarseGrid),
+    AselfInvEven(_cbgrid),
+    AselfInvOdd(_cbgrid),
+    dag_factor(nbasis*nbasis)
+  {
+    fillFactor();
+  };
+
+  CoarsenedMatrix(GridCartesian &CoarseGrid, GridRedBlackCartesian &CoarseRBGrid, int hermitian_=0) 	:
+
+    _grid(&CoarseGrid),
+    _cbgrid(&CoarseRBGrid),
+    geom(CoarseGrid._ndimension),
+    hermitian(hermitian_),
+    Stencil(&CoarseGrid,geom.npoint,Even,geom.directions,geom.displacements),
+    StencilEven(&CoarseRBGrid,geom.npoint,Even,geom.directions,geom.displacements),
+    StencilOdd(&CoarseRBGrid,geom.npoint,Odd,geom.directions,geom.displacements),
+    A(geom.npoint,&CoarseGrid),
+    Aeven(geom.npoint,&CoarseRBGrid),
+    Aodd(geom.npoint,&CoarseRBGrid),
+    AselfInv(&CoarseGrid),
+    AselfInvEven(&CoarseRBGrid),
+    AselfInvOdd(&CoarseRBGrid),
+    dag_factor(nbasis*nbasis)
+  {
+    fillFactor();
+  };
+
+  void fillFactor() {
+    Eigen::MatrixXd dag_factor_eigen = Eigen::MatrixXd::Ones(nbasis, nbasis);
+    if(!hermitian) {
+      const int nb = nbasis/2;
+      dag_factor_eigen.block(0,nb,nb,nb) *= -1.0;
+      dag_factor_eigen.block(nb,0,nb,nb) *= -1.0;
+    }
+
+    // GPU readable prefactor
+    thread_for(i, nbasis*nbasis, {
+      int j = i/nbasis;
+      int k = i%nbasis;
+      dag_factor[i] = dag_factor_eigen(j, k);
+    });
+  }
+
+  void CoarsenOperator(GridBase *FineGrid,LinearOperatorBase<Lattice<Fobj> > &linop,
+		       Aggregation<Fobj,CComplex,nbasis> & Subspace)
+  {
+    typedef Lattice<typename Fobj::tensor_reduced> FineComplexField;
+    typedef typename Fobj::scalar_type scalar_type;
+
+    std::cout << GridLogMessage<< "CoarsenMatrix "<< std::endl;
+
+    FineComplexField one(FineGrid); one=scalar_type(1.0,0.0);
+    FineComplexField zero(FineGrid); zero=scalar_type(0.0,0.0);
+
+    std::vector<FineComplexField> masks(geom.npoint,FineGrid);
+    FineComplexField imask(FineGrid); // contributions from within this block
+    FineComplexField omask(FineGrid); // contributions from outwith this block
+
+    FineComplexField evenmask(FineGrid);
+    FineComplexField oddmask(FineGrid); 
+
+    FineField     phi(FineGrid);
+    FineField     tmp(FineGrid);
+    FineField     zz(FineGrid); zz=Zero();
+    FineField    Mphi(FineGrid);
+    FineField    Mphie(FineGrid);
+    FineField    Mphio(FineGrid);
+    std::vector<FineField>     Mphi_p(geom.npoint,FineGrid);
+
+    Lattice<iScalar<vInteger> > coor (FineGrid);
+    Lattice<iScalar<vInteger> > bcoor(FineGrid);
+    Lattice<iScalar<vInteger> > bcb  (FineGrid); bcb = Zero();
+
+    CoarseVector iProj(Grid()); 
+    CoarseVector oProj(Grid()); 
+    CoarseVector SelfProj(Grid()); 
+    CoarseComplexField iZProj(Grid()); 
+    CoarseComplexField oZProj(Grid()); 
+
+    CoarseScalar InnerProd(Grid()); 
+
+    std::cout << GridLogMessage<< "CoarsenMatrix Orthog "<< std::endl;
+    // Orthogonalise the subblocks over the basis
+    blockOrthogonalise(InnerProd,Subspace.subspace);
+
+    // Compute the matrix elements of linop between this orthonormal
+    // set of vectors.
+    std::cout << GridLogMessage<< "CoarsenMatrix masks "<< std::endl;
+    int self_stencil=-1;
+    for(int p=0;p<geom.npoint;p++)
+    { 
+      int dir   = geom.directions[p];
+      int disp  = geom.displacements[p];
+      A[p]=Zero();
+      if( geom.displacements[p]==0){
+	self_stencil=p;
+      }
+
+      Integer block=(FineGrid->_rdimensions[dir])/(Grid()->_rdimensions[dir]);
+
+      LatticeCoordinate(coor,dir);
+
+      ///////////////////////////////////////////////////////
+      // Work out even and odd block checkerboarding for fast diagonal term
+      ///////////////////////////////////////////////////////
+      if ( disp==1 ) {
+	bcb   = bcb + div(coor,block);
+      }
+	
+      if ( disp==0 ) {
+	  masks[p]= Zero();
+      } else if ( disp==1 ) {
+	masks[p] = where(mod(coor,block)==(block-1),one,zero);
+      } else if ( disp==-1 ) {
+	masks[p] = where(mod(coor,block)==(Integer)0,one,zero);
+      }
+    }
+    evenmask = where(mod(bcb,2)==(Integer)0,one,zero);
+    oddmask  = one-evenmask;
+
+    assert(self_stencil!=-1);
+
+    for(int i=0;i<nbasis;i++){
+
+      phi=Subspace.subspace[i];
+
+      std::cout << GridLogMessage<< "CoarsenMatrix vector "<<i << std::endl;
+      linop.OpDirAll(phi,Mphi_p);
+      linop.OpDiag  (phi,Mphi_p[geom.npoint-1]);
+
+      for(int p=0;p<geom.npoint;p++){ 
+
+	Mphi = Mphi_p[p];
+
+	int dir   = geom.directions[p];
+	int disp  = geom.displacements[p];
+
+	if ( (disp==-1) || (!hermitian ) ) {
+
+	  ////////////////////////////////////////////////////////////////////////
+	  // Pick out contributions coming from this cell and neighbour cell
+	  ////////////////////////////////////////////////////////////////////////
+	  omask = masks[p];
+	  imask = one-omask;
+	
+	  for(int j=0;j<nbasis;j++){
+	    
+	    blockMaskedInnerProduct(oZProj,omask,Subspace.subspace[j],Mphi);
+	    
+	    autoView( iZProj_v , iZProj, AcceleratorRead) ;
+	    autoView( oZProj_v , oZProj, AcceleratorRead) ;
+	    autoView( A_p     ,  A[p], AcceleratorWrite);
+	    autoView( A_self  , A[self_stencil], AcceleratorWrite);
+
+	    accelerator_for(ss, Grid()->oSites(), Fobj::Nsimd(),{ coalescedWrite(A_p[ss](j,i),oZProj_v(ss)); });
+	    if ( hermitian && (disp==-1) ) {
+	      for(int pp=0;pp<geom.npoint;pp++){// Find the opposite link and set <j|A|i> = <i|A|j>*
+		int dirp   = geom.directions[pp];
+		int dispp  = geom.displacements[pp];
+		if ( (dirp==dir) && (dispp==1) ){
+		  auto sft = conjugate(Cshift(oZProj,dir,1));
+		  autoView( sft_v    ,  sft  , AcceleratorWrite);
+		  autoView( A_pp     ,  A[pp], AcceleratorWrite);
+		  accelerator_for(ss, Grid()->oSites(), Fobj::Nsimd(),{ coalescedWrite(A_pp[ss](i,j),sft_v(ss)); });
+		}
+	      }
+	    }
+
+	  }
+	}
+      }
+
+      ///////////////////////////////////////////
+      // Faster alternate self coupling.. use hermiticity to save 2x
+      ///////////////////////////////////////////
+      {
+	mult(tmp,phi,evenmask);  linop.Op(tmp,Mphie);
+	mult(tmp,phi,oddmask );  linop.Op(tmp,Mphio);
+
+	{
+	  autoView( tmp_      , tmp, AcceleratorWrite);
+	  autoView( evenmask_ , evenmask, AcceleratorRead);
+	  autoView( oddmask_  ,  oddmask, AcceleratorRead);
+	  autoView( Mphie_    ,  Mphie, AcceleratorRead);
+	  autoView( Mphio_    ,  Mphio, AcceleratorRead);
+	  accelerator_for(ss, FineGrid->oSites(), Fobj::Nsimd(),{ 
+	      coalescedWrite(tmp_[ss],evenmask_(ss)*Mphie_(ss) + oddmask_(ss)*Mphio_(ss));
+	    });
+	}
+
+	blockProject(SelfProj,tmp,Subspace.subspace);
+
+	autoView( SelfProj_ , SelfProj, AcceleratorRead);
+	autoView( A_self  , A[self_stencil], AcceleratorWrite);
+
+	accelerator_for(ss, Grid()->oSites(), Fobj::Nsimd(),{
+	  for(int j=0;j<nbasis;j++){
+	    coalescedWrite(A_self[ss](j,i), SelfProj_(ss)(j));
+	  }
+	});
+
+      }
+    }
+    if(hermitian) {
+      std::cout << GridLogMessage << " ForceHermitian, new code "<<std::endl;
+    }
+
+    InvertSelfStencilLink(); std::cout << GridLogMessage << "Coarse self link inverted" << std::endl;
+    FillHalfCbs(); std::cout << GridLogMessage << "Coarse half checkerboards filled" << std::endl;
+  }
+
+  void InvertSelfStencilLink() {
+    std::cout << GridLogDebug << "CoarsenedMatrix::InvertSelfStencilLink" << std::endl;
+    int localVolume = Grid()->lSites();
+
+    typedef typename Cobj::scalar_object scalar_object;
+
+    autoView(Aself_v,    A[geom.npoint-1], CpuRead);
+    autoView(AselfInv_v, AselfInv,         CpuWrite);
+    thread_for(site, localVolume, { // NOTE: Not able to bring this to GPU because of Eigen + peek/poke
+      Eigen::MatrixXcd selfLinkEigen    = Eigen::MatrixXcd::Zero(nbasis, nbasis);
+      Eigen::MatrixXcd selfLinkInvEigen = Eigen::MatrixXcd::Zero(nbasis, nbasis);
+
+      scalar_object selfLink    = Zero();
+      scalar_object selfLinkInv = Zero();
+
+      Coordinate lcoor;
+
+      Grid()->LocalIndexToLocalCoor(site, lcoor);
+      peekLocalSite(selfLink, Aself_v, lcoor);
+
+      for (int i = 0; i < nbasis; ++i)
+        for (int j = 0; j < nbasis; ++j)
+          selfLinkEigen(i, j) = static_cast<ComplexD>(TensorRemove(selfLink(i, j)));
+
+      selfLinkInvEigen = selfLinkEigen.inverse();
+
+      for(int i = 0; i < nbasis; ++i)
+        for(int j = 0; j < nbasis; ++j)
+          selfLinkInv(i, j) = selfLinkInvEigen(i, j);
+
+      pokeLocalSite(selfLinkInv, AselfInv_v, lcoor);
+    });
+  }
+
+  void FillHalfCbs() {
+    std::cout << GridLogDebug << "CoarsenedMatrix::FillHalfCbs" << std::endl;
+    for(int p = 0; p < geom.npoint; ++p) {
+      pickCheckerboard(Even, Aeven[p], A[p]);
+      pickCheckerboard(Odd, Aodd[p], A[p]);
+    }
+    pickCheckerboard(Even, AselfInvEven, AselfInv);
+    pickCheckerboard(Odd, AselfInvOdd, AselfInv);
+  }
+};
+
+NAMESPACE_END(Grid);
+#endif
diff --git a/Grid/algorithms/multigrid/GeneralCoarsenedMatrix.h b/Grid/algorithms/multigrid/GeneralCoarsenedMatrix.h
new file mode 100644
index 00000000..c15a2a63
--- /dev/null
+++ b/Grid/algorithms/multigrid/GeneralCoarsenedMatrix.h
@@ -0,0 +1,418 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/GeneralCoarsenedMatrix.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+#include <Grid/qcd/QCD.h> // needed for Dagger(Yes|No), Inverse(Yes|No)
+
+#include <Grid/lattice/PaddedCell.h>
+#include <Grid/stencil/GeneralLocalStencil.h>
+
+NAMESPACE_BEGIN(Grid);
+
+// Fine Object == (per site) type of fine field
+// nbasis      == number of deflation vectors
+template<class Fobj,class CComplex,int nbasis>
+class GeneralCoarsenedMatrix : public SparseMatrixBase<Lattice<iVector<CComplex,nbasis > > >  {
+public:
+
+  typedef GeneralCoarsenedMatrix<Fobj,CComplex,nbasis> GeneralCoarseOp;
+  typedef iVector<CComplex,nbasis >           siteVector;
+  typedef iMatrix<CComplex,nbasis >           siteMatrix;
+  typedef Lattice<iScalar<CComplex> >         CoarseComplexField;
+  typedef Lattice<siteVector>                 CoarseVector;
+  typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
+  typedef iMatrix<CComplex,nbasis >  Cobj;
+  typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
+  typedef Lattice<Fobj >        FineField;
+  typedef CoarseVector Field;
+  ////////////////////
+  // Data members
+  ////////////////////
+  int hermitian;
+  GridBase      *       _FineGrid; 
+  GridCartesian *       _CoarseGrid; 
+  NonLocalStencilGeometry &geom;
+  PaddedCell Cell;
+  GeneralLocalStencil Stencil;
+  
+  std::vector<CoarseMatrix> _A;
+  std::vector<CoarseMatrix> _Adag;
+
+  ///////////////////////
+  // Interface
+  ///////////////////////
+  GridBase      * Grid(void)           { return _FineGrid; };   // this is all the linalg routines need to know
+  GridBase      * FineGrid(void)       { return _FineGrid; };   // this is all the linalg routines need to know
+  GridCartesian * CoarseGrid(void)     { return _CoarseGrid; };   // this is all the linalg routines need to know
+
+  void ProjectNearestNeighbour(RealD shift, GeneralCoarseOp &CopyMe)
+  {
+    int nfound=0;
+    std::cout << GridLogMessage <<"GeneralCoarsenedMatrix::ProjectNearestNeighbour "<< CopyMe._A[0].Grid()<<std::endl;
+    for(int p=0;p<geom.npoint;p++){
+      for(int pp=0;pp<CopyMe.geom.npoint;pp++){
+ 	// Search for the same relative shift
+	// Avoids brutal handling of Grid pointers
+	if ( CopyMe.geom.shifts[pp]==geom.shifts[p] ) {
+	  _A[p] = CopyMe.Cell.Extract(CopyMe._A[pp]);
+	  _Adag[p] = CopyMe.Cell.Extract(CopyMe._Adag[pp]);
+	  nfound++;
+	}
+      }
+    }
+    assert(nfound==geom.npoint);
+    ExchangeCoarseLinks();
+  }
+  
+  GeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridBase *FineGrid, GridCartesian * CoarseGrid)
+    : geom(_geom),
+      _FineGrid(FineGrid),
+      _CoarseGrid(CoarseGrid),
+      hermitian(1),
+      Cell(_geom.Depth(),_CoarseGrid),
+      Stencil(Cell.grids.back(),geom.shifts)
+  {
+    {
+      int npoint = _geom.npoint;
+      autoView( Stencil_v  , Stencil, AcceleratorRead);
+      int osites=Stencil.Grid()->oSites();
+      for(int ss=0;ss<osites;ss++){
+	for(int point=0;point<npoint;point++){
+	  auto SE = Stencil_v.GetEntry(point,ss);
+	  int o = SE->_offset;
+	  assert( o< osites);
+	}
+      }    
+    }
+
+    _A.resize(geom.npoint,CoarseGrid);
+    _Adag.resize(geom.npoint,CoarseGrid);
+  }
+  void M (const CoarseVector &in, CoarseVector &out)
+  {
+    Mult(_A,in,out);
+  }
+  void Mdag (const CoarseVector &in, CoarseVector &out)
+  {
+    if ( hermitian ) M(in,out);
+    else Mult(_Adag,in,out);
+  }
+  void Mult (std::vector<CoarseMatrix> &A,const CoarseVector &in, CoarseVector &out)
+  {
+    RealD tviews=0;
+    RealD ttot=0;
+    RealD tmult=0;
+    RealD texch=0;
+    RealD text=0;
+    ttot=-usecond();
+    conformable(CoarseGrid(),in.Grid());
+    conformable(in.Grid(),out.Grid());
+    out.Checkerboard() = in.Checkerboard();
+    CoarseVector tin=in;
+
+    texch-=usecond();
+    CoarseVector pin  = Cell.Exchange(tin);
+    texch+=usecond();
+
+    CoarseVector pout(pin.Grid()); pout=Zero();
+
+    int npoint = geom.npoint;
+    typedef LatticeView<Cobj> Aview;
+      
+    const int Nsimd = CComplex::Nsimd();
+    
+    int osites=pin.Grid()->oSites();
+    //    int gsites=pin.Grid()->gSites();
+
+    RealD flops = 1.0* npoint * nbasis * nbasis * 8 * osites;
+    RealD bytes = (1.0*osites*sizeof(siteMatrix)*npoint+2.0*osites*sizeof(siteVector))*npoint;
+      
+    //    for(int point=0;point<npoint;point++){
+    //      conformable(A[point],pin);
+    //    }
+
+    {
+      tviews-=usecond();
+      autoView( in_v , pin, AcceleratorRead);
+      autoView( out_v , pout, AcceleratorWrite);
+      autoView( Stencil_v  , Stencil, AcceleratorRead);
+      tviews+=usecond();
+      
+      for(int point=0;point<npoint;point++){
+	tviews-=usecond();
+	autoView( A_v, A[point],AcceleratorRead);
+	tviews+=usecond();
+	tmult-=usecond();
+	accelerator_for(sss, osites*nbasis, Nsimd, {
+
+	    typedef decltype(coalescedRead(in_v[0]))    calcVector;
+
+	    int ss = sss/nbasis;
+	    int b  = sss%nbasis;
+
+	    auto SE  = Stencil_v.GetEntry(point,ss);
+	    auto nbr = coalescedReadGeneralPermute(in_v[SE->_offset],SE->_permute,Nd);
+	    auto res = out_v(ss)(b);
+	    for(int bb=0;bb<nbasis;bb++) {
+	      res = res + coalescedRead(A_v[ss](b,bb))*nbr(bb);
+	    }
+	    coalescedWrite(out_v[ss](b),res);
+	});
+
+	tmult+=usecond();
+      }
+    }
+    text-=usecond();
+    out = Cell.Extract(pout);
+    text+=usecond();
+    ttot+=usecond();
+
+    std::cout << GridLogPerformance<<"Coarse Mult Aviews "<<tviews<<" us"<<std::endl;
+    std::cout << GridLogPerformance<<"Coarse Mult exch "<<texch<<" us"<<std::endl;
+    std::cout << GridLogPerformance<<"Coarse Mult mult "<<tmult<<" us"<<std::endl;
+    std::cout << GridLogPerformance<<"Coarse Mult ext  "<<text<<" us"<<std::endl;
+    std::cout << GridLogPerformance<<"Coarse Mult tot  "<<ttot<<" us"<<std::endl;
+    std::cout << GridLogPerformance<<"Coarse Kernel flop/s "<< flops/tmult<<" mflop/s"<<std::endl;
+    std::cout << GridLogPerformance<<"Coarse Kernel bytes/s"<< bytes/tmult<<" MB/s"<<std::endl;
+    std::cout << GridLogPerformance<<"Coarse overall flops/s "<< flops/ttot<<" mflop/s"<<std::endl;
+    std::cout << GridLogPerformance<<"Coarse total bytes   "<< bytes/1e6<<" MB"<<std::endl;
+  };
+
+  void PopulateAdag(void)
+  {
+    for(int64_t bidx=0;bidx<CoarseGrid()->gSites() ;bidx++){
+      Coordinate bcoor;
+      CoarseGrid()->GlobalIndexToGlobalCoor(bidx,bcoor);
+      
+      for(int p=0;p<geom.npoint;p++){
+	Coordinate scoor = bcoor;
+	for(int mu=0;mu<bcoor.size();mu++){
+	  int L = CoarseGrid()->GlobalDimensions()[mu];
+	  scoor[mu] = (bcoor[mu] - geom.shifts[p][mu] + L) % L; // Modulo arithmetic
+	}
+	// Flip to poke/peekLocalSite and not too bad
+	auto link = peekSite(_A[p],scoor);
+	int pp = geom.Reverse(p);
+	pokeSite(adj(link),_Adag[pp],bcoor);
+      }
+    }
+  }
+  /////////////////////////////////////////////////////////////
+  // 
+  // A) Only reduced flops option is to use a padded cell of depth 4
+  // and apply MpcDagMpc in the padded cell.
+  //
+  // Makes for ONE application of MpcDagMpc per vector instead of 30 or 80.
+  // With the effective cell size around (B+8)^4 perhaps 12^4/4^4 ratio
+  // Cost is 81x more, same as stencil size.
+  //
+  // But: can eliminate comms and do as local dirichlet.
+  //
+  // Local exchange gauge field once.
+  // Apply to all vectors, local only computation.
+  // Must exchange ghost subcells in reverse process of PaddedCell to take inner products
+  //
+  // B) Can reduce cost: pad by 1, apply Deo      (4^4+6^4+8^4+8^4 )/ (4x 4^4)
+  //                     pad by 2, apply Doe
+  //                     pad by 3, apply Deo
+  //                     then break out 8x directions; cost is ~10x MpcDagMpc per vector
+  //
+  // => almost factor of 10 in setup cost, excluding data rearrangement
+  //
+  // Intermediates -- ignore the corner terms, leave approximate and force Hermitian
+  // Intermediates -- pad by 2 and apply 1+8+24 = 33 times.
+  /////////////////////////////////////////////////////////////
+
+    //////////////////////////////////////////////////////////
+    // BFM HDCG style approach: Solve a system of equations to get Aij
+    //////////////////////////////////////////////////////////
+    /*
+     *     Here, k,l index which possible shift within the 3^Nd "ball" connected by MdagM.
+     *
+     *     conj(phases[block]) proj[k][ block*Nvec+j ] =  \sum_ball  e^{i q_k . delta} < phi_{block,j} | MdagM | phi_{(block+delta),i} > 
+     *                                                 =  \sum_ball e^{iqk.delta} A_ji
+     *
+     *     Must invert matrix M_k,l = e^[i q_k . delta_l]
+     *
+     *     Where q_k = delta_k . (2*M_PI/global_nb[mu])
+     */
+  void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
+		       Aggregation<Fobj,CComplex,nbasis> & Subspace)
+  {
+    std::cout << GridLogMessage<< "GeneralCoarsenMatrix "<< std::endl;
+    GridBase *grid = FineGrid();
+
+    RealD tproj=0.0;
+    RealD teigen=0.0;
+    RealD tmat=0.0;
+    RealD tphase=0.0;
+    RealD tinv=0.0;
+
+    /////////////////////////////////////////////////////////////
+    // Orthogonalise the subblocks over the basis
+    /////////////////////////////////////////////////////////////
+    CoarseScalar InnerProd(CoarseGrid()); 
+    blockOrthogonalise(InnerProd,Subspace.subspace);
+
+    const int npoint = geom.npoint;
+      
+    Coordinate clatt = CoarseGrid()->GlobalDimensions();
+    int Nd = CoarseGrid()->Nd();
+
+      /*
+       *     Here, k,l index which possible momentum/shift within the N-points connected by MdagM.
+       *     Matrix index i is mapped to this shift via 
+       *               geom.shifts[i]
+       *
+       *     conj(pha[block]) proj[k (which mom)][j (basis vec cpt)][block] 
+       *       =  \sum_{l in ball}  e^{i q_k . delta_l} < phi_{block,j} | MdagM | phi_{(block+delta_l),i} > 
+       *       =  \sum_{l in ball} e^{iqk.delta_l} A_ji^{b.b+l}
+       *       = M_{kl} A_ji^{b.b+l}
+       *
+       *     Must assemble and invert matrix M_k,l = e^[i q_k . delta_l]
+       *  
+       *     Where q_k = delta_k . (2*M_PI/global_nb[mu])
+       *
+       *     Then A{ji}^{b,b+l} = M^{-1}_{lm} ComputeProj_{m,b,i,j}
+       */
+    teigen-=usecond();
+    Eigen::MatrixXcd Mkl    = Eigen::MatrixXcd::Zero(npoint,npoint);
+    Eigen::MatrixXcd invMkl = Eigen::MatrixXcd::Zero(npoint,npoint);
+    ComplexD ci(0.0,1.0);
+    for(int k=0;k<npoint;k++){ // Loop over momenta
+
+      for(int l=0;l<npoint;l++){ // Loop over nbr relative
+	ComplexD phase(0.0,0.0);
+	for(int mu=0;mu<Nd;mu++){
+	  RealD TwoPiL =  M_PI * 2.0/ clatt[mu];
+	  phase=phase+TwoPiL*geom.shifts[k][mu]*geom.shifts[l][mu];
+	}
+	phase=exp(phase*ci);
+	Mkl(k,l) = phase;
+      }
+    }
+    invMkl = Mkl.inverse();
+    teigen+=usecond();
+
+    ///////////////////////////////////////////////////////////////////////
+    // Now compute the matrix elements of linop between the orthonormal
+    // set of vectors.
+    ///////////////////////////////////////////////////////////////////////
+    FineField phaV(grid); // Phased block basis vector
+    FineField MphaV(grid);// Matrix applied
+    CoarseVector coarseInner(CoarseGrid());
+
+    std::vector<CoarseVector> ComputeProj(npoint,CoarseGrid());
+    std::vector<CoarseVector>          FT(npoint,CoarseGrid());
+    for(int i=0;i<nbasis;i++){// Loop over basis vectors
+      std::cout << GridLogMessage<< "CoarsenMatrixColoured vec "<<i<<"/"<<nbasis<< std::endl;
+      for(int p=0;p<npoint;p++){ // Loop over momenta in npoint
+	/////////////////////////////////////////////////////
+	// Stick a phase on every block
+	/////////////////////////////////////////////////////
+	tphase-=usecond();
+	CoarseComplexField coor(CoarseGrid());
+	CoarseComplexField pha(CoarseGrid());	pha=Zero();
+	for(int mu=0;mu<Nd;mu++){
+	  LatticeCoordinate(coor,mu);
+	  RealD TwoPiL =  M_PI * 2.0/ clatt[mu];
+	  pha = pha + (TwoPiL * geom.shifts[p][mu]) * coor;
+	}
+	pha  =exp(pha*ci);
+	phaV=Zero();
+	blockZAXPY(phaV,pha,Subspace.subspace[i],phaV);
+	tphase+=usecond();
+
+	/////////////////////////////////////////////////////////////////////
+	// Multiple phased subspace vector by matrix and project to subspace
+	// Remove local bulk phase to leave relative phases
+	/////////////////////////////////////////////////////////////////////
+	tmat-=usecond();
+	linop.Op(phaV,MphaV);
+	tmat+=usecond();
+
+	tproj-=usecond();
+	blockProject(coarseInner,MphaV,Subspace.subspace);
+	coarseInner = conjugate(pha) * coarseInner;
+
+	ComputeProj[p] = coarseInner;
+	tproj+=usecond();
+
+      }
+
+      tinv-=usecond();
+      for(int k=0;k<npoint;k++){
+	FT[k] = Zero();
+	for(int l=0;l<npoint;l++){
+	  FT[k]= FT[k]+ invMkl(l,k)*ComputeProj[l];
+	}
+      
+	int osites=CoarseGrid()->oSites();
+	autoView( A_v  , _A[k], AcceleratorWrite);
+	autoView( FT_v  , FT[k], AcceleratorRead);
+	accelerator_for(sss, osites, 1, {
+	    for(int j=0;j<nbasis;j++){
+	      A_v[sss](j,i) = FT_v[sss](j);
+	    }
+        });
+      }
+      tinv+=usecond();
+    }
+
+    for(int p=0;p<geom.npoint;p++){
+      Coordinate coor({0,0,0,0,0});
+      auto sval = peekSite(_A[p],coor);
+    }
+
+    // Only needed if nonhermitian
+    if ( ! hermitian ) {
+      std::cout << GridLogMessage<<"PopulateAdag  "<<std::endl;
+      PopulateAdag();
+    }
+
+    // Need to write something to populate Adag from A
+    ExchangeCoarseLinks();
+    std::cout << GridLogMessage<<"CoarsenOperator eigen  "<<teigen<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator phase  "<<tphase<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator mat    "<<tmat <<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator proj   "<<tproj<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator inv    "<<tinv<<" us"<<std::endl;
+  }
+  void ExchangeCoarseLinks(void){
+    for(int p=0;p<geom.npoint;p++){
+      _A[p] = Cell.Exchange(_A[p]);
+      _Adag[p]= Cell.Exchange(_Adag[p]);
+    }
+  }
+  virtual  void Mdiag    (const Field &in, Field &out){ assert(0);};
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);};
+  virtual  void MdirAll  (const Field &in, std::vector<Field> &out){assert(0);};
+};
+
+
+NAMESPACE_END(Grid);
diff --git a/Grid/algorithms/multigrid/Geometry.h b/Grid/algorithms/multigrid/Geometry.h
new file mode 100644
index 00000000..74db4417
--- /dev/null
+++ b/Grid/algorithms/multigrid/Geometry.h
@@ -0,0 +1,243 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/GeneralCoarsenedMatrix.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+
+/////////////////////////////////////////////////////////////////
+// Geometry class in cartesian case
+/////////////////////////////////////////////////////////////////
+
+class Geometry {
+public:
+  int npoint;
+  int base;
+  std::vector<int> directions   ;
+  std::vector<int> displacements;
+  std::vector<int> points_dagger;
+
+  Geometry(int _d)  {
+    
+    base = (_d==5) ? 1:0;
+
+    // make coarse grid stencil for 4d , not 5d
+    if ( _d==5 ) _d=4;
+
+    npoint = 2*_d+1;
+    directions.resize(npoint);
+    displacements.resize(npoint);
+    points_dagger.resize(npoint);
+    for(int d=0;d<_d;d++){
+      directions[d   ] = d+base;
+      directions[d+_d] = d+base;
+      displacements[d  ] = +1;
+      displacements[d+_d]= -1;
+      points_dagger[d   ] = d+_d;
+      points_dagger[d+_d] = d;
+    }
+    directions   [2*_d]=0;
+    displacements[2*_d]=0;
+    points_dagger[2*_d]=2*_d;
+  }
+
+  int point(int dir, int disp) {
+    assert(disp == -1 || disp == 0 || disp == 1);
+    assert(base+0 <= dir && dir < base+4);
+
+    // directions faster index = new indexing
+    // 4d (base = 0):
+    // point 0  1  2  3  4  5  6  7  8
+    // dir   0  1  2  3  0  1  2  3  0
+    // disp +1 +1 +1 +1 -1 -1 -1 -1  0
+    // 5d (base = 1):
+    // point 0  1  2  3  4  5  6  7  8
+    // dir   1  2  3  4  1  2  3  4  0
+    // disp +1 +1 +1 +1 -1 -1 -1 -1  0
+
+    // displacements faster index = old indexing
+    // 4d (base = 0):
+    // point 0  1  2  3  4  5  6  7  8
+    // dir   0  0  1  1  2  2  3  3  0
+    // disp +1 -1 +1 -1 +1 -1 +1 -1  0
+    // 5d (base = 1):
+    // point 0  1  2  3  4  5  6  7  8
+    // dir   1  1  2  2  3  3  4  4  0
+    // disp +1 -1 +1 -1 +1 -1 +1 -1  0
+
+    if(dir == 0 and disp == 0)
+      return 8;
+    else // New indexing
+      return (1 - disp) / 2 * 4 + dir - base;
+    // else // Old indexing
+    //   return (4 * (dir - base) + 1 - disp) / 2;
+  }
+};
+
+/////////////////////////////////////////////////////////////////
+// Less local equivalent of Geometry class in cartesian case
+/////////////////////////////////////////////////////////////////
+class NonLocalStencilGeometry {
+public:
+  int depth;
+  int hops;
+  int npoint;
+  std::vector<Coordinate> shifts;
+  Coordinate stencil_size;
+  Coordinate stencil_lo;
+  Coordinate stencil_hi;
+  GridCartesian *grid;
+  GridCartesian *Grid() {return grid;};
+  int Depth(void){return 1;};   // Ghost zone depth
+  int Hops(void){return hops;}; // # of hops=> level of corner fill in in stencil
+
+  virtual int DimSkip(void) =0;
+
+  virtual ~NonLocalStencilGeometry() {};
+
+  int  Reverse(int point)
+  {
+    int Nd = Grid()->Nd();
+    Coordinate shft = shifts[point];
+    Coordinate rev(Nd);
+    for(int mu=0;mu<Nd;mu++) rev[mu]= -shft[mu];
+    for(int p=0;p<npoint;p++){
+      if(rev==shifts[p]){
+	return p;
+      }
+    }
+    assert(0);
+    return -1;
+  }
+  void BuildShifts(void)
+  {
+    this->shifts.resize(0);
+    int Nd = this->grid->Nd();
+
+    int dd = this->DimSkip();
+    for(int s0=this->stencil_lo[dd+0];s0<=this->stencil_hi[dd+0];s0++){
+    for(int s1=this->stencil_lo[dd+1];s1<=this->stencil_hi[dd+1];s1++){
+    for(int s2=this->stencil_lo[dd+2];s2<=this->stencil_hi[dd+2];s2++){
+    for(int s3=this->stencil_lo[dd+3];s3<=this->stencil_hi[dd+3];s3++){
+      Coordinate sft(Nd,0);
+      sft[dd+0] = s0;
+      sft[dd+1] = s1;
+      sft[dd+2] = s2;
+      sft[dd+3] = s3;
+      int nhops = abs(s0)+abs(s1)+abs(s2)+abs(s3);
+      if(nhops<=this->hops) this->shifts.push_back(sft);
+    }}}}
+    this->npoint = this->shifts.size();
+    std::cout << GridLogMessage << "NonLocalStencilGeometry has "<< this->npoint << " terms in stencil "<<std::endl;
+  }
+  
+  NonLocalStencilGeometry(GridCartesian *_coarse_grid,int _hops) : grid(_coarse_grid), hops(_hops)
+  {
+    Coordinate latt = grid->GlobalDimensions();
+    stencil_size.resize(grid->Nd());
+    stencil_lo.resize(grid->Nd());
+    stencil_hi.resize(grid->Nd());
+    for(int d=0;d<grid->Nd();d++){
+     if ( latt[d] == 1 ) {
+      stencil_lo[d] = 0;
+      stencil_hi[d] = 0;
+      stencil_size[d]= 1;
+     } else if ( latt[d] == 2 ) {
+      stencil_lo[d] = -1;
+      stencil_hi[d] = 0;
+      stencil_size[d]= 2;
+     } else if ( latt[d] > 2 ) {
+       stencil_lo[d] = -1;
+       stencil_hi[d] =  1;
+       stencil_size[d]= 3;
+     }
+    }
+  };
+
+};
+
+// Need to worry about red-black now
+class NonLocalStencilGeometry4D : public NonLocalStencilGeometry {
+public:
+  virtual int DimSkip(void) { return 0;};
+  NonLocalStencilGeometry4D(GridCartesian *Coarse,int _hops) : NonLocalStencilGeometry(Coarse,_hops) { };
+  virtual ~NonLocalStencilGeometry4D() {};
+};
+class NonLocalStencilGeometry5D : public NonLocalStencilGeometry {
+public:
+  virtual int DimSkip(void) { return 1; }; 
+  NonLocalStencilGeometry5D(GridCartesian *Coarse,int _hops) : NonLocalStencilGeometry(Coarse,_hops)  { };
+  virtual ~NonLocalStencilGeometry5D() {};
+};
+/*
+ * Bunch of different options classes
+ */
+class NextToNextToNextToNearestStencilGeometry4D : public NonLocalStencilGeometry4D {
+public:
+  NextToNextToNextToNearestStencilGeometry4D(GridCartesian *Coarse) :  NonLocalStencilGeometry4D(Coarse,4)
+  {
+    this->BuildShifts();
+  };
+};
+class NextToNextToNextToNearestStencilGeometry5D : public  NonLocalStencilGeometry5D {
+public:
+  NextToNextToNextToNearestStencilGeometry5D(GridCartesian *Coarse) :  NonLocalStencilGeometry5D(Coarse,4)
+  {
+    this->BuildShifts();
+  };
+};
+class NextToNearestStencilGeometry4D : public  NonLocalStencilGeometry4D {
+public:
+  NextToNearestStencilGeometry4D(GridCartesian *Coarse) :  NonLocalStencilGeometry4D(Coarse,2)
+  {
+    this->BuildShifts();
+  };
+};
+class NextToNearestStencilGeometry5D : public  NonLocalStencilGeometry5D {
+public:
+  NextToNearestStencilGeometry5D(GridCartesian *Coarse) :  NonLocalStencilGeometry5D(Coarse,2)
+  {
+    this->BuildShifts();
+  };
+};
+class NearestStencilGeometry4D : public  NonLocalStencilGeometry4D {
+public:
+  NearestStencilGeometry4D(GridCartesian *Coarse) :  NonLocalStencilGeometry4D(Coarse,1)
+  {
+    this->BuildShifts();
+  };
+};
+class NearestStencilGeometry5D : public  NonLocalStencilGeometry5D {
+public:
+  NearestStencilGeometry5D(GridCartesian *Coarse) :  NonLocalStencilGeometry5D(Coarse,1)
+  {
+    this->BuildShifts();
+  };
+};
+
+NAMESPACE_END(Grid);
diff --git a/Grid/algorithms/multigrid/Multigrid.h b/Grid/algorithms/multigrid/Multigrid.h
new file mode 100644
index 00000000..a1bf400c
--- /dev/null
+++ b/Grid/algorithms/multigrid/Multigrid.h
@@ -0,0 +1,33 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: Grid/algorithms/multigrid/MultiGrid.h
+
+    Copyright (C) 2023
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#pragma once
+
+#include <Grid/algorithms/multigrid/Aggregates.h>
+#include <Grid/algorithms/multigrid/Geometry.h>
+#include <Grid/algorithms/multigrid/CoarsenedMatrix.h>
+#include <Grid/algorithms/multigrid/GeneralCoarsenedMatrix.h>