Clean organisation

Grid/algorithms/multigrid/HDCGMultiGrid.h

@@ -0,0 +1,632 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/HDCGMultiGrid.h
+
+    Copyright (C) 2026
+
+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.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h>
+#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h>
+#include <Grid/algorithms/multigrid/PVdagMMultiGrid.h>
+#include <Grid/algorithms/multigrid/HDCGMultiGridParams.h>
+
+NAMESPACE_BEGIN(Grid);
+
+//////////////////////////////////////////////////////////////////////////////////////
+// The two-level mrhs HDCG (arXiv:1402.2585, mrhs form arXiv:2409.03904) as
+// objects, the shape of PVdagMMultiGrid.h:
+//
+//   MGCoarseGrids / MGFineGridsF   shared with the PVdagM chain
+//   HDCGCoarsening                 everything that is a function of the gauge
+//                                  field: the raw near-null basis (and its
+//                                  refinement), the transfer operators, the
+//                                  Galerkin coarse operator (V2, Hermitian),
+//                                  the coarse eigenvectors for deflation
+//   HDCGSolver                     the solve chain on a borrowed coarsening:
+//                                  smoother, deflated coarse solve, ADEF-2
+//                                  preconditioner (fp64 or fp32 behind the
+//                                  seam), outer mrhs fPcg / BlockCGrQ
+//
+// The Hermitian operator is the Schur-preconditioned M^dag M on one
+// checkerboard; the near-null basis and the fine fields live on the
+// checkerboarded grid.  The coarse space is unvectorised and ALWAYS the D+1
+// multiRHS field; Nrhs 1 is the extent-1 case.  Grids outlive the
+// coarsening outlives the solver.
+//////////////////////////////////////////////////////////////////////////////////////
+template<class Fobj,class CComplex,int nbasis>
+class HDCGCoarsening {
+public:
+  typedef Lattice<Fobj>                                          FineField;
+  typedef MultiGeneralCoarsenedOperatorV2<Fobj,CComplex,nbasis>  CoarseOperator;
+  typedef typename CoarseOperator::CoarseVector                  CoarseVector;
+  typedef typename GridTypeMapper<Fobj>::SinglePrecision         FobjF;
+  typedef Lattice<FobjF>                                         FineFieldF;
+
+  MGCoarseGrids                    &Grids;     // borrowed
+  MGFineGridsF                     &GridsF;    // borrowed
+  GridBase                         *FrbGrid;   // borrowed: the checkerboarded fine grid
+  int                               cb;
+  HDCGMultiGridParams               Params;
+  NextToNextToNextToNearestStencilGeometry5D geom;   // the coarse stencil of the Schur operator
+  CoarseOperator                    CoarseOp;
+  //////////////////////////////////////////////////////////////////////
+  // ONE fine transfer operator.  Its STORE follows the coarse sector --
+  // the coarse space is what it feeds -- while its import and export
+  // accept EITHER fine precision, because they are already a layout
+  // transformation and a scalar conversion inside one costs nothing.
+  //////////////////////////////////////////////////////////////////////
+  static const bool CoarseIsSingle =
+    ( sizeof(typename GridTypeMapper<CComplex>::scalar_type) == sizeof(ComplexF) );
+  typedef typename std::conditional<CoarseIsSingle,FineFieldF,FineField>::type ProjectorField;
+  typedef MultiRHSBlockProject<ProjectorField>                        ProjectorL1_t;
+
+  ProjectorL1_t                     MrhsProjectorL1;
+  std::vector<FineField>            rawNull;         // RAW fine near-null basis
+  std::vector<CoarseVector>         evec;            // coarse eigenvectors, on Coarse5d
+  std::vector<RealD>                eval;
+  MultiRHSDeflation<CoarseVector>   Deflator;        // nev==0 until CoarseLanczos
+  GridCartesian                    *CMrhs;           // transient solve grid, owned
+  int                               nrhs;
+
+  HDCGCoarsening(MGCoarseGrids &_Grids, MGFineGridsF &_GridsF, GridBase *_FrbGrid, int _cb,
+                 const HDCGMultiGridParams &P)
+    : Grids(_Grids), GridsF(_GridsF), FrbGrid(_FrbGrid), cb(_cb), Params(P),
+      geom(_Grids.Coarse5d),
+      CoarseOp(geom,_Grids.Coarse5d),
+      CMrhs(nullptr),
+      nrhs(-1)
+  {
+    Deflator.Deallocate();
+    // The coarse stencil is a max-norm-1 box (Geometry.h) and the Schur
+    // operator M^dag M reaches four fine hops, so a block shorter than 4
+    // leaves couplings outside the stencil: the Fourier extraction then
+    // folds them into the retained points and the Galerkin certificate
+    // reads it.  Warn here; the certificate is the guard.
+    for(int d=0;d<4;d++){
+      if ( Params.Setup.Block1[d] < 4 )
+        std::cout << GridLogWarning << "HDCGCoarsening: Block1[" << d << "] = " << Params.Setup.Block1[d]
+                  << " < 4: the Schur operator's range exceeds the coarse stencil; "
+                  << "expect a Galerkin defect" << std::endl;
+    }
+  };
+  ~HDCGCoarsening()
+  {
+    CoarseOp.ReleaseGrid();       // borrowers let go before their grids die
+    if ( CMrhs ) delete CMrhs;
+  }
+
+  int FileExists(const std::string &fname)
+  {
+    uint64_t exists=0;
+    if ( fname.length() ) {
+      if ( Grids.FGrid->IsBoss() ){ std::ifstream f(fname); exists=f.good()?1:0; }
+      Grids.FGrid->GlobalSum(exists);
+    }
+    return (int)exists;
+  }
+  FineField FineNoise(GridParallelRNG &RNG)
+  {
+    FineField v(FrbGrid); v.Checkerboard()=cb; gaussian(RNG,v);
+    return v;
+  }
+  // Chebyshev upper bound: the parameter, or 1.1 x a power-method estimate
+  RealD UpperBound(RealD param, LinearOperatorBase<FineField> &HermOp, GridParallelRNG &RNG)
+  {
+    if ( param > 0.0 ) return param;
+    FineField src = FineNoise(RNG);
+    PowerMethod<FineField> PM;
+    RealD lam = PM(HermOp,src);
+    std::cout << GridLogMessage << "HDCGCoarsening: fine power method " << lam
+              << ", Chebyshev upper bound " << 1.1*lam << std::endl;
+    return 1.1*lam;
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // The RAW basis: a refined basis on disk wins, then a raw basis on
+  // disk, else built by the chosen method (and saved if named).  CG
+  // refinement happens here; HDCG refinement needs a solver, see
+  // HDCGRefineSubspace below.  The Aggregation is scaffolding for the
+  // Chebyshev filters only: they run entirely on the fine grid.
+  ////////////////////////////////////////////////////////////////////
+  void GetSubspace(GridParallelRNG &RNG, LinearOperatorBase<FineField> &HermOp)
+  {
+    const HDCGSubspaceParams &S = Params.Subspace;
+    rawNull.clear();
+    rawNull.reserve(nbasis);
+    for(int k=0;k<nbasis;k++){ rawNull.push_back(FineField(FrbGrid)); rawNull[k].Checkerboard()=cb; }
+
+    if ( FileExists(S.RefinedSubspaceFile) ) {
+      std::cout << GridLogMessage << "HDCGCoarsening: loading REFINED subspace " << S.RefinedSubspaceFile << std::endl;
+      loadSubspace(rawNull, S.RefinedSubspaceFile, S.FileControl);
+      return;
+    }
+    if ( FileExists(Params.Setup.SubspaceFile) ) {
+      std::cout << GridLogMessage << "HDCGCoarsening: loading subspace " << Params.Setup.SubspaceFile << std::endl;
+      loadSubspace(rawNull, Params.Setup.SubspaceFile, S.FileControl);
+    } else {
+      Aggregation<Fobj,CComplex,nbasis> Agg(Grids.Coarse5d,FrbGrid,cb);
+      if ( S.Method == HDCGSubspaceMethod::ChebyshevNew ) {
+        RealD hi = UpperBound(S.ChebyHi,HermOp,RNG);
+        std::cout << GridLogMessage << "HDCGCoarsening: CreateSubspaceChebyshevNew hi " << hi << std::endl;
+        Agg.CreateSubspaceChebyshevNew(RNG,HermOp,hi);
+      } else if ( S.Method == HDCGSubspaceMethod::Chebyshev ) {
+        RealD hi = UpperBound(S.ChebyHi,HermOp,RNG);
+        std::cout << GridLogMessage << "HDCGCoarsening: CreateSubspaceChebyshev hi " << hi
+                  << " lo " << S.ChebyLo << " order " << S.ChebyOrder << std::endl;
+        Agg.CreateSubspaceChebyshev(RNG,HermOp,nbasis,hi,S.ChebyLo,S.ChebyOrder);
+      } else {
+        std::cout << GridLogMessage << "HDCGCoarsening: subspace from fine eigenvectors " << S.FineEvecFile
+                  << " summed " << S.FineEvecSample << " per vector" << std::endl;
+        LoadFineEvecSum(Agg.subspace,S.FineEvecFile,S.FineEvecSample,S.FileControl);
+      }
+      for(int k=0;k<nbasis;k++) rawNull[k]=Agg.subspace[k];
+      if ( Params.Setup.SubspaceFile.length() )
+        saveSubspace(rawNull, Params.Setup.SubspaceFile, S.FileControl);
+    }
+    if ( S.Refine == HDCGRefineMethod::CG ) {
+      std::cout << GridLogMessage << "HDCGCoarsening: CG refinement shift " << S.RefineShift
+                << " tol " << S.RefineTol << " maxit " << S.RefineMaxIt << std::endl;
+      Aggregation<Fobj,CComplex,nbasis> Agg(Grids.Coarse5d,FrbGrid,cb);
+      for(int k=0;k<nbasis;k++) Agg.subspace[k]=rawNull[k];
+      Agg.RefineSubspace(HermOp,S.RefineShift,S.RefineTol,S.RefineMaxIt);
+      for(int k=0;k<nbasis;k++) rawNull[k]=Agg.subspace[k];
+      if ( S.RefinedSubspaceFile.length() )
+        saveSubspace(rawNull, S.RefinedSubspaceFile, S.FileControl);
+    }
+  }
+  // Basis vector b = precisionChange of the sum of `sample` consecutive fp32 records
+  void LoadFineEvecSum(std::vector<FineField> &sub, const std::string &file, int sample, int control)
+  {
+#ifdef HAVE_LIME
+    emptyUserRecord record;
+    FineFieldF tmp(GridsF.FrbGridF), sum(GridsF.FrbGridF);
+    tmp.Checkerboard()=cb; sum.Checkerboard()=cb;
+    ScidacReader RD;
+    RD.open(file);
+    for(int b=0;b<nbasis;b++){
+      sum = Zero();
+      for(int n=0;n<sample;n++){
+        RD.readScidacFieldRecord(tmp,record,control);
+        sum = sum + tmp;
+      }
+      precisionChange(sub[b],sum);
+    }
+    RD.close();
+#else
+    GRID_ASSERT(0);
+#endif
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // HDCG refinement of the raw basis (the paper's loop): blocks of
+  // nrhs vectors, normalised, solved with the refining solver, and
+  // renormalised.  The solver was built on the SHIFTED operator and
+  // the coarsening of it; the caller re-coarsens afterwards.
+  ////////////////////////////////////////////////////////////////////
+  void RefineSubspaceHDCG(LinearOperatorBase<FineField> &HermOp, TwoLevelCGmrhs<FineField> &Refiner, int nr)
+  {
+    std::vector<FineField> src(nr,FrbGrid), res(nr,FrbGrid);
+    FineField tmp(FrbGrid); tmp.Checkerboard()=cb;
+    for(int r=0;r<nr;r++){ src[r].Checkerboard()=cb; res[r].Checkerboard()=cb; }
+    for(int b=0;b<nbasis;b+=nr){
+      int nb = std::min(nbasis-b,nr);
+      for(int r=0;r<nb;r++){
+        RealD scale = std::pow(norm2(rawNull[b+r]),-0.5);
+        src[r] = rawNull[b+r]*scale;
+      }
+      for(int r=nb;r<nr;r++) src[r] = src[0];   // pad the block; results discarded
+      for(int r=0;r<nr;r++)  res[r] = Zero();
+      HermOp.HermOp(src[0],tmp);
+      std::cout << GridLogMessage << "HDCGCoarsening: refine block " << b << " before <n|A|n> " << norm2(tmp) << std::endl;
+      Refiner(src,res);
+      for(int r=0;r<nb;r++){
+        RealD scale = std::pow(norm2(res[r]),-0.5);
+        rawNull[b+r] = res[r]*scale;
+      }
+      HermOp.HermOp(rawNull[b],tmp);
+      std::cout << GridLogMessage << "HDCGCoarsening: refine block " << b << " after  <n|A|n> " << norm2(tmp) << std::endl;
+    }
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // The Galerkin coarse operator and transfer operators from the RAW
+  // basis.  CoarsenOperator block-orthonormalises its input in place,
+  // so a working copy is taken and the RAW vectors survive.  The
+  // operator handed in must apply the Hermitian operator through Op
+  // (HermOpAdaptor).
+  //
+  // The RAW basis has done its work once this returns, so it is freed --
+  // it is nbasis fine vectors, the single largest resident block of the
+  // setup.  retain_basis keeps it, for the two callers that coarsen more
+  // than once: the HDCG subspace refinement below, and a fixed-basis
+  // rebuild on a changed gauge field (HMC).  Setup.RetainSubspace keeps
+  // it for the whole life of the object.
+  ////////////////////////////////////////////////////////////////////
+  void Coarsen(LinearOperatorBase<FineField> &HermOp, int retain_basis=0)
+  {
+    GRID_ASSERT( rawNull.size() == nbasis );
+    std::vector<FineField> sub(rawNull);
+
+    CoarseOp.SetGrid(Grids.CoarseBatch);
+    std::cout << GridLogMessage << "HDCGCoarsening: CoarsenOperator, batch " << Grids.batch << std::endl;
+    // CoarsenOperator orthonormalises sub in place and imports it into OUR
+    // transfer operator, which it then uses.  One basis store, imported once.
+    // The projector's grid carries the STORE's SIMD layout: the fp32 fine
+    // red-black grid when the coarse sector is fp32.
+    MrhsProjectorL1.Allocate(nbasis,
+                             CoarseIsSingle ? (GridBase *)GridsF.FrbGridF : (GridBase *)FrbGrid,
+                             Grids.Coarse5d);
+    CoarseOp.CoarsenOperator(HermOp,sub,Grids.Coarse5d,Grids.batch,MrhsProjectorL1);
+    if ( !retain_basis && !Params.Setup.RetainSubspace ) DiscardBasis();
+    MrhsProjectorL1.ReleaseScratch();
+    nrhs = -1;                     // the batch grid is set; SetNrhs retargets
+    ReportDeviceFootprint("after Coarsen");
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // Retarget the coarse operator to Nrhs: the D+1 grid is owned here
+  // and every borrower re-captures it after this call.
+  ////////////////////////////////////////////////////////////////////
+  void SetNrhs(int nr)
+  {
+    if ( nr == nrhs ) return;
+    CoarseOp.ReleaseGrid();
+    if ( CMrhs ) { delete CMrhs; CMrhs=nullptr; }
+    Coordinate cml({nr,1,Grids.clatt[0],Grids.clatt[1],Grids.clatt[2],Grids.clatt[3]});
+    CMrhs = new GridCartesian(cml,Grids.cmsimd,Grids.cmmpi);
+    CoarseOp.SetGrid(CMrhs);
+    nrhs = nr;
+    std::cout << GridLogMessage << "HDCGCoarsening: coarse operator at Nrhs " << nr << std::endl;
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // Coarse eigenvectors for the deflated coarse solve.  Block Lanczos
+  // runs at the solve's Nrhs on the D+1 operator; plain Lanczos at
+  // Nrhs 1 (the extent-1 D+1 grid), its vectors then sliced to the D
+  // grid so both feed the deflator the same way.  Loaded from EvecFile
+  // when present, saved there when named.
+  ////////////////////////////////////////////////////////////////////
+  void CoarseLanczos(int nr)
+  {
+    const HDCGLanczosParams &L = Params.Lanczos;
+    evec.clear(); eval.clear();
+    Deflator.Deallocate();
+    if ( L.Method == HDCGLanczosMethod::None ) {
+      std::cout << GridLogMessage << "HDCGCoarsening: no coarse deflation" << std::endl;
+      return;
+    }
+    if ( FileExists(L.EvecFile) ) {
+      std::cout << GridLogMessage << "HDCGCoarsening: loading coarse eigenpairs " << L.EvecFile << std::endl;
+      loadEigenpairs(evec,eval,L.EvecFile,Grids.Coarse5d,Params.Subspace.FileControl);
+    } else {
+      int Nconv=0;
+      SetNrhs( L.Method == HDCGLanczosMethod::Block ? nr : 1 );
+      HermitianLinearOperator<CoarseOperator,CoarseVector> HermOpC(CoarseOp);
+      GridParallelRNG cRNG(CMrhs); cRNG.SeedFixedIntegers(std::vector<int>({31,32,33,34}));
+      RealD hi = L.ChebyHi;
+      if ( hi <= 0.0 ) {
+        CoarseVector s(CMrhs); random(cRNG,s);
+        PowerMethod<CoarseVector> PM;
+        hi = 1.1*PM(HermOpC,s);
+        std::cout << GridLogMessage << "HDCGCoarsening: coarse Chebyshev upper bound " << hi << std::endl;
+      }
+      Chebyshev<CoarseVector> Cheby(L.ChebyLo,hi,L.ChebyOrder);
+
+      if ( L.Method == HDCGLanczosMethod::Block ) {
+        std::cout << GridLogMessage << "HDCGCoarsening: block Lanczos, Nrhs " << nr
+                  << " Nstop " << L.Nstop << " Nk " << L.Nk << " Nm " << L.Nm << std::endl;
+        ImplicitlyRestartedBlockLanczosCoarse<CoarseVector>
+          IRL(HermOpC,Grids.Coarse5d,CMrhs,nr,Cheby,L.Nstop,1,nr,L.Nk,L.Nm,L.Tol,L.MaxIt);
+        evec.resize(L.Nm,Grids.Coarse5d);
+        eval.resize(L.Nm);
+        std::vector<CoarseVector> src(nr,Grids.Coarse5d);
+        GridParallelRNG dRNG(Grids.Coarse5d); dRNG.SeedFixedIntegers(std::vector<int>({41,42,43,44}));
+        for(int r=0;r<nr;r++) random(dRNG,src[r]);
+        IRL.calc(eval,evec,src,Nconv,LanczosType::irbl);
+        evec.resize(Nconv,Grids.Coarse5d);                 // only the converged pairs
+        eval.resize(Nconv);
+      } else {
+        std::cout << GridLogMessage << "HDCGCoarsening: plain Lanczos at Nrhs 1, Nstop " << L.Nstop
+                  << " Nk " << L.Nk << " Nm " << L.Nm << std::endl;
+        PlainHermOp<CoarseVector>    Op(HermOpC);
+        FunctionHermOp<CoarseVector> OpCheby(Cheby,HermOpC);
+        ImplicitlyRestartedLanczos<CoarseVector> IRL(OpCheby,Op,L.Nstop,L.Nk,L.Nm,L.Tol,L.MaxIt);
+        std::vector<CoarseVector> evec1(L.Nm,CMrhs);
+        std::vector<RealD>        eval1(L.Nm);
+        CoarseVector src(CMrhs); random(cRNG,src);
+        IRL.calc(eval1,evec1,src,Nconv);
+        evec.resize(Nconv,Grids.Coarse5d);
+        eval.resize(Nconv);
+        for(int e=0;e<Nconv;e++){ ExtractSliceFast(evec[e],evec1[e],0,0); eval[e]=eval1[e]; }
+      }
+      // The Lanczos routines set their outputs only on reaching Nstop;
+      // a partial result carries no usable eigenpairs.
+      if ( Nconv < L.Nstop ) {
+        std::cout << GridLogWarning << "HDCGCoarsening: coarse Lanczos converged " << Nconv << " < Nstop " << L.Nstop
+                  << " modes (Chebyshev window [" << L.ChebyLo << "," << hi << "] must hold at least Nstop eigenvalues"
+                  << " and Nm must resolve them); coarse solve undeflated" << std::endl;
+        evec.clear(); eval.clear();
+        return;
+      }
+      if ( L.EvecFile.length() ) saveEigenpairs(evec,eval,L.EvecFile,Params.Subspace.FileControl);
+    }
+    const int nev = evec.size();
+    Deflator.ImportEigenBasis(evec,eval);
+    // The Lattice-layout copies have been imported into the BLAS store and
+    // saved if they were going to be; the solve reads only the BLAS store.
+    evec.clear(); evec.shrink_to_fit();
+    std::cout << GridLogMessage << "HDCGCoarsening: " << nev << " coarse eigenvectors imported for deflation"
+              << (nev ? " (lowest "+std::to_string(eval[0])+", highest "+std::to_string(eval[nev-1])+")" : "")
+              << std::endl;
+    ReportDeviceFootprint("after CoarseLanczos");
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // Certificates at Nrhs 1 through the solve's own transfer operators:
+  // Galerkin ||A_c x - P^dag A P x|| / ||A_c x||  (the rounding level of
+  // the coarsening: fp64 coarse ~1e-13, fp32 ~1e-6), and Hermiticity
+  // |<x,A_c y> - <A_c x,y>| / |<x,A_c y>|, which CG relies on.
+  ////////////////////////////////////////////////////////////////////
+  RealD CertifyCoarsening(LinearOperatorBase<FineField> &HermOp)
+  {
+    SetNrhs(1);
+    CoarseVector xc(CMrhs), yc(CMrhs), Acx(CMrhs), Acy(CMrhs), PtAPx(CMrhs);
+    GridParallelRNG cRNG(CMrhs); cRNG.SeedFixedIntegers(std::vector<int>({21,22,23,24}));
+    random(cRNG,xc);
+    random(cRNG,yc);
+
+    CoarseOp.M(xc,Acx);                                     // A_c x
+    CoarseOp.M(yc,Acy);                                     // A_c y
+
+    std::vector<FineField> Px(1,FrbGrid), APx(1,FrbGrid);
+    Px[0].Checkerboard()=cb; APx[0].Checkerboard()=cb;
+    MrhsProjectorL1.blockPromote(Px,xc);                      // P x
+    HermOp.HermOp(Px[0],APx[0]);                            // A P x
+    MrhsProjectorL1.blockProject(APx,PtAPx);                  // P^dag A P x
+
+    CoarseVector d(CMrhs); d = Acx - PtAPx;
+    RealD rel = std::sqrt(norm2(d)/norm2(Acx));
+    ComplexD xAy = TensorRemove(innerProduct(xc,Acy));
+    ComplexD Axy = TensorRemove(innerProduct(Acx,yc));
+    RealD herm = abs(xAy-Axy)/abs(xAy);
+    std::cout << GridLogMessage << "HDCGCoarsening: GALERKIN CERTIFICATE ||A_c x - P^dag A P x||/||A_c x|| = " << rel << std::endl;
+    std::cout << GridLogMessage << "HDCGCoarsening: HERMITICITY CERTIFICATE |<x,A_c y> - <A_c x,y>|/|<x,A_c y>| = " << herm << std::endl;
+    return rel;
+  }
+
+
+  //////////////////////////////////////////////////////////////////////
+  // The device memory this coarsening holds that the memory manager
+  // CANNOT evict: BLAS stores, which are plain device allocations.
+  // Lattice fields are omitted on purpose -- they are evictable, so they
+  // cost eviction traffic, not failure.  What is reported is what has to
+  // fit alongside the manager's cache, the comms buffers and the shm
+  // segment, which is the budget an out-of-memory run actually breaks.
+  //////////////////////////////////////////////////////////////////////
+  void ReportDeviceFootprint(const std::string &when)
+  {
+    uint64_t p1 = MrhsProjectorL1.DeviceBytes();
+    uint64_t o1 = CoarseOp.DeviceBytes();
+    uint64_t df = Deflator.DeviceBytes();
+    const double G = 1024.*1024.*1024.;
+    std::cout << GridLogMessage << "HDCGCoarsening (" << when << "): resident device memory/rank "
+              << (p1+o1+df)/G << " GB = transfer L1 " << p1/G
+              << " + operator L1 " << o1/G
+              << " + deflation " << df/G << " GB" << std::endl;
+    std::cout << GridLogMessage << "HDCGCoarsening (" << when << "): not evictable; must fit alongside the "
+              << MemoryManager::DeviceMaxBytes/G
+              << " GB manager cache, the comms buffers and the shm segment" << std::endl;
+  }
+
+  void DiscardBasis(void)
+  {
+    rawNull.clear(); rawNull.shrink_to_fit();
+  }
+};
+
+//////////////////////////////////////////////////////////////////////
+// The solve chain on a borrowed coarsening.  Both fine precisions are
+// built; the outer picks the fp32 preconditioner through the seam when
+// Setup.FinePrecision is fp32.  OperatorShift is 0 in production; the
+// HDCG-refinement solver runs on the shifted operator.
+//////////////////////////////////////////////////////////////////////
+template<class Matrix,class MatrixF,class Coarsening>
+class HDCGSolver {
+public:
+  typedef typename Coarsening::FineField      FineField;
+  typedef typename Coarsening::FineFieldF     FineFieldF;
+  typedef typename Coarsening::CoarseOperator CoarseOperator;
+  typedef typename Coarsening::CoarseVector   CoarseVector;
+
+  Coarsening          &C;
+  HDCGMultiGridParams  Params;
+  int                  nrhs;
+  int                  _regrid;     // FIRST member-like init: SetNrhs before members capture grids
+  RealD                OperatorShift;
+
+  // fp64 fine chain
+  SchurDiagMooeeOperator<Matrix,FineField>    HermOpEO;
+  ShiftedHermOpLinearOperator<FineField>      HermOpShifted;
+  HermOpAdaptor<FineField>                    FineOp;         // Op = HermOp: the operator of the solve
+  ShiftedHermOpLinearOperator<FineField>      SmootherOp;
+  CGSmoother<FineField>                       SmootherCG;
+  ChebyshevSmoother<FineField>                SmootherCheby;
+  // coarse
+  HermitianLinearOperator<CoarseOperator,CoarseVector> LinOpC;
+  ConjugateGradient<CoarseVector>             CoarseCG;
+  DoNothingGuesser<CoarseVector>              NoGuess;
+  HPDSolver<CoarseVector>                     CoarseDeflatedCG;
+  ChebyshevInverter<CoarseVector>             CoarseCheby;
+  MrhsADEF2Preconditioner<FineField,CoarseVector,typename Coarsening::ProjectorL1_t> ADEF2;
+  // fp32 fine chain, behind the seam; the coarse solve is shared
+  SchurDiagMooeeOperator<MatrixF,FineFieldF>  HermOpEOF;
+  ShiftedHermOpLinearOperator<FineFieldF>     HermOpShiftedF;
+  HermOpAdaptor<FineFieldF>                   FineOpF;
+  ShiftedHermOpLinearOperator<FineFieldF>     SmootherOpF;
+  CGSmoother<FineFieldF>                      SmootherCGF;
+  ChebyshevSmoother<FineFieldF>               SmootherChebyF;
+  MrhsADEF2Preconditioner<FineFieldF,CoarseVector,typename Coarsening::ProjectorL1_t> ADEF2F;
+  MrhsMixedPrecPreconditioner<FineField,FineFieldF> PreconSeam;
+  TwoLevelCGmrhs<FineField>                   Outer;
+
+  LinearFunction<FineField>  &PickSmoother(void) {
+    return (Params.Smoother.Type==HDCGSmootherType::CG)
+      ? static_cast<LinearFunction<FineField>&>(SmootherCG)
+      : static_cast<LinearFunction<FineField>&>(SmootherCheby);
+  }
+  LinearFunction<FineFieldF> &PickSmootherF(void) {
+    return (Params.Smoother.Type==HDCGSmootherType::CG)
+      ? static_cast<LinearFunction<FineFieldF>&>(SmootherCGF)
+      : static_cast<LinearFunction<FineFieldF>&>(SmootherChebyF);
+  }
+  LinearFunction<CoarseVector> &PickCoarseSolve(void) {
+    return (Params.CoarseSolver.Method==HDCGCoarseMethod::DeflatedCG)
+      ? static_cast<LinearFunction<CoarseVector>&>(CoarseDeflatedCG)
+      : static_cast<LinearFunction<CoarseVector>&>(CoarseCheby);
+  }
+  template<class F, class Op, class ShiftedOp>
+  static LinearOperatorBase<F> &Unshifted(RealD shift, Op &op, ShiftedOp &shifted) {
+    return (shift==0.0) ? static_cast<LinearOperatorBase<F>&>(op) : static_cast<LinearOperatorBase<F>&>(shifted);
+  }
+
+  HDCGSolver(Matrix &Ddwf, MatrixF &DdwfF, Coarsening &_C, const HDCGMultiGridParams &P, int nr, RealD shift=0.0)
+    : C(_C), Params(P), nrhs(nr),
+      _regrid((C.SetNrhs(nr),0)),
+      OperatorShift(shift),
+      HermOpEO(Ddwf),
+      HermOpShifted(HermOpEO,shift),
+      FineOp(Unshifted<FineField>(shift,HermOpEO,HermOpShifted)),
+      SmootherOp(FineOp,P.Smoother.Shift),
+      SmootherCG(P.Smoother.Order,SmootherOp),
+      SmootherCheby(P.Smoother.ChebyLo,P.Smoother.ChebyHi,P.Smoother.Order,FineOp),
+      LinOpC(C.CoarseOp),
+      CoarseCG(P.CoarseSolver.Tol,P.CoarseSolver.MaxIt,false),
+      CoarseDeflatedCG(LinOpC,CoarseCG,NoGuess),
+      CoarseCheby(P.CoarseSolver.ChebyLo,P.CoarseSolver.ChebyHi,P.CoarseSolver.ChebyOrder,LinOpC),
+      ADEF2(FineOp,PickSmoother(),PickCoarseSolve(),PickCoarseSolve(),C.MrhsProjectorL1,C.Deflator,C.CMrhs,C.FrbGrid),
+      HermOpEOF(DdwfF),
+      HermOpShiftedF(HermOpEOF,shift),
+      FineOpF(Unshifted<FineFieldF>(shift,HermOpEOF,HermOpShiftedF)),
+      SmootherOpF(FineOpF,P.Smoother.Shift),
+      SmootherCGF(P.Smoother.Order,SmootherOpF),
+      SmootherChebyF(P.Smoother.ChebyLo,P.Smoother.ChebyHi,P.Smoother.Order,FineOpF),
+      ADEF2F(FineOpF,PickSmootherF(),PickCoarseSolve(),PickCoarseSolve(),C.MrhsProjectorL1,C.Deflator,C.CMrhs,C.GridsF.FrbGridF),
+      PreconSeam(ADEF2F,C.FrbGrid,C.GridsF.FrbGridF,nr),
+      Outer(P.Outer.Tol,P.Outer.MaxIt,FineOp,
+            (P.Setup.FinePrecision==MGPrecision::fp32)
+              ? static_cast<MrhsPreconditioner<FineField>&>(PreconSeam)
+              : static_cast<MrhsPreconditioner<FineField>&>(ADEF2),
+            C.FrbGrid,
+            (P.Outer.Algorithm==HDCGOuterAlgorithm::PrecBlockCGrQ) ? MrhsCGAlgorithm::PrecBlockCGrQ : MrhsCGAlgorithm::fPcg)
+  {
+    // BlockCGrQ assumes a stationary preconditioner: the fixed-iteration CG
+    // smoother is fixed WORK but a nonlinear map, and a coarse CG solved to a
+    // tolerance is not stationary either.  Chebyshev on both is.
+    if ( Params.Outer.Algorithm==HDCGOuterAlgorithm::PrecBlockCGrQ &&
+         ( Params.CoarseSolver.Method!=HDCGCoarseMethod::Chebyshev ||
+           Params.Smoother.Type     !=HDCGSmootherType::Chebyshev ) )
+      std::cout << GridLogWarning << "HDCGSolver: PrecBlockCGrQ with a non-stationary preconditioner "
+                << "(needs Smoother.Type Chebyshev AND CoarseSolver.Method Chebyshev); expect stagnation"
+                << std::endl;
+    std::cout << GridLogMessage << "HDCGSolver: Nrhs " << nr
+              << ", operator shift " << shift
+              << ", smoother " << (Params.Smoother.Type==HDCGSmootherType::CG ? "CG" : "Chebyshev")
+              << " order " << Params.Smoother.Order
+              << ", coarse " << (Params.CoarseSolver.Method==HDCGCoarseMethod::DeflatedCG ? "deflated CG" : "Chebyshev")
+              << " (deflation vectors " << C.Deflator.nev << ")"
+              << ", preconditioner fine level "
+              << (Params.Setup.FinePrecision==MGPrecision::fp32 ? "fp32 (seam at the outer CG)" : "fp64")
+              << ", outer " << (Params.Outer.Algorithm==HDCGOuterAlgorithm::PrecBlockCGrQ ? "PrecBlockCGrQ" : "fPcg")
+              << " EXACT fp64" << std::endl;
+  }
+
+  // One right-hand side through the LinearFunction interface (nrhs must be 1)
+  void Solve(const FineField &src, FineField &sol)
+  {
+    GRID_ASSERT( nrhs == 1 );
+    sol.Checkerboard() = src.Checkerboard();
+    GridStopWatch w; w.Start();
+    Outer(src,sol);
+    w.Stop();
+    FineField Ax(src.Grid()); Ax.Checkerboard()=src.Checkerboard();
+    FineOp.HermOp(sol,Ax);
+    Ax=Ax-src;
+    std::cout << GridLogMessage << "HDCGSolver: single rhs total " << w.Elapsed()
+              << "; FINAL Nrhs 1 (LinearFunction): true residual = " << std::sqrt(norm2(Ax)/norm2(src)) << std::endl;
+  }
+  void Solve(std::vector<FineField> &src, std::vector<FineField> &sol)
+  {
+    GRID_ASSERT( (int)src.size() == nrhs );
+    for(int r=0;r<nrhs;r++) sol[r].Checkerboard() = src[r].Checkerboard();
+    MemoryManager::EvictAll();
+    MemoryManager::DropCache();
+
+    GridStopWatch w; w.Start();
+    Outer(src,sol);
+    w.Stop();
+    std::cout << GridLogMessage << "HDCGSolver: Nrhs "<<nrhs<<" total " << w.Elapsed()
+              << "  (per RHS: " << w.useconds()/1.0e6/nrhs << " s)" << std::endl;
+    {
+      FineField Ax(src[0].Grid()); Ax.Checkerboard()=src[0].Checkerboard();
+      RealD worst=0.0;
+      for(int r=0;r<nrhs;r++){
+        FineOp.HermOp(sol[r],Ax);
+        Ax=Ax-src[r];
+        RealD rn=std::sqrt(norm2(Ax)/norm2(src[r]));
+        std::cout << GridLogMessage << "FINAL Nrhs "<<nrhs<<": rhs["<<r<<"] true residual = " << rn << std::endl;
+        worst=std::max(worst,rn);
+      }
+      std::cout << GridLogMessage << "FINAL Nrhs "<<nrhs<<": worst-case residual = " << worst << std::endl;
+    }
+  }
+};
+
+//////////////////////////////////////////////////////////////////////
+// HDCG refinement of the basis (the paper's loop): coarsen the SHIFTED
+// operator, deflate it, solve each raw basis vector on it with a
+// short-tolerance HDCG, keep the normalised solutions.  A no-op unless
+// Subspace.Refine is HDCG and no refined basis was loaded.  The caller
+// then coarsens the unshifted operator as usual.
+//////////////////////////////////////////////////////////////////////
+template<class Matrix,class MatrixF,class Coarsening>
+void HDCGRefineSubspace(Coarsening &C, Matrix &Ddwf, MatrixF &DdwfF,
+                        LinearOperatorBase<typename Coarsening::FineField> &HermOp, int nr)
+{
+  typedef typename Coarsening::FineField FineField;
+  const HDCGSubspaceParams &S = C.Params.Subspace;
+  if ( S.Refine != HDCGRefineMethod::HDCG ) return;
+  if ( C.FileExists(S.RefinedSubspaceFile) ) return;
+
+  std::cout << GridLogMessage << "HDCGRefineSubspace: shift " << S.RefineShift << " tol " << S.RefineTol
+            << " maxit " << S.RefineMaxIt << " smoother order " << S.RefineSmootherOrder << std::endl;
+  ShiftedHermOpLinearOperator<FineField> Shifted(HermOp,S.RefineShift);
+  HermOpAdaptor<FineField>               ShiftedOp(Shifted);
+  C.Coarsen(ShiftedOp,1);   // refinement below still needs the RAW basis
+  C.CoarseLanczos(nr);
+
+  HDCGMultiGridParams PR = C.Params;
+  PR.Outer.Tol      = S.RefineTol;
+  PR.Outer.MaxIt    = S.RefineMaxIt;
+  PR.Smoother.Order = S.RefineSmootherOrder;
+  HDCGSolver<Matrix,MatrixF,Coarsening> Refiner(Ddwf,DdwfF,C,PR,nr,S.RefineShift);
+  C.RefineSubspaceHDCG(HermOp,Refiner.Outer,nr);
+
+  if ( S.RefinedSubspaceFile.length() )
+    saveSubspace(C.rawNull, S.RefinedSubspaceFile, S.FileControl);
+}
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/HDCGMultiGridParams.h

@@ -0,0 +1,177 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/HDCGMultiGridParams.h
+
+    Copyright (C) 2026
+
+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.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+//////////////////////////////////////////////////////////////////////////////////////
+// Parameters of the two-level mrhs HDCG (Hermitian, Schur-preconditioned)
+// multigrid, HDCGMultiGrid.h.  Serialisable: read from XML, template written
+// on a missing file.  Setup (blocking, batch, basis file, fine precision) is
+// the same struct as the PVdagM chain's; Block2 is only used when a
+// coarse-coarse level exists but must still divide the coarse lattice.
+//////////////////////////////////////////////////////////////////////////////////////
+
+GRID_SERIALIZABLE_ENUM(HDCGSubspaceMethod, undef, ChebyshevNew, 1, Chebyshev, 2, FineEvecs, 3);
+GRID_SERIALIZABLE_ENUM(HDCGRefineMethod,   undef, None, 1, CG, 2, HDCG, 3);
+GRID_SERIALIZABLE_ENUM(HDCGLanczosMethod,  undef, None, 1, Block, 2, Plain, 3);
+GRID_SERIALIZABLE_ENUM(HDCGCoarseMethod,   undef, DeflatedCG, 1, Chebyshev, 2);
+GRID_SERIALIZABLE_ENUM(HDCGSmootherType,   undef, CG, 1, Chebyshev, 2);
+GRID_SERIALIZABLE_ENUM(HDCGOuterAlgorithm, undef, fPcg, 1, PrecBlockCGrQ, 2);
+
+// Near-null basis.  ChebyshevNew filters noise through the Hermitian operator
+// with the upper bound only; Chebyshev takes lo/hi/order; FineEvecs sums
+// FineEvecSample consecutive fp32 fine eigenvector records per basis vector.
+// A ChebyHi <= 0 means 1.1 x the power-method estimate.  FileControl is the
+// scidac binary-IO control word of the basis files (0 for pre-2026 HDCG files).
+struct HDCGSubspaceParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGSubspaceParams,
+                                  HDCGSubspaceMethod, Method,
+                                  RealD,       ChebyHi,
+                                  RealD,       ChebyLo,
+                                  int,         ChebyOrder,
+                                  std::string, FineEvecFile,
+                                  int,         FineEvecSample,
+                                  int,         FileControl,
+                                  HDCGRefineMethod, Refine,
+                                  RealD,       RefineShift,         // CG: the shift of the refining CG; HDCG: of the refining operator
+                                  RealD,       RefineTol,
+                                  int,         RefineMaxIt,
+                                  int,         RefineSmootherOrder, // HDCG refine only
+                                  std::string, RefinedSubspaceFile);
+  HDCGSubspaceParams()
+    : Method(HDCGSubspaceMethod::ChebyshevNew), ChebyHi(0.0), ChebyLo(0.01), ChebyOrder(500),
+      FineEvecFile(""), FineEvecSample(1),
+      FileControl(BinaryIO::BINARYIO_LEXICOGRAPHIC|BinaryIO::BINARYIO_AGGREGATE),
+      Refine(HDCGRefineMethod::None), RefineShift(1.0e-3), RefineTol(1.0e-3), RefineMaxIt(500),
+      RefineSmootherOrder(11), RefinedSubspaceFile("") {};
+};
+
+// Coarse eigenvectors for the deflated coarse solve: block Lanczos at the
+// solve's Nrhs, or plain Lanczos at Nrhs 1, on the same coarse operator.
+// ChebyHi <= 0 means 1.1 x the power-method estimate on the coarse operator.
+struct HDCGLanczosParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGLanczosParams,
+                                  HDCGLanczosMethod, Method,
+                                  int,   Nstop,
+                                  int,   Nk,
+                                  int,   Nm,
+                                  RealD, ChebyLo,
+                                  RealD, ChebyHi,
+                                  int,   ChebyOrder,
+                                  RealD, Tol,
+                                  int,   MaxIt,
+                                  std::string, EvecFile);
+  HDCGLanczosParams()
+    : Method(HDCGLanczosMethod::Block), Nstop(32), Nk(32), Nm(128),
+      ChebyLo(0.01), ChebyHi(0.0), ChebyOrder(201), Tol(1.0e-4), MaxIt(100), EvecFile("") {};
+};
+
+struct HDCGCoarseSolverParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGCoarseSolverParams,
+                                  HDCGCoarseMethod, Method,
+                                  RealD, Tol,            // DeflatedCG
+                                  int,   MaxIt,
+                                  RealD, ChebyLo,        // Chebyshev polynomial inverse
+                                  RealD, ChebyHi,
+                                  int,   ChebyOrder);
+  HDCGCoarseSolverParams()
+    : Method(HDCGCoarseMethod::DeflatedCG), Tol(5.0e-2), MaxIt(3000),
+      ChebyLo(1.0e-2), ChebyHi(40.0), ChebyOrder(120) {};
+};
+
+// CG: Order iterations of CG on the operator shifted by Shift.
+// Chebyshev: order Order polynomial on [ChebyLo,ChebyHi] of the unshifted operator.
+struct HDCGSmootherParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGSmootherParams,
+                                  HDCGSmootherType, Type,
+                                  RealD, Shift,
+                                  int,   Order,
+                                  RealD, ChebyLo,
+                                  RealD, ChebyHi);
+  HDCGSmootherParams() : Type(HDCGSmootherType::CG), Shift(2.0), Order(7), ChebyLo(2.0), ChebyHi(92.0) {};
+};
+
+struct HDCGOuterParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGOuterParams,
+                                  RealD, Tol,
+                                  int,   MaxIt,
+                                  HDCGOuterAlgorithm, Algorithm);
+  HDCGOuterParams() : Tol(1.0e-8), MaxIt(500), Algorithm(HDCGOuterAlgorithm::fPcg) {};
+};
+
+struct HDCGMultiGridParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGMultiGridParams,
+                                  MGSetupParams,          Setup,
+                                  HDCGSubspaceParams,     Subspace,
+                                  HDCGLanczosParams,      Lanczos,
+                                  HDCGCoarseSolverParams, CoarseSolver,
+                                  HDCGSmootherParams,     Smoother,
+                                  HDCGOuterParams,        Outer);
+  HDCGMultiGridParams() {};
+};
+
+inline void CheckValidity(const HDCGMultiGridParams &P)
+{
+  GRID_ASSERT( P.Setup.Block1.size() == 4 );
+  GRID_ASSERT( P.Setup.Block2.size() == 4 );
+  GRID_ASSERT( P.Setup.CoarsenBatch  >= 1 );
+  GRID_ASSERT( P.Setup.FinePrecision != MGPrecision::undef );
+  GRID_ASSERT( P.Subspace.Method     != HDCGSubspaceMethod::undef );
+  GRID_ASSERT( P.Subspace.Refine     != HDCGRefineMethod::undef );
+  GRID_ASSERT( P.Subspace.FineEvecSample >= 1 );
+  GRID_ASSERT( P.Lanczos.Method      != HDCGLanczosMethod::undef );
+  GRID_ASSERT( P.Lanczos.Nstop <= P.Lanczos.Nk && P.Lanczos.Nk <= P.Lanczos.Nm );
+  GRID_ASSERT( P.CoarseSolver.Method != HDCGCoarseMethod::undef );
+  GRID_ASSERT( P.CoarseSolver.Tol     > 0.0 );
+  GRID_ASSERT( P.Smoother.Type       != HDCGSmootherType::undef );
+  GRID_ASSERT( P.Smoother.Order       > 0 );
+  GRID_ASSERT( P.Outer.Tol            > 0.0 );
+  GRID_ASSERT( P.Outer.Algorithm     != HDCGOuterAlgorithm::undef );
+}
+
+// As ReadPVdagMMultiGridParams: a missing file gets a template written next
+// to it and false returned; an empty filename means defaults.
+inline bool ReadHDCGMultiGridParams(HDCGMultiGridParams &P, const std::string &xmlfile)
+{
+  if ( xmlfile.length() ) {
+    bool good;
+    { std::ifstream f(xmlfile); good = f.good(); }
+    if ( !good ) {
+      std::cout << GridLogMessage << "HDCGMultiGridParams: " << xmlfile << " does not exist" << std::endl;
+      if ( GlobalSharedMemory::WorldRank == 0 ) {
+        XmlWriter WR(xmlfile+".templ");
+        write(WR, "HDCGMultiGridParams", P);
+        std::cout << GridLogMessage << "HDCGMultiGridParams: template written to "
+                  << xmlfile << ".templ" << std::endl;
+      }
+      return false;
+    }
+    XmlReader RD(xmlfile);
+    read(RD, "HDCGMultiGridParams", P);
+  }
+  CheckValidity(P);
+  std::cout << GridLogMessage << "HDCGMultiGridParams ("
+            << (xmlfile.length() ? xmlfile : std::string("defaults")) << "):" << std::endl;
+  std::cout << P << std::endl;
+  return true;
+}
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/MrhsPreconditioner.h

