diff --git a/Grid/algorithms/multigrid/deprecated/GeneralCoarsenedMatrix.h b/Grid/algorithms/multigrid/deprecated/GeneralCoarsenedMatrix.h
new file mode 100644
index 000000000..2b988b64c
--- /dev/null
+++ b/Grid/algorithms/multigrid/deprecated/GeneralCoarsenedMatrix.h
@@ -0,0 +1,704 @@
+/*************************************************************************************
+
+    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>
+#include <Grid/algorithms/deflation/MultiRHSBlockProject.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 iVector<CComplex,nbasis >  Cvec;
+  typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
+  typedef Lattice<Fobj >        FineField;
+  typedef Lattice<CComplex >    FineComplexField;
+  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;
+  std::vector<CoarseVector> MultTemporaries;
+
+  int64_t MultCalls;
+  double  MultFlopsAccum;
+  double  MultUsecAccum;
+
+  ///////////////////////
+  // Interface
+  ///////////////////////
+  //////////////////////////////////////////////////////////////////////////
+  // Bilingual accessors: everything a consumer needs to read the operator
+  // without knowing which of the three coarse classes it holds. The D
+  // dimensional grid the elements live on, the geometry they are indexed by,
+  // and one unpadded point at a time (a whole npoint vector is too much
+  // memory at production nbasis).
+  //////////////////////////////////////////////////////////////////////////
+  GridCartesian * CoarseGridD(void)             { return _CoarseGrid; };
+  NonLocalStencilGeometry & Geometry(void)      { return geom; };
+  void ExtractMatrix(int p,CoarseMatrix &A)     { A = Cell.Extract(_A[p]); };
+
+  GridBase      * Grid(void)           { return _CoarseGrid; };   // 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 ShiftMatrix(RealD shift)
+  {
+    int Nd=_FineGrid->Nd(); 
+    Coordinate zero_shift(Nd,0);
+    for(int p=0;p<geom.npoint;p++){
+      if ( zero_shift==geom.shifts[p] ) {
+	_A[p] = _A[p]+shift;
+	//	_Adag[p] = _Adag[p]+shift;
+      }
+    }    
+  }
+  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++;
+	}
+      }
+    }
+    GRID_ASSERT(nfound==geom.npoint);
+    ExchangeCoarseLinks();
+  }
+  */
+  
+  GeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridBase *FineGrid, GridCartesian * CoarseGrid,int _herm=1)
+    : geom(_geom),
+      _FineGrid(FineGrid),
+      _CoarseGrid(CoarseGrid),
+      hermitian(_herm),
+      Cell(_geom.Depth(),_CoarseGrid),
+      Stencil(Cell.grids.back(),geom.shifts),
+      MultCalls(0), MultFlopsAccum(0.0), MultUsecAccum(0.0)
+  {
+    {
+      int npoint = _geom.npoint;
+    }
+    _A.resize(geom.npoint,CoarseGrid);
+    if ( !hermitian ) _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)
+      Mult(_A,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; RealD ttemps=0; RealD tcopy=0;
+    RealD tmult2=0;
+
+    ttot=-usecond();
+    conformable(CoarseGrid(),in.Grid());
+    conformable(in.Grid(),out.Grid());
+    out.Checkerboard() = in.Checkerboard();
+    CoarseVector tin=in;
+
+    texch-=usecond();
+    CoarseVector pin = Cell.ExchangePeriodic(tin);
+    texch+=usecond();
+
+    CoarseVector pout(pin.Grid());
+
+    int npoint = geom.npoint;
+    typedef LatticeView<Cobj> Aview;
+    typedef LatticeView<Cvec> Vview;
+      
+    const int Nsimd = CComplex::Nsimd();
+    
+    int64_t osites=pin.Grid()->oSites();
+
+    RealD flops = 1.0* npoint * nbasis * nbasis * 8.0 * osites * CComplex::Nsimd();
+    RealD bytes = 1.0*osites*sizeof(siteMatrix)*npoint
+                + 2.0*osites*sizeof(siteVector)*npoint;
+      
+    {
+      tviews-=usecond();
+      autoView( in_v , pin, AcceleratorRead);
+      autoView( out_v , pout, AcceleratorWriteDiscard);
+      autoView( Stencil_v  , Stencil, AcceleratorRead);
+      tviews+=usecond();
+
+      // Static and prereserve to keep UVM region live and not resized across multiple calls
+      ttemps-=usecond();
+      MultTemporaries.resize(npoint,pin.Grid());       
+      ttemps+=usecond();
+      std::vector<Aview> AcceleratorViewContainer_h;
+      std::vector<Vview> AcceleratorVecViewContainer_h; 
+
+      tviews-=usecond();
+      for(int p=0;p<npoint;p++) {
+	AcceleratorViewContainer_h.push_back(      A[p].View(AcceleratorRead));
+	AcceleratorVecViewContainer_h.push_back(MultTemporaries[p].View(AcceleratorWrite));
+      }
+      tviews+=usecond();
+
+      static deviceVector<Aview> AcceleratorViewContainer; AcceleratorViewContainer.resize(npoint);
+      static deviceVector<Vview> AcceleratorVecViewContainer; AcceleratorVecViewContainer.resize(npoint); 
+      
+      auto Aview_p = &AcceleratorViewContainer[0];
+      auto Vview_p = &AcceleratorVecViewContainer[0];
+      tcopy-=usecond();
+      acceleratorCopyToDevice(&AcceleratorViewContainer_h[0],&AcceleratorViewContainer[0],npoint *sizeof(Aview));
+      acceleratorCopyToDevice(&AcceleratorVecViewContainer_h[0],&AcceleratorVecViewContainer[0],npoint *sizeof(Vview));
+      tcopy+=usecond();
+
+      tmult-=usecond();
+      accelerator_for(spb, osites*nbasis*npoint, Nsimd, {
+	  typedef decltype(coalescedRead(in_v[0](0))) calcComplex;
+	  int32_t ss   = spb/(nbasis*npoint);
+	  int32_t bp   = spb%(nbasis*npoint);
+	  int32_t point= bp/nbasis;
+	  int32_t b    = bp%nbasis;
+	  auto SE  = Stencil_v.GetEntry(point,ss);
+	  auto nbr = coalescedReadGeneralPermute(in_v[SE->_offset],SE->_permute,Nd);
+	  auto res = coalescedRead(Aview_p[point][ss](0,b))*nbr(0);
+	  for(int bb=1;bb<nbasis;bb++) {
+	    res = res + coalescedRead(Aview_p[point][ss](bb,b))*nbr(bb);
+	  }
+	  coalescedWrite(Vview_p[point][ss](b),res);
+      });
+      tmult2-=usecond();
+      accelerator_for(sb, osites*nbasis, Nsimd, {
+	  int ss = sb/nbasis;
+	  int b  = sb%nbasis;
+	  auto res = coalescedRead(Vview_p[0][ss](b));
+	  for(int point=1;point<npoint;point++){
+	    res = res + coalescedRead(Vview_p[point][ss](b));
+	  }
+	  coalescedWrite(out_v[ss](b),res);
+      });
+      tmult2+=usecond();
+      tmult+=usecond();
+      for(int p=0;p<npoint;p++) {
+	AcceleratorViewContainer_h[p].ViewClose();
+	AcceleratorVecViewContainer_h[p].ViewClose();
+      }
+    }
+
+    text-=usecond();
+    out = Cell.Extract(pout);
+    text+=usecond();
+    ttot+=usecond();
+    
+    MultCalls++;
+    MultFlopsAccum += flops;
+    MultUsecAccum  += ttot;
+    std::cout << GridLogPerformance
+	      << "Coarse Mult call " << MultCalls
+	      << "  tot " << ttot << " us"
+	      << "  kernel " << tmult << " us"
+	      << "  kernel " << flops/tmult*1e-3 << " GFlop/s"
+	      << "  overall " << MultFlopsAccum/MultUsecAccum*1e-3 << " GFlop/s (cumul)"
+	      << "  bw " << bytes/tmult*1e-3 << " GB/s"
+	      << std::endl;
+
+  };
+  
+  void PopulateAdag(void)
+  {
+#if 0
+    // Serial global peek/poke reference implementation
+    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);
+      }
+    }
+#else
+    // Parallel: _Adag[pp](x) = adj( _A[p](x + s_pp) ),  pp = Reverse(p), s_pp = -s_p.
+    // The neighbour fetch reuses the same padded-cell + stencil machinery as Mult,
+    // reading one matrix element per coalesced access so no whole site matrix
+    // (230KB at nbasis=60) ever lands on a GPU thread stack (HIP limit 128KB).
+    // Halo sites compute garbage neighbours; Cell.Extract discards them.
+    // Must run on the unpadded _A, i.e. before ExchangeCoarseLinks.
+    const int Nsimd = CComplex::Nsimd();
+    for(int p=0;p<geom.npoint;p++){
+      int pp = geom.Reverse(p);
+      CoarseMatrix Apad = Cell.ExchangePeriodic(_A[p]);
+      CoarseMatrix Dpad(Apad.Grid());
+      int64_t osites = Apad.Grid()->oSites();
+      {
+	autoView( Apad_v   , Apad, AcceleratorRead);
+	autoView( Dpad_v   , Dpad, AcceleratorWriteDiscard);
+	autoView( Stencil_v, Stencil, AcceleratorRead);
+	accelerator_for(sj, osites*nbasis, Nsimd, {
+	    int32_t ss = sj/nbasis;
+	    int32_t j  = sj%nbasis;
+	    auto SE = Stencil_v.GetEntry(pp,ss);
+	    for(int i=0;i<nbasis;i++){
+	      auto z = coalescedReadGeneralPermute(Apad_v[SE->_offset](i,j),SE->_permute,Nd);
+	      coalescedWrite(Dpad_v[ss](j,i),conjugate(z));
+	    }
+	});
+      }
+      _Adag[pp] = Cell.Extract(Dpad);
+    }
+#endif
+  }
+  /////////////////////////////////////////////////////////////
+  // 
+  // 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])
+     */
+#if 0
+  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, nbasis, {
+#ifdef GRID_SIMT
+	    int j = acceleratorSIMTlane(nbasis);
+	    A_v[sss](i,j) = FT_v[sss](j);
+#else
+	    // CPU build: acceleratorSIMTlane()==0 -- an un-looped SIMT tensor
+	    // index writes ONLY j=0 and silently drops the other nbasis-1
+	    // columns (caught by Test_schur_dense_coarse import certificate,
+	    // 2026-08-14).  Loop explicitly.
+	    for(int j=0;j<nbasis;j++){
+	      A_v[sss](i,j) = FT_v[sss](j);
+	    }
+#endif
+        });
+      }
+      tinv+=usecond();
+    }
+
+    // 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;
+  }
+#else
+  //////////////////////////////////////////////////////////////////////
+  // Galerkin projection of matrix
+  //////////////////////////////////////////////////////////////////////
+  void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
+		       Aggregation<Fobj,CComplex,nbasis> & Subspace)
+  {
+    CoarsenOperator(linop,Subspace,Subspace);
+  }
+  //////////////////////////////////////////////////////////////////////
+  // Petrov - Galerkin projection of matrix
+  //////////////////////////////////////////////////////////////////////
+  void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
+		       Aggregation<Fobj,CComplex,nbasis> & U,
+		       Aggregation<Fobj,CComplex,nbasis> & V)
+  {
+    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 tphaseBZ=0.0;
+    RealD tinv=0.0;
+
+    /////////////////////////////////////////////////////////////
+    // Orthogonalise the subblocks over the basis
+    /////////////////////////////////////////////////////////////
+    CoarseScalar InnerProd(CoarseGrid()); 
+    blockOrthogonalise(InnerProd,V.subspace);
+    blockOrthogonalise(InnerProd,U.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.
+    ///////////////////////////////////////////////////////////////////////
+    std::vector<FineComplexField> phaF(npoint,grid);
+    std::vector<CoarseComplexField> pha(npoint,CoarseGrid());
+
+    typedef typename CComplex::scalar_type SComplex;
+    FineComplexField one(grid); one=SComplex(1.0);
+    FineComplexField zz(grid); zz = Zero();
+    tphase=-usecond();
+    for(int p=0;p<npoint;p++){ // Loop over momenta in npoint
+      /////////////////////////////////////////////////////
+      // Stick a phase on every block
+      /////////////////////////////////////////////////////
+      CoarseComplexField coor(CoarseGrid());
+      pha[p]=Zero();
+      for(int mu=0;mu<Nd;mu++){
+	LatticeCoordinate(coor,mu);
+	RealD TwoPiL =  M_PI * 2.0/ clatt[mu];
+	pha[p] = pha[p] + (TwoPiL * geom.shifts[p][mu]) * coor;
+      }
+      pha[p]  =exp(pha[p]*ci);
+
+      blockZAXPY(phaF[p],pha[p],one,zz);
+
+    }
+    tphase+=usecond();
+
+    // Import basis into BLAS layout once; blockProject then reads it once per
+    // basis vector rather than once per (i,p) as in scalar blockProject.
+    // Process all npoint in a single batch.
+    MultiRHSBlockProject<FineField> Projector;
+    Projector.Allocate(nbasis, grid, CoarseGrid());
+    Projector.ImportBasis(U.subspace);
+
+    std::vector<FineField>    phaV_batch(npoint, grid);
+    std::vector<FineField>    MphaV_batch(npoint, grid);
+    std::vector<CoarseVector> proj_batch(npoint, CoarseGrid());
+    std::vector<CoarseVector> ComputeProj(npoint, CoarseGrid());
+    std::vector<CoarseVector>          FT(npoint, CoarseGrid());
+
+    // Pre-allocate BLAS_F and BLAS_C to avoid repeated hipMalloc/hipFree of
+    // ~5.6 GB per blockProject call, which hangs on ROCm for large allocations.
+    Projector.BLAS_F.resize(Projector.fine_vol * Projector.words * npoint);
+    Projector.BLAS_C.resize(Projector.coarse_vol * nbasis * npoint);
+
+    for(int i=0;i<nbasis;i++){// Loop over basis vectors
+      accelerator_barrier(); // ensure prior iteration's async writes are retired
+      std::cout << GridLogMessage<< "CoarsenMatrixColoured vec "<<i<<"/"<<nbasis<< std::endl;
+
+      tphaseBZ-=usecond();
+      for(int p=0;p<npoint;p++)
+	phaV_batch[p] = phaF[p] * V.subspace[i];
+      tphaseBZ+=usecond();
+      std::cout << GridLogMessage<< "CoarsenMatrixColoured vec "<<i<<" phaseBZ done"<< std::endl;
+
+      tmat-=usecond();
+      for(int p=0;p<npoint;p++)
+	linop.Op(phaV_batch[p], MphaV_batch[p]);
+      tmat+=usecond();
+      std::cout << GridLogMessage<< "CoarsenMatrixColoured vec "<<i<<" mat done"<< std::endl;
+
+      // One batched GEMM reads BLAS_V once for all npoint vectors.
+      tproj-=usecond();
+      Projector.blockProject(MphaV_batch, proj_batch);
+      std::cout << GridLogMessage<< "CoarsenMatrixColoured vec "<<i<<" blockProject done"<< std::endl;
+      for(int p=0;p<npoint;p++)
+	ComputeProj[p] = conjugate(pha[p]) * proj_batch[p];
+      tproj+=usecond();
+      std::cout << GridLogMessage<< "CoarsenMatrixColoured vec "<<i<<" proj done"<< std::endl;
+
+      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, nbasis, {
+#ifdef GRID_SIMT
+	    int j = acceleratorSIMTlane(nbasis);
+	    A_v[sss](i,j) = FT_v[sss](j);
+#else
+	    // CPU build: acceleratorSIMTlane()==0 -- an un-looped SIMT tensor
+	    // index writes ONLY j=0 and silently drops the other nbasis-1
+	    // columns (caught by Test_schur_dense_coarse import certificate,
+	    // 2026-08-14).  Loop explicitly.
+	    for(int j=0;j<nbasis;j++){
+	      A_v[sss](i,j) = FT_v[sss](j);
+	    }
+#endif
+        });
+      }
+      tinv+=usecond();
+    }
+
+    // Only needed if nonhermitian
+    if ( ! hermitian ) {
+      std::cout << GridLogMessage<<"PopulateAdag  "<<std::endl;
+      PopulateAdag();
+    }
+
+    //    for(int p=0;p<geom.npoint;p++){
+    //      std::cout << " _A["<<p<<"] "<<norm2(_A[p])<<std::endl;
+    //    }
+
+    // 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 phaseBZ "<<tphaseBZ<<" 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;
+  }
+#endif  
+  void ExchangeCoarseLinks(void){
+    for(int p=0;p<geom.npoint;p++){
+      _A[p] = Cell.ExchangePeriodic(_A[p]);
+      if ( !hermitian ) _Adag[p]= Cell.ExchangePeriodic(_Adag[p]);
+    }
+  }
+  virtual  void Mdiag    (const Field &in, Field &out){ GRID_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/deprecated/GeneralCoarsenedMatrixMultiRHS.h b/Grid/algorithms/multigrid/deprecated/GeneralCoarsenedMatrixMultiRHS.h
new file mode 100644
index 000000000..292f6f035
--- /dev/null
+++ b/Grid/algorithms/multigrid/deprecated/GeneralCoarsenedMatrixMultiRHS.h
@@ -0,0 +1,738 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/GeneralCoarsenedMatrixMultiRHS.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);
+
+
+// Fine Object == (per site) type of fine field
+// nbasis      == number of deflation vectors
+template<class Fobj,class CComplex,int nbasis>
+class MultiGeneralCoarsenedMatrix : public SparseMatrixBase<Lattice<iVector<CComplex,nbasis > > >  {
+public:
+  typedef typename CComplex::scalar_object SComplex;
+  typedef GeneralCoarsenedMatrix<Fobj,CComplex,nbasis> GeneralCoarseOp;
+  typedef MultiGeneralCoarsenedMatrix<Fobj,CComplex,nbasis> MultiGeneralCoarseOp;
+
+  typedef iVector<CComplex,nbasis >           siteVector;
+  typedef iMatrix<CComplex,nbasis >           siteMatrix;
+  typedef iVector<SComplex,nbasis >           calcVector;
+  typedef iMatrix<SComplex,nbasis >           calcMatrix;
+  typedef Lattice<iScalar<CComplex> >         CoarseComplexField;
+  typedef Lattice<siteVector>                 CoarseVector;
+  typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
+  typedef iMatrix<CComplex,nbasis >  Cobj;
+  typedef iVector<CComplex,nbasis >  Cvec;
+  typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
+  typedef Lattice<Fobj >        FineField;
+  typedef Lattice<CComplex >    FineComplexField;
+  typedef CoarseVector Field;
+
+  ////////////////////
+  // Data members
+  ////////////////////
+  GridCartesian *       _CoarseGridMulti; 
+  NonLocalStencilGeometry geom;
+  NonLocalStencilGeometry geom_srhs;
+  PaddedCell Cell;
+  GeneralLocalStencil Stencil;
+
+  deviceVector<calcVector> BLAS_B;
+  deviceVector<calcVector> BLAS_C;
+  std::vector<deviceVector<calcMatrix> > BLAS_A;
+
+  std::vector<deviceVector<ComplexD *> > BLAS_AP;
+  std::vector<deviceVector<ComplexD *> > BLAS_BP;
+  deviceVector<ComplexD *>               BLAS_CP;
+
+  ///////////////////////
+  // Interface
+  ///////////////////////
+  GridBase      * Grid(void)           { return _CoarseGridMulti; };   // this is all the linalg routines need to know
+  GridCartesian * CoarseGrid(void)     { return _CoarseGridMulti; };   // this is all the linalg routines need to know
+
+  //////////////////////////////////////////////////////////////////////////
+  // Bilingual accessors, matching GeneralCoarsenedMatrix. Grid() here is the
+  // D+1 multiRHS grid and this class never holds the D dimensional one, so
+  // ExtractMatrix writes into whatever grid the caller's lattice is on.
+  //////////////////////////////////////////////////////////////////////////
+  NonLocalStencilGeometry & Geometry(void)      { return geom_srhs; };
+  void ExtractMatrix(int p,CoarseMatrix &A)     { BLAStoGrid(A,BLAS_A[p]); };
+
+  // I/O on the operator matrices, via the BLAS layout array. The parameter is
+  // a vector over the geometry points; the body indexes A[p].
+  void SetMatrix (int p,std::vector<CoarseMatrix> & A)
+  {
+    GRID_ASSERT(A.size()==geom_srhs.npoint);
+    GridtoBLAS(A[p],BLAS_A[p]);
+  }
+  void GetMatrix (int p,std::vector<CoarseMatrix> & A)
+  {
+    GRID_ASSERT(A.size()==geom_srhs.npoint);
+    BLAStoGrid(A[p],BLAS_A[p]);
+  }
+  void CopyMatrix (GeneralCoarseOp &_Op)
+  {
+    for(int p=0;p<geom.npoint;p++){
+      auto Aup = _Op.Cell.Extract(_Op._A[p]);
+      //Unpadded
+      GridtoBLAS(Aup,BLAS_A[p]);
+    }
+  }
+  /*
+  void CheckMatrix (GeneralCoarseOp &_Op)
+  {
+    std::cout <<"************* Checking the little direc operator mRHS"<<std::endl;
+    for(int p=0;p<geom.npoint;p++){
+      //Unpadded
+      auto Aup = _Op.Cell.Extract(_Op._A[p]);
+      auto Ack = Aup;
+      BLAStoGrid(Ack,BLAS_A[p]);
+      std::cout << p<<" Ack "<<norm2(Ack)<<std::endl;
+      std::cout << p<<" Aup "<<norm2(Aup)<<std::endl;
+    }
+    std::cout <<"************* "<<std::endl;
+  }
+  */
+  
+  MultiGeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridCartesian *CoarseGridMulti) :
+    _CoarseGridMulti(CoarseGridMulti),
+    geom_srhs(_geom),
+    geom(_CoarseGridMulti,_geom.hops,_geom.skip+1),
+    Cell(geom.Depth(),_CoarseGridMulti),
+    Stencil(Cell.grids.back(),geom.shifts) // padded cell stencil
+  {
+    int32_t padded_sites   = Cell.grids.back()->lSites();
+    int32_t unpadded_sites = CoarseGridMulti->lSites();
+    
+    int32_t nrhs  = CoarseGridMulti->FullDimensions()[0];  // # RHS
+    int32_t orhs  = nrhs/CComplex::Nsimd();
+
+    padded_sites   = padded_sites/nrhs;
+    unpadded_sites = unpadded_sites/nrhs;
+    
+    /////////////////////////////////////////////////
+    // Device data vector storage
+    /////////////////////////////////////////////////
+    BLAS_A.resize(geom.npoint);
+    for(int p=0;p<geom.npoint;p++){
+      BLAS_A[p].resize (unpadded_sites); // no ghost zone, npoint elements
+    }
+    
+    BLAS_B.resize(nrhs *padded_sites);   // includes ghost zone
+    BLAS_C.resize(nrhs *unpadded_sites); // no ghost zone
+    BLAS_AP.resize(geom.npoint);
+    BLAS_BP.resize(geom.npoint);
+    for(int p=0;p<geom.npoint;p++){
+      BLAS_AP[p].resize(unpadded_sites);
+      BLAS_BP[p].resize(unpadded_sites);
+    }
+    BLAS_CP.resize(unpadded_sites);
+
+    /////////////////////////////////////////////////
+    // Pointers to data
+    /////////////////////////////////////////////////
+
+    // Site identity mapping for A
+    for(int p=0;p<geom.npoint;p++){
+      for(int ss=0;ss<unpadded_sites;ss++){
+	ComplexD *ptr = (ComplexD *)&BLAS_A[p][ss];
+	acceleratorPut(BLAS_AP[p][ss],ptr);
+      }
+    }
+    // Site identity mapping for C
+    for(int ss=0;ss<unpadded_sites;ss++){
+      ComplexD *ptr = (ComplexD *)&BLAS_C[ss*nrhs];
+      acceleratorPut(BLAS_CP[ss],ptr);
+    }
+
+    // Neighbour table is more complicated
+    int32_t j=0; // Interior point counter (unpadded)
+    for(int32_t s=0;s<padded_sites;s++){ // 4 volume, padded
+      int ghost_zone=0;
+      for(int32_t point = 0 ; point < geom.npoint; point++){
+	int i=s*orhs*geom.npoint+point;
+	if( Stencil._entries[i]._wrap ) { // stencil is indexed by the oSite of the CoarseGridMulti, hence orhs factor
+	  ghost_zone=1; // If general stencil wrapped in any direction, wrap=1
+	}
+      }
+
+      if( ghost_zone==0) {
+	for(int32_t point = 0 ; point < geom.npoint; point++){
+	  int i=s*orhs*geom.npoint+point;
+ 	  int32_t nbr = Stencil._entries[i]._offset*CComplex::Nsimd(); // oSite -> lSite
+	  GRID_ASSERT(nbr<BLAS_B.size());
+	  ComplexD * ptr = (ComplexD *)&BLAS_B[nbr];
+	  acceleratorPut(BLAS_BP[point][j],ptr); // neighbour indexing in ghost zone volume
+	}
+	j++;
+      }
+    }
+    GRID_ASSERT(j==unpadded_sites);
+  }
+  template<class vobj> void GridtoBLAS(const Lattice<vobj> &from,deviceVector<typename vobj::scalar_object> &to)
+  {
+  typedef typename vobj::scalar_object sobj;
+  typedef typename vobj::scalar_type scalar_type;
+  typedef typename vobj::vector_type vector_type;
+
+  GridBase *Fg = from.Grid();
+  GRID_ASSERT(!Fg->_isCheckerBoarded);
+  int nd = Fg->_ndimension;
+
+  to.resize(Fg->lSites());
+
+  Coordinate LocalLatt = Fg->LocalDimensions();
+  size_t nsite = 1;
+  for(int i=0;i<nd;i++) nsite *= LocalLatt[i];
+
+  ////////////////////////////////////////////////////////////////////////////////////////////////
+  // do the index calc on the GPU
+  ////////////////////////////////////////////////////////////////////////////////////////////////
+  Coordinate f_ostride = Fg->_ostride;
+  Coordinate f_istride = Fg->_istride;
+  Coordinate f_rdimensions = Fg->_rdimensions;
+
+  autoView(from_v,from,AcceleratorRead);
+  auto to_v = &to[0];
+
+  const int words=sizeof(vobj)/sizeof(vector_type);
+  accelerator_for(idx,nsite,1,{
+      
+      Coordinate from_coor, base;
+      Lexicographic::CoorFromIndex(base,idx,LocalLatt);
+      for(int i=0;i<nd;i++){
+	from_coor[i] = base[i];
+      }
+      int from_oidx = 0; for(int d=0;d<nd;d++) from_oidx+=f_ostride[d]*(from_coor[d]%f_rdimensions[d]);
+      int from_lane = 0; for(int d=0;d<nd;d++) from_lane+=f_istride[d]*(from_coor[d]/f_rdimensions[d]);
+
+      const vector_type* from = (const vector_type *)&from_v[from_oidx];
+      scalar_type* to = (scalar_type *)&to_v[idx];
+      
+      scalar_type stmp;
+      for(int w=0;w<words;w++){
+	stmp = getlane(from[w], from_lane);
+	to[w] = stmp;
+      }
+    });
+  }    
+  template<class vobj> void BLAStoGrid(Lattice<vobj> &grid,deviceVector<typename vobj::scalar_object> &in)
+  {
+  typedef typename vobj::scalar_object sobj;
+  typedef typename vobj::scalar_type scalar_type;
+  typedef typename vobj::vector_type vector_type;
+
+  GridBase *Tg = grid.Grid();
+  GRID_ASSERT(!Tg->_isCheckerBoarded);
+  int nd = Tg->_ndimension;
+  
+  GRID_ASSERT(in.size()==Tg->lSites());
+
+  Coordinate LocalLatt = Tg->LocalDimensions();
+  size_t nsite = 1;
+  for(int i=0;i<nd;i++) nsite *= LocalLatt[i];
+
+  ////////////////////////////////////////////////////////////////////////////////////////////////
+  // do the index calc on the GPU
+  ////////////////////////////////////////////////////////////////////////////////////////////////
+  Coordinate t_ostride = Tg->_ostride;
+  Coordinate t_istride = Tg->_istride;
+  Coordinate t_rdimensions = Tg->_rdimensions;
+
+  autoView(to_v,grid,AcceleratorWrite);
+  auto from_v = &in[0];
+
+  const int words=sizeof(vobj)/sizeof(vector_type);
+  accelerator_for(idx,nsite,1,{
+      
+      Coordinate to_coor, base;
+      Lexicographic::CoorFromIndex(base,idx,LocalLatt);
+      for(int i=0;i<nd;i++){
+	to_coor[i] = base[i];
+      }
+      int to_oidx = 0; for(int d=0;d<nd;d++) to_oidx+=t_ostride[d]*(to_coor[d]%t_rdimensions[d]);
+      int to_lane = 0; for(int d=0;d<nd;d++) to_lane+=t_istride[d]*(to_coor[d]/t_rdimensions[d]);
+
+      vector_type* to = (vector_type *)&to_v[to_oidx];
+      scalar_type* from = (scalar_type *)&from_v[idx];
+      
+      scalar_type stmp;
+      for(int w=0;w<words;w++){
+	stmp=from[w];
+	putlane(to[w], stmp, to_lane);
+      }
+    });
+  }
+  void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
+		       Aggregation<Fobj,CComplex,nbasis> & Subspace,
+		       GridBase *CoarseGrid)
+  {
+#if 0
+    std::cout << GridLogMessage<< "GeneralCoarsenMatrixMrhs "<< std::endl;
+
+    GridBase *grid = Subspace.FineGrid;
+
+    /////////////////////////////////////////////////////////////
+    // Orthogonalise the subblocks over the basis
+    /////////////////////////////////////////////////////////////
+    CoarseScalar InnerProd(CoarseGrid); 
+    blockOrthogonalise(InnerProd,Subspace.subspace);
+
+    const int npoint = geom_srhs.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}
+       */
+    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_srhs.shifts[k][mu]*geom_srhs.shifts[l][mu];
+	}
+	phase=exp(phase*ci);
+	Mkl(k,l) = phase;
+      }
+    }
+    invMkl = Mkl.inverse();
+
+    ///////////////////////////////////////////////////////////////////////
+    // 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
+    std::vector<FineComplexField> phaF(npoint,grid);
+    std::vector<CoarseComplexField> pha(npoint,CoarseGrid);
+    
+    CoarseVector coarseInner(CoarseGrid);
+    
+    typedef typename CComplex::scalar_type SComplex;
+    FineComplexField one(grid); one=SComplex(1.0);
+    FineComplexField zz(grid); zz = Zero();
+    for(int p=0;p<npoint;p++){ // Loop over momenta in npoint
+      /////////////////////////////////////////////////////
+      // Stick a phase on every block
+      /////////////////////////////////////////////////////
+      CoarseComplexField coor(CoarseGrid);
+      pha[p]=Zero();
+      for(int mu=0;mu<Nd;mu++){
+	LatticeCoordinate(coor,mu);
+	RealD TwoPiL =  M_PI * 2.0/ clatt[mu];
+	pha[p] = pha[p] + (TwoPiL * geom_srhs.shifts[p][mu]) * coor;
+      }
+      pha[p]  =exp(pha[p]*ci);	
+
+      blockZAXPY(phaF[p],pha[p],one,zz);
+    }
+
+    // Could save on temporary storage here
+    std::vector<CoarseMatrix> _A;
+    _A.resize(geom_srhs.npoint,CoarseGrid);
+
+    std::vector<CoarseVector> ComputeProj(npoint,CoarseGrid);
+    CoarseVector          FT(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
+
+	phaV = phaF[p]*Subspace.subspace[i];
+
+	/////////////////////////////////////////////////////////////////////
+	// Multiple phased subspace vector by matrix and project to subspace
+	// Remove local bulk phase to leave relative phases
+	/////////////////////////////////////////////////////////////////////
+	linop.Op(phaV,MphaV);
+
+	// Fixme, could use batched block projector here
+	blockProject(coarseInner,MphaV,Subspace.subspace);
+
+	coarseInner = conjugate(pha[p]) * coarseInner;
+
+	ComputeProj[p] = coarseInner;
+      }
+
+      // Could do this with a block promote or similar BLAS call via the MultiRHSBlockProjector with a const matrix.
+      for(int k=0;k<npoint;k++){
+
+	FT = Zero();
+	for(int l=0;l<npoint;l++){
+	  FT= FT+ invMkl(l,k)*ComputeProj[l];
+	}
+      
+	int osites=CoarseGrid->oSites();
+	autoView( A_v  , _A[k], AcceleratorWrite);
+	autoView( FT_v  , FT, AcceleratorRead);
+	accelerator_for(sss, osites, 1, {
+	    for(int j=0;j<nbasis;j++){
+	      A_v[sss](i,j) = FT_v[sss](j);
+	    }
+        });
+      }
+    }
+
+    // Only needed if nonhermitian
+    //    if ( ! hermitian ) {
+    //      std::cout << GridLogMessage<<"PopulateAdag  "<<std::endl;
+    //      PopulateAdag();
+    //    }
+    // Need to write something to populate Adag from A
+
+    for(int p=0;p<geom_srhs.npoint;p++){
+      GridtoBLAS(_A[p],BLAS_A[p]);
+    }
+    /*
+Grid : Message : 11698.730546 s : CoarsenOperator eigen  1334 us
+Grid : Message : 11698.730563 s : CoarsenOperator phase  34729 us
+Grid : Message : 11698.730565 s : CoarsenOperator phaseBZ 2423814 us
+Grid : Message : 11698.730566 s : CoarsenOperator mat    127890998 us
+Grid : Message : 11698.730567 s : CoarsenOperator proj   515840840 us
+Grid : Message : 11698.730568 s : CoarsenOperator inv    103948313 us
+Takes 600s to compute matrix elements, DOMINATED by the block project.
+Easy to speed up with the batched block project.
+Store npoint vectors, get npoint x Nbasis block projection, and 81 fold faster.
+
+// Block project below taks to 240s
+Grid : Message : 328.193418 s : CoarsenOperator phase      38338 us
+Grid : Message : 328.193434 s : CoarsenOperator phaseBZ  1711226 us
+Grid : Message : 328.193436 s : CoarsenOperator mat    122213270 us
+//Grid : Message : 328.193438 s : CoarsenOperator proj   1181154 us <-- this is mistimed
+//Grid : Message : 11698.730568 s : CoarsenOperator inv  103948313 us <-- Cut this ~10x if lucky by loop fusion
+     */
+#else
+    RealD tproj=0.0;
+    RealD tmat=0.0;
+    RealD tphase=0.0;
+    RealD tphaseBZ=0.0;
+    RealD tinv=0.0;
+
+    std::cout << GridLogMessage<< "GeneralCoarsenMatrixMrhs "<< std::endl;
+
+    GridBase *grid = Subspace.FineGrid;
+
+    /////////////////////////////////////////////////////////////
+    // Orthogonalise the subblocks over the basis
+    /////////////////////////////////////////////////////////////
+    CoarseScalar InnerProd(CoarseGrid); 
+    blockOrthogonalise(InnerProd,Subspace.subspace);
+
+
+    MultiRHSBlockProject<Lattice<Fobj> >    Projector;
+    Projector.Allocate(nbasis,grid,CoarseGrid);
+    Projector.ImportBasis(Subspace.subspace);
+    
+    const int npoint = geom_srhs.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}
+       */
+    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_srhs.shifts[k][mu]*geom_srhs.shifts[l][mu];
+	}
+	phase=exp(phase*ci);
+	Mkl(k,l) = phase;
+      }
+    }
+    invMkl = Mkl.inverse();
+
+    ///////////////////////////////////////////////////////////////////////
+    // 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
+    std::vector<FineComplexField> phaF(npoint,grid);
+    std::vector<CoarseComplexField> pha(npoint,CoarseGrid);
+    
+    CoarseVector coarseInner(CoarseGrid);
+    
+    tphase=-usecond();
+    typedef typename CComplex::scalar_type SComplex;
+    FineComplexField one(grid); one=SComplex(1.0);
+    FineComplexField zz(grid); zz = Zero();
+    for(int p=0;p<npoint;p++){ // Loop over momenta in npoint
+      /////////////////////////////////////////////////////
+      // Stick a phase on every block
+      /////////////////////////////////////////////////////
+      CoarseComplexField coor(CoarseGrid);
+      pha[p]=Zero();
+      for(int mu=0;mu<Nd;mu++){
+	LatticeCoordinate(coor,mu);
+	RealD TwoPiL =  M_PI * 2.0/ clatt[mu];
+	pha[p] = pha[p] + (TwoPiL * geom_srhs.shifts[p][mu]) * coor;
+      }
+      pha[p]  =exp(pha[p]*ci);	
+
+      blockZAXPY(phaF[p],pha[p],one,zz);
+    }
+    tphase+=usecond();
+
+    // Could save on temporary storage here
+    std::vector<CoarseMatrix> _A;
+    _A.resize(geom_srhs.npoint,CoarseGrid);
+
+    // Count use small chunks than npoint == 81 and save memory
+    int batch = 9;
+    std::vector<FineField>    _MphaV(batch,grid);
+    std::vector<CoarseVector> TmpProj(batch,CoarseGrid);
+
+    std::vector<CoarseVector> ComputeProj(npoint,CoarseGrid);
+    CoarseVector          FT(CoarseGrid);
+    for(int i=0;i<nbasis;i++){// Loop over basis vectors
+      std::cout << GridLogMessage<< "CoarsenMatrixColoured vec "<<i<<"/"<<nbasis<< std::endl;
+
+      //      std::cout << GridLogMessage << " phasing the fine vector "<<std::endl;
+      // Fixme : do this in batches
+      for(int p=0;p<npoint;p+=batch){ // Loop over momenta in npoint
+
+	for(int b=0;b<MIN(batch,npoint-p);b++){
+	  tphaseBZ-=usecond();
+	  phaV = phaF[p+b]*Subspace.subspace[i];
+	  tphaseBZ+=usecond();
+
+	  /////////////////////////////////////////////////////////////////////
+	  // Multiple phased subspace vector by matrix and project to subspace
+	  // Remove local bulk phase to leave relative phases
+	  /////////////////////////////////////////////////////////////////////
+	  // Memory footprint was an issue
+	  tmat-=usecond();
+	  linop.Op(phaV,MphaV);
+	  _MphaV[b] = MphaV;
+	  tmat+=usecond();
+	}      
+
+	//	std::cout << GridLogMessage << " Calling block project "<<std::endl;
+	tproj-=usecond();
+	Projector.blockProject(_MphaV,TmpProj);
+	tproj+=usecond();
+	
+	//	std::cout << GridLogMessage << " conj phasing the coarse vectors "<<std::endl;
+	for(int b=0;b<MIN(batch,npoint-p);b++){
+	  ComputeProj[p+b] = conjugate(pha[p+b])*TmpProj[b];
+	}
+      }
+
+      // Could do this with a block promote or similar BLAS call via the MultiRHSBlockProjector with a const matrix.
+      
+      // std::cout << GridLogMessage << " Starting FT inv "<<std::endl;
+      tinv-=usecond();
+      for(int k=0;k<npoint;k++){
+	FT = Zero();
+	// 81 kernel calls as many ComputeProj vectors
+	// Could fuse with a vector of views, but ugly
+	// Could unroll the expression and run fewer kernels -- much more attractive
+	// Could also do non blocking.
+#if 0	
+	for(int l=0;l<npoint;l++){
+	  FT= FT+ invMkl(l,k)*ComputeProj[l];
+	}
+#else
+	const int radix = 9;
+	int ll;
+	for(ll=0;ll+radix-1<npoint;ll+=radix){
+	  // When ll = npoint-radix, ll+radix-1 = npoint-1, and we do it all.
+	  FT = FT 
+	    + invMkl(ll+0,k)*ComputeProj[ll+0]
+	    + invMkl(ll+1,k)*ComputeProj[ll+1]
+	    + invMkl(ll+2,k)*ComputeProj[ll+2]
+	    + invMkl(ll+3,k)*ComputeProj[ll+3]
+	    + invMkl(ll+4,k)*ComputeProj[ll+4]
+	    + invMkl(ll+5,k)*ComputeProj[ll+5]
+	    + invMkl(ll+6,k)*ComputeProj[ll+6]
+	    + invMkl(ll+7,k)*ComputeProj[ll+7]
+	    + invMkl(ll+8,k)*ComputeProj[ll+8];
+	}
+	for(int l=ll;l<npoint;l++){
+	  FT= FT+ invMkl(l,k)*ComputeProj[l];
+	}
+#endif
+      
+	// 1 kernel call -- must be cheaper
+	int osites=CoarseGrid->oSites();
+	autoView( A_v  , _A[k], AcceleratorWrite);
+	autoView( FT_v  , FT, AcceleratorRead);
+	accelerator_for(sss, osites, 1, {
+	    for(int j=0;j<nbasis;j++){
+	      A_v[sss](i,j) = FT_v[sss](j);
+	    }
+        });
+      }
+      tinv+=usecond();
+    }
+
+    // Only needed if nonhermitian
+    //    if ( ! hermitian ) {
+    //      std::cout << GridLogMessage<<"PopulateAdag  "<<std::endl;
+    //      PopulateAdag();
+    //    }
+    // Need to write something to populate Adag from A
+    //    std::cout << GridLogMessage << " Calling GridtoBLAS "<<std::endl;
+    for(int p=0;p<geom_srhs.npoint;p++){
+      GridtoBLAS(_A[p],BLAS_A[p]);
+    }
+    std::cout << GridLogMessage<<"CoarsenOperator phase  "<<tphase<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator phaseBZ "<<tphaseBZ<<" 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;
+#endif
+  }
+  void Mdag(const CoarseVector &in, CoarseVector &out)
+  {
+    this->M(in,out);
+  }
+  void M (const CoarseVector &in, CoarseVector &out)
+  {
+    //    std::cout << GridLogMessage << "New Mrhs coarse"<<std::endl;
+    conformable(CoarseGrid(),in.Grid());
+    conformable(in.Grid(),out.Grid());
+    out.Checkerboard() = in.Checkerboard();
+
+    RealD t_tot;
+    RealD t_exch;
+    RealD t_GtoB;
+    RealD t_BtoG;
+    RealD t_mult;
+
+    t_tot=-usecond();
+    CoarseVector tin=in;
+    t_exch=-usecond();
+    CoarseVector pin = Cell.ExchangePeriodic(tin); //padded input
+    t_exch+=usecond();
+
+    CoarseVector pout(pin.Grid());
+
+    int npoint = geom.npoint;
+    typedef calcMatrix* Aview;
+    typedef LatticeView<Cvec> Vview;
+      
+    const int Nsimd = CComplex::Nsimd();
+
+    int64_t nrhs  =pin.Grid()->GlobalDimensions()[0];
+    GRID_ASSERT(nrhs>=1);
+
+    RealD flops,bytes;
+    int64_t osites=in.Grid()->oSites(); // unpadded
+    int64_t unpadded_vol = CoarseGrid()->lSites()/nrhs;
+    
+    flops = 1.0* npoint * nbasis * nbasis * 8.0 * osites * CComplex::Nsimd();
+    bytes = 1.0*osites*sizeof(siteMatrix)*npoint/pin.Grid()->GlobalDimensions()[0]
+          + 2.0*osites*sizeof(siteVector)*npoint;
+    
+
+    t_GtoB=-usecond();
+    GridtoBLAS(pin,BLAS_B);
+    t_GtoB+=usecond();
+
+    GridBLAS BLAS;
+
+    t_mult=-usecond();
+    for(int p=0;p<geom.npoint;p++){
+      RealD c = 1.0;
+      if (p==0) c = 0.0;
+      ComplexD beta(c);
+
+      BLAS.gemmBatched(nbasis,nrhs,nbasis,
+		       ComplexD(1.0),
+		       BLAS_AP[p], 
+		       BLAS_BP[p], 
+		       ComplexD(c), 
+		       BLAS_CP);
+    }
+    BLAS.synchronise();
+    t_mult+=usecond();
+
+    t_BtoG=-usecond();
+    BLAStoGrid(out,BLAS_C);
+    t_BtoG+=usecond();
+    t_tot+=usecond();
+    /*
+    std::cout << GridLogMessage << "New Mrhs coarse DONE "<<std::endl;
+    std::cout << GridLogMessage<<"Coarse Mult exch "<<t_exch<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"Coarse Mult mult "<<t_mult<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"Coarse Mult GtoB  "<<t_GtoB<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"Coarse Mult BtoG  "<<t_BtoG<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"Coarse Mult tot  "<<t_tot<<" us"<<std::endl;
+    */
+    //    std::cout << GridLogMessage<<std::endl;
+    //    std::cout << GridLogMessage<<"Coarse Kernel flops "<< flops<<std::endl;
+    //    std::cout << GridLogMessage<<"Coarse Kernel flop/s "<< flops/t_mult<<" mflop/s"<<std::endl;
+    //    std::cout << GridLogMessage<<"Coarse Kernel bytes/s "<< bytes/t_mult/1000<<" GB/s"<<std::endl;
+    //    std::cout << GridLogMessage<<"Coarse overall flops/s "<< flops/t_tot<<" mflop/s"<<std::endl;
+    //    std::cout << GridLogMessage<<"Coarse total bytes   "<< bytes/1e6<<" MB"<<std::endl;
+  };
+  virtual  void Mdiag    (const Field &in, Field &out){ GRID_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/deprecated/TwoLevelADEF2.h b/Grid/algorithms/multigrid/deprecated/TwoLevelADEF2.h
new file mode 100644
index 000000000..52ff3e851
--- /dev/null
+++ b/Grid/algorithms/multigrid/deprecated/TwoLevelADEF2.h
@@ -0,0 +1,136 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./Grid/algorithms/multigrid/deprecated/TwoLevelADEF2.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <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
+
+// DEPRECATED with the V1 coarse operators: the single-RHS ADEF-2 on an
+// Aggregation (ProjectToSubspace/PromoteFromSubspace) and a D-dimensional
+// coarse operator.  The mrhs solver in AdefMrhs.h covers one right-hand
+// side through its LinearFunction interface, on the V2 coarse operator.
+
+NAMESPACE_BEGIN(Grid);
+
+template<class Field, class CoarseField, class Aggregation>
+class TwoLevelADEF2 : public TwoLevelCG<Field>
+{
+ public:
+  ///////////////////////////////////////////////////////////////////////////////////
+  // Need something that knows how to get from Coarse to fine and back again
+  //  void ProjectToSubspace(CoarseVector &CoarseVec,const FineField &FineVec){
+  //  void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){
+  ///////////////////////////////////////////////////////////////////////////////////
+  GridBase *coarsegrid;
+  Aggregation &_Aggregates;                    
+  LinearFunction<CoarseField> &_CoarseSolver;
+  LinearFunction<CoarseField> &_CoarseSolverPrecise;
+  ///////////////////////////////////////////////////////////////////////////////////
+  
+  // more most opertor functions
+  TwoLevelADEF2(RealD tol,
+		Integer maxit,
+		LinearOperatorBase<Field>    &FineLinop,
+		LinearFunction<Field>        &Smoother,
+		LinearFunction<CoarseField>  &CoarseSolver,
+		LinearFunction<CoarseField>  &CoarseSolverPrecise,
+		Aggregation &Aggregates
+		) :
+      TwoLevelCG<Field>(tol,maxit,FineLinop,Smoother,Aggregates.FineGrid),
+      _CoarseSolver(CoarseSolver),
+      _CoarseSolverPrecise(CoarseSolverPrecise),
+      _Aggregates(Aggregates)
+  {
+    coarsegrid = Aggregates.CoarseGrid;
+  };
+
+  virtual void PcgM1(Field & in, Field & out)
+  {
+    GRID_TRACE("MultiGridPreconditioner ");
+    // [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
+
+    Field tmp(this->grid);
+    Field Min(this->grid);
+    CoarseField PleftProj(this->coarsegrid);
+    CoarseField PleftMss_proj(this->coarsegrid);
+
+    this->SmoothTimer.Start();
+    this->_Smoother(in,Min);
+    this->SmoothTimer.Stop();
+    this->SmoothCalls++;
+
+    this->MatrixTimer.Start();
+    this->_FineLinop.HermOp(Min,out);
+    this->MatrixTimer.Stop();
+    this->MatrixCalls++;
+    axpy(tmp,-1.0,out,in);          // tmp  = in - A Min
+
+    this->ProjectTimer.Start();
+    this->_Aggregates.ProjectToSubspace(PleftProj,tmp);
+    this->ProjectTimer.Stop();
+    this->ProjectCalls++;
+    this->CoarseTimer.Start();
+    this->_CoarseSolver(PleftProj,PleftMss_proj); // Ass^{-1} [in - A Min]_s
+    this->CoarseTimer.Stop();
+    this->CoarseCalls++;
+    this->PromoteTimer.Start();
+    this->_Aggregates.PromoteFromSubspace(PleftMss_proj,tmp);// tmp = Q[in - A Min]
+    this->PromoteTimer.Stop();
+    this->PromoteCalls++;
+
+    axpy(out,1.0,Min,tmp); // Min+tmp
+  }
+
+  virtual void Vstart(Field & x,const Field & src)
+  {
+    std::cout << GridLogMessage<<"HDCG: fPcg Vstart "<<std::endl;
+    ///////////////////////////////////
+    // Choose x_0 such that 
+    // x_0 = guess +  (A_ss^inv) r_s = guess + Ass_inv [src -Aguess]
+    //                               = [1 - Ass_inv A] Guess + Assinv src
+    //                               = P^T guess + Assinv src 
+    //                               = Vstart  [Tang notation]
+    // This gives:
+    // W^T (src - A x_0) = src_s - A guess_s - r_s
+    //                   = src_s - (A guess)_s - src_s  + (A guess)_s 
+    //                   = 0 
+    ///////////////////////////////////
+    Field r(this->grid);
+    Field mmp(this->grid);
+    CoarseField PleftProj(this->coarsegrid);
+    CoarseField PleftMss_proj(this->coarsegrid);
+
+    std::cout << GridLogMessage<<"HDCG: fPcg Vstart projecting "<<std::endl;
+    this->_Aggregates.ProjectToSubspace(PleftProj,src);     
+    std::cout << GridLogMessage<<"HDCG: fPcg Vstart coarse solve "<<std::endl;
+    this->_CoarseSolverPrecise(PleftProj,PleftMss_proj); // Ass^{-1} r_s
+    std::cout << GridLogMessage<<"HDCG: fPcg Vstart promote "<<std::endl;
+    this->_Aggregates.PromoteFromSubspace(PleftMss_proj,x);  
+
+  }
+
+};
+
+NAMESPACE_END(Grid);