@@ -0,0 +1,73 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/MrhsPreconditioner.h
+
+    Copyright (C) 2026
+
+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.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+//////////////////////////////////////////////////////////////////////
+// The mrhs interfaces.  This header is in the Algorithms.h include
+// chain (AdefMrhs.h needs it); the concrete mrhs classes are in
+// MrhsMultiGrid.h, reached through PVdagMMultiGrid.h / HDCGMultiGrid.h.
+//
+// MrhsLinearFunction: vector-of-fields in, vector out.  The outer level
+// carries mrhs as std::vector<Field>; below it mrhs is PACKED into a
+// single D+1 field (rhs = dim 0) and the coarse classes are plain
+// LinearFunctions on that.
+//////////////////////////////////////////////////////////////////////
+template<class Field>
+class MrhsLinearFunction {
+public:
+  virtual void operator()(std::vector<Field> &in, std::vector<Field> &out) = 0;
+  virtual ~MrhsLinearFunction(){};
+};
+
+//////////////////////////////////////////////////////////////////////
+// An mrhs preconditioner for the two-level CG family (AdefMrhs.h):
+// operator() is M1; Vstart chooses x_0 from the source (Tang's ADEF-2
+// start x_0 = Q b; the default is the zero start).  The timers are the
+// M1 breakdown the outer solver reports.
+//////////////////////////////////////////////////////////////////////
+template<class Field>
+class MrhsPreconditioner : public MrhsLinearFunction<Field> {
+public:
+  GridStopWatch ProjectTimer, PromoteTimer, DeflateTimer, CoarseTimer, FineTimer, SmoothTimer;
+  virtual void Vstart(std::vector<Field> &x, std::vector<Field> &src)
+  {
+    for(int r=0;r<(int)x.size();r++) x[r] = Zero();
+  }
+  virtual void ResetTimers(void)
+  {
+    ProjectTimer.Reset(); PromoteTimer.Reset(); DeflateTimer.Reset();
+    CoarseTimer.Reset();  FineTimer.Reset();    SmoothTimer.Reset();
+  }
+  virtual void ReportTimers(const std::string &prefix)
+  {
+    std::cout<<GridLogMessage<<"**** M1 breakdown:"<<std::endl;
+    std::cout<<GridLogMessage<<prefix<<"Project "<<ProjectTimer.Elapsed()<<std::endl;
+    std::cout<<GridLogMessage<<prefix<<"Promote "<<PromoteTimer.Elapsed()<<std::endl;
+    std::cout<<GridLogMessage<<prefix<<"Deflate "<<DeflateTimer.Elapsed()<<std::endl;
+    std::cout<<GridLogMessage<<prefix<<"Coarse  "<<CoarseTimer.Elapsed()<<std::endl;
+    std::cout<<GridLogMessage<<prefix<<"Fine    "<<FineTimer.Elapsed()<<std::endl;
+    std::cout<<GridLogMessage<<prefix<<"Smooth  "<<SmoothTimer.Elapsed()<<std::endl;
+  }
+  virtual ~MrhsPreconditioner(){};
+};
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/MultiGridIO.h

@@ -0,0 +1,94 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/MultiGridIO.h
+
+    Copyright (C) 2026
+
+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.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+//////////////////////////////////////////////////////////////////////
+// Multigrid setup I/O shared by the PVdagM and HDCG chains: a set of
+// bare-vector scidac records (a near-null basis, or a set of
+// eigenvectors with their eigenvalues beside them in XML).
+//
+// The coarse OPERATOR is not written: coarsening is cheap relative
+// to producing the basis, and is redone.
+//
+// `control` is the scidac binary-IO control word.  The library default
+// is lexicographic + aggregated; pass 0 to read files written by the
+// pre-2026 HDCG drivers, which used control 0.
+//////////////////////////////////////////////////////////////////////
+static const int MultiGridIODefaultControl = BinaryIO::BINARYIO_LEXICOGRAPHIC|BinaryIO::BINARYIO_AGGREGATE;
+
+template <class Field>
+void saveSubspace(std::vector<Field> &subspace, std::string const fname,
+		  int control = MultiGridIODefaultControl)
+{
+#ifdef HAVE_LIME
+  Grid::emptyUserRecord record;
+  Grid::ScidacWriter SW(subspace[0].Grid()->IsBoss());
+  SW.open(fname);
+  for (int k = 0; k < (int)subspace.size(); k++) {
+    SW.writeScidacFieldRecord(subspace[k], record, 0, control);
+  }
+  SW.close();
+#endif
+}
+template <class Field>
+void loadSubspace(std::vector<Field> &subspace, std::string const fname,
+		  int control = MultiGridIODefaultControl)
+{
+#ifdef HAVE_LIME
+  Grid::emptyUserRecord record;
+  Grid::ScidacReader SR;
+  SR.open(fname);
+  for (int k = 0; k < (int)subspace.size(); k++) {
+    SR.readScidacFieldRecord(subspace[k], record, control);
+  }
+  SR.close();
+#endif
+}
+
+//////////////////////////////////////////////////////////////////////
+// Eigenpairs: the vectors as a subspace file, the values in
+// <fname>.evals.xml.  evec.size() fixes how many are read.
+//////////////////////////////////////////////////////////////////////
+template <class Field>
+void saveEigenpairs(std::vector<Field> &evec, std::vector<RealD> &eval, std::string const fname,
+		    int control = MultiGridIODefaultControl)
+{
+  GRID_ASSERT(evec.size()==eval.size());
+  saveSubspace(evec, fname, control);
+  if ( evec[0].Grid()->IsBoss() ) {
+    XmlWriter WR(fname+".evals.xml");
+    write(WR,"evals",eval);
+  }
+}
+// The eigenvalue file sets the count; the vectors are allocated on grid.
+template <class Field>
+void loadEigenpairs(std::vector<Field> &evec, std::vector<RealD> &eval, std::string const fname,
+		    GridBase *grid, int control = MultiGridIODefaultControl)
+{
+  XmlReader RD(fname+".evals.xml");
+  read(RD,"evals",eval);
+  evec.clear();
+  evec.resize(eval.size(), Field(grid));
+  loadSubspace(evec, fname, control);
+}
+
+NAMESPACE_END(Grid);