From 3a3a20b8c9cfbc65f34a43912a622870256ecf01 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 29 Sep 2026 13:46:09 -0400
Subject: [PATCH] Updates

---
 Grid/algorithms/multigrid/BlockCyclic.h       |    3 +-
 .../GeneralCoarsenedMatrixMultiRHS.h          | 1272 +++++++++++++++++
 MPI_benchmark/simple_reproducer_noncache.cc   |  214 +++
 ...Example_pvdagm_3level_DenseCoarseMatrix.cc |  990 +++++++++++++
 .../Example_pvdagm_3level_FourierSmoother.cc  | 1145 +++++++++++++++
 systems/WorkArounds.txt                       |   85 +-
 tests/debug/Test_blockcyclic.cc               |    2 +-
 tests/debug/Test_coarse.cc                    |  261 ++++
 tests/debug/Test_coarse_coarsen.cc            |  279 ++++
 tests/debug/Test_mrhs_deflation.cc            |   89 ++
 tests/debug/Test_sloppy_dagger.cc             |  207 +++
 11 files changed, 4506 insertions(+), 41 deletions(-)
 create mode 100644 Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHS.h
 create mode 100644 MPI_benchmark/simple_reproducer_noncache.cc
 create mode 100644 examples/Example_pvdagm_3level_DenseCoarseMatrix.cc
 create mode 100644 examples/Example_pvdagm_3level_FourierSmoother.cc
 create mode 100644 tests/debug/Test_coarse.cc
 create mode 100644 tests/debug/Test_coarse_coarsen.cc
 create mode 100644 tests/debug/Test_mrhs_deflation.cc
 create mode 100644 tests/debug/Test_sloppy_dagger.cc

diff --git a/Grid/algorithms/multigrid/BlockCyclic.h b/Grid/algorithms/multigrid/BlockCyclic.h
index c27f53768..eda021195 100644
--- a/Grid/algorithms/multigrid/BlockCyclic.h
+++ b/Grid/algorithms/multigrid/BlockCyclic.h
@@ -39,8 +39,7 @@ typedef ComplexD DenseInverseScalar;
 // BlockCyclicLayout: the index arithmetic of a 2D block-cyclic distribution
 // of an N x N matrix over a Pr x Pc logical process grid with block size nb.
 //
-// This is stage 1 of the 2D distributed dense inverse
-// (documentation/DistributedDenseInverse2D.tex).  It is deliberately
+// This is stage 1 of the 2D distributed dense inverse.  It is deliberately
 // COMMUNICATOR-FREE: every mapping is a static pure function of
 // (N, nb, Pr, Pc), so the whole layout is exhaustively unit-testable on one
 // rank with no MPI in the loop (Test_blockcyclic).  A thin instance layer
diff --git a/Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHS.h b/Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHS.h
new file mode 100644
index 000000000..54c0d67b5
--- /dev/null
+++ b/Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHS.h
@@ -0,0 +1,1272 @@
+/*************************************************************************************
+
+    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 MultiGeneralCoarsenedOperator : public SparseMatrixBase<Lattice<iVector<CComplex,nbasis > > >  {
+public:
+  typedef typename CComplex::scalar_object SComplex;
+  // The BLAS scalar of this level follows the coefficient precision:
+  // ComplexF for an fp32 coarse space, ComplexD for fp64.
+  typedef typename GridTypeMapper<CComplex>::scalar_type CoarseBLASScalar;
+  typedef MultiGeneralCoarsenedOperator<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 CoarseVector Field;
+
+  // Block operations on the fine vectors carry the fine layout, which need
+  // not be the coarse one
+  typedef decltype(innerProduct(Fobj(),Fobj())) FineInner;
+  typedef Lattice<FineInner>                    FineComplexField;
+  typedef Lattice<FineInner>                    BlockComplexField;
+
+  ////////////////////
+  // Data members
+  //
+  // Nrhs independent: the D dimensional coarse grid, the geometry, the padded
+  // cell that supplies the stencil grid, the stencil, and the matrix elements.
+  //
+  // Nrhs dependent: the D+1 grid, its padded cell, and the BLAS B/C buffers
+  // with their pointer tables. Owned by SetNRHS().
+  ////////////////////
+  GridCartesian *       _CoarseGrid;        // D dimensional
+  NonLocalStencilGeometry geom;
+  NonLocalStencilGeometry geom_srhs;
+  PaddedCell CellD;                         // D dimensional, supplies stencil grid
+  GeneralLocalStencil Stencil;              // D dimensional
+
+  int                   _Nrhs;
+  GridCartesian *       _CoarseGridMulti;   // D+1 dimensional, SetNRHS
+  PaddedCell *          CellMulti;          // D+1 dimensional, SetNRHS
+
+  deviceVector<calcVector> BLAS_B;
+  deviceVector<calcVector> BLAS_C;
+  // The matrix, ONE buffer, site major: a site's npoint blocks are contiguous.
+  // Both organisations read it through their pointer tables, the folded one a
+  // site at a time and the unfolded one a point at a time with a stride of the
+  // whole stencil, so nothing is duplicated.  Writing or reading one point
+  // goes through a single-point scratch.
+  deviceVector<calcMatrix>               BLAS_Amat;   // [site][point]
+  deviceVector<calcMatrix>               BLAS_Ascr;   // one point, for I/O
+
+  std::vector<deviceVector<CoarseBLASScalar *> > BLAS_AP;
+  std::vector<deviceVector<CoarseBLASScalar *> > BLAS_BP;
+  deviceVector<CoarseBLASScalar *>               BLAS_CP;
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Stencil legs carried in the BATCH dimension, LegGroup at a time.
+  //
+  // One call per stencil point batches over the local coarse volume alone,
+  // which at the production point is 1024 sites.  That is too few for the
+  // single-precision kernel: it then runs at half the bandwidth the double
+  // one reaches and takes the same time, so an fp32 coarse space buys
+  // nothing in the operator.  The same shape at four times the batch runs
+  // 1.7x faster in fp32 than fp64, so the fix is more batch, not fewer
+  // bytes.  Grouping G legs multiplies the batch by G and divides the launch
+  // count by G, at the price of G partial outputs summed at the end.
+  //
+  // G must divide npoint; 3 divides both the 33 point (next-to-nearest) and
+  // 81 point (full box) stencils.  G=1 is the ungrouped form.  The only
+  // numerical difference is the ORDER of the sum over stencil points.
+  ///////////////////////////////////////////////////////////////////////////
+  // The group need not divide npoint: the last call carries the remainder,
+  // with its own smaller tables.  The batch a call presents is LegGroup times
+  // the local coarse volume, so the group that saturates the kernel depends
+  // on the decomposition, and the default is chosen from it rather than
+  // fixed.  TargetBatch is where the measured single-precision curve flattens
+  // (530 GB/s at 1024, 695 at 9216, 728 at 33792 for 60x12x60 on MI250X).
+  //
+  // 9216 is nine times the production local coarse volume, so the default
+  // lands on nine legs per call: 81 = 9x9 for the full-box stencil and
+  // 33 = 9+9+9+6 for the next-to-nearest one, which is the grouping the
+  // earlier HDCG coarse operator used.
+  static const int TargetBatch = 9216;
+
+  int LegGroup = 1;
+  std::vector<deviceVector<CoarseBLASScalar *> > BLAS_APg; // per group
+  std::vector<deviceVector<CoarseBLASScalar *> > BLAS_BPg;
+  deviceVector<CoarseBLASScalar *>               BLAS_CPg; // LegGroup*sites
+  deviceVector<CoarseBLASScalar *>               BLAS_CPr; // remainder*sites
+  deviceVector<calcVector>                       BLAS_Cg;  // LegGroup partials
+
+  int LegGroups(void)    const { return geom.npoint/LegGroup; }          // full
+  int LegRemainder(void) const { return geom.npoint%LegGroup; }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Matrix pointer table for the grouped call.  Group g holds legs
+  // g*LegGroup .. and the batch index runs site fastest within a leg.  The
+  // last group is short when LegGroup does not divide npoint.
+  ///////////////////////////////////////////////////////////////////////////
+  void BuildGroupedA(void)
+  {
+    int32_t sites = _CoarseGrid->lSites();
+    int ng = LegGroups() + (LegRemainder()?1:0);
+    BLAS_APg.resize(ng);
+    for(int g=0;g<ng;g++){
+      int legs = (g<LegGroups()) ? LegGroup : LegRemainder();
+      BLAS_APg[g].resize(legs*sites);
+      for(int l=0;l<legs;l++){
+	int p = g*LegGroup+l;
+	for(int32_t ss=0;ss<sites;ss++){
+	  CoarseBLASScalar *ptr = (CoarseBLASScalar *)&BLAS_Amat[(uint64_t)ss*geom.npoint+p];
+	  acceleratorPut(BLAS_APg[g][l*sites+ss],ptr);
+	}
+      }
+    }
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Stencil legs folded along the CONTRACTION instead of the batch.
+  //
+  // One GEMM per site of [nbasis x (KFold*nbasis)] . [(KFold*nbasis) x nrhs],
+  // accumulating straight into C, so there are no partial results to combine.
+  // It needs the site's matrices contiguous, which is a repack done once, and
+  // the neighbour vectors gathered into one block per site, which replaces the
+  // copy the ungrouped path already does.
+  //
+  // Measured on MI250X at nbasis 60, batch 1024, single precision: the short
+  // form (K=60) reaches 487 GB/s at one right-hand side and 625 at twelve,
+  // the folded form (K=1980) 912 and 823.  Beyond about K=2000 it collapses:
+  // K=4860, the full 81 point box, gives 304 GB/s.  So the fold is chunked to
+  // keep K near 2000, and the chunks accumulate.
+  //
+  // The gather costs one extra write of nbasis*nrhs per leg per site, which is
+  // 3% of the traffic at one right-hand side and 20% at twelve; the fold wins
+  // at both, and by most at one, which is the HMC case.
+  ///////////////////////////////////////////////////////////////////////////
+  static const int TargetK = 2048;
+
+  int KFold = 0;                                     // legs per GEMM; 0 = off
+  deviceVector<CoarseBLASScalar>                 BLAS_BK;   // [site][rhs][leg][nbasis]
+  std::vector<deviceVector<CoarseBLASScalar *> > BLAS_AKP;  // per chunk, per site
+  std::vector<deviceVector<CoarseBLASScalar *> > BLAS_BKP;  // per chunk, per site
+  deviceVector<CoarseBLASScalar *>               BLAS_BPflat; // [point][site] neighbours
+
+  int KChunks(void) const { return (geom.npoint + KFold - 1)/KFold; }
+  int KLegs(int g)  const { return std::min(KFold, geom.npoint - g*KFold); }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Repack the matrix so a site's stencil points are contiguous: the folded
+  // GEMM reads each site's block as one column-major nbasis x (legs*nbasis)
+  // matrix, which is also a single sequential stream rather than npoint
+  // interleaved ones.  Costs one pass over the matrix, once per solve setup.
+  ///////////////////////////////////////////////////////////////////////////
+  void PackFoldedMatrix(void)
+  {
+    // Nothing to pack: the matrix is already site major.  Only the pointer
+    // table is built, one entry per site per chunk, at the chunk's first leg.
+    int32_t sites = _CoarseGrid->lSites();
+    const int np = geom.npoint;
+    BLAS_AKP.resize(KChunks());
+    for(int g=0;g<KChunks();g++){
+      BLAS_AKP[g].resize(sites);
+      for(int32_t ss=0;ss<sites;ss++){
+	CoarseBLASScalar *ptr = (CoarseBLASScalar *)&BLAS_Amat[(uint64_t)ss*np + (uint64_t)g*KFold];
+	acceleratorPut(BLAS_AKP[g][ss],ptr);
+      }
+    }
+  }
+
+  // Legs per GEMM that keeps K near where the kernel runs best.
+  int AutoKFold(void) const
+  {
+    int K = TargetK/nbasis;
+    if ( K < 1 )             K = 1;
+    if ( K > geom.npoint )   K = geom.npoint;
+    return K;
+  }
+  // Legs per call that brings the batch up to where the kernel saturates.
+  int AutoLegGroup(void) const
+  {
+    int32_t sites = _CoarseGrid->lSites();
+    int G = (TargetBatch + sites - 1)/sites;
+    if ( G < 1 )             G = 1;
+    if ( G > geom.npoint )   G = geom.npoint;
+    return G;
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // The fold along K gathers the neighbours, which costs nbasis*nrhs per leg
+  // per site: 16 MB of traffic per apply at one right hand side and 195 at
+  // twelve, against a matrix of 973.  Measured at the production point it is
+  // worth 7.6% at one right hand side and costs 7.6% at twelve, so the fold
+  // is used only for a small batch, and a large one carries its legs in the
+  // batch dimension instead.
+  ///////////////////////////////////////////////////////////////////////////
+  static const int FoldNrhsMax = 4;
+
+  void ChooseOrganisation(int nrhs)
+  {
+    if ( nrhs <= FoldNrhsMax ) {
+      KFold = AutoKFold();
+    } else {
+      KFold = 0;
+    }
+    LegGroup = AutoLegGroup();
+    if ( KFold ) {
+      PackFoldedMatrix();
+    } else {
+      BuildGroupedA();
+    }
+  }
+
+  ///////////////////////
+  // Interface
+  ///////////////////////
+  GridBase      * Grid(void)           { CheckGridSet(); return _CoarseGridMulti; };
+  GridCartesian * CoarseGrid(void)     { CheckGridSet(); return _CoarseGridMulti; };
+  GridCartesian * CoarseGridD(void)    { return _CoarseGrid; };        // lower dimensional grid
+  int             Nrhs(void)           { CheckGridSet(); return _Nrhs; };
+
+  void CheckGridSet(void)
+  {
+    if ( _CoarseGridMulti == nullptr ) {
+      std::cout << GridLogError
+		<< "MultiGeneralCoarsenedOperator: the multiRHS grid has not been set."
+		<< std::endl;
+      std::cout << GridLogError
+		<< "  Call SetGrid(CoarseGridMulti) with the D+1 dimensional grid your"
+		<< std::endl;
+      std::cout << GridLogError
+		<< "  coarse vectors live on, before Grid(), Nrhs() or M()."
+		<< std::endl;
+      GRID_ASSERT(_CoarseGridMulti != nullptr);
+    }
+  }
+
+  //////////////////////////////////////////////////////////////////////////
+  // Accessors.  Note Grid() is the D+1 multiRHS grid, so a consumer wanting
+  // the space the matrix elements live on must ask for CoarseGridD().
+  //////////////////////////////////////////////////////////////////////////
+  // Resident device memory this operator holds (matrix elements and the
+  // Nrhs-dependent BLAS buffers).  Not evictable.
+  uint64_t DeviceBytes(void)
+  {
+    uint64_t b = (uint64_t)(BLAS_B.capacity()+BLAS_C.capacity()+BLAS_Cg.capacity())*sizeof(calcVector)
+               + (uint64_t)BLAS_BK.capacity()*sizeof(CoarseBLASScalar);
+    b += (uint64_t)(BLAS_Amat.capacity()+BLAS_Ascr.capacity())*sizeof(calcMatrix);
+    return b;
+  }
+
+  NonLocalStencilGeometry & Geometry(void)      { return geom_srhs; };
+
+  // One stencil point in or out of the site-major buffer.
+  void MatrixPointIn(int p,const deviceVector<calcMatrix> &src)
+  {
+    int32_t sites = _CoarseGrid->lSites();
+    const int np = geom.npoint;
+    calcMatrix       *dst = &BLAS_Amat[0];
+    const calcMatrix *sp  = &src[0];
+    accelerator_for(ss,sites,1,{ dst[(uint64_t)ss*np+p] = sp[ss]; });
+  }
+  void MatrixPointOut(int p,deviceVector<calcMatrix> &dst)
+  {
+    int32_t sites = _CoarseGrid->lSites();
+    const int np = geom.npoint;
+    const calcMatrix *sp = &BLAS_Amat[0];
+    dst.resize(sites);
+    calcMatrix *dp = &dst[0];
+    accelerator_for(ss,sites,1,{ dp[ss] = sp[(uint64_t)ss*np+p]; });
+  }
+  void ExtractMatrix(int p,CoarseMatrix &A)
+  {
+    MatrixPointOut(p,BLAS_Ascr);
+    BLAStoGrid(A,BLAS_Ascr);
+  }
+
+  // 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_Ascr);
+    MatrixPointIn(p,BLAS_Ascr);
+  }
+  void GetMatrix (int p,std::vector<CoarseMatrix> & A)
+  {
+    GRID_ASSERT(A.size()==geom_srhs.npoint);
+    MatrixPointOut(p,BLAS_Ascr);
+    BLAStoGrid(A[p],BLAS_Ascr);
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Constructor takes the D dimensional coarse grid. Everything built here
+  // is independent of Nrhs, in particular the matrix elements, which must
+  // survive a change of Nrhs untouched.
+  ///////////////////////////////////////////////////////////////////////////
+  MultiGeneralCoarsenedOperator(NonLocalStencilGeometry &_geom,GridCartesian *CoarseGrid) :
+    _CoarseGrid(CoarseGrid),
+    geom_srhs(_geom),
+    geom(CoarseGrid,_geom.hops,_geom.skip),
+    CellD(geom.Depth(),CoarseGrid),
+    Stencil(CellD.grids.back(),geom.shifts), // D dimensional padded cell stencil
+    _Nrhs(-1),
+    _CoarseGridMulti(nullptr),
+    CellMulti(nullptr)
+  {
+    int32_t unpadded_sites = _CoarseGrid->lSites();
+
+    /////////////////////////////////////////////////
+    // Matrix elements and their pointer table
+    /////////////////////////////////////////////////
+    BLAS_Amat.resize((uint64_t)unpadded_sites*geom.npoint);
+    BLAS_Ascr.resize(unpadded_sites);
+    BLAS_AP.resize(geom.npoint);
+    for(int p=0;p<geom.npoint;p++){
+      BLAS_AP[p].resize(unpadded_sites);
+    }
+
+    // Point p at site ss: the unfolded path strides by the whole stencil.
+    for(int p=0;p<geom.npoint;p++){
+      for(int ss=0;ss<unpadded_sites;ss++){
+	CoarseBLASScalar *ptr = (CoarseBLASScalar *)&BLAS_Amat[(uint64_t)ss*geom.npoint+p];
+	acceleratorPut(BLAS_AP[p][ss],ptr);
+      }
+    }
+
+    // The organisation is chosen in SetGrid, where Nrhs is known.
+    KFold    = 0;
+    LegGroup = AutoLegGroup();
+    BuildGroupedA();
+    std::cout << GridLogMessage << "MultiGeneralCoarsenedOperator: stencil "
+	      << geom.npoint << " points, local coarse volume " << unpadded_sites
+	      << ", legs per GEMM " << LegGroup
+	      << " (batch " << LegGroup*unpadded_sites << ", "
+	      << LegGroups() << " full call(s)"
+	      << (LegRemainder() ? " + a remainder of "+std::to_string(LegRemainder()) : "")
+	      << ")" << std::endl;
+  }
+
+  virtual ~MultiGeneralCoarsenedOperator()
+  {
+    ReleaseGrid();
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Free everything SetGrid allocated. The D+1 grid is borrowed from the
+  // caller and is never deleted here. Safe to call repeatedly and before
+  // the destructor.
+  ///////////////////////////////////////////////////////////////////////////
+  void ReleaseGrid(void)
+  {
+    if ( CellMulti != nullptr ) { delete CellMulti; CellMulti = nullptr; }
+
+    _CoarseGridMulti = nullptr;  // borrowed, not owned
+    _Nrhs            = -1;
+
+    BLAS_B.resize(0);
+    BLAS_C.resize(0);
+    BLAS_Cg.resize(0);
+    BLAS_CPg.resize(0);
+    BLAS_CPr.resize(0);
+    BLAS_BK.resize(0);
+    BLAS_BPflat.resize(0);
+    for(int g=0;g<(int)BLAS_BKP.size();g++){ BLAS_BKP[g].resize(0); }
+    BLAS_BKP.resize(0);
+    for(int g=0;g<(int)BLAS_BPg.size();g++){ BLAS_BPg[g].resize(0); }
+    BLAS_BPg.resize(0);
+    for(int p=0;p<BLAS_BP.size();p++){
+      BLAS_BP[p].resize(0);
+    }
+    BLAS_BP.resize(0);
+    BLAS_CP.resize(0);
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Everything that depends on Nrhs. Idempotent; lazy called on demand.
+  //
+  // The stencil lives on the D dimensional padded grid. PaddedCell only pads
+  // a dimension when it is distributed, and the rhs direction never is, so
+  // the D+1 padded grid is exactly Nrhs copies of the D dimensional padded
+  // grid with rhs innermost. The neighbour offset therefore carries an Nrhs
+  // factor, in the same way the Nsimd factor is carried.
+  ///////////////////////////////////////////////////////////////////////////
+  void SetGrid(GridCartesian *CoarseGridMulti)
+  {
+    GRID_ASSERT(CoarseGridMulti != nullptr);
+
+    if ( CoarseGridMulti == _CoarseGridMulti ) return; // idempotent on identity
+
+    ReleaseGrid();
+
+    /////////////////////////////////////////////////
+    // The D+1 grid is supplied and owned by the caller. Two operators over
+    // the same coarse space must share one grid object or their fields will
+    // not conform, so this is never manufactured internally.
+    /////////////////////////////////////////////////
+    int nd = _CoarseGrid->_ndimension;
+
+    GRID_ASSERT(CoarseGridMulti->_ndimension == nd+1);
+    GRID_ASSERT(CoarseGridMulti->_processors[0] == 1);   // rhs is not distributed
+    for(int d=0;d<nd;d++){
+      GRID_ASSERT(CoarseGridMulti->_fdimensions[d+1] == _CoarseGrid->_fdimensions[d]);
+      GRID_ASSERT(CoarseGridMulti->_processors [d+1] == _CoarseGrid->_processors [d]);
+      GRID_ASSERT(CoarseGridMulti->_simd_layout[d+1] == _CoarseGrid->_simd_layout[d]);
+    }
+
+    _CoarseGridMulti = CoarseGridMulti;
+    _Nrhs            = CoarseGridMulti->_fdimensions[0];
+    GRID_ASSERT(_Nrhs>=1);
+
+    int nrhs = _Nrhs;
+
+    CellMulti = new PaddedCell(geom.Depth(),_CoarseGridMulti);
+
+    int32_t padded_sites   = CellD.grids.back()->lSites(); // D dimensional
+    int32_t unpadded_sites = _CoarseGrid->lSites();        // D dimensional
+
+    // The neighbour offset multiplication by nrhs is exact only if the D+1
+    // padded volume is nrhs copies of the D dimensional one. Check it.
+    GRID_ASSERT(CellMulti->grids.back()->lSites() == nrhs*padded_sites);
+    GRID_ASSERT(_CoarseGridMulti->lSites()        == nrhs*unpadded_sites);
+
+    /////////////////////////////////////////////////
+    // Choose the organisation FIRST: every table below is sized from
+    // LegGroup or KFold, and is indexed with them in the apply, so the two
+    // must be decided before anything is allocated.
+    /////////////////////////////////////////////////
+    ChooseOrganisation(nrhs);
+    std::cout << GridLogMessage << "MultiGeneralCoarsenedOperator: Nrhs " << nrhs
+	      << ", coarse apply "
+	      << ( KFold ? "folds "+std::to_string(KFold)+" legs along K (K = "
+		           +std::to_string(KFold*nbasis)+", "+std::to_string(KChunks())+" chunk(s))"
+		         : "carries "+std::to_string(LegGroup)+" legs in the batch" )
+	      << std::endl;
+
+    /////////////////////////////////////////////////
+    // Device data vector storage
+    /////////////////////////////////////////////////
+    BLAS_B.resize(nrhs *padded_sites);   // includes ghost zone
+    BLAS_C.resize(nrhs *unpadded_sites); // no ghost zone
+    BLAS_BP.resize(geom.npoint);
+    for(int p=0;p<geom.npoint;p++){
+      BLAS_BP[p].resize(unpadded_sites);
+    }
+    BLAS_CP.resize(unpadded_sites);
+
+    // Site identity mapping for C
+    for(int ss=0;ss<unpadded_sites;ss++){
+      CoarseBLASScalar *ptr = (CoarseBLASScalar *)&BLAS_C[ss*nrhs];
+      acceleratorPut(BLAS_CP[ss],ptr);
+    }
+
+    // Grouped output: LegGroup partial results, each the full C volume, and
+    // their pointer table.  The batch index runs site fastest within a leg,
+    // matching the grouped A and B tables.  The remainder call writes into
+    // the first few partials, so it needs a truncated copy of the table.
+    int nfull = LegGroups();
+    int rem   = LegRemainder();
+    int ng    = nfull + (rem?1:0);
+    BLAS_Cg.resize (LegGroup*nrhs*unpadded_sites);
+    BLAS_CPg.resize(LegGroup*unpadded_sites);
+    BLAS_CPr.resize(rem*unpadded_sites);
+    for(int l=0;l<LegGroup;l++){
+      for(int ss=0;ss<unpadded_sites;ss++){
+	CoarseBLASScalar *ptr = (CoarseBLASScalar *)&BLAS_Cg[(l*unpadded_sites+ss)*nrhs];
+	acceleratorPut(BLAS_CPg[l*unpadded_sites+ss],ptr);
+	if ( l<rem ) acceleratorPut(BLAS_CPr[l*unpadded_sites+ss],ptr);
+      }
+    }
+    BLAS_BPg.resize(ng);
+    for(int g=0;g<ng;g++){
+      int legs = (g<nfull) ? LegGroup : rem;
+      BLAS_BPg[g].resize(legs*unpadded_sites);
+    }
+
+    // Folded-along-K working set: one chunk of gathered neighbour vectors,
+    // laid out per site as a column-major (legs*nbasis) x nrhs block, plus a
+    // pointer table per chunk because a short last chunk has a shorter site
+    // stride.  BLAS_BPflat is the neighbour table the gather kernel reads.
+    if ( KFold ) {
+      BLAS_BK.resize((uint64_t)unpadded_sites*nrhs*KFold*nbasis);
+      BLAS_BKP.resize(KChunks());
+      for(int g=0;g<KChunks();g++){
+	int legs = KLegs(g);
+	BLAS_BKP[g].resize(unpadded_sites);
+	for(int32_t ss=0;ss<unpadded_sites;ss++){
+	  CoarseBLASScalar *ptr = &BLAS_BK[(uint64_t)ss*nrhs*legs*nbasis];
+	  acceleratorPut(BLAS_BKP[g][ss],ptr);
+	}
+      }
+      BLAS_BPflat.resize((uint64_t)geom.npoint*unpadded_sites);
+    }
+
+    // Neighbour table is more complicated
+    int32_t j=0; // Interior point counter (unpadded)
+    for(int32_t s=0;s<padded_sites;s++){ // D volume, padded
+      int ghost_zone=0;
+      for(int32_t point = 0 ; point < geom.npoint; point++){
+	int i=s*geom.npoint+point;
+	if( Stencil._entries[i]._wrap ) { // stencil is indexed by the oSite of the D dim grid
+	  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*geom.npoint+point;
+ 	  int32_t nbr = Stencil._entries[i]._offset*CComplex::Nsimd(); // oSite -> lSite, D dim
+	  nbr = nbr*nrhs;                                              // D -> D+1, rhs innermost
+	  GRID_ASSERT(nbr<BLAS_B.size());
+	  CoarseBLASScalar * ptr = (CoarseBLASScalar *)&BLAS_B[nbr];
+	  acceleratorPut(BLAS_BP[point][j],ptr); // neighbour indexing in ghost zone volume
+	  acceleratorPut(BLAS_BPg[point/LegGroup][(point%LegGroup)*unpadded_sites+j],ptr);
+	  if ( KFold ) acceleratorPut(BLAS_BPflat[(uint64_t)point*unpadded_sites+j],ptr);
+	}
+	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);
+      }
+    });
+  }
+  ///////////////////////////////////////////////////////////////////////////
+  // Shared by both CoarsenOperator variants
+  //
+  //     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}
+  //
+  //     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}
+  ///////////////////////////////////////////////////////////////////////////
+  ///////////////////////////////////////////////////////////////////////////
+  // Probe momenta.  The stencil DISPLACEMENTS are the +-1 box; the momenta
+  // used to separate them are free, and are chosen here to span the
+  // Brillouin zone: k_mu = K_mu * shift_mu with K_mu ~ L_mu/3, so a phase
+  // is ~2pi/3 per unit displacement rather than 2pi/L_mu.
+  //
+  // This is what conditions the extraction.  With K_mu = 1 every entry of
+  // the phase matrix tends to 1 as the coarse lattice grows, the matrix
+  // tends to rank one, and its inverse amplifies any error in the measured
+  // projections: at 24.24.16.32 the condition number is 2.7e5, so fp32
+  // projections give coarse matrix elements wrong by ~1e-3.  Spread momenta
+  // bring it to ~20, independent of the lattice size.
+  //
+  // K_mu is stepped away from L_mu/2, where +k and -k alias into the same
+  // momentum and the matrix is singular.
+  ///////////////////////////////////////////////////////////////////////////
+  void CoarsenMomenta(GridBase *CoarseGrid,Coordinate &K)
+  {
+    Coordinate clatt = CoarseGrid->GlobalDimensions();
+    int Nd = CoarseGrid->Nd();
+    K.resize(Nd);
+    for(int mu=0;mu<Nd;mu++){
+      int L = clatt[mu];
+      int k = (L>=3) ? (int)std::lround(L/3.0) : 1;
+      if ( k < 1 ) k = 1;
+      if ( (L>2) && ((2*k)%L == 0) ) k = k-1;     // +k and -k must differ
+      if ( k < 1 ) k = 1;
+      K[mu] = k;
+    }
+  }
+
+  void CoarsenFourierMatrix(GridBase *CoarseGrid,Eigen::MatrixXcd &invMkl)
+  {
+    const int npoint = geom_srhs.npoint;
+    Coordinate clatt = CoarseGrid->GlobalDimensions();
+    int Nd = CoarseGrid->Nd();
+    Coordinate K;  CoarsenMomenta(CoarseGrid,K);
+
+    Eigen::MatrixXcd Mkl = 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*K[mu]*geom_srhs.shifts[k][mu]*geom_srhs.shifts[l][mu];
+	}
+	phase=exp(phase*ci);
+	Mkl(k,l) = phase;
+      }
+    }
+    invMkl = Mkl.inverse();
+
+    // The extraction's error amplification, printed because it is the thing
+    // that decides whether an fp32 coarse sector is usable.
+    Eigen::JacobiSVD<Eigen::MatrixXcd> svd(Mkl);
+    RealD cond = svd.singularValues()(0)/svd.singularValues()(npoint-1);
+    std::cout << GridLogMessage << "CoarsenOperator: probe momenta "<<K
+	      <<" on coarse lattice "<<clatt
+	      <<" ; Fourier matrix condition number "<<cond<<std::endl;
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // blockOrthogonalise and blockZAXPY are block operations on the fine
+  // vectors, using a coarse shaped field only as an index set. They need a
+  // grid carrying the fine SIMD layout, which the coarse space no longer
+  // does. Constructed local to the caller so it cannot be mistaken for the
+  // coarse grid.
+  ///////////////////////////////////////////////////////////////////////////
+  void CoarsenBlockGridLayout(GridBase *grid,GridBase *CoarseGrid,
+			      Coordinate &latt,Coordinate &simd,Coordinate &mpi)
+  {
+    int nd = CoarseGrid->_ndimension;
+    latt.resize(nd); simd.resize(nd); mpi.resize(nd);
+    for(int d=0;d<nd;d++){
+      latt[d] = CoarseGrid->_fdimensions[d];
+      simd[d] = grid->_simd_layout[d];
+      mpi [d] = CoarseGrid->_processors[d];
+    }
+  }
+
+  // D+1 coarse grid holding the batch, rhs innermost and unvectorised
+  void CoarsenBatchGridLayout(GridBase *CoarseGrid,int batch,
+			      Coordinate &latt,Coordinate &simd,Coordinate &mpi)
+  {
+    latt.resize(1,batch); simd.resize(1,1); mpi.resize(1,1);
+    latt[0]=batch; simd[0]=1; mpi[0]=1;
+    for(int d=0;d<CoarseGrid->_ndimension;d++){
+      latt.push_back(CoarseGrid->_fdimensions[d]);
+      simd.push_back(CoarseGrid->_simd_layout[d]);
+      mpi .push_back(CoarseGrid->_processors[d]);
+    }
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // The Fourier inverse needs the phase in the coarse layout and the basis
+  // phasing needs it in the fine layout; each is built from its own
+  // coordinates rather than transferred.
+  ///////////////////////////////////////////////////////////////////////////
+  void CoarsenPhases(GridBase *grid,GridBase *CoarseGrid,GridCartesian *BlockGrid,
+		     std::vector<CoarseComplexField> &pha,
+		     std::vector<FineComplexField> &phaF)
+  {
+    const int npoint = geom_srhs.npoint;
+    Coordinate clatt = CoarseGrid->GlobalDimensions();
+    int Nd = CoarseGrid->Nd();
+    Coordinate K;  CoarsenMomenta(CoarseGrid,K);
+    ComplexD ci(0.0,1.0);
+
+    // The fine-side scratch carries the FINE level's precision, which need
+    // not be the coarse one (fp64 fine, fp32 coarse at L1).
+    typedef typename GridTypeMapper<FineInner>::scalar_type FineScalar;
+    FineComplexField one(grid); one=FineScalar(1.0);
+    FineComplexField zz(grid);  zz = Zero();
+    BlockComplexField pha_blk (BlockGrid);
+    BlockComplexField blk_coor(BlockGrid);
+
+    for(int p=0;p<npoint;p++){ // Loop over momenta in npoint
+      CoarseComplexField coor(CoarseGrid);
+      pha[p] =Zero();
+      pha_blk=Zero();
+      for(int mu=0;mu<Nd;mu++){
+	RealD TwoPiL =  M_PI * 2.0/ clatt[mu];
+	LatticeCoordinate(coor,mu);
+	pha[p]  = pha[p]  + (TwoPiL * K[mu] * geom_srhs.shifts[p][mu]) * coor;
+	LatticeCoordinate(blk_coor,mu);
+	pha_blk = pha_blk + (TwoPiL * K[mu] * geom_srhs.shifts[p][mu]) * blk_coor;
+      }
+      pha[p] =exp(pha[p] *ci);
+      pha_blk=exp(pha_blk*ci);
+
+      blockZAXPY(phaF[p],pha_blk,one,zz);
+    }
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Remove the bulk phase from the batch of coarse projections and
+  // accumulate the Fourier inverse into A. Both variants reach here with
+  // TmpProj in the same batch coarse order, so this is shared verbatim.
+  ///////////////////////////////////////////////////////////////////////////
+  void CoarsenAccumulate(int p,int i0,int nbv,int batch,
+			 Eigen::MatrixXcd &invMkl,
+			 std::vector<CoarseComplexField> &pha,
+			 CoarseComplexField &phaB,
+			 CoarseVector &TmpProj,
+			 std::vector<CoarseMatrix> &_A,
+			 GridBase *CoarseGrid)
+  {
+    typedef typename CComplex::scalar_type SComplex;
+    const int npoint = geom_srhs.npoint;
+
+    for(int b=0;b<batch;b++) InsertSliceFast(pha[p],phaB,b,0);
+    TmpProj = conjugate(phaB)*TmpProj;
+
+    int osites=CoarseGrid->oSites();
+    for(int k=0;k<npoint;k++){
+      SComplex sc(invMkl(p,k).real(),invMkl(p,k).imag());
+      CComplex coef(sc);
+      autoView( A_v  , _A[k], AcceleratorWrite);
+      autoView( TP_v , TmpProj, AcceleratorRead);
+      accelerator_for(sss, osites, 1, {
+	  for(int b=0;b<nbv;b++){
+	    for(int j=0;j<nbasis;j++){
+	      A_v[sss](i0+b,j) = A_v[sss](i0+b,j) + coef*TP_v[b+batch*sss](j);
+	    }
+	  }
+      });
+    }
+  }
+
+  void CoarsenReport(RealD tphase,RealD tphaseBZ,RealD tslice,
+		     RealD tmat,RealD tproj,RealD tinv)
+  {
+    std::cout << GridLogMessage<<"CoarsenOperator phase  "<<tphase<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator phaseBZ "<<tphaseBZ<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator slice  "<<tslice <<" 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;
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Coarsen a NATIVELY multiRHS fine operator.
+  //
+  // linop acts on the D+1 dimensional fine grid FineGridMulti, with the batch
+  // of phased basis vectors carried in the rhs direction. A single RHS
+  // operator can be promoted with MrhsPromotedOperator, but that pays an
+  // ExtractSlice/InsertSlice pair per rhs; prefer the single RHS variant
+  // below in that case.
+  ///////////////////////////////////////////////////////////////////////////
+  // Owning the transfer operator: the caller passes one, this imports the
+  // orthonormalised basis into it, and the caller keeps it for the solve.
+  // The overload below makes its own and is kept for the pre-2026 drivers;
+  // it costs a second full basis store, which is why nothing new should use
+  // it.
+  void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
+		       GridCartesian *FineGridMulti,
+		       std::vector<FineField> &Subspace,
+		       GridBase *CoarseGrid)
+  {
+    MultiRHSBlockProject<Lattice<Fobj> > Projector;
+    CoarsenOperator(linop,FineGridMulti,Subspace,CoarseGrid,Projector);
+  }
+  template<class Projector_t>
+  void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
+		       GridCartesian *FineGridMulti,
+		       std::vector<FineField> &Subspace,
+		       GridBase *CoarseGrid,
+		       Projector_t &Projector)
+  {
+    RealD tproj=0.0, tmat=0.0, tphase=0.0, tphaseBZ=0.0, tslice=0.0, tinv=0.0;
+
+    std::cout << GridLogMessage<< "GeneralCoarsenMatrixMrhs (multiRHS fine operator)"<< std::endl;
+
+    GRID_ASSERT(Subspace.size()==nbasis);
+    GridBase *grid = Subspace[0].Grid();
+
+    GRID_ASSERT(FineGridMulti->_ndimension == grid->_ndimension+1);
+    GRID_ASSERT(FineGridMulti->_processors[0] == 1);
+    for(int d=0;d<grid->_ndimension;d++){
+      GRID_ASSERT(FineGridMulti->_fdimensions[d+1] == grid->_fdimensions[d]);
+      GRID_ASSERT(FineGridMulti->_processors [d+1] == grid->_processors [d]);
+    }
+    int batch = FineGridMulti->_fdimensions[0];
+
+    Coordinate blatt,bsimd,bmpi;
+    CoarsenBlockGridLayout(grid,CoarseGrid,blatt,bsimd,bmpi);
+    GridCartesian BlockGrid(blatt,bsimd,bmpi);
+
+    BlockComplexField InnerProd(&BlockGrid);
+    blockOrthogonalise(InnerProd,Subspace);
+
+    // the caller owns it; import the basis we have just orthonormalised
+    Projector.Allocate(nbasis,grid,CoarseGrid);
+    Projector.ImportBasis(Subspace);
+
+    const int npoint = geom_srhs.npoint;
+
+    Eigen::MatrixXcd invMkl;
+    CoarsenFourierMatrix(CoarseGrid,invMkl);
+
+    FineField phaV(grid);
+    std::vector<FineComplexField>   phaF(npoint,grid);
+    std::vector<CoarseComplexField> pha (npoint,CoarseGrid);
+
+    tphase=-usecond();
+    CoarsenPhases(grid,CoarseGrid,&BlockGrid,pha,phaF);
+    tphase+=usecond();
+
+    std::vector<CoarseMatrix> _A;
+    _A.resize(npoint,CoarseGrid);
+    for(int k=0;k<npoint;k++) _A[k] = Zero();
+
+    Coordinate cmlatt,cmsimd,cmmpi;
+    CoarsenBatchGridLayout(CoarseGrid,batch,cmlatt,cmsimd,cmmpi);
+    GridCartesian CoarseBatchGrid(cmlatt,cmsimd,cmmpi);
+
+    CoarseVector       TmpProj(&CoarseBatchGrid);
+    CoarseComplexField phaB(&CoarseBatchGrid);
+
+    FineField hi_in (FineGridMulti);
+    FineField hi_out(FineGridMulti);
+    FineField zzF(grid); zzF = Zero();
+
+    for(int i0=0;i0<nbasis;i0+=batch){ // Loop over batches of basis vectors
+
+      int nbv = MIN(batch,nbasis-i0);
+      std::cout << GridLogMessage<< "CoarsenMatrixColoured vec "<<i0<<"/"<<nbasis<< std::endl;
+
+      for(int p=0;p<npoint;p++){ // Loop over momenta
+
+	// One phase, applied to the whole batch. Tail slices are zeroed so
+	// the operator never sees undefined data.
+	for(int b=0;b<nbv;b++){
+	  tphaseBZ-=usecond();
+	  phaV = phaF[p]*Subspace[i0+b];
+	  tphaseBZ+=usecond();
+	  tslice-=usecond();
+	  InsertSliceFast(phaV,hi_in,b,0);
+	  tslice+=usecond();
+	}
+	tslice-=usecond();
+	for(int b=nbv;b<batch;b++){
+	  InsertSliceFast(zzF,hi_in,b,0);
+	}
+	tslice+=usecond();
+
+	tmat-=usecond();
+	linop.Op(hi_in,hi_out);
+	tmat+=usecond();
+
+	tproj-=usecond();
+	Projector.blockProject(hi_out,TmpProj);
+	tproj+=usecond();
+
+	tinv-=usecond();
+	CoarsenAccumulate(p,i0,nbv,batch,invMkl,pha,phaB,TmpProj,_A,CoarseGrid);
+	tinv+=usecond();
+      }
+    }
+
+    for(int p=0;p<npoint;p++){
+      GridtoBLAS(_A[p],BLAS_Ascr);
+      MatrixPointIn(p,BLAS_Ascr);
+    }
+    if ( KFold ) PackFoldedMatrix();   // the matrix is only valid from here
+    CoarsenReport(tphase,tphaseBZ,tslice,tmat,tproj,tinv);
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Coarsen a SINGLE RHS fine operator.
+  //
+  // No multiRHS packing: the operator is applied once per phased basis
+  // vector and the batch is assembled on the coarse side by the mixed
+  // blockProject, which takes a vector of fine fields and writes the batch
+  // coarse field the accumulate expects. Only nbv applications per momentum,
+  // so a batch that does not divide nbasis wastes nothing, and the live fine
+  // storage is batch fields rather than two D+1 fields of extent batch.
+  ///////////////////////////////////////////////////////////////////////////
+  void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
+		       std::vector<FineField> &Subspace,
+		       GridBase *CoarseGrid,
+		       int batch)
+  {
+    MultiRHSBlockProject<Lattice<Fobj> > Projector;
+    CoarsenOperator(linop,Subspace,CoarseGrid,batch,Projector);
+  }
+  template<class Projector_t>
+  void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
+		       std::vector<FineField> &Subspace,
+		       GridBase *CoarseGrid,
+		       int batch,
+		       Projector_t &Projector)
+  {
+    RealD tproj=0.0, tmat=0.0, tphase=0.0, tphaseBZ=0.0, tslice=0.0, tinv=0.0;
+
+    std::cout << GridLogMessage<< "GeneralCoarsenMatrixMrhs (single RHS fine operator)"<< std::endl;
+
+    GRID_ASSERT(Subspace.size()==nbasis);
+    GRID_ASSERT(batch>=1);
+    GridBase *grid = Subspace[0].Grid();
+
+    Coordinate blatt,bsimd,bmpi;
+    CoarsenBlockGridLayout(grid,CoarseGrid,blatt,bsimd,bmpi);
+    GridCartesian BlockGrid(blatt,bsimd,bmpi);
+
+    BlockComplexField InnerProd(&BlockGrid);
+    blockOrthogonalise(InnerProd,Subspace);
+
+    // the caller owns it; import the basis we have just orthonormalised
+    Projector.Allocate(nbasis,grid,CoarseGrid);
+    Projector.ImportBasis(Subspace);
+
+    const int npoint = geom_srhs.npoint;
+
+    Eigen::MatrixXcd invMkl;
+    CoarsenFourierMatrix(CoarseGrid,invMkl);
+
+    FineField phaV(grid);
+    std::vector<FineComplexField>   phaF(npoint,grid);
+    std::vector<CoarseComplexField> pha (npoint,CoarseGrid);
+
+    tphase=-usecond();
+    CoarsenPhases(grid,CoarseGrid,&BlockGrid,pha,phaF);
+    tphase+=usecond();
+
+    std::vector<CoarseMatrix> _A;
+    _A.resize(npoint,CoarseGrid);
+    for(int k=0;k<npoint;k++) _A[k] = Zero();
+
+    Coordinate cmlatt,cmsimd,cmmpi;
+    CoarsenBatchGridLayout(CoarseGrid,batch,cmlatt,cmsimd,cmmpi);
+    GridCartesian CoarseBatchGrid(cmlatt,cmsimd,cmmpi);
+
+    CoarseVector       TmpProj(&CoarseBatchGrid);
+    CoarseComplexField phaB(&CoarseBatchGrid);
+
+    std::vector<FineField> MphaV(batch,grid);
+
+    for(int i0=0;i0<nbasis;i0+=batch){ // Loop over batches of basis vectors
+
+      int nbv = MIN(batch,nbasis-i0);
+      std::cout << GridLogMessage<< "CoarsenMatrixColoured vec "<<i0<<"/"<<nbasis<< std::endl;
+
+      for(int p=0;p<npoint;p++){ // Loop over momenta
+
+	for(int b=0;b<nbv;b++){
+	  tphaseBZ-=usecond();
+	  phaV = phaF[p]*Subspace[i0+b];
+	  tphaseBZ+=usecond();
+	  tmat-=usecond();
+	  linop.Op(phaV,MphaV[b]);
+	  tmat+=usecond();
+	}
+	// The accumulate reads only the first nbv slices, but the projector
+	// sees the whole vector, so the tail must not be undefined.
+	for(int b=nbv;b<batch;b++) MphaV[b] = Zero();
+
+	tproj-=usecond();
+	Projector.blockProject(MphaV,TmpProj);
+	tproj+=usecond();
+
+	tinv-=usecond();
+	CoarsenAccumulate(p,i0,nbv,batch,invMkl,pha,phaB,TmpProj,_A,CoarseGrid);
+	tinv+=usecond();
+      }
+    }
+
+    for(int p=0;p<npoint;p++){
+      GridtoBLAS(_A[p],BLAS_Ascr);
+      MatrixPointIn(p,BLAS_Ascr);
+    }
+    if ( KFold ) PackFoldedMatrix();   // the matrix is only valid from here
+    CoarsenReport(tphase,tphaseBZ,tslice,tmat,tproj,tinv);
+  }
+  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;
+
+    CheckGridSet();
+    if ( in.Grid() != _CoarseGridMulti ) {
+      std::cout << GridLogError
+		<< "MultiGeneralCoarsenedOperator::M called with a field on a"
+		<< std::endl;
+      std::cout << GridLogError
+		<< "  different grid object from the one given to SetGrid(). Two"
+		<< std::endl;
+      std::cout << GridLogError
+		<< "  grids of identical shape do not conform; share one object."
+		<< std::endl;
+      GRID_ASSERT(in.Grid() == _CoarseGridMulti);
+    }
+
+    GRID_TRACE("CoarseMult");
+    t_tot=-usecond();
+    t_exch=-usecond();
+    // The exchange takes its input by const reference, so it reads the
+    // caller's field directly.  lambda scope so the roctx range covers
+    // exactly the exchange; the PaddedCellFwd/BwdMPI markers inside it then
+    // nest properly.
+    CoarseVector pin = [&](){ GRID_TRACE("CoarseExchange");
+                              return CellMulti->ExchangePeriodic(in); }(); //padded input
+    t_exch+=usecond();
+
+    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();
+    { GRID_TRACE("CoarseGridToBLAS");
+      GridtoBLAS(pin,BLAS_B);
+    }
+    t_GtoB+=usecond();
+
+    GridBLAS BLAS;
+
+    t_mult=-usecond();
+    { GRID_TRACE("CoarseStencilGEMM");
+    // The scalar type selects the GEMM: Cgemm for an fp32 coarse space,
+    // Zgemm for fp64.
+    if ( KFold ) {
+      // Gather one chunk of neighbour vectors per site, then one GEMM whose
+      // contraction runs over the chunk's legs, accumulating into C.  No
+      // partial results, so nothing to combine afterwards.
+      const int    nb  = nbasis;
+      const int64_t ns = (int64_t)_CoarseGrid->lSites();  // D dimensional: CoarseGrid() is D+1
+      GRID_ASSERT( (int64_t)BLAS_BPflat.size() == (int64_t)geom.npoint*ns );
+      CoarseBLASScalar * const bk = &BLAS_BK[0];
+      CoarseBLASScalar * const * const bp = &BLAS_BPflat[0];
+      for(int g=0;g<KChunks();g++){
+	const int64_t legs = KLegs(g);
+	const int64_t leg0 = (int64_t)g*KFold;
+	{ GRID_TRACE("CoarseLegGather");
+	  accelerator_for(idx, ns*legs*nrhs, 1, {
+	      int64_t j  = idx % nrhs;          // right hand side  (column)
+	      int64_t t  = idx / nrhs;
+	      int64_t l  = t % legs;            // leg within the chunk
+	      int64_t ss = t / legs;            // site
+	      const CoarseBLASScalar *src = bp[(leg0+l)*ns + ss] + j*nb;
+	      CoarseBLASScalar       *dst = bk + ((ss*nrhs + j)*legs + l)*nb;
+	      for(int i=0;i<nb;i++) dst[i] = src[i];
+	    });
+	}
+	RealD c = (g==0) ? 0.0 : 1.0;
+	BLAS.gemmBatched(nbasis,nrhs,legs*nbasis,
+			 CoarseBLASScalar(1.0),
+			 BLAS_AKP[g],
+			 BLAS_BKP[g],
+			 CoarseBLASScalar(c),
+			 BLAS_CP);
+	// One gather buffer serves every chunk and the multiply is issued
+	// asynchronously, so the next chunk's gather must not start until this
+	// multiply has read it.  One synchronise per chunk, none when the whole
+	// stencil folds into a single call.
+	if ( g+1 < KChunks() ) BLAS.synchronise();
+      }
+    } else if ( LegGroup == 1 ) {
+      for(int p=0;p<geom.npoint;p++){
+	RealD c = (p==0) ? 0.0 : 1.0;
+	BLAS.gemmBatched(nbasis,nrhs,nbasis,
+			 CoarseBLASScalar(1.0),
+			 BLAS_AP[p],
+			 BLAS_BP[p],
+			 CoarseBLASScalar(c),
+			 BLAS_CP);
+      }
+    } else {
+      // LegGroup legs per call: the batch is LegGroup times the local volume
+      // and there are LegGroup partial outputs, accumulated across the
+      // groups and summed below.  A final short call carries the remainder
+      // when the group does not divide the stencil; it writes into the first
+      // few partials, which the full calls have already initialised.
+      int nfull = LegGroups();
+      int rem   = LegRemainder();
+      for(int g=0;g<nfull;g++){
+	RealD c = (g==0) ? 0.0 : 1.0;
+	BLAS.gemmBatched(nbasis,nrhs,nbasis,
+			 CoarseBLASScalar(1.0),
+			 BLAS_APg[g],
+			 BLAS_BPg[g],
+			 CoarseBLASScalar(c),
+			 BLAS_CPg);
+      }
+      if ( rem ) {
+	RealD c = (nfull==0) ? 0.0 : 1.0;
+	BLAS.gemmBatched(nbasis,nrhs,nbasis,
+			 CoarseBLASScalar(1.0),
+			 BLAS_APg[nfull],
+			 BLAS_BPg[nfull],
+			 CoarseBLASScalar(c),
+			 BLAS_CPr);
+      }
+    }
+    BLAS.synchronise();
+    }
+    t_mult+=usecond();
+
+    if ( !KFold && LegGroup > 1 ) { GRID_TRACE("CoarseLegSum");
+      // Sum the LegGroup partial results into the single output buffer.
+      // Flat over scalars: a calcVector is nbasis of them and the partials
+      // are contiguous, one whole C volume after another.
+      int64_t nscalar = (int64_t)BLAS_C.size()*nbasis;   // one whole C volume
+      const int G = LegGroup;
+      CoarseBLASScalar *dst = (CoarseBLASScalar *)&BLAS_C[0];
+      CoarseBLASScalar *src = (CoarseBLASScalar *)&BLAS_Cg[0];
+      accelerator_for(i,nscalar,1,{
+	  CoarseBLASScalar sum = src[i];
+	  for(int l=1;l<G;l++) sum = sum + src[(int64_t)l*nscalar + i];
+	  dst[i] = sum;
+	});
+    }
+
+    t_BtoG=-usecond();
+    { GRID_TRACE("CoarseBLASToGrid");
+      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/MPI_benchmark/simple_reproducer_noncache.cc b/MPI_benchmark/simple_reproducer_noncache.cc
new file mode 100644
index 000000000..5c3aa37ad
--- /dev/null
+++ b/MPI_benchmark/simple_reproducer_noncache.cc
@@ -0,0 +1,214 @@
+// Minimal reproducer for the CXI uncached-registration defect.
+//
+// Hypothesis: with the libfabric MR cache disabled, a device pointer at a non-zero offset
+// into a hipMalloc'd allocation is registered as the whole allocation and the offset is
+// discarded, so a send transmits from the base and a receive lands at the base.
+//
+// Each test therefore prints, for both the data and its landing place, what a correct
+// implementation must give, what the hypothesis predicts instead, and what was observed.
+//
+//   hipcc -O2 -o simple_reproducer_noncache simple_reproducer_noncache.cc \
+//         -I$MPICH_DIR/include -L$MPICH_DIR/lib -lmpi -L$MPICH_DIR/gtl/lib -lmpi_gtl_hsa
+//
+//   MPICH_GPU_SUPPORT_ENABLED=1 FI_MR_CACHE_MAX_COUNT=0 \
+//     srun -N2 -n2 --ntasks-per-node=1 ./simple_reproducer_noncache
+
+#include <mpi.h>
+#include <hip/hip_runtime.h>
+#include <cstdio>
+#include <cstdint>
+#include <cstdlib>
+#include <vector>
+
+// The message must be large enough to go rendezvous: an eager message is staged through a host
+// bounce buffer with hipMemcpy, which honours the offset, and never registers the user pointer.
+static const int NWORDS = 262144; // 2 MB buffers
+static int       NMSG   = 16384;  // 128 KB messages          (argv[1], in words)
+static int       OFF    = 16384;  // offending offset in words (argv[2]); clearest at >= NMSG,
+                                  // which keeps the intended and predicted windows disjoint
+
+static int rank, size, peer;
+static uint64_t *A, *B;
+static std::vector<uint64_t> h(NWORDS);
+
+#define HIP(cmd) do { hipError_t e=(cmd); if ( e != hipSuccess ) {			\
+      printf("hip error %s at line %d\n",hipGetErrorString(e),__LINE__);		\
+      MPI_Abort(MPI_COMM_WORLD,1); } } while(0)
+
+static uint64_t word(uint64_t tag,int r,int n)
+{
+  return (tag<<60) | ((uint64_t)r<<32) | (uint64_t)n;
+}
+
+static const char *verdict(int observed,int correct,int predicted)
+{
+  if ( observed == correct   ) return "CORRECT";
+  if ( observed == predicted ) return "WRONG, as predicted";
+  return "WRONG, and not as predicted";
+}
+
+static void refresh(void)
+{
+  for(int n=0;n<NWORDS;n++) h[n] = word(0xa,rank,n);
+  HIP(hipMemcpy(A,h.data(),NWORDS*sizeof(uint64_t),hipMemcpyHostToDevice));
+  for(int n=0;n<NWORDS;n++) h[n] = word(0xb,rank,n);
+  HIP(hipMemcpy(B,h.data(),NWORDS*sizeof(uint64_t),hipMemcpyHostToDevice));
+  MPI_Barrier(MPI_COMM_WORLD);
+}
+
+static void exchange(int sendoff,int recvoff)
+{
+  MPI_Request req[2];
+  MPI_Irecv(B+recvoff,NMSG*sizeof(uint64_t),MPI_BYTE,peer,0,MPI_COMM_WORLD,&req[0]);
+  MPI_Isend(A+sendoff,NMSG*sizeof(uint64_t),MPI_BYTE,peer,0,MPI_COMM_WORLD,&req[1]);
+  MPI_Waitall(2,req,MPI_STATUSES_IGNORE);
+}
+
+// Report the state of one NMSG-word window of B: how much of it was overwritten, and by what.
+static void window(const char *label,int at)
+{
+  int mod=0, m0=-1;
+  for(int i=0;i<NMSG && at+i<NWORDS;i++) {
+    if ( h[at+i] != word(0xb,rank,at+i) ) { if ( m0 < 0 ) m0 = at+i; mod++; }
+  }
+  if ( mod == 0 ) {
+    printf("rank %d    %-9s B[%d..%d] : untouched\n",rank,label,at,at+NMSG-1);
+    return;
+  }
+  uint64_t w = h[m0];
+  printf("rank %d    %-9s B[%d..%d] : %d of %d words overwritten, B[%d] holds A[%d] from rank %d\n",
+	 rank,label,at,at+NMSG-1,mod,NMSG,m0,(int)(w & 0xffffffff),(int)((w>>32) & 0xfffffff));
+}
+
+// Name a word: still the receive pattern, or a word of the neighbour's send buffer.
+static const char *describe(uint64_t w,int idx)
+{
+  static char s[64];
+  if      ( w == word(0xb,rank,idx) ) snprintf(s,sizeof(s),"untouched");
+  else if ( (w>>60) == 0xa )          snprintf(s,sizeof(s),"A[%d] of rank %d",
+					       (int)(w & 0xffffffff),(int)((w>>32) & 0xfffffff));
+  else                                snprintf(s,sizeof(s),"unrecognised");
+  return s;
+}
+
+static void report(const char *name,int sendoff,int recvoff)
+{
+  HIP(hipMemcpy(h.data(),B,NWORDS*sizeof(uint64_t),hipMemcpyDeviceToHost));
+
+  int first=-1, last=-1;
+  for(int n=0;n<NWORDS;n++) {
+    if ( h[n] != word(0xb,rank,n) ) { if ( first < 0 ) first = n; last = n; }
+  }
+
+  printf("rank %d  %s : send &A[%d] -> recv &B[%d]\n",rank,name,sendoff,recvoff);
+  if ( first < 0 ) { printf("rank %d    nothing was written\n",rank); return; }
+
+  uint64_t w   = h[first];
+  int      src = (int)(w & 0xffffffff);
+  int      sr  = (int)((w>>32) & 0xfffffff);
+
+  // The hypothesis predicts both offsets are discarded: the bytes at A[0], landing at B[0].
+  printf("rank %d    data      : want A[%d]  predict A[0]  got A[%d] from rank %d   %s\n",
+	 rank,sendoff,src,sr,verdict(src,sendoff,0));
+  printf("rank %d    place     : want B[%d]  predict B[0]  got B[%d]                %s\n",
+	 rank,recvoff,first,verdict(first,recvoff,0));
+
+  if ( src == sendoff && first == recvoff ) return;
+
+  // Both candidate landing sites in our receive buffer, and both candidate source words in the
+  // neighbour's send buffer, which the fill pattern fixes by construction.
+  printf("rank %d    Boff %d\n",rank,recvoff);
+  printf("rank %d      B[0]      = 0x%016llx  %s\n",
+	 rank,(unsigned long long)h[0],describe(h[0],0));
+  if ( recvoff != 0 )
+    printf("rank %d      B[%d]\t= 0x%016llx  %s\n",
+	   rank,recvoff,(unsigned long long)h[recvoff],describe(h[recvoff],recvoff));
+  printf("rank %d    Neighbour rank %d send buffer holds\n",rank,peer);
+  printf("rank %d      A[0]\t= 0x%016llx\n",
+	 rank,(unsigned long long)word(0xa,peer,0));
+  if ( sendoff != 0 )
+    printf("rank %d      A[%d]\t= 0x%016llx\n",
+	   rank,sendoff,(unsigned long long)word(0xa,peer,sendoff));
+
+  window("intended",recvoff);
+  if ( recvoff != 0 ) window("predicted",0);
+  if ( first < recvoff || last >= recvoff+NMSG ) {
+    if ( first != 0 || last != NMSG-1 ) {
+      printf("rank %d    stray     : modified words span B[%d..%d], outside both windows\n",
+	     rank,first,last);
+    }
+  }
+
+  int bad=0;
+  for(int i=0;i<NMSG && first+i<NWORDS;i++) {
+    if ( h[first+i] != word(0xa,peer,src+i) ) bad++;
+  }
+  if ( bad ) printf("rank %d    body      : %d of %d words are not a contiguous run from A[%d]\n",
+		    rank,bad,NMSG,src);
+}
+
+// Let each rank print in turn: flush, then a barrier so rank r is on the page before r+1 starts.
+static void ordered(void (*print)(void))
+{
+  for(int r=0;r<size;r++) {
+    if ( r == rank ) { print(); fflush(stdout); }
+    MPI_Barrier(MPI_COMM_WORLD);
+  }
+}
+
+static const char *test_name;
+static int         test_sendoff, test_recvoff;
+static void        print_report(void) { report(test_name,test_sendoff,test_recvoff); }
+static void        print_setup (void)
+{
+  printf("rank %d  A %p  B %p  message %d words (%d bytes)  offset %d words (%d bytes)\n",
+	 rank,(void *)A,(void *)B,NMSG,(int)(NMSG*sizeof(uint64_t)),OFF,(int)(OFF*sizeof(uint64_t)));
+}
+
+int main(int argc,char **argv)
+{
+  MPI_Init(&argc,&argv);
+  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+  MPI_Comm_size(MPI_COMM_WORLD,&size);
+  if ( size != 2 ) {
+    if ( rank == 0 ) printf("run with two ranks, one per node\n");
+    MPI_Finalize();
+    return 1;
+  }
+  peer = 1-rank;
+
+  if ( argc > 1 ) NMSG = atoi(argv[1]);
+  if ( argc > 2 ) OFF  = atoi(argv[2]);
+  if ( NMSG < 1 || OFF < 0 || OFF+NMSG > NWORDS ) {
+    if ( rank == 0 ) printf("message %d and offset %d words do not fit in %d\n",NMSG,OFF,NWORDS);
+    MPI_Finalize();
+    return 1;
+  }
+
+  int ndev;
+  HIP(hipGetDeviceCount(&ndev));
+  if ( ndev < 1 ) {
+    printf("rank %d sees no GPU\n",rank);
+    MPI_Abort(MPI_COMM_WORLD,1);
+  }
+  HIP(hipSetDevice(rank%ndev));
+  HIP(hipMalloc(&A,NWORDS*sizeof(uint64_t)));
+  HIP(hipMalloc(&B,NWORDS*sizeof(uint64_t)));
+
+  ordered(print_setup);
+
+  const int  offsets[4][2] = { {0,0}, {OFF,0}, {0,OFF}, {OFF,OFF} };
+  const char *names[4]     = { "test 1","test 2","test 3","test 4" };
+
+  for(int t=0;t<4;t++) {
+    refresh();
+    exchange(offsets[t][0],offsets[t][1]);
+    test_name = names[t]; test_sendoff = offsets[t][0]; test_recvoff = offsets[t][1];
+    ordered(print_report);
+  }
+
+  HIP(hipFree(A));
+  HIP(hipFree(B));
+  MPI_Finalize();
+  return 0;
+}
diff --git a/examples/Example_pvdagm_3level_DenseCoarseMatrix.cc b/examples/Example_pvdagm_3level_DenseCoarseMatrix.cc
new file mode 100644
index 000000000..7b6ce5b18
--- /dev/null
+++ b/examples/Example_pvdagm_3level_DenseCoarseMatrix.cc
@@ -0,0 +1,990 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./examples/Example_pvdagm_3level_DenseCoarseMatrix.cc
+
+    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.
+
+    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 */
+
+//
+// PVdagM three level multigrid on MultiGeneralCoarsenedOperator.
+//
+// STAGES ONE AND TWO: grids, types, subspace, and the L1 and L2 coarsenings.
+// The dense bottom and the solves are not here yet.
+//
+// Differences from Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc:
+//
+//  * The coarse space is UNVECTORISED (sComplexD). The fine space stays
+//    vectorised. MultiRHSBlockProject carries the mixed layout.
+//
+//  * One operator, not two. The deprecated path needed
+//    DeprecatedGeneralCoarsenedMatrix to coarsen and
+//    DeprecatedMultiGeneralCoarsenedMatrix to apply, bridged by CopyMatrix.
+//    This operator does both, and single versus multiRHS is SetGrid on the
+//    same object with the matrix elements built once.
+//
+//  * Nrhs is unconstrained. The deprecated path required nrhs % vComplex::Nsimd() == 0 because
+//    its multiRHS grid carried the SIMD in the rhs direction.
+//
+//  * CoarsenOperator takes the subspace vectors, not an Aggregation. It block
+//    orthonormalises them IN PLACE -- the vectors are far too large to copy
+//    defensively -- so rawNull is taken first and the RAW vectors are what
+//    define the L2 null space. Do not insert an Orthogonalise() anywhere:
+//    projecting a block-orthonormal vector onto its own block-orthonormalised
+//    aggregation gives e_k, and the near null content is silently gone. The
+//    ||<psi|psi> - I||_F guard below is what catches that.
+//
+// Env: LATT LS MASS NBASIS(compile time) NRHS BLOCK BLOCK2 COARSEN_BATCH
+//      HOT_START CONFIG SUBSPACE_FILE DEPRECATED_CHECK MRHS_COARSEN
+//
+
+#include <Grid/Grid.h>
+#include <algorithm>              // std::sort, for the runtime-environment dump in ParseEnvironment
+#include <Grid/lattice/PaddedCell.h>
+#include <Grid/stencil/GeneralLocalStencil.h>
+#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
+#include <Grid/algorithms/multigrid/MrhsMultiGrid.h>
+#include <Grid/algorithms/multigrid/DenseCoarseMatrix.h>
+
+#include <memory>
+
+using namespace std;
+using namespace Grid;
+
+// Compile time so it can be cut down for laptop runs: -DNBASIS=8
+#ifndef NBASIS
+#define NBASIS 60
+#endif
+
+RealD mass      = 0.00078;
+int   Nrhs      = 12;
+int   Ls        = 24;
+int   CoarsenBatch = 9;
+std::vector<int> lat_size({48,48,48,96});
+
+// Solver tuning.  PRINCIPLE (PB, 2026-08-24): the defaults ARE the current
+// optimum, so an unset environment reproduces the best banked result; they
+// are updated as and when a better point is found, and every change is
+// dated here.  Environment variables of the same names override for sweeps.
+//
+// Current optimum: 2026-08-24, slurm-5335492 F4, 48^3x96 Ls=24 on 288 GCDs,
+//   17.2 s/RHS at Nrhs=4, 32.2 s at Nrhs=1 (exact-halo FINAL ~1e-8 pending
+//   the exact-outer rerun).  Smoother mmax == order (full GCR history);
+//   PB's mmax=1 trial gave 72 vs ~60 outer iterations and was slower.
+RealD FineSmootherShift    = 0.1;
+int   FineSmootherOrder    = 6;
+int   FineSmootherMmax     = 6;
+RealD CoarseSmootherShift  = 0.1;
+int   PowerIterations      = 0;     // >0: power-iterate the smoother operators before the solves (spectral edge)
+// Smoother implementation per level (Smoothers.h):
+//   gcr    : the adaptive PGCR (default, as always)
+//   replay : run the PGCR with a coefficient recorder for the first
+//            PolyRecordIters outer steps, then switch to GCRReplaySmoother
+//            (same polynomial, no inner products) -- 1402.2585 p.13 revisited
+//   cheb   : ChebyshevNonHermitianSmoother, 1/x on [ChebLo,ChebHi], order
+//            = the GCR step count of that level
+std::string FineSmootherMode   = "gcr";
+std::string CoarseSmootherMode = "gcr";
+int   PolyRecordIters      = 4;
+int   PolyRecordStart      = 0;     // outer step at which recording begins (0: from the first step)
+std::string PolyRecordSelect = "last"; // which recorded call to replay: last|first|mean (mean is the bad one)
+int   PolyRefresh          = 0;     // >0: every PolyRefresh outer steps, one adaptive step re-records the polynomial (HDCG: every 10)
+int   PolyVerbose          = 0;     // 1: fixed-polynomial smoothers print |r_m|/|r_0| per call
+RealD FineChebLo   = 3.0,  FineChebHi   = 137.0;   // from the harvested polynomial and PowerIterations edge
+RealD CoarseChebLo = 8.0,  CoarseChebHi = 45.0;
+int   CoarseSmootherNstep  = 2;
+int   CoarseSmootherMmax   = 2;
+RealD CoarseSolverTol      = 0.05;
+int   CoarseSolverOrder    = 200;
+int   CoarseSolverMmax     = 16;
+RealD OuterTol             = 1.0e-8;
+int   OuterMmax            = 4;
+int   OuterNstep           = 8;
+
+// "It's legal to get the same answer faster, not to get a less correct
+// answer."  (PB, 2026-08-24)
+//
+// Halo-precision POLICY: reduced-precision (fp32 wire)
+// halos belong in the PRECONDITIONER -- the smoother, the V-cycle's own
+// residuals, and the coarsening -- and NEVER in the outer Krylov.  The
+// outer operator's applications define what "converged" means; making them
+// sloppy turns the stopping criterion into a statement about the wrong
+// operator (measured: solver stops at computed 9.8e-9 while the true
+// residual is 3.4e-8).  Exactness costs one exact fine matvec per outer
+// iteration, a few percent of the solve.
+//
+// Stencil::SloppyComms is a free runtime setter (Stencil.h:303), so the
+// policy is implemented by SCOPED toggling: SetFineSloppy(1) on entering
+// the preconditioner / coarsening, SetFineSloppy(0) on leaving.  The
+// operators default to EXACT.  FineSloppyComms therefore now means
+// "sloppy inside the preconditioner"; =0 makes everything exact.
+int   FineSloppyComms      = 1;
+// SmootherCoeffLog=1 : print the GCR step lengths a_k and orthogonalisation
+// coefficients b_kj of BOTH smoothers every call -- the harvest for a fixed
+// polynomial smoother (stable coefficients => stationary p(A), no reductions).
+int   SmootherCoeffLog     = 0;
+std::function<void(int)> SetFineSloppy = [](int){};
+
+void ParseEnvironment(void)
+{
+  if(getenv("MASS"))          mass        = atof(getenv("MASS"));
+  if(getenv("NRHS"))          Nrhs        = atoi(getenv("NRHS"));
+  if(getenv("LS"))            Ls          = atoi(getenv("LS"));
+  if(getenv("COARSEN_BATCH")) CoarsenBatch= atoi(getenv("COARSEN_BATCH"));
+  if(getenv("FineSmootherShift"))  FineSmootherShift  = atof(getenv("FineSmootherShift"));
+  if(getenv("FineSmootherOrder"))  FineSmootherOrder  = atoi(getenv("FineSmootherOrder"));
+  if(getenv("FineSmootherMmax"))   FineSmootherMmax   = atoi(getenv("FineSmootherMmax"));
+  if(getenv("CoarseSmootherShift"))CoarseSmootherShift= atof(getenv("CoarseSmootherShift"));
+  if(getenv("PowerIterations"))    PowerIterations    = atoi(getenv("PowerIterations"));
+  if(getenv("FineSmootherMode"))   FineSmootherMode   = getenv("FineSmootherMode");
+  if(getenv("CoarseSmootherMode")) CoarseSmootherMode = getenv("CoarseSmootherMode");
+  if(getenv("PolyRecordIters"))    PolyRecordIters    = atoi(getenv("PolyRecordIters"));
+  if(getenv("PolyRecordStart"))    PolyRecordStart    = atoi(getenv("PolyRecordStart"));
+  if(getenv("PolyRecordSelect"))   PolyRecordSelect   = getenv("PolyRecordSelect");
+  if(getenv("PolyRefresh"))        PolyRefresh        = atoi(getenv("PolyRefresh"));
+  if(getenv("PolyVerbose"))        PolyVerbose        = atoi(getenv("PolyVerbose"));
+  if(getenv("FineChebLo"))         FineChebLo         = atof(getenv("FineChebLo"));
+  if(getenv("FineChebHi"))         FineChebHi         = atof(getenv("FineChebHi"));
+  if(getenv("CoarseChebLo"))       CoarseChebLo       = atof(getenv("CoarseChebLo"));
+  if(getenv("CoarseChebHi"))       CoarseChebHi       = atof(getenv("CoarseChebHi"));
+  if(getenv("CoarseSmootherNstep"))CoarseSmootherNstep= atoi(getenv("CoarseSmootherNstep"));
+  if(getenv("CoarseSmootherMmax")) CoarseSmootherMmax = atoi(getenv("CoarseSmootherMmax"));
+  if(getenv("CoarseSolverTol"))    CoarseSolverTol    = atof(getenv("CoarseSolverTol"));
+  if(getenv("CoarseSolverOrder"))  CoarseSolverOrder  = atoi(getenv("CoarseSolverOrder"));
+  if(getenv("CoarseSolverMmax"))   CoarseSolverMmax   = atoi(getenv("CoarseSolverMmax"));
+  if(getenv("OuterTol"))           OuterTol           = atof(getenv("OuterTol"));
+  if(getenv("OuterMmax"))          OuterMmax          = atoi(getenv("OuterMmax"));
+  if(getenv("FineSloppyComms"))    FineSloppyComms    = atoi(getenv("FineSloppyComms"));
+  if(getenv("SmootherCoeffLog"))   SmootherCoeffLog   = atoi(getenv("SmootherCoeffLog"));
+  if(getenv("OuterNstep"))         OuterNstep         = atoi(getenv("OuterNstep"));
+  if(getenv("LATT")){
+    Coordinate l;
+    GridCmdOptionIntVector(std::string(getenv("LATT")),l);
+    GRID_ASSERT(l.size()==4);
+    for(int d=0;d<4;d++) lat_size[d]=l[d];
+  }
+
+  std::cout << GridLogMessage << "PARAM: LATT               "
+	    << lat_size[0]<<"."<<lat_size[1]<<"."<<lat_size[2]<<"."<<lat_size[3] << std::endl;
+  std::cout << GridLogMessage << "PARAM: LS                 " << Ls           << std::endl;
+  std::cout << GridLogMessage << "PARAM: MASS               " << mass         << std::endl;
+  std::cout << GridLogMessage << "PARAM: NBASIS             " << NBASIS       << std::endl;
+  std::cout << GridLogMessage << "PARAM: NRHS               " << Nrhs         << std::endl;
+  std::cout << GridLogMessage << "PARAM: COARSEN_BATCH      " << CoarsenBatch << std::endl;
+  // EVERY knob this programme parsed, so a log identifies its own run (2026-08-28:
+  // four jobs in the queue differing in FineSloppyComms/NRHS/mmax and none of it
+  // printed).  Interim until the Serializable parameter struct replaces all of this.
+  auto P = [](const char *n, auto v){ std::cout << GridLogMessage << "PARAM: " << std::left << std::setw(20) << n << std::right << " " << v << std::endl; };
+  P("FineSmootherShift",FineSmootherShift);   P("FineSmootherOrder",FineSmootherOrder);   P("FineSmootherMmax",FineSmootherMmax);
+  P("FineSmootherMode",FineSmootherMode);     P("FineChebLo",FineChebLo);                 P("FineChebHi",FineChebHi);
+  P("FineSloppyComms",FineSloppyComms);
+  P("CoarseSmootherShift",CoarseSmootherShift); P("CoarseSmootherNstep",CoarseSmootherNstep); P("CoarseSmootherMmax",CoarseSmootherMmax);
+  P("CoarseSmootherMode",CoarseSmootherMode); P("CoarseChebLo",CoarseChebLo);             P("CoarseChebHi",CoarseChebHi);
+  P("CoarseSolverTol",CoarseSolverTol);       P("CoarseSolverOrder",CoarseSolverOrder);   P("CoarseSolverMmax",CoarseSolverMmax);
+  P("OuterTol",OuterTol);                     P("OuterMmax",OuterMmax);                   P("OuterNstep",OuterNstep);
+  P("PowerIterations",PowerIterations);       P("PolyRecordIters",PolyRecordIters);       P("PolyRecordStart",PolyRecordStart);
+  P("PolyRecordSelect",PolyRecordSelect);     P("PolyRefresh",PolyRefresh);               P("PolyVerbose",PolyVerbose);
+  P("SmootherCoeffLog",SmootherCoeffLog);
+  // library-side knobs, read straight from the environment as the library will
+  const char *envs[] = {"BLOCK","BLOCK2","SUBSPACE_FILE","CONFIG","HOT_START","MRHS_COARSEN","DEPRECATED_CHECK",
+                        "DENSE_CC","DENSE_SCHUR","DENSE_SCHUR2D","DENSE_DEVICE_SUM","DENSE_SPLITK","DENSE_NB",
+                        "SCHUR2D_LEAF_SPAN","SCHUR2D_LEAF_LU","SCHUR2D_PROBE","SUMMA_HANDSHAKE","DENSE_APPLY_PROFILE","SLAB_FILE",
+                        "SOLVE_SRHS"};
+  for(auto e : envs) P(e, getenv(e) ? std::string(getenv(e)) : std::string("(unset)"));
+
+  ////////////////////////////////////////////////////////////////////////////////////////
+  // The COMPLETE fabric/runtime environment, scanned from environ rather than a curated
+  // list.  2026-08-30: we could not tell from a failing log whether FI_CXI_ATS was set,
+  // because it was not in the list -- and a module or the submitting shell can set
+  // anything.  A log must describe its own run without reference to the job script.
+  ////////////////////////////////////////////////////////////////////////////////////////
+  {
+    extern char **environ;
+    const char *prefixes[] = {"FI_","OFI_","MPICH_","MPIR_","PMI_","CXI_","HSA_","HIP_","ROCR_",
+                              "ROCM_","AMD_","GPU_","NUMA_","OMP_","GOMP_","GRID_","CRAY_","LIBFABRIC"};
+    std::vector<std::string> hits;
+    for(char **e = environ; e && *e; e++){
+      std::string s(*e);
+      for(auto p : prefixes){
+        if( s.compare(0,strlen(p),p)==0 ){
+          if(s.size()>200) s = s.substr(0,197)+"...";     // paths can be enormous
+          hits.push_back(s);
+          break;
+        }
+      }
+    }
+    std::sort(hits.begin(),hits.end());
+    std::cout << GridLogMessage << "PARAM: ---- runtime environment: "<<hits.size()<<" variables set ----"<<std::endl;
+    for(auto &s : hits) std::cout << GridLogMessage << "PARAM: ENV " << s << std::endl;
+    // Variables whose ABSENCE is as informative as their value (defaults bite: an unset
+    // FI_MR_CACHE_MONITOR means libfabric's defective memhooks, see systems/WorkArounds.txt).
+    const char *notable[] = {"FI_MR_CACHE_MONITOR","FI_MR_CACHE_MAX_COUNT","FI_MR_ROCR_CACHE_MONITOR_ENABLED",
+                             "FI_CXI_ATS","FI_CXI_RDZV_THRESHOLD","FI_CXI_DISABLE_HMEM_DEV_REGISTER",
+                             "FI_HMEM_ROCR_USE_DMABUF","MPICH_GPU_SUPPORT_ENABLED","MPICH_OFI_NIC_POLICY",
+                             "MPICH_SMP_SINGLE_COPY_MODE","OMP_NUM_THREADS","GRID_ALLOC_NCACHE_LARGE"};
+    for(auto n : notable) if(!getenv(n))
+      std::cout << GridLogMessage << "PARAM: ENV " << n << " (UNSET - provider/runtime default applies)" << std::endl;
+  }
+}
+
+template <class Field>
+void saveSubspace(std::vector<Field> &subspace, std::string const fname){
+#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);
+  SW.close();
+#endif
+}
+template <class Field>
+void loadSubspace(std::vector<Field> &subspace, std::string const fname){
+#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);
+  SR.close();
+#endif
+}
+
+//////////////////////////////////////////////////////////////////////
+// A = PV^dag M (non-Hermitian)
+//////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field>
+class PVdagMLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat; Matrix &_PV;
+public:
+  PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV) {};
+  void OpDiag (const Field &in, Field &out) { assert(0); }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) { assert(0); }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){ assert(0); };
+  void Op     (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); }
+  void AdjOp  (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(in,tmp); _Mat.Mdag(tmp,out); }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ HermOp(in,out); ComplexD d=innerProduct(in,out); n1=real(d); n2=norm2(out); }
+  void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+//////////////////////////////////////////////////////////////////////
+// ||<v|v> - I||_F over a set of coarse vectors. Small means the raw near null
+// content survived the projection; see GramGuard for where a leak lands.
+//////////////////////////////////////////////////////////////////////
+template<class CoarseField>
+RealD GramDefect(std::vector<CoarseField> &v)
+{
+  RealD s2=0.0;
+  for(int i=0;i<(int)v.size();i++){
+    for(int j=0;j<(int)v.size();j++){
+      ComplexD sij=TensorRemove(innerProduct(v[i],v[j]));
+      ComplexD d=sij-(i==j?ComplexD(1.0):ComplexD(0.0));
+      s2+=real(d)*real(d)+imag(d)*imag(d);
+    }
+  }
+  return std::sqrt(s2);
+}
+
+// On a leak every image collapses to the block unit e_k, the Gram becomes
+// N*I, and the defect lands at (N-1)*sqrt(nbasis) -- orders above the ~0.2
+// of a content preserving projection. Trip well below that so a mis-set
+// threshold costs a log line rather than the run.
+template<class CoarseField>
+void GramGuard(const std::string &name,std::vector<CoarseField> &v,GridBase *grid)
+{
+  RealD defect = GramDefect(v);
+  RealD N      = (RealD)grid->gSites();
+  RealD leak   = (N-1.0)*std::sqrt((RealD)v.size());
+  RealD trip   = std::sqrt(N);
+  std::cout << GridLogMessage << "GUARD: ||<"<<name<<"|"<<name<<"> - I||_F = " << defect
+	    << "   (e_k leak would be " << leak << ", trip at " << trip << ")" << std::endl;
+  GRID_ASSERT( defect < trip );
+}
+
+
+//////////////////////////////////////////////////////////////////////
+// Shifted variants for the smoothers
+//////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field>
+class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat; Matrix &_PV;
+public:
+  RealD shift;
+  ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
+  void OpDiag (const Field &in, Field &out) { assert(0); }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) { assert(0); }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){ assert(0); };
+  void Op     (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); out = out + shift*in; }
+  void AdjOp  (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(tmp,out); _Mat.Mdag(in,tmp); out = out + shift*in; }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
+  void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+template<class Field>
+class ShiftedLinearOperator : public LinearOperatorBase<Field> {
+  LinearOperatorBase<Field> &_Op; RealD shift;
+public:
+  ShiftedLinearOperator(RealD _shift, LinearOperatorBase<Field> &Op) : _Op(Op), shift(_shift) {}
+  void OpDiag  (const Field &in, Field &out) { assert(0); }
+  void OpDir   (const Field &in, Field &out,int dir,int disp) { assert(0); }
+  void OpDirAll (const Field &in, std::vector<Field> &out) { assert(0); }
+  void Op      (const Field &in, Field &out) { _Op.Op(in,out);    out = out + shift*in; }
+  void AdjOp   (const Field &in, Field &out) { _Op.AdjOp(in,out); out = out + shift*in; }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
+  void HermOp  (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+//////////////////////////////////////////////////////////////////////
+// Power iteration on a (non-Hermitian) operator: the spectral edge the
+// smoother polynomial must not exceed.  Reports
+//   step 0 : |A v|/|v| on a RANDOM unit v -- a one-sample lower bound on
+//            sigma_max(A).  If this and the converged value agree, the
+//            operator is near-normal and the spectral picture (R_m(lambda)
+//            on the spectrum) is trustworthy; if not, the field of values
+//            sets the safe interval and the spectrum understates it.
+//   step k : |A v_k|/|v_k| -> |lambda_max| as v_k -> the dominant
+//            eigenvector; the complex Rayleigh quotient <v,Av> gives its
+//            phase (real => on the axis).  A non-converging oscillation
+//            means a complex-conjugate pair of equal modulus at the top.
+// Uses Op(), not HermOp(): this is the operator the smoother sees.
+//////////////////////////////////////////////////////////////////////
+template<class Field>
+void PowerIteration(const std::string &name, LinearOperatorBase<Field> &Op, GridBase *grid, int iters)
+{
+  GRID_TRACE("PowerIteration");
+  GridParallelRNG RNG(grid); RNG.SeedFixedIntegers(std::vector<int>({7,11,13,17}));
+  Field v(grid), Av(grid);
+  gaussian(RNG,v);
+  RealD nv = std::sqrt(norm2(v)); v = v*(1.0/nv);
+  RealD ratio=0.0, ratio0=0.0; ComplexD rq(0.0);
+  for(int i=0;i<iters;i++){
+    Op.Op(v,Av);
+    RealD nAv = std::sqrt(norm2(Av));
+    rq = innerProduct(v,Av);                 // v is unit
+    ratio = nAv;
+    if ( i==0 ) {
+      ratio0 = ratio;
+      std::cout << GridLogMessage << "PowerIteration " << name << " step 0 (random v): |Av|/|v| = " << ratio
+                << "   [lower bound on sigma_max]" << std::endl;
+    }
+    if ( (i%10==0) || (i==iters-1) )
+      std::cout << GridLogMessage << "PowerIteration " << name << " step " << i << " |Av|/|v| = " << ratio
+                << "  Rayleigh = (" << real(rq) << "," << imag(rq) << ")" << std::endl;
+    v = Av*(1.0/nAv);
+  }
+  std::cout << GridLogMessage << "PowerIteration " << name << " SUMMARY: |lambda_max| ~ " << ratio
+            << "  Rayleigh (" << real(rq) << "," << imag(rq) << ")"
+            << "  phase " << std::atan2(imag(rq),real(rq)) << " rad"
+            << "  step-0 ratio / converged = " << ratio0/ratio
+            << (ratio0/ratio > 1.2 ? "   ** non-normal: sigma_max well above |lambda_max| **" : "   (near-normal)")
+            << std::endl;
+}
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+  ParseEnvironment();
+
+  RealD M5=1.8, b=1.5, c=0.5;
+  const int nbasis=NBASIS;
+  const int nrhs=Nrhs;
+  const int batch=CoarsenBatch;
+
+  Coordinate mpi = GridDefaultMpi();
+  Coordinate fsimd= GridDefaultSimd(Nd,vComplex::Nsimd());
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(lat_size,fsimd,mpi);
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  // Level 1 blocking (default 2^4)
+  Coordinate clatt = lat_size;
+  Coordinate Block({2,2,3,3});         // L1->L2 blocking; banked optimum 2026-08-24 (env BLOCK overrides)
+  if ( getenv("BLOCK") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK")),Block); GRID_ASSERT(Block.size()==4); }
+  for(int d=0;d<4;d++){ GRID_ASSERT(lat_size[d]%Block[d]==0); clatt[d]=lat_size[d]/Block[d]; }
+  std::cout << GridLogMessage << "Block  " << Block  << "  coarse lattice " << clatt << std::endl;
+
+  //////////////////////////////////////////////////////////////////////
+  // The coarse space is unvectorised. The 5D coarse grid is built here
+  // rather than through SpaceTimeGrid so the SIMD layout is ours.
+  //////////////////////////////////////////////////////////////////////
+  Coordinate c5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]});
+  Coordinate c5simd({1,1,1,1,1});
+  Coordinate c5mpi ({1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  GridCartesian *Coarse5d = new GridCartesian(c5latt,c5simd,c5mpi);
+
+  // 6D coarse multiRHS grid: rhs is dim 0, undistributed and unvectorised.
+  // No divisibility constraint on nrhs, unlike the deprecated operators.
+  Coordinate cmlatt({nrhs,1,clatt[0],clatt[1],clatt[2],clatt[3]});
+  Coordinate cmsimd({1,1,1,1,1,1});
+  Coordinate cmmpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  GridCartesian *CoarseMrhs = new GridCartesian(cmlatt,cmsimd,cmmpi);
+
+  // 6D coarse grid at the coarsening batch, used only while CoarsenOperator
+  // runs. The matrix elements survive the change back to nrhs.
+  Coordinate cblatt({batch,1,clatt[0],clatt[1],clatt[2],clatt[3]});
+  GridCartesian *CoarseBatch = new GridCartesian(cblatt,cmsimd,cmmpi);
+
+  // 6D fine grid carrying the coarsening batch: fine SIMD layout preserved
+  Coordinate fmlatt({batch,Ls,lat_size[0],lat_size[1],lat_size[2],lat_size[3]});
+  Coordinate fmsimd({1,1,fsimd[0],fsimd[1],fsimd[2],fsimd[3]});
+  Coordinate fmmpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  GridCartesian *FineMrhs = new GridCartesian(fmlatt,fmsimd,fmmpi);
+
+  std::cout << GridLogMessage << "Nsimd  fine " << FGrid->Nsimd()
+	    << "   coarse " << Coarse5d->Nsimd() << std::endl;
+
+  GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers({1,2,3,4});
+  GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8});
+
+  //////////////////////////////////////////////////////////////////////
+  // Gauge field
+  //////////////////////////////////////////////////////////////////////
+  LatticeGaugeField Umu(UGrid);
+  if ( getenv("HOT_START") ) {
+    std::cout << GridLogMessage << "Hot start gauge field" << std::endl;
+    SU<Nc>::HotConfiguration(RNG4,Umu);
+  } else {
+    std::string file("/ccs/home/poare/ckpoint_lat.1000");
+    if ( getenv("CONFIG") ) file = std::string(getenv("CONFIG"));
+    std::cout << GridLogMessage << "Reading gauge field " << file << std::endl;
+    FieldMetaData header;
+    NerscIO::readConfiguration(Umu,header,file);
+  }
+
+  MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
+  MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c);
+
+  // PVdagM and ShiftedPVdagM are thin wrappers over these same two objects,
+  // so this one callback controls every fine halo in the program.  Default
+  // EXACT; the preconditioner and the coarsening turn sloppiness on for
+  // their own scope only (policy note at FineSloppyComms).
+  SetFineSloppy = [&Ddwf,&Dpv](int sloppy){
+    Ddwf.SloppyComms(sloppy);
+    Dpv .SloppyComms(sloppy);
+  };
+  SetFineSloppy(0);
+  std::cout << GridLogMessage << "Fine halo policy: preconditioner+coarsening "
+            << (FineSloppyComms ? "SLOPPY (fp32 wire)" : "exact")
+            << ", outer Krylov EXACT" << std::endl;
+
+  typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD>        PVdagM_t;
+  typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
+  PVdagM_t PVdagM(Ddwf,Dpv);
+
+  //////////////////////////////////////////////////////////////////////
+  // Level 1 types: unvectorised coarse scalar
+  //////////////////////////////////////////////////////////////////////
+  typedef sTComplexD                                                          CComplexS;
+  typedef MultiGeneralCoarsenedOperator<vSpinColourVector,CComplexS,nbasis> CoarseOperator;
+  typedef CoarseOperator::CoarseVector                                        CoarseVector;
+  typedef Aggregation<vSpinColourVector,CComplexS,nbasis>                     Subspace;
+
+  NextToNearestStencilGeometry5D geom(Coarse5d);
+
+  //////////////////////////////////////////////////////////////////////
+  // Subspace: load RAW (no Orthogonalise!), or generate.
+  //
+  // The Aggregation is scaffolding for CreateSubspaceGCR only. That runs
+  // entirely on the fine grid and ends in GlobalOrthonormalise, which is a
+  // whole-lattice Gram-Schmidt, so the coarse grid it holds is never
+  // dereferenced and may be the unvectorised one.
+  //////////////////////////////////////////////////////////////////////
+  std::string subspace_file = "subspace_nb" + std::to_string(nbasis) + ".scidac";
+  if ( getenv("SUBSPACE_FILE") ) subspace_file = std::string(getenv("SUBSPACE_FILE"));
+  uint64_t file_exists=0;
+  if ( UGrid->IsBoss() ){ std::ifstream f(subspace_file); file_exists=f.good()?1:0; }
+  UGrid->GlobalSum(file_exists);
+
+  const int cb=0;
+  Subspace AggregatesGCR(Coarse5d,FGrid,cb);
+  if ( file_exists ){
+    std::cout << GridLogMessage << "*** Loading subspace from disk (kept RAW) ***" << std::endl;
+    loadSubspace(AggregatesGCR.subspace, subspace_file);
+  } else {
+    std::cout << GridLogMessage << "*** GCR subspace generation ***" << std::endl;
+    AggregatesGCR.CreateSubspaceGCR(RNG5,PVdagM,nbasis);
+    saveSubspace(AggregatesGCR.subspace, subspace_file);
+  }
+
+  // RAW copy BEFORE CoarsenOperator block-orthonormalises in place.
+  std::vector<LatticeFermionD> rawNull(nbasis,FGrid);
+  for(int k=0;k<nbasis;k++) rawNull[k]=AggregatesGCR.subspace[k];
+
+  //////////////////////////////////////////////////////////////////////
+  // L1 coarsening. The fine operator is single RHS, so it is promoted to
+  // the 6D batch grid; a natively multiRHS fine operator would substitute
+  // here with no other change.
+  //////////////////////////////////////////////////////////////////////
+  CoarseOperator CoarseOpPV(geom,Coarse5d);
+  CoarseOpPV.SetGrid(CoarseBatch);
+
+  std::cout << GridLogMessage << "*** L1 CoarsenOperator, batch "<<batch<<" ***" << std::endl;
+  SetFineSloppy(FineSloppyComms);   // coarsening builds the PRECONDITIONER
+  if ( getenv("MRHS_COARSEN") ) {
+    // Promote the single RHS operator and pack the batch: costs an
+    // ExtractSlice/InsertSlice pair per rhs. Here for the A/B only; this is
+    // the path a natively multiRHS fine operator would take.
+    MrhsPromotedOperator<LatticeFermionD> MrhsPVdagM(PVdagM,FGrid,batch);
+    CoarseOpPV.CoarsenOperator(MrhsPVdagM,FineMrhs,AggregatesGCR.subspace,Coarse5d);
+  } else {
+    // PVdagM is single RHS: apply it directly, batch on the coarse side.
+    CoarseOpPV.CoarsenOperator(PVdagM,AggregatesGCR.subspace,Coarse5d,batch);
+  }
+  SetFineSloppy(0);
+
+  // (moved here 2026-08-28: it needs AggregatesGCR.subspace, which is freed right after
+  //  the projector imports it below -- the 60 fine vectors are 20 GB of host memory per
+  //  rank and 60 LRU-eligible device fields competing with the solver's working set)
+  //////////////////////////////////////////////////////////////////////
+  // Optional cross check of the coarse matrix elements against the
+  // deprecated path, which needs a vectorised coarse space. Block Gram-Schmidt
+  // is idempotent, so the deprecated path may re-orthonormalise the same vectors in place
+  // without a second copy of the subspace.
+  //////////////////////////////////////////////////////////////////////
+  if ( getenv("DEPRECATED_CHECK") ) {
+
+    typedef DeprecatedGeneralCoarsenedMatrix     <vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
+    typedef DeprecatedMultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MrhsLittleDiracOperator;
+    typedef Aggregation                <vSpinColourVector,vTComplex,nbasis> SubspaceV;
+
+    Coordinate v5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]});
+    Coordinate v5simd({1,fsimd[0],fsimd[1],fsimd[2],fsimd[3]});
+    GridCartesian *Coarse5dV = new GridCartesian(v5latt,v5simd,c5mpi);
+
+    int nrhs_dep = vComplex::Nsimd();
+    Coordinate vmlatt({nrhs_dep,1,clatt[0],clatt[1],clatt[2],clatt[3]});
+    Coordinate vmsimd({vComplex::Nsimd(),1,1,1,1,1});
+    GridCartesian *CoarseMrhsV = new GridCartesian(vmlatt,vmsimd,cmmpi);
+
+    NextToNearestStencilGeometry5D geomV(Coarse5dV);
+
+    SubspaceV AggV(Coarse5dV,FGrid,cb);
+    for(int k=0;k<nbasis;k++) AggV.subspace[k]=AggregatesGCR.subspace[k];
+
+    LittleDiracOperator LittleDiracOpPV(geomV,FGrid,Coarse5dV);
+    std::cout << GridLogMessage << "*** deprecated CoarsenOperator (cross check) ***" << std::endl;
+    SetFineSloppy(FineSloppyComms);
+    LittleDiracOpPV.CoarsenOperator(PVdagM,AggV);
+    SetFineSloppy(0);
+
+    MrhsLittleDiracOperator mrhsDep(geomV,CoarseMrhsV);
+    mrhsDep.CopyMatrix(LittleDiracOpPV);
+
+    // BLAS_A is written by GridtoBLAS in lSite order and both sides carry
+    // the same scalar_object, so the two are directly comparable.
+    typedef MrhsLittleDiracOperator::calcMatrix calcMatrix;
+    int npoint = geom.npoint;
+    RealD num=0.0, den=0.0;
+    // The operator holds ONE site-major matrix buffer; MatrixPointOut hands back
+    // stencil point p in the same lSite order as the deprecated per-point
+    // array, and sizes a2.
+    deviceVector<calcMatrix> a2;
+    for(int p=0;p<npoint;p++){
+      int64_t sites = mrhsDep.BLAS_A[p].size();
+      CoarseOpPV.MatrixPointOut(p,a2);
+      GRID_ASSERT(sites == (int64_t)a2.size());
+      std::vector<calcMatrix> h1(sites),h2(sites);
+      acceleratorCopyFromDevice(&mrhsDep.BLAS_A[p][0],&h1[0],sites*sizeof(calcMatrix));
+      acceleratorCopyFromDevice(&a2[0],              &h2[0],sites*sizeof(calcMatrix));
+      ComplexD *w1=(ComplexD *)&h1[0];
+      ComplexD *w2=(ComplexD *)&h2[0];
+      int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD);
+      for(int64_t i=0;i<words;i++){
+        ComplexD d=w1[i]-w2[i];
+        num += real(d)*real(d)+imag(d)*imag(d);
+        den += real(w1[i])*real(w1[i])+imag(w1[i])*imag(w1[i]);
+      }
+    }
+    std::cout << GridLogMessage << "DEPRECATED_CHECK: |A_dep|^2 = " << den << std::endl;
+    std::cout << GridLogMessage << "DEPRECATED_CHECK: |A_dep - A|^2 / |A_dep|^2 = " << num/den << std::endl;
+    GRID_ASSERT( den > 0.0 );
+    GRID_ASSERT( num/den < 1.0e-18 );
+  }
+
+  // Stay on the batch grid: the L2 coarsening drives this operator at the
+  // batch. It is switched to the solve Nrhs once L2 is built.
+
+  //////////////////////////////////////////////////////////////////////
+  // psi_coarse = P^dag (RAW fine null) -> Galerkin images that carry the
+  // near null content, and are free: A_c (P psi) = P A psi.
+  //////////////////////////////////////////////////////////////////////
+  MultiRHSBlockProject<LatticeFermionD> MrhsProjector;
+  MrhsProjector.Allocate(nbasis,FGrid,Coarse5d);
+  MrhsProjector.ImportBasis(AggregatesGCR.subspace);   // block orthonormal basis
+  // The basis now lives in the projector's BLAS_V (device); nothing after this
+  // point reads AggregatesGCR.subspace (DEPRECATED_CHECK moved above).  Free it: 60 fine
+  // fields = 20 GB host per rank, and 60 LRU-eligible device copies.
+  AggregatesGCR.subspace.clear(); AggregatesGCR.subspace.shrink_to_fit();
+
+  std::vector<CoarseVector> psi_coarse(nbasis,Coarse5d);
+  MrhsProjector.blockProject(rawNull,psi_coarse);      // RAW vectors in
+  rawNull.clear(); rawNull.shrink_to_fit();
+
+  GramGuard("psi_coarse",psi_coarse,Coarse5d);
+
+  //////////////////////////////////////////////////////////////////////
+  // STAGE TWO: L2 -> L3.
+  //
+  // The fine operator here is the L1 operator, which is natively multiRHS, so the
+  // multiRHS driver applies with no promotion adapter: its D+1 grid IS the
+  // batch grid the L1 operator is currently set to.
+  //////////////////////////////////////////////////////////////////////
+  Coordinate cclatt = clatt;
+  Coordinate Block2({4,4,2,4});        // L2->L3 blocking; banked optimum 2026-08-24 (env BLOCK2 overrides)
+  if ( getenv("BLOCK2") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK2")),Block2); GRID_ASSERT(Block2.size()==4); }
+  for(int d=0;d<4;d++){ GRID_ASSERT(clatt[d]%Block2[d]==0); cclatt[d]=clatt[d]/Block2[d]; }
+  std::cout << GridLogMessage << "Block2 " << Block2 << "  coarse-coarse lattice " << cclatt << std::endl;
+
+  Coordinate cc5latt({1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+  GridCartesian *CoarseCoarse5d = new GridCartesian(cc5latt,c5simd,c5mpi);
+
+  Coordinate ccmlatt({nrhs,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+  GridCartesian *CoarseCoarseMrhs = new GridCartesian(ccmlatt,cmsimd,cmmpi);
+
+  Coordinate ccblatt({batch,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+  GridCartesian *CoarseCoarseBatch = new GridCartesian(ccblatt,cmsimd,cmmpi);
+
+  // Coarsening deepens the tensor nest by one iScalar
+  typedef CoarseVector::vector_object                                         CoarseSiteObj;
+  typedef iScalar<CComplexS>                                                  CComplexS2;
+  typedef MultiGeneralCoarsenedOperator<CoarseSiteObj,CComplexS2,nbasis>    CoarseCoarseOperator;
+  typedef CoarseCoarseOperator::CoarseVector                                  CoarseCoarseVector;
+
+  NextToNearestStencilGeometry5D geom2(CoarseCoarse5d);
+
+  // RAW copy of the coarse null vectors, for the same reason as rawNull:
+  // the L2 CoarsenOperator block-orthonormalises its subspace in place, and
+  // the L3 basis must be defined by the vectors that still carry content.
+  std::vector<CoarseVector> rawPsi(nbasis,Coarse5d);
+  for(int k=0;k<nbasis;k++) rawPsi[k]=psi_coarse[k];
+
+  CoarseCoarseOperator CoarseOpL2(geom2,CoarseCoarse5d);
+  CoarseOpL2.SetGrid(CoarseCoarseBatch);
+
+  NonHermitianLinearOperator<CoarseOperator,CoarseVector> LinOpCoarse(CoarseOpPV);
+
+  std::cout << GridLogMessage << "*** L2 CoarsenOperator, batch "<<batch<<" ***" << std::endl;
+  CoarseOpL2.CoarsenOperator(LinOpCoarse,CoarseBatch,psi_coarse,CoarseCoarse5d);
+
+  //////////////////////////////////////////////////////////////////////
+  // Both operators to the solve Nrhs. The matrix elements are Nrhs
+  // independent and survive the change.
+  //////////////////////////////////////////////////////////////////////
+  CoarseOpPV.SetGrid(CoarseMrhs);
+  CoarseOpL2.SetGrid(CoarseCoarseMrhs);
+  std::cout << GridLogMessage << "L1 operator at Nrhs " << CoarseOpPV.Nrhs()
+	    << ",  L2 operator at Nrhs " << CoarseOpL2.Nrhs() << std::endl;
+
+  //////////////////////////////////////////////////////////////////////
+  // psi_cc from the RAW coarse null vectors, and the same guard
+  //////////////////////////////////////////////////////////////////////
+  MultiRHSBlockProject<CoarseVector> MrhsProjectorL2;
+  MrhsProjectorL2.Allocate(nbasis,Coarse5d,CoarseCoarse5d);
+  MrhsProjectorL2.ImportBasis(psi_coarse);            // block orthonormal basis
+
+  {
+    std::vector<CoarseCoarseVector> psi_cc(nbasis,CoarseCoarse5d);
+    MrhsProjectorL2.blockProject(rawPsi,psi_cc);      // RAW vectors in
+
+    GramGuard("psi_cc",psi_cc,CoarseCoarse5d);
+  }
+  rawPsi.clear(); rawPsi.shrink_to_fit();
+
+  //////////////////////////////////////////////////////////////////////
+  // Both coarse operators apply on their solve grids
+  //////////////////////////////////////////////////////////////////////
+  {
+    GridParallelRNG cRNG(Coarse5d); cRNG.SeedFixedIntegers({3,4,5,6});
+    CoarseVector cin(CoarseMrhs), cout_(CoarseMrhs);
+    random(cRNG,cin);
+    CoarseOpPV.M(cin,cout_);
+    std::cout << GridLogMessage << "L1 apply |in|^2 = " << norm2(cin)
+	      << "   |M in|^2 = " << norm2(cout_) << std::endl;
+    GRID_ASSERT( norm2(cout_) > 0.0 );
+
+    GridParallelRNG ccRNG(CoarseCoarse5d); ccRNG.SeedFixedIntegers({7,8,9,10});
+    CoarseCoarseVector ccin(CoarseCoarseMrhs), ccout(CoarseCoarseMrhs);
+    random(ccRNG,ccin);
+    CoarseOpL2.M(ccin,ccout);
+    std::cout << GridLogMessage << "L2 apply |in|^2 = " << norm2(ccin)
+	      << "   |M in|^2 = " << norm2(ccout) << std::endl;
+    GRID_ASSERT( norm2(ccout) > 0.0 );
+  }
+
+  //////////////////////////////////////////////////////////////////////
+  // STAGE THREE (part one): the dense bottom on L2.
+  //
+  // DenseCoarseMatrix is bilingual: it takes the elements through
+  // Geometry()/ExtractMatrix(), so the L2 operator serves directly. It does
+  // detect that a multiRHS op cannot apply on the D dimensional grid and
+  // skips its own certificate and VERIFY, so the equivalent check is done
+  // here instead, driving the L2 operator at Nrhs 1 through a slice.
+  //////////////////////////////////////////////////////////////////////
+  typedef DenseCoarseMatrix<CComplexS2,nbasis> DenseCC_t;
+  std::unique_ptr<DenseCC_t> DenseCC;
+
+  if ( getenv("DENSE_CC")==nullptr || atoi(getenv("DENSE_CC")) ) {
+
+    std::cout << GridLogMessage << "*** L3 dense bottom: import from the L2 operator ***" << std::endl;
+    DenseCC.reset(new DenseCC_t(CoarseCoarse5d));
+    DenseCC->Import(CoarseOpL2);
+
+    ////////////////////////////////////////////////////////////////////
+    // ||A Ainv x - x|| / ||x||, the check Import could not run itself
+    ////////////////////////////////////////////////////////////////////
+    Coordinate cc1latt({1,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+    GridCartesian *CoarseCoarseOne = new GridCartesian(cc1latt,cmsimd,cmmpi);
+
+    CoarseOpL2.SetGrid(CoarseCoarseOne);
+
+    CoarseCoarseVector x(CoarseCoarse5d),y(CoarseCoarse5d),z(CoarseCoarse5d);
+    GridParallelRNG dRNG(CoarseCoarse5d); dRNG.SeedFixedIntegers({11,12,13,14});
+    random(dRNG,x);
+
+    (*DenseCC)(x,y);                    // y = Ainv x
+
+    CoarseCoarseVector y1(CoarseCoarseOne),z1(CoarseCoarseOne);
+    InsertSliceFast(y,y1,0,0);
+    CoarseOpL2.M(y1,z1);                // z = A y
+    ExtractSliceFast(z,z1,0,0);
+
+    z = z - x;
+    RealD rel = std::sqrt(norm2(z)/norm2(x));
+    std::cout << GridLogMessage << "L3 dense: ||A Ainv x - x||/||x|| = " << rel << std::endl;
+    GRID_ASSERT( rel < 1.0e-2 );
+
+    CoarseOpL2.SetGrid(CoarseCoarseMrhs);
+    delete CoarseCoarseOne;
+  }
+
+  //////////////////////////////////////////////////////////////////////
+  // STAGE THREE (part two): the solves.
+  //
+  // Both operators are driven from the SAME objects at whatever Nrhs is
+  // asked for -- the matrix elements were built once and survive SetGrid --
+  // so single RHS and multiRHS are the same code path with a different grid.
+  //////////////////////////////////////////////////////////////////////
+  GRID_ASSERT(DenseCC != nullptr);   // the PGCR bottom is not ported yet
+
+  typedef PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> FineSmoother_t;
+
+  ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv);
+  TrivialPrecon<LatticeFermionD> simple_fine;
+  TrivialPrecon<CoarseVector>    simpleC;
+
+  auto RunSolve = [&](int nr)
+  {
+    std::cout << GridLogMessage << "**********************************************" << std::endl;
+    std::cout << GridLogMessage << " THREE-level solve, Nrhs = " << nr << std::endl;
+    // Device-memory budget BEFORE the solve (2026-08-28: NRHS=6 died in hipMalloc at the
+    // first fine-smoother history allocation -- reported as an asynchronous "memory
+    // access fault" unless AMD_SERIALIZE_KERNEL/COPY made it a clean OOM).  The outer
+    // mRHS GCR holds src, sol, r, Az and OuterMmax x (p,q) fine fields PER RHS; the fine
+    // smoother adds FineSmootherMmax x (p,q) once.  The MemoryManager LRU cap
+    // (--device-mem) must be BELOW what is physically left after the non-LRU allocations
+    // (comms buffers, dense slab, stencil buffers), or the device fills before anything
+    // is evicted.  Print the estimate, the LRU state, and the device's own free count.
+    {
+      uint64_t fieldBytes = (uint64_t)FGrid->lSites()*sizeof(typename LatticeFermionD::scalar_object);
+      double outerGB = (double)nr*(4 + 2*OuterMmax)*fieldBytes/1.0e9;
+      double smthGB  = (double)(2*FineSmootherMmax + 4)*fieldBytes/1.0e9;
+      std::cout << GridLogMessage << "Device budget: fine field " << fieldBytes/1.0e6 << " MB; outer GCR history "
+                << nr << " x (4 + 2 x " << OuterMmax << ") fields = " << outerGB << " GB; fine smoother history + temps ~ "
+                << smthGB << " GB; MemoryManager device LRU " << MemoryManager::DeviceCacheBytes()/1.0e9 << " GB now, cap "
+                << MemoryManager::DeviceMaxBytes/1.0e9 << " GB" << std::endl;
+      // Empty the device LRU: every setup-era Lattice copy (coarse null vectors,
+      // coarsening temporaries) goes back to host, so the solve's working set
+      // starts from a clean device and the cap applies to it alone.
+      MemoryManager::EvictAll();
+      // ...and release the allocation caches' held blocks (setup-era deviceVector
+      // scratch that is "free" to the caller but not to hipMalloc).
+      MemoryManager::DropCache();
+      MemoryManager::PrintBytes();
+      acceleratorMem();
+    }
+    std::cout << GridLogMessage << "**********************************************" << std::endl;
+
+    Coordinate cml({nr,1,clatt[0],clatt[1],clatt[2],clatt[3]});
+    Coordinate ccml({nr,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+    GridCartesian *CMrhs  = new GridCartesian(cml, cmsimd,cmmpi);
+    GridCartesian *CCMrhs = new GridCartesian(ccml,cmsimd,cmmpi);
+
+    CoarseOpPV.SetGrid(CMrhs);
+    CoarseOpL2.SetGrid(CCMrhs);
+
+    NonHermitianLinearOperator<CoarseOperator,CoarseVector>            LinOpC (CoarseOpPV);
+    NonHermitianLinearOperator<CoarseCoarseOperator,CoarseCoarseVector> LinOpCC(CoarseOpL2);
+
+    MrhsDenseCCSolve<DenseCC_t,CoarseCoarseVector> ccSolve(*DenseCC,nr);
+
+    ShiftedLinearOperator<CoarseVector> ShiftedC(CoarseSmootherShift, LinOpC);
+    if ( PowerIterations > 0 ) {
+      // Spectral edges the smoother polynomials must respect (see
+      // scripts/gcr_polynomial.py: |R_m|>1 beyond the edge = amplification).
+      PowerIteration<CoarseVector>   ("CoarseSmootherOp(shift="+std::to_string(CoarseSmootherShift)+")", ShiftedC, CMrhs, PowerIterations);
+      PowerIteration<CoarseVector>   ("CoarseOp(unshifted)", LinOpC, CMrhs, PowerIterations);
+      PowerIteration<LatticeFermionD>("FineSmootherOp(shift="+std::to_string(FineSmootherShift)+")", ShiftedPVdagM, Ddwf.FermionGrid(), PowerIterations);
+    }
+    PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
+      CoarseSmootherGCR(0.01,1,ShiftedC,simpleC,CoarseSmootherMmax,CoarseSmootherNstep);
+    CoarseSmootherGCR.Level(2); CoarseSmootherGCR.Name("Csmoother"); CoarseSmootherGCR.SetZeroGuess(1);
+
+    SwitchableSmoother<CoarseVector> CoarseSmootherSlot(CoarseSmootherGCR,"Csmoother GCR");
+    MrhsCoarseThreeLevelPrec<CoarseVector,CoarseCoarseVector>
+      L2to3Precon(LinOpC, CoarseSmootherSlot, MrhsProjectorL2, ccSolve,
+                  Coarse5d, CoarseCoarse5d, CCMrhs, nr);
+
+    PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
+      L2PGCR(CoarseSolverTol, CoarseSolverOrder/16, LinOpC, L2to3Precon, CoarseSolverMmax, 16);
+    L2PGCR.Level(2); L2PGCR.Name("Couter"); L2PGCR.SetZeroGuess(1);
+
+    FineSmoother_t SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,FineSmootherMmax,FineSmootherOrder);
+    SmootherGCR.Level(1); SmootherGCR.Name("Fsmoother"); SmootherGCR.SetZeroGuess(1);
+    SmootherGCR.LogCoefficients(SmootherCoeffLog);
+    CoarseSmootherGCR.LogCoefficients(SmootherCoeffLog);
+
+    SwitchableSmoother<LatticeFermionD> FineSmootherSlot(SmootherGCR,"Fsmoother GCR");
+    MrhsTwoLevelMG<LatticeFermionD,CoarseVector,SwitchableSmoother<LatticeFermionD> >
+      ThreeLevelPrecon(PVdagM, FineSmootherSlot, MrhsProjector, L2PGCR, Coarse5d, CMrhs);
+    ThreeLevelPrecon.SetSloppy   = SetFineSloppy;
+    ThreeLevelPrecon.SloppyComms = FineSloppyComms;
+
+    MrhsPGCRNonHermitian<LatticeFermionD>
+      L1PGCR(OuterTol,1000,PVdagM,ThreeLevelPrecon,OuterMmax,OuterNstep);
+    L1PGCR.Level(1); L1PGCR.Name("Fouter"); L1PGCR.SetZeroGuess(1);
+
+    //////////////////////////////////////////////////////////////////////
+    // Smoother modes.  Objects live for the duration of this RunSolve.
+    //////////////////////////////////////////////////////////////////////
+    std::unique_ptr<ChebyshevNonHermitianSmoother<LatticeFermionD> > FineCheb;
+    std::unique_ptr<ChebyshevNonHermitianSmoother<CoarseVector> >    CoarseCheb;
+    std::unique_ptr<GCRReplaySmoother<LatticeFermionD> > FineReplay;
+    std::unique_ptr<GCRReplaySmoother<CoarseVector> >    CoarseReplay;
+    GCRCoefficients recF, recC;
+    {
+      GCRCoefficients::Select sel = GCRCoefficients::Last;
+      if ( PolyRecordSelect=="first" ) sel = GCRCoefficients::First;
+      if ( PolyRecordSelect=="mean" )  sel = GCRCoefficients::Mean;
+      recF.select = sel; recC.select = sel;
+    }
+    if ( FineSmootherMode == "cheb" ) {
+      FineCheb.reset(new ChebyshevNonHermitianSmoother<LatticeFermionD>(FineChebLo,FineChebHi,FineSmootherOrder,ShiftedPVdagM));
+      FineCheb->Verbose = PolyVerbose; FineCheb->name = "Fsmoother";
+      FineSmootherSlot.Set(*FineCheb,"Fsmoother Chebyshev");
+    }
+    if ( CoarseSmootherMode == "cheb" ) {
+      CoarseCheb.reset(new ChebyshevNonHermitianSmoother<CoarseVector>(CoarseChebLo,CoarseChebHi,CoarseSmootherNstep,ShiftedC));
+      CoarseCheb->Verbose = PolyVerbose; CoarseCheb->name = "Csmoother";
+      CoarseSmootherSlot.Set(*CoarseCheb,"Csmoother Chebyshev");
+    }
+    // Recording window [PolyRecordStart, PolyRecordStart+PolyRecordIters).
+    // M3 (2026-08-26): the GCR polynomial changes fast over the first outer
+    // steps (per-call |r|/|r0| 0.0034 -> 0.017 over steps 1-4) and the MEAN
+    // of those is a poor smoother (replay 0.033-0.049 per call); record a
+    // settled window instead.
+    if ( PolyRecordStart == 0 ) {
+      if ( FineSmootherMode   == "replay" ) SmootherGCR.SetCoefficientRecorder(&recF);
+      if ( CoarseSmootherMode == "replay" ) CoarseSmootherGCR.SetCoefficientRecorder(&recC);
+    }
+    // Record -> replay, with optional periodic re-recording ("re-record, not
+    // fade away": HDCG refreshed its polynomial every 10 steps, tracking the
+    // evolving spectral content of the residual).  Schedule on outer steps:
+    //   [PolyRecordStart, +PolyRecordIters)          adaptive GCR, recording
+    //   then replay of the selected recorded call;
+    //   if PolyRefresh>0: every PolyRefresh steps, ONE adaptive recording
+    //   step, then replay of that call.
+    auto BuildReplays = [&](void){
+      if ( FineSmootherMode == "replay" ) {
+        SmootherGCR.SetCoefficientRecorder(nullptr);
+        recF.Flush(); recF.Report("Fsmoother");
+        FineReplay.reset(new GCRReplaySmoother<LatticeFermionD>(ShiftedPVdagM,recF));
+        FineReplay->Verbose = PolyVerbose; FineReplay->name = "Fsmoother";
+        FineSmootherSlot.Set(*FineReplay,"Fsmoother replay");
+        SmootherGCR.ReleaseHistory();
+      }
+      if ( CoarseSmootherMode == "replay" ) {
+        CoarseSmootherGCR.SetCoefficientRecorder(nullptr);
+        recC.Flush(); recC.Report("Csmoother");
+        CoarseReplay.reset(new GCRReplaySmoother<CoarseVector>(ShiftedC,recC));
+        CoarseReplay->Verbose = PolyVerbose; CoarseReplay->name = "Csmoother";
+        CoarseSmootherSlot.Set(*CoarseReplay,"Csmoother replay");
+        CoarseSmootherGCR.ReleaseHistory();
+      }
+    };
+    auto StartRecording = [&](int step){
+      if ( FineSmootherMode   == "replay" ) { { auto sel=recF.select; recF = GCRCoefficients(); recF.select=sel; } SmootherGCR.SetCoefficientRecorder(&recF);       FineSmootherSlot.Set(SmootherGCR,"Fsmoother GCR (recording)"); }
+      if ( CoarseSmootherMode == "replay" ) { { auto sel=recC.select; recC = GCRCoefficients(); recC.select=sel; } CoarseSmootherGCR.SetCoefficientRecorder(&recC); CoarseSmootherSlot.Set(CoarseSmootherGCR,"Csmoother GCR (recording)"); }
+      std::cout << GridLogMessage << "Smoother coefficient recording starts at outer step " << step << std::endl;
+    };
+    int switchStep = PolyRecordStart + PolyRecordIters;
+    L1PGCR.OnStep = [&](int step){
+      if ( FineSmootherMode != "replay" && CoarseSmootherMode != "replay" ) return;
+      if ( step == PolyRecordStart && PolyRecordStart > 0 ) StartRecording(step);
+      if ( step == switchStep ) { BuildReplays(); return; }
+      if ( PolyRefresh > 0 && step > switchStep ) {
+        int since = step - switchStep;
+        if ( since % PolyRefresh == 0 )     { StartRecording(step); }       // one adaptive, recorded step
+        if ( since % PolyRefresh == 1 )     { BuildReplays(); }             // then replay it
+      }
+    };
+    std::cout << GridLogMessage << "Smoother modes: fine " << FineSmootherMode << "  coarse " << CoarseSmootherMode
+              << (FineSmootherMode=="replay"||CoarseSmootherMode=="replay" ? "  (record outer steps "+std::to_string(PolyRecordStart)+".."+std::to_string(PolyRecordStart+PolyRecordIters)+")" : "")
+              << std::endl;
+
+    std::vector<LatticeFermionD> src(nr,FGrid), sol(nr,FGrid);
+    for(int r=0;r<nr;r++){ gaussian(RNG5,src[r]); sol[r]=Zero(); }
+
+    GridStopWatch w; w.Start();
+    L1PGCR(src,sol);
+    w.Stop();
+    std::cout << GridLogMessage << "3-level solve Nrhs "<<nr<<" total " << w.Elapsed()
+              << "  (per RHS: " << w.useconds()/1.0e6/nr << " s)" << std::endl;
+
+    // The outer operator is exact by policy; assert the state rather than
+    // trust it -- a preconditioner that failed to restore would surface here.
+    SetFineSloppy(0);
+    { LatticeFermionD Ax(FGrid); RealD worst=0.0;
+      for(int r=0;r<nr;r++){ PVdagM.Op(sol[r],Ax); Ax=Ax-src[r];
+        RealD rn=std::sqrt(norm2(Ax)/norm2(src[r]));
+        std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": rhs["<<r<<"] true residual = " << rn << std::endl;
+        worst=std::max(worst,rn); }
+      std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": worst-case residual = " << worst
+                << "   (exact-halo verification)" << std::endl;
+    }
+
+
+    // The operators borrow these grids and build a PaddedCell on them, so
+    // they must let go before the grids are destroyed.
+    CoarseOpPV.ReleaseGrid();
+    CoarseOpL2.ReleaseGrid();
+    delete CMrhs; delete CCMrhs;
+  };
+
+  RunSolve(nrhs);
+  if ( getenv("SOLVE_SRHS")==nullptr || atoi(getenv("SOLVE_SRHS")) ) RunSolve(1);
+
+  std::cout << GridLogMessage << "*** stage three complete: solves done ***" << std::endl;
+
+  Grid_finalize();
+}
diff --git a/examples/Example_pvdagm_3level_FourierSmoother.cc b/examples/Example_pvdagm_3level_FourierSmoother.cc
new file mode 100644
index 000000000..d0351adde
--- /dev/null
+++ b/examples/Example_pvdagm_3level_FourierSmoother.cc
@@ -0,0 +1,1145 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./examples/Example_pvdagm_3level_FourierSmoother.cc
+
+    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.
+
+    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 */
+
+//
+// PVdagM three level multigrid on MultiGeneralCoarsenedOperator.
+//
+// STAGES ONE AND TWO: grids, types, subspace, and the L1 and L2 coarsenings.
+// The dense bottom and the solves are not here yet.
+//
+// Differences from Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc:
+//
+//  * The coarse space is UNVECTORISED (sComplexD). The fine space stays
+//    vectorised. MultiRHSBlockProject carries the mixed layout.
+//
+//  * One operator, not two. The deprecated path needed
+//    DeprecatedGeneralCoarsenedMatrix to coarsen and
+//    DeprecatedMultiGeneralCoarsenedMatrix to apply, bridged by CopyMatrix.
+//    This operator does both, and single versus multiRHS is SetGrid on the
+//    same object with the matrix elements built once.
+//
+//  * Nrhs is unconstrained. The deprecated path required nrhs % vComplex::Nsimd() == 0 because
+//    its multiRHS grid carried the SIMD in the rhs direction.
+//
+//  * CoarsenOperator takes the subspace vectors, not an Aggregation. It block
+//    orthonormalises them IN PLACE -- the vectors are far too large to copy
+//    defensively -- so rawNull is taken first and the RAW vectors are what
+//    define the L2 null space. Do not insert an Orthogonalise() anywhere:
+//    projecting a block-orthonormal vector onto its own block-orthonormalised
+//    aggregation gives e_k, and the near null content is silently gone. The
+//    ||<psi|psi> - I||_F guard below is what catches that.
+//
+// Env: LATT LS MASS NBASIS(compile time) NRHS BLOCK BLOCK2 COARSEN_BATCH
+//      HOT_START CONFIG SUBSPACE_FILE DEPRECATED_CHECK MRHS_COARSEN
+//
+
+#include <Grid/Grid.h>
+#include <Grid/qcd/utils/FreeMobius5D_claude.h>
+#include <algorithm>              // std::sort, for the runtime-environment dump in ParseEnvironment
+#include <Grid/lattice/PaddedCell.h>
+#include <Grid/stencil/GeneralLocalStencil.h>
+#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
+#include <Grid/algorithms/multigrid/MrhsMultiGrid.h>
+#include <Grid/algorithms/multigrid/DenseCoarseMatrix.h>
+
+#include <memory>
+
+using namespace std;
+using namespace Grid;
+
+// Compile time so it can be cut down for laptop runs: -DNBASIS=8
+#ifndef NBASIS
+#define NBASIS 60
+#endif
+
+RealD mass      = 0.00078;
+int   Nrhs      = 12;
+int   Ls        = 24;
+int   CoarsenBatch = 9;
+std::vector<int> lat_size({48,48,48,96});
+
+// Solver tuning.  PRINCIPLE (PB, 2026-08-24): the defaults ARE the current
+// optimum, so an unset environment reproduces the best banked result; they
+// are updated as and when a better point is found, and every change is
+// dated here.  Environment variables of the same names override for sweeps.
+//
+// Current optimum: 2026-08-24, slurm-5335492 F4, 48^3x96 Ls=24 on 288 GCDs,
+//   17.2 s/RHS at Nrhs=4, 32.2 s at Nrhs=1 (exact-halo FINAL ~1e-8 pending
+//   the exact-outer rerun).  Smoother mmax == order (full GCR history);
+//   PB's mmax=1 trial gave 72 vs ~60 outer iterations and was slower.
+RealD FineSmootherShift    = 0.1;
+int   FineSmootherOrder    = 6;
+int   FineSmootherMmax     = 6;
+RealD CoarseSmootherShift  = 0.1;
+int   PowerIterations      = 0;     // >0: power-iterate the smoother operators before the solves (spectral edge)
+// Smoother implementation per level (Smoothers.h):
+//   gcr    : the adaptive PGCR (default, as always)
+//   replay : run the PGCR with a coefficient recorder for the first
+//            PolyRecordIters outer steps, then switch to GCRReplaySmoother
+//            (same polynomial, no inner products) -- 1402.2585 p.13 revisited
+//   cheb   : ChebyshevNonHermitianSmoother, 1/x on [ChebLo,ChebHi], order
+//            = the GCR step count of that level
+std::string FineSmootherMode   = "gcr";
+std::string CoarseSmootherMode = "gcr";
+int   PolyRecordIters      = 4;
+int   PolyRecordStart      = 0;     // outer step at which recording begins (0: from the first step)
+std::string PolyRecordSelect = "last"; // which recorded call to replay: last|first|mean (mean is the bad one)
+int   PolyRefresh          = 0;     // >0: every PolyRefresh outer steps, one adaptive step re-records the polynomial (HDCG: every 10)
+int   PolyVerbose          = 0;     // 1: fixed-polynomial smoothers print |r_m|/|r_0| per call
+RealD FineChebLo   = 3.0,  FineChebHi   = 137.0;   // from the harvested polynomial and PowerIterations edge
+RealD CoarseChebLo = 8.0,  CoarseChebHi = 45.0;
+// Landau frame + Fourier fine smoother (FineSmootherMode "fourier"):
+// the configuration is rotated to the Landau frame at load; the smoother is the
+// free Mobius PVdagM inverse, F(1)^-dag then F(m_s), by 4D FFT (FreeMobius5D).
+// FourierSmootherMass caps the free IR gain (structural at large Ls: the free
+// operator's own e^{-alpha Ls} cap vanishes, so m_s must supply it).
+RealD GaugeFixAlpha       = 0.05;   // stable window 0.02-0.05 (Grid FA normalisation)
+RealD GaugeFixTol         = 1.0e-5; // smoothness saturates here (tol scan 2026-09-01)
+RealD FourierSmootherMass = 0.1;
+RealD SmootherGainCap     = 0.0;    // >0: cap free-block singular values (IR gain limiter); 0 = off
+int   RotateSubspace      = 1;      // rotate a LOADED subspace into the frame (legacy files)
+int   CoarseSmootherNstep  = 2;
+int   CoarseSmootherMmax   = 2;
+RealD CoarseSolverTol      = 0.05;
+int   CoarseSolverOrder    = 200;
+int   CoarseSolverMmax     = 16;
+RealD OuterTol             = 1.0e-8;
+int   OuterMmax            = 4;
+int   OuterNstep           = 8;
+
+// "It's legal to get the same answer faster, not to get a less correct
+// answer."  (PB, 2026-08-24)
+//
+// Halo-precision POLICY: reduced-precision (fp32 wire)
+// halos belong in the PRECONDITIONER -- the smoother, the V-cycle's own
+// residuals, and the coarsening -- and NEVER in the outer Krylov.  The
+// outer operator's applications define what "converged" means; making them
+// sloppy turns the stopping criterion into a statement about the wrong
+// operator (measured: solver stops at computed 9.8e-9 while the true
+// residual is 3.4e-8).  Exactness costs one exact fine matvec per outer
+// iteration, a few percent of the solve.
+//
+// Stencil::SloppyComms is a free runtime setter (Stencil.h:303), so the
+// policy is implemented by SCOPED toggling: SetFineSloppy(1) on entering
+// the preconditioner / coarsening, SetFineSloppy(0) on leaving.  The
+// operators default to EXACT.  FineSloppyComms therefore now means
+// "sloppy inside the preconditioner"; =0 makes everything exact.
+int   FineSloppyComms      = 1;
+// SmootherCoeffLog=1 : print the GCR step lengths a_k and orthogonalisation
+// coefficients b_kj of BOTH smoothers every call -- the harvest for a fixed
+// polynomial smoother (stable coefficients => stationary p(A), no reductions).
+int   SmootherCoeffLog     = 0;
+std::function<void(int)> SetFineSloppy = [](int){};
+
+void ParseEnvironment(void)
+{
+  if(getenv("MASS"))          mass        = atof(getenv("MASS"));
+  if(getenv("NRHS"))          Nrhs        = atoi(getenv("NRHS"));
+  if(getenv("LS"))            Ls          = atoi(getenv("LS"));
+  if(getenv("COARSEN_BATCH")) CoarsenBatch= atoi(getenv("COARSEN_BATCH"));
+  if(getenv("FineSmootherShift"))  FineSmootherShift  = atof(getenv("FineSmootherShift"));
+  if(getenv("FineSmootherOrder"))  FineSmootherOrder  = atoi(getenv("FineSmootherOrder"));
+  if(getenv("FineSmootherMmax"))   FineSmootherMmax   = atoi(getenv("FineSmootherMmax"));
+  if(getenv("CoarseSmootherShift"))CoarseSmootherShift= atof(getenv("CoarseSmootherShift"));
+  if(getenv("PowerIterations"))    PowerIterations    = atoi(getenv("PowerIterations"));
+  if(getenv("FineSmootherMode"))   FineSmootherMode   = getenv("FineSmootherMode");
+  if(getenv("CoarseSmootherMode")) CoarseSmootherMode = getenv("CoarseSmootherMode");
+  if(getenv("PolyRecordIters"))    PolyRecordIters    = atoi(getenv("PolyRecordIters"));
+  if(getenv("PolyRecordStart"))    PolyRecordStart    = atoi(getenv("PolyRecordStart"));
+  if(getenv("PolyRecordSelect"))   PolyRecordSelect   = getenv("PolyRecordSelect");
+  if(getenv("PolyRefresh"))        PolyRefresh        = atoi(getenv("PolyRefresh"));
+  if(getenv("PolyVerbose"))        PolyVerbose        = atoi(getenv("PolyVerbose"));
+  if(getenv("GaugeFixAlpha"))       GaugeFixAlpha       = atof(getenv("GaugeFixAlpha"));
+  if(getenv("GaugeFixTol"))         GaugeFixTol         = atof(getenv("GaugeFixTol"));
+  if(getenv("FourierSmootherMass")) FourierSmootherMass = atof(getenv("FourierSmootherMass"));
+  if(getenv("SmootherGainCap"))     SmootherGainCap     = atof(getenv("SmootherGainCap"));
+  if(getenv("RotateSubspace"))      RotateSubspace      = atoi(getenv("RotateSubspace"));
+  if(getenv("FineChebLo"))         FineChebLo         = atof(getenv("FineChebLo"));
+  if(getenv("FineChebHi"))         FineChebHi         = atof(getenv("FineChebHi"));
+  if(getenv("CoarseChebLo"))       CoarseChebLo       = atof(getenv("CoarseChebLo"));
+  if(getenv("CoarseChebHi"))       CoarseChebHi       = atof(getenv("CoarseChebHi"));
+  if(getenv("CoarseSmootherNstep"))CoarseSmootherNstep= atoi(getenv("CoarseSmootherNstep"));
+  if(getenv("CoarseSmootherMmax")) CoarseSmootherMmax = atoi(getenv("CoarseSmootherMmax"));
+  if(getenv("CoarseSolverTol"))    CoarseSolverTol    = atof(getenv("CoarseSolverTol"));
+  if(getenv("CoarseSolverOrder"))  CoarseSolverOrder  = atoi(getenv("CoarseSolverOrder"));
+  if(getenv("CoarseSolverMmax"))   CoarseSolverMmax   = atoi(getenv("CoarseSolverMmax"));
+  if(getenv("OuterTol"))           OuterTol           = atof(getenv("OuterTol"));
+  if(getenv("OuterMmax"))          OuterMmax          = atoi(getenv("OuterMmax"));
+  if(getenv("FineSloppyComms"))    FineSloppyComms    = atoi(getenv("FineSloppyComms"));
+  if(getenv("SmootherCoeffLog"))   SmootherCoeffLog   = atoi(getenv("SmootherCoeffLog"));
+  if(getenv("OuterNstep"))         OuterNstep         = atoi(getenv("OuterNstep"));
+  if(getenv("LATT")){
+    Coordinate l;
+    GridCmdOptionIntVector(std::string(getenv("LATT")),l);
+    GRID_ASSERT(l.size()==4);
+    for(int d=0;d<4;d++) lat_size[d]=l[d];
+  }
+
+  std::cout << GridLogMessage << "PARAM: LATT               "
+	    << lat_size[0]<<"."<<lat_size[1]<<"."<<lat_size[2]<<"."<<lat_size[3] << std::endl;
+  std::cout << GridLogMessage << "PARAM: LS                 " << Ls           << std::endl;
+  std::cout << GridLogMessage << "PARAM: MASS               " << mass         << std::endl;
+  std::cout << GridLogMessage << "PARAM: NBASIS             " << NBASIS       << std::endl;
+  std::cout << GridLogMessage << "PARAM: NRHS               " << Nrhs         << std::endl;
+  std::cout << GridLogMessage << "PARAM: COARSEN_BATCH      " << CoarsenBatch << std::endl;
+  // EVERY knob this programme parsed, so a log identifies its own run (2026-08-28:
+  // four jobs in the queue differing in FineSloppyComms/NRHS/mmax and none of it
+  // printed).  Interim until the Serializable parameter struct replaces all of this.
+  auto P = [](const char *n, auto v){ std::cout << GridLogMessage << "PARAM: " << std::left << std::setw(20) << n << std::right << " " << v << std::endl; };
+  P("FineSmootherShift",FineSmootherShift);   P("FineSmootherOrder",FineSmootherOrder);   P("FineSmootherMmax",FineSmootherMmax);
+  P("FineSmootherMode",FineSmootherMode);     P("FineChebLo",FineChebLo);                 P("FineChebHi",FineChebHi);
+  P("FineSloppyComms",FineSloppyComms);
+  P("CoarseSmootherShift",CoarseSmootherShift); P("CoarseSmootherNstep",CoarseSmootherNstep); P("CoarseSmootherMmax",CoarseSmootherMmax);
+  P("CoarseSmootherMode",CoarseSmootherMode); P("CoarseChebLo",CoarseChebLo);             P("CoarseChebHi",CoarseChebHi);
+  P("CoarseSolverTol",CoarseSolverTol);       P("CoarseSolverOrder",CoarseSolverOrder);   P("CoarseSolverMmax",CoarseSolverMmax);
+  P("OuterTol",OuterTol);                     P("OuterMmax",OuterMmax);                   P("OuterNstep",OuterNstep);
+  P("PowerIterations",PowerIterations);       P("PolyRecordIters",PolyRecordIters);       P("PolyRecordStart",PolyRecordStart);
+  P("PolyRecordSelect",PolyRecordSelect);     P("PolyRefresh",PolyRefresh);               P("PolyVerbose",PolyVerbose);
+  P("SmootherCoeffLog",SmootherCoeffLog);
+  P("GaugeFixAlpha",GaugeFixAlpha);           P("GaugeFixTol",GaugeFixTol);
+  P("FourierSmootherMass",FourierSmootherMass); P("SmootherGainCap",SmootherGainCap); P("RotateSubspace",RotateSubspace);
+  // library-side knobs, read straight from the environment as the library will
+  const char *envs[] = {"BLOCK","BLOCK2","SUBSPACE_FILE","CONFIG","HOT_START","MRHS_COARSEN","DEPRECATED_CHECK",
+                        "DENSE_CC","DENSE_SCHUR","DENSE_SCHUR2D","DENSE_DEVICE_SUM","DENSE_SPLITK","DENSE_NB",
+                        "SCHUR2D_LEAF_SPAN","SCHUR2D_LEAF_LU","SCHUR2D_PROBE","SUMMA_HANDSHAKE","DENSE_APPLY_PROFILE","SLAB_FILE",
+                        "SOLVE_SRHS"};
+  for(auto e : envs) P(e, getenv(e) ? std::string(getenv(e)) : std::string("(unset)"));
+
+  ////////////////////////////////////////////////////////////////////////////////////////
+  // The COMPLETE fabric/runtime environment, scanned from environ rather than a curated
+  // list.  2026-08-30: we could not tell from a failing log whether FI_CXI_ATS was set,
+  // because it was not in the list -- and a module or the submitting shell can set
+  // anything.  A log must describe its own run without reference to the job script.
+  ////////////////////////////////////////////////////////////////////////////////////////
+  {
+    extern char **environ;
+    const char *prefixes[] = {"FI_","OFI_","MPICH_","MPIR_","PMI_","CXI_","HSA_","HIP_","ROCR_",
+                              "ROCM_","AMD_","GPU_","NUMA_","OMP_","GOMP_","GRID_","CRAY_","LIBFABRIC"};
+    std::vector<std::string> hits;
+    for(char **e = environ; e && *e; e++){
+      std::string s(*e);
+      for(auto p : prefixes){
+        if( s.compare(0,strlen(p),p)==0 ){
+          if(s.size()>200) s = s.substr(0,197)+"...";     // paths can be enormous
+          hits.push_back(s);
+          break;
+        }
+      }
+    }
+    std::sort(hits.begin(),hits.end());
+    std::cout << GridLogMessage << "PARAM: ---- runtime environment: "<<hits.size()<<" variables set ----"<<std::endl;
+    for(auto &s : hits) std::cout << GridLogMessage << "PARAM: ENV " << s << std::endl;
+    // Variables whose ABSENCE is as informative as their value (defaults bite: an unset
+    // FI_MR_CACHE_MONITOR means libfabric's defective memhooks, see systems/WorkArounds.txt).
+    const char *notable[] = {"FI_MR_CACHE_MONITOR","FI_MR_CACHE_MAX_COUNT","FI_MR_ROCR_CACHE_MONITOR_ENABLED",
+                             "FI_CXI_ATS","FI_CXI_RDZV_THRESHOLD","FI_CXI_DISABLE_HMEM_DEV_REGISTER",
+                             "FI_HMEM_ROCR_USE_DMABUF","MPICH_GPU_SUPPORT_ENABLED","MPICH_OFI_NIC_POLICY",
+                             "MPICH_SMP_SINGLE_COPY_MODE","OMP_NUM_THREADS","GRID_ALLOC_NCACHE_LARGE"};
+    for(auto n : notable) if(!getenv(n))
+      std::cout << GridLogMessage << "PARAM: ENV " << n << " (UNSET - provider/runtime default applies)" << std::endl;
+  }
+}
+
+// Apply the 4D gauge transformation g to every s slice of a 5D fermion.
+void RotateBasis5D(LatticeFermionD &ferm5,const LatticeColourMatrixD &g)
+{
+  GridBase *FGrid = ferm5.Grid();
+  GridBase *UGrid = g.Grid();
+  int Ls5 = FGrid->_fdimensions[0];
+  LatticeFermionD ferm4((GridCartesian *)UGrid);
+  for(int s=0;s<Ls5;s++){
+    ExtractSlice(ferm4,ferm5,s,0);
+    ferm4 = g*ferm4;
+    InsertSlice(ferm4,ferm5,s,0);
+  }
+}
+
+// Fourier fine smoother: the free Mobius PVdagM inverse in the Landau frame,
+//   (Dpv^dag Dm)^-1_free = F(m_s) . F(1)^-dag .
+// F(1)^-dag has the same phase/FFT structure as F(1) with every momentum block
+// conjugate-transposed; DaggerBlocks() post-processes a built F in place.
+// Cap the singular values of every per-momentum free block inverse at Gmax.
+// The smoother is a hyperparameter, not a physical action, so S(p) is unconstrained
+// by Nielsen-Ninomiya: this bounds the free-IR gain (kills the amplifying band that
+// destabilises a stationary smoother) while leaving the UV -- where the block gain is
+// already < Gmax -- untouched.  One scalar, self-targeting to the IR pole.  The
+// interacting IR is Q's job (the nbasis subspace); the cap only stops the smoother
+// harming it.  <=0 disables.
+void GainCap(FreeMobius5DInverse<WilsonImplD> &F, double Gmax)
+{
+  if ( Gmax <= 0.0 ) return;
+  int n5 = 4*F.Ls;
+  size_t nslot = F.Minv_dev.size()/((size_t)n5*n5);
+  for(size_t slot=0;slot<nslot;slot++){
+    size_t base = slot*(size_t)n5*n5;
+    Eigen::MatrixXcd B(n5,n5);
+    for(int r=0;r<n5;r++) for(int cc=0;cc<n5;cc++){
+      Grid::Complex z = F.Minv_dev[base+(size_t)r*n5+cc];
+      B(r,cc) = std::complex<double>(z.real(),z.imag());
+    }
+    Eigen::JacobiSVD<Eigen::MatrixXcd> svd(B,Eigen::ComputeFullU|Eigen::ComputeFullV);
+    Eigen::VectorXd s = svd.singularValues();
+    for(int i=0;i<s.size();i++) if(s(i)>Gmax) s(i)=Gmax;
+    Eigen::MatrixXcd Bc = svd.matrixU()*s.asDiagonal()*svd.matrixV().adjoint();
+    for(int r=0;r<n5;r++) for(int cc=0;cc<n5;cc++){
+      std::complex<double> z = Bc(r,cc);
+      F.Minv_dev[base+(size_t)r*n5+cc] = Grid::Complex(z.real(),z.imag());
+    }
+  }
+}
+
+void DaggerBlocks(FreeMobius5DInverse<WilsonImplD> &F)
+{
+  int n5 = 4*F.Ls;
+  for(int i=0;i<(int)F.Minv.size();i++) F.Minv[i] = F.Minv[i].adjoint().eval();
+  size_t nslot = F.Minv_dev.size()/((size_t)n5*n5);
+  for(size_t slot=0;slot<nslot;slot++){
+    size_t base = slot*(size_t)n5*n5;
+    for(int r=0;r<n5;r++){
+      for(int cc=r+1;cc<n5;cc++){
+        Grid::Complex a = F.Minv_dev[base+(size_t)r*n5+cc];
+        F.Minv_dev[base+(size_t)r*n5+cc] = conjugate(F.Minv_dev[base+(size_t)cc*n5+r]);
+        F.Minv_dev[base+(size_t)cc*n5+r] = conjugate(a);
+      }
+      F.Minv_dev[base+(size_t)r*n5+r] = conjugate(F.Minv_dev[base+(size_t)r*n5+r]);
+    }
+  }
+}
+
+class FourierPVdagMSmoother : public LinearFunction<LatticeFermionD> {
+public:
+  using LinearFunction<LatticeFermionD>::operator();
+  FreeMobius5DInverse<WilsonImplD> &Fm;     // F(m_s)
+  FreeMobius5DInverse<WilsonImplD> &F1dag;  // F(1) with daggered blocks
+  FourierPVdagMSmoother(FreeMobius5DInverse<WilsonImplD> &_Fm,
+                        FreeMobius5DInverse<WilsonImplD> &_F1dag) : Fm(_Fm), F1dag(_F1dag) {}
+  void operator()(const LatticeFermionD &in,LatticeFermionD &out){
+    // Outer operator A = PV^dag M (PVdagMLinearOperator::Op = _Mat.M then _PV.Mdag).
+    // A^-1 = M^-1 (PV^dag)^-1 = F(m) . F(1)^-dag:  apply F(1)^-dag first, then F(m).
+    LatticeFermionD tmp(in.Grid());
+    F1dag(in,tmp);
+    Fm(tmp,out);
+  }
+};
+
+template <class Field>
+void saveSubspace(std::vector<Field> &subspace, std::string const fname){
+#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);
+  SW.close();
+#endif
+}
+template <class Field>
+void loadSubspace(std::vector<Field> &subspace, std::string const fname){
+#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);
+  SR.close();
+#endif
+}
+
+//////////////////////////////////////////////////////////////////////
+// A = PV^dag M (non-Hermitian)
+//////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field>
+class PVdagMLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat; Matrix &_PV;
+public:
+  PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV) {};
+  void OpDiag (const Field &in, Field &out) { assert(0); }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) { assert(0); }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){ assert(0); };
+  void Op     (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); }
+  void AdjOp  (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(in,tmp); _Mat.Mdag(tmp,out); }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ HermOp(in,out); ComplexD d=innerProduct(in,out); n1=real(d); n2=norm2(out); }
+  void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+//////////////////////////////////////////////////////////////////////
+// ||<v|v> - I||_F over a set of coarse vectors. Small means the raw near null
+// content survived the projection; see GramGuard for where a leak lands.
+//////////////////////////////////////////////////////////////////////
+template<class CoarseField>
+RealD GramDefect(std::vector<CoarseField> &v)
+{
+  RealD s2=0.0;
+  for(int i=0;i<(int)v.size();i++){
+    for(int j=0;j<(int)v.size();j++){
+      ComplexD sij=TensorRemove(innerProduct(v[i],v[j]));
+      ComplexD d=sij-(i==j?ComplexD(1.0):ComplexD(0.0));
+      s2+=real(d)*real(d)+imag(d)*imag(d);
+    }
+  }
+  return std::sqrt(s2);
+}
+
+// On a leak every image collapses to the block unit e_k, the Gram becomes
+// N*I, and the defect lands at (N-1)*sqrt(nbasis) -- orders above the ~0.2
+// of a content preserving projection. Trip well below that so a mis-set
+// threshold costs a log line rather than the run.
+template<class CoarseField>
+void GramGuard(const std::string &name,std::vector<CoarseField> &v,GridBase *grid)
+{
+  RealD defect = GramDefect(v);
+  RealD N      = (RealD)grid->gSites();
+  RealD leak   = (N-1.0)*std::sqrt((RealD)v.size());
+  RealD trip   = std::sqrt(N);
+  std::cout << GridLogMessage << "GUARD: ||<"<<name<<"|"<<name<<"> - I||_F = " << defect
+	    << "   (e_k leak would be " << leak << ", trip at " << trip << ")" << std::endl;
+  GRID_ASSERT( defect < trip );
+}
+
+
+//////////////////////////////////////////////////////////////////////
+// Shifted variants for the smoothers
+//////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field>
+class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat; Matrix &_PV;
+public:
+  RealD shift;
+  ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
+  void OpDiag (const Field &in, Field &out) { assert(0); }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) { assert(0); }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){ assert(0); };
+  void Op     (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); out = out + shift*in; }
+  void AdjOp  (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(tmp,out); _Mat.Mdag(in,tmp); out = out + shift*in; }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
+  void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+template<class Field>
+class ShiftedLinearOperator : public LinearOperatorBase<Field> {
+  LinearOperatorBase<Field> &_Op; RealD shift;
+public:
+  ShiftedLinearOperator(RealD _shift, LinearOperatorBase<Field> &Op) : _Op(Op), shift(_shift) {}
+  void OpDiag  (const Field &in, Field &out) { assert(0); }
+  void OpDir   (const Field &in, Field &out,int dir,int disp) { assert(0); }
+  void OpDirAll (const Field &in, std::vector<Field> &out) { assert(0); }
+  void Op      (const Field &in, Field &out) { _Op.Op(in,out);    out = out + shift*in; }
+  void AdjOp   (const Field &in, Field &out) { _Op.AdjOp(in,out); out = out + shift*in; }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
+  void HermOp  (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+//////////////////////////////////////////////////////////////////////
+// Power iteration on a (non-Hermitian) operator: the spectral edge the
+// smoother polynomial must not exceed.  Reports
+//   step 0 : |A v|/|v| on a RANDOM unit v -- a one-sample lower bound on
+//            sigma_max(A).  If this and the converged value agree, the
+//            operator is near-normal and the spectral picture (R_m(lambda)
+//            on the spectrum) is trustworthy; if not, the field of values
+//            sets the safe interval and the spectrum understates it.
+//   step k : |A v_k|/|v_k| -> |lambda_max| as v_k -> the dominant
+//            eigenvector; the complex Rayleigh quotient <v,Av> gives its
+//            phase (real => on the axis).  A non-converging oscillation
+//            means a complex-conjugate pair of equal modulus at the top.
+// Uses Op(), not HermOp(): this is the operator the smoother sees.
+//////////////////////////////////////////////////////////////////////
+template<class Field>
+void PowerIteration(const std::string &name, LinearOperatorBase<Field> &Op, GridBase *grid, int iters)
+{
+  GRID_TRACE("PowerIteration");
+  GridParallelRNG RNG(grid); RNG.SeedFixedIntegers(std::vector<int>({7,11,13,17}));
+  Field v(grid), Av(grid);
+  gaussian(RNG,v);
+  RealD nv = std::sqrt(norm2(v)); v = v*(1.0/nv);
+  RealD ratio=0.0, ratio0=0.0; ComplexD rq(0.0);
+  for(int i=0;i<iters;i++){
+    Op.Op(v,Av);
+    RealD nAv = std::sqrt(norm2(Av));
+    rq = innerProduct(v,Av);                 // v is unit
+    ratio = nAv;
+    if ( i==0 ) {
+      ratio0 = ratio;
+      std::cout << GridLogMessage << "PowerIteration " << name << " step 0 (random v): |Av|/|v| = " << ratio
+                << "   [lower bound on sigma_max]" << std::endl;
+    }
+    if ( (i%10==0) || (i==iters-1) )
+      std::cout << GridLogMessage << "PowerIteration " << name << " step " << i << " |Av|/|v| = " << ratio
+                << "  Rayleigh = (" << real(rq) << "," << imag(rq) << ")" << std::endl;
+    v = Av*(1.0/nAv);
+  }
+  std::cout << GridLogMessage << "PowerIteration " << name << " SUMMARY: |lambda_max| ~ " << ratio
+            << "  Rayleigh (" << real(rq) << "," << imag(rq) << ")"
+            << "  phase " << std::atan2(imag(rq),real(rq)) << " rad"
+            << "  step-0 ratio / converged = " << ratio0/ratio
+            << (ratio0/ratio > 1.2 ? "   ** non-normal: sigma_max well above |lambda_max| **" : "   (near-normal)")
+            << std::endl;
+}
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+  ParseEnvironment();
+
+  RealD M5=1.8, b=1.5, c=0.5;
+  const int nbasis=NBASIS;
+  const int nrhs=Nrhs;
+  const int batch=CoarsenBatch;
+
+  Coordinate mpi = GridDefaultMpi();
+  Coordinate fsimd= GridDefaultSimd(Nd,vComplex::Nsimd());
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(lat_size,fsimd,mpi);
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  // Level 1 blocking (default 2^4)
+  Coordinate clatt = lat_size;
+  Coordinate Block({2,2,3,3});         // L1->L2 blocking; banked optimum 2026-08-24 (env BLOCK overrides)
+  if ( getenv("BLOCK") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK")),Block); GRID_ASSERT(Block.size()==4); }
+  for(int d=0;d<4;d++){ GRID_ASSERT(lat_size[d]%Block[d]==0); clatt[d]=lat_size[d]/Block[d]; }
+  std::cout << GridLogMessage << "Block  " << Block  << "  coarse lattice " << clatt << std::endl;
+
+  //////////////////////////////////////////////////////////////////////
+  // The coarse space is unvectorised. The 5D coarse grid is built here
+  // rather than through SpaceTimeGrid so the SIMD layout is ours.
+  //////////////////////////////////////////////////////////////////////
+  Coordinate c5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]});
+  Coordinate c5simd({1,1,1,1,1});
+  Coordinate c5mpi ({1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  GridCartesian *Coarse5d = new GridCartesian(c5latt,c5simd,c5mpi);
+
+  // 6D coarse multiRHS grid: rhs is dim 0, undistributed and unvectorised.
+  // No divisibility constraint on nrhs, unlike the deprecated operators.
+  Coordinate cmlatt({nrhs,1,clatt[0],clatt[1],clatt[2],clatt[3]});
+  Coordinate cmsimd({1,1,1,1,1,1});
+  Coordinate cmmpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  GridCartesian *CoarseMrhs = new GridCartesian(cmlatt,cmsimd,cmmpi);
+
+  // 6D coarse grid at the coarsening batch, used only while CoarsenOperator
+  // runs. The matrix elements survive the change back to nrhs.
+  Coordinate cblatt({batch,1,clatt[0],clatt[1],clatt[2],clatt[3]});
+  GridCartesian *CoarseBatch = new GridCartesian(cblatt,cmsimd,cmmpi);
+
+  // 6D fine grid carrying the coarsening batch: fine SIMD layout preserved
+  Coordinate fmlatt({batch,Ls,lat_size[0],lat_size[1],lat_size[2],lat_size[3]});
+  Coordinate fmsimd({1,1,fsimd[0],fsimd[1],fsimd[2],fsimd[3]});
+  Coordinate fmmpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  GridCartesian *FineMrhs = new GridCartesian(fmlatt,fmsimd,fmmpi);
+
+  std::cout << GridLogMessage << "Nsimd  fine " << FGrid->Nsimd()
+	    << "   coarse " << Coarse5d->Nsimd() << std::endl;
+
+  GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers({1,2,3,4});
+  GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8});
+
+  //////////////////////////////////////////////////////////////////////
+  // Gauge field
+  //////////////////////////////////////////////////////////////////////
+  LatticeGaugeField Umu(UGrid);
+  if ( getenv("HOT_START") ) {
+    std::cout << GridLogMessage << "Hot start gauge field" << std::endl;
+    SU<Nc>::HotConfiguration(RNG4,Umu);
+  } else {
+    std::string file("/ccs/home/poare/ckpoint_lat.1000");
+    if ( getenv("CONFIG") ) file = std::string(getenv("CONFIG"));
+    std::cout << GridLogMessage << "Reading gauge field " << file << std::endl;
+    FieldMetaData header;
+    NerscIO::readConfiguration(Umu,header,file);
+  }
+
+  //////////////////////////////////////////////////////////////////////
+  // Landau frame.  Fix a copy, project Omega back onto the group, rotate
+  // Umu in place: everything downstream runs in the smooth frame.  The
+  // plaquette is gauge invariant (checked); Omega is retained to rotate a
+  // subspace that was generated in the original frame.
+  //////////////////////////////////////////////////////////////////////
+  LatticeColourMatrixD Omega(UGrid);
+  {
+    LatticeGaugeField Ufix(UGrid); Ufix = Umu;
+    FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(
+        Ufix,Omega,GaugeFixAlpha,10000,GaugeFixTol,GaugeFixTol,true,-1,false);
+    Omega = ProjectOnGroup(Omega);
+    RealD plaq0 = WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
+    SU<Nc>::GaugeTransform<PeriodicGimplR>(Umu,Omega);
+    RealD plaq1 = WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
+    std::cout << GridLogMessage << "Landau frame: plaquette " << plaq0 << " -> " << plaq1
+              << " (invariant); 1-linkTrace/Nc = "
+              << 1.0-WilsonLoops<PeriodicGimplR>::linkTrace(Umu) << std::endl;
+    GRID_ASSERT(fabs(plaq1-plaq0) < 1.0e-10);
+  }
+
+  MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
+  MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c);
+
+  // PVdagM and ShiftedPVdagM are thin wrappers over these same two objects,
+  // so this one callback controls every fine halo in the program.  Default
+  // EXACT; the preconditioner and the coarsening turn sloppiness on for
+  // their own scope only (policy note at FineSloppyComms).
+  SetFineSloppy = [&Ddwf,&Dpv](int sloppy){
+    Ddwf.SloppyComms(sloppy);
+    Dpv .SloppyComms(sloppy);
+  };
+  SetFineSloppy(0);
+  std::cout << GridLogMessage << "Fine halo policy: preconditioner+coarsening "
+            << (FineSloppyComms ? "SLOPPY (fp32 wire)" : "exact")
+            << ", outer Krylov EXACT" << std::endl;
+
+  typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD>        PVdagM_t;
+  typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
+  PVdagM_t PVdagM(Ddwf,Dpv);
+
+  //////////////////////////////////////////////////////////////////////
+  // Level 1 types: unvectorised coarse scalar
+  //////////////////////////////////////////////////////////////////////
+  typedef sTComplexD                                                          CComplexS;
+  typedef MultiGeneralCoarsenedOperator<vSpinColourVector,CComplexS,nbasis> CoarseOperator;
+  typedef CoarseOperator::CoarseVector                                        CoarseVector;
+  typedef Aggregation<vSpinColourVector,CComplexS,nbasis>                     Subspace;
+
+  NextToNearestStencilGeometry5D geom(Coarse5d);
+
+  //////////////////////////////////////////////////////////////////////
+  // Subspace: load RAW (no Orthogonalise!), or generate.
+  //
+  // The Aggregation is scaffolding for CreateSubspaceGCR only. That runs
+  // entirely on the fine grid and ends in GlobalOrthonormalise, which is a
+  // whole-lattice Gram-Schmidt, so the coarse grid it holds is never
+  // dereferenced and may be the unvectorised one.
+  //////////////////////////////////////////////////////////////////////
+  std::string subspace_file = "subspace_nb" + std::to_string(nbasis) + ".scidac";
+  if ( getenv("SUBSPACE_FILE") ) subspace_file = std::string(getenv("SUBSPACE_FILE"));
+  uint64_t file_exists=0;
+  if ( UGrid->IsBoss() ){ std::ifstream f(subspace_file); file_exists=f.good()?1:0; }
+  UGrid->GlobalSum(file_exists);
+
+  const int cb=0;
+  Subspace AggregatesGCR(Coarse5d,FGrid,cb);
+  if ( file_exists ){
+    std::cout << GridLogMessage << "*** Loading subspace from disk (kept RAW) ***" << std::endl;
+    loadSubspace(AggregatesGCR.subspace, subspace_file);
+    if ( RotateSubspace ) {
+      // Legacy file, generated in the original frame: the near-null space is
+      // gauge covariant, so rotating by Omega gives the exact frame basis.
+      std::cout << GridLogMessage << "*** Rotating loaded subspace into the Landau frame ***" << std::endl;
+      for(int k=0;k<nbasis;k++) RotateBasis5D(AggregatesGCR.subspace[k],Omega);
+    }
+  } else {
+    std::cout << GridLogMessage << "*** GCR subspace generation ***" << std::endl;
+    AggregatesGCR.CreateSubspaceGCR(RNG5,PVdagM,nbasis);
+    saveSubspace(AggregatesGCR.subspace, subspace_file);
+  }
+
+  // RAW copy BEFORE CoarsenOperator block-orthonormalises in place.
+  std::vector<LatticeFermionD> rawNull(nbasis,FGrid);
+  for(int k=0;k<nbasis;k++) rawNull[k]=AggregatesGCR.subspace[k];
+
+  //////////////////////////////////////////////////////////////////////
+  // L1 coarsening. The fine operator is single RHS, so it is promoted to
+  // the 6D batch grid; a natively multiRHS fine operator would substitute
+  // here with no other change.
+  //////////////////////////////////////////////////////////////////////
+  CoarseOperator CoarseOpPV(geom,Coarse5d);
+  CoarseOpPV.SetGrid(CoarseBatch);
+
+  std::cout << GridLogMessage << "*** L1 CoarsenOperator, batch "<<batch<<" ***" << std::endl;
+  SetFineSloppy(FineSloppyComms);   // coarsening builds the PRECONDITIONER
+  if ( getenv("MRHS_COARSEN") ) {
+    // Promote the single RHS operator and pack the batch: costs an
+    // ExtractSlice/InsertSlice pair per rhs. Here for the A/B only; this is
+    // the path a natively multiRHS fine operator would take.
+    MrhsPromotedOperator<LatticeFermionD> MrhsPVdagM(PVdagM,FGrid,batch);
+    CoarseOpPV.CoarsenOperator(MrhsPVdagM,FineMrhs,AggregatesGCR.subspace,Coarse5d);
+  } else {
+    // PVdagM is single RHS: apply it directly, batch on the coarse side.
+    CoarseOpPV.CoarsenOperator(PVdagM,AggregatesGCR.subspace,Coarse5d,batch);
+  }
+  SetFineSloppy(0);
+
+  // (moved here 2026-08-28: it needs AggregatesGCR.subspace, which is freed right after
+  //  the projector imports it below -- the 60 fine vectors are 20 GB of host memory per
+  //  rank and 60 LRU-eligible device fields competing with the solver's working set)
+  //////////////////////////////////////////////////////////////////////
+  // Optional cross check of the coarse matrix elements against the
+  // deprecated path, which needs a vectorised coarse space. Block Gram-Schmidt
+  // is idempotent, so the deprecated path may re-orthonormalise the same vectors in place
+  // without a second copy of the subspace.
+  //////////////////////////////////////////////////////////////////////
+  if ( getenv("DEPRECATED_CHECK") ) {
+
+    typedef DeprecatedGeneralCoarsenedMatrix     <vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
+    typedef DeprecatedMultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MrhsLittleDiracOperator;
+    typedef Aggregation                <vSpinColourVector,vTComplex,nbasis> SubspaceV;
+
+    Coordinate v5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]});
+    Coordinate v5simd({1,fsimd[0],fsimd[1],fsimd[2],fsimd[3]});
+    GridCartesian *Coarse5dV = new GridCartesian(v5latt,v5simd,c5mpi);
+
+    int nrhs_dep = vComplex::Nsimd();
+    Coordinate vmlatt({nrhs_dep,1,clatt[0],clatt[1],clatt[2],clatt[3]});
+    Coordinate vmsimd({vComplex::Nsimd(),1,1,1,1,1});
+    GridCartesian *CoarseMrhsV = new GridCartesian(vmlatt,vmsimd,cmmpi);
+
+    NextToNearestStencilGeometry5D geomV(Coarse5dV);
+
+    SubspaceV AggV(Coarse5dV,FGrid,cb);
+    for(int k=0;k<nbasis;k++) AggV.subspace[k]=AggregatesGCR.subspace[k];
+
+    LittleDiracOperator LittleDiracOpPV(geomV,FGrid,Coarse5dV);
+    std::cout << GridLogMessage << "*** deprecated CoarsenOperator (cross check) ***" << std::endl;
+    SetFineSloppy(FineSloppyComms);
+    LittleDiracOpPV.CoarsenOperator(PVdagM,AggV);
+    SetFineSloppy(0);
+
+    MrhsLittleDiracOperator mrhsDep(geomV,CoarseMrhsV);
+    mrhsDep.CopyMatrix(LittleDiracOpPV);
+
+    // BLAS_A is written by GridtoBLAS in lSite order and both sides carry
+    // the same scalar_object, so the two are directly comparable.
+    typedef MrhsLittleDiracOperator::calcMatrix calcMatrix;
+    int npoint = geom.npoint;
+    RealD num=0.0, den=0.0;
+    // The operator holds ONE site-major matrix buffer; MatrixPointOut hands back
+    // stencil point p in the same lSite order as the deprecated per-point
+    // array, and sizes a2.
+    deviceVector<calcMatrix> a2;
+    for(int p=0;p<npoint;p++){
+      int64_t sites = mrhsDep.BLAS_A[p].size();
+      CoarseOpPV.MatrixPointOut(p,a2);
+      GRID_ASSERT(sites == (int64_t)a2.size());
+      std::vector<calcMatrix> h1(sites),h2(sites);
+      acceleratorCopyFromDevice(&mrhsDep.BLAS_A[p][0],&h1[0],sites*sizeof(calcMatrix));
+      acceleratorCopyFromDevice(&a2[0],              &h2[0],sites*sizeof(calcMatrix));
+      ComplexD *w1=(ComplexD *)&h1[0];
+      ComplexD *w2=(ComplexD *)&h2[0];
+      int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD);
+      for(int64_t i=0;i<words;i++){
+        ComplexD d=w1[i]-w2[i];
+        num += real(d)*real(d)+imag(d)*imag(d);
+        den += real(w1[i])*real(w1[i])+imag(w1[i])*imag(w1[i]);
+      }
+    }
+    std::cout << GridLogMessage << "DEPRECATED_CHECK: |A_dep|^2 = " << den << std::endl;
+    std::cout << GridLogMessage << "DEPRECATED_CHECK: |A_dep - A|^2 / |A_dep|^2 = " << num/den << std::endl;
+    GRID_ASSERT( den > 0.0 );
+    GRID_ASSERT( num/den < 1.0e-18 );
+  }
+
+  // Stay on the batch grid: the L2 coarsening drives this operator at the
+  // batch. It is switched to the solve Nrhs once L2 is built.
+
+  //////////////////////////////////////////////////////////////////////
+  // psi_coarse = P^dag (RAW fine null) -> Galerkin images that carry the
+  // near null content, and are free: A_c (P psi) = P A psi.
+  //////////////////////////////////////////////////////////////////////
+  MultiRHSBlockProject<LatticeFermionD> MrhsProjector;
+  MrhsProjector.Allocate(nbasis,FGrid,Coarse5d);
+  MrhsProjector.ImportBasis(AggregatesGCR.subspace);   // block orthonormal basis
+  // The basis now lives in the projector's BLAS_V (device); nothing after this
+  // point reads AggregatesGCR.subspace (DEPRECATED_CHECK moved above).  Free it: 60 fine
+  // fields = 20 GB host per rank, and 60 LRU-eligible device copies.
+  AggregatesGCR.subspace.clear(); AggregatesGCR.subspace.shrink_to_fit();
+
+  std::vector<CoarseVector> psi_coarse(nbasis,Coarse5d);
+  MrhsProjector.blockProject(rawNull,psi_coarse);      // RAW vectors in
+  rawNull.clear(); rawNull.shrink_to_fit();
+
+  GramGuard("psi_coarse",psi_coarse,Coarse5d);
+
+  //////////////////////////////////////////////////////////////////////
+  // STAGE TWO: L2 -> L3.
+  //
+  // The fine operator here is the L1 operator, which is natively multiRHS, so the
+  // multiRHS driver applies with no promotion adapter: its D+1 grid IS the
+  // batch grid the L1 operator is currently set to.
+  //////////////////////////////////////////////////////////////////////
+  Coordinate cclatt = clatt;
+  Coordinate Block2({4,4,2,4});        // L2->L3 blocking; banked optimum 2026-08-24 (env BLOCK2 overrides)
+  if ( getenv("BLOCK2") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK2")),Block2); GRID_ASSERT(Block2.size()==4); }
+  for(int d=0;d<4;d++){ GRID_ASSERT(clatt[d]%Block2[d]==0); cclatt[d]=clatt[d]/Block2[d]; }
+  std::cout << GridLogMessage << "Block2 " << Block2 << "  coarse-coarse lattice " << cclatt << std::endl;
+
+  Coordinate cc5latt({1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+  GridCartesian *CoarseCoarse5d = new GridCartesian(cc5latt,c5simd,c5mpi);
+
+  Coordinate ccmlatt({nrhs,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+  GridCartesian *CoarseCoarseMrhs = new GridCartesian(ccmlatt,cmsimd,cmmpi);
+
+  Coordinate ccblatt({batch,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+  GridCartesian *CoarseCoarseBatch = new GridCartesian(ccblatt,cmsimd,cmmpi);
+
+  // Coarsening deepens the tensor nest by one iScalar
+  typedef CoarseVector::vector_object                                         CoarseSiteObj;
+  typedef iScalar<CComplexS>                                                  CComplexS2;
+  typedef MultiGeneralCoarsenedOperator<CoarseSiteObj,CComplexS2,nbasis>    CoarseCoarseOperator;
+  typedef CoarseCoarseOperator::CoarseVector                                  CoarseCoarseVector;
+
+  NextToNearestStencilGeometry5D geom2(CoarseCoarse5d);
+
+  // RAW copy of the coarse null vectors, for the same reason as rawNull:
+  // the L2 CoarsenOperator block-orthonormalises its subspace in place, and
+  // the L3 basis must be defined by the vectors that still carry content.
+  std::vector<CoarseVector> rawPsi(nbasis,Coarse5d);
+  for(int k=0;k<nbasis;k++) rawPsi[k]=psi_coarse[k];
+
+  CoarseCoarseOperator CoarseOpL2(geom2,CoarseCoarse5d);
+  CoarseOpL2.SetGrid(CoarseCoarseBatch);
+
+  NonHermitianLinearOperator<CoarseOperator,CoarseVector> LinOpCoarse(CoarseOpPV);
+
+  std::cout << GridLogMessage << "*** L2 CoarsenOperator, batch "<<batch<<" ***" << std::endl;
+  CoarseOpL2.CoarsenOperator(LinOpCoarse,CoarseBatch,psi_coarse,CoarseCoarse5d);
+
+  //////////////////////////////////////////////////////////////////////
+  // Both operators to the solve Nrhs. The matrix elements are Nrhs
+  // independent and survive the change.
+  //////////////////////////////////////////////////////////////////////
+  CoarseOpPV.SetGrid(CoarseMrhs);
+  CoarseOpL2.SetGrid(CoarseCoarseMrhs);
+  std::cout << GridLogMessage << "L1 operator at Nrhs " << CoarseOpPV.Nrhs()
+	    << ",  L2 operator at Nrhs " << CoarseOpL2.Nrhs() << std::endl;
+
+  //////////////////////////////////////////////////////////////////////
+  // psi_cc from the RAW coarse null vectors, and the same guard
+  //////////////////////////////////////////////////////////////////////
+  MultiRHSBlockProject<CoarseVector> MrhsProjectorL2;
+  MrhsProjectorL2.Allocate(nbasis,Coarse5d,CoarseCoarse5d);
+  MrhsProjectorL2.ImportBasis(psi_coarse);            // block orthonormal basis
+
+  {
+    std::vector<CoarseCoarseVector> psi_cc(nbasis,CoarseCoarse5d);
+    MrhsProjectorL2.blockProject(rawPsi,psi_cc);      // RAW vectors in
+
+    GramGuard("psi_cc",psi_cc,CoarseCoarse5d);
+  }
+  rawPsi.clear(); rawPsi.shrink_to_fit();
+
+  //////////////////////////////////////////////////////////////////////
+  // Both coarse operators apply on their solve grids
+  //////////////////////////////////////////////////////////////////////
+  {
+    GridParallelRNG cRNG(Coarse5d); cRNG.SeedFixedIntegers({3,4,5,6});
+    CoarseVector cin(CoarseMrhs), cout_(CoarseMrhs);
+    random(cRNG,cin);
+    CoarseOpPV.M(cin,cout_);
+    std::cout << GridLogMessage << "L1 apply |in|^2 = " << norm2(cin)
+	      << "   |M in|^2 = " << norm2(cout_) << std::endl;
+    GRID_ASSERT( norm2(cout_) > 0.0 );
+
+    GridParallelRNG ccRNG(CoarseCoarse5d); ccRNG.SeedFixedIntegers({7,8,9,10});
+    CoarseCoarseVector ccin(CoarseCoarseMrhs), ccout(CoarseCoarseMrhs);
+    random(ccRNG,ccin);
+    CoarseOpL2.M(ccin,ccout);
+    std::cout << GridLogMessage << "L2 apply |in|^2 = " << norm2(ccin)
+	      << "   |M in|^2 = " << norm2(ccout) << std::endl;
+    GRID_ASSERT( norm2(ccout) > 0.0 );
+  }
+
+  //////////////////////////////////////////////////////////////////////
+  // STAGE THREE (part one): the dense bottom on L2.
+  //
+  // DenseCoarseMatrix is bilingual: it takes the elements through
+  // Geometry()/ExtractMatrix(), so the L2 operator serves directly. It does
+  // detect that a multiRHS op cannot apply on the D dimensional grid and
+  // skips its own certificate and VERIFY, so the equivalent check is done
+  // here instead, driving the L2 operator at Nrhs 1 through a slice.
+  //////////////////////////////////////////////////////////////////////
+  typedef DenseCoarseMatrix<CComplexS2,nbasis> DenseCC_t;
+  std::unique_ptr<DenseCC_t> DenseCC;
+
+  if ( getenv("DENSE_CC")==nullptr || atoi(getenv("DENSE_CC")) ) {
+
+    std::cout << GridLogMessage << "*** L3 dense bottom: import from the L2 operator ***" << std::endl;
+    DenseCC.reset(new DenseCC_t(CoarseCoarse5d));
+    DenseCC->Import(CoarseOpL2);
+
+    ////////////////////////////////////////////////////////////////////
+    // ||A Ainv x - x|| / ||x||, the check Import could not run itself
+    ////////////////////////////////////////////////////////////////////
+    Coordinate cc1latt({1,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+    GridCartesian *CoarseCoarseOne = new GridCartesian(cc1latt,cmsimd,cmmpi);
+
+    CoarseOpL2.SetGrid(CoarseCoarseOne);
+
+    CoarseCoarseVector x(CoarseCoarse5d),y(CoarseCoarse5d),z(CoarseCoarse5d);
+    GridParallelRNG dRNG(CoarseCoarse5d); dRNG.SeedFixedIntegers({11,12,13,14});
+    random(dRNG,x);
+
+    (*DenseCC)(x,y);                    // y = Ainv x
+
+    CoarseCoarseVector y1(CoarseCoarseOne),z1(CoarseCoarseOne);
+    InsertSliceFast(y,y1,0,0);
+    CoarseOpL2.M(y1,z1);                // z = A y
+    ExtractSliceFast(z,z1,0,0);
+
+    z = z - x;
+    RealD rel = std::sqrt(norm2(z)/norm2(x));
+    std::cout << GridLogMessage << "L3 dense: ||A Ainv x - x||/||x|| = " << rel << std::endl;
+    GRID_ASSERT( rel < 1.0e-2 );
+
+    CoarseOpL2.SetGrid(CoarseCoarseMrhs);
+    delete CoarseCoarseOne;
+  }
+
+  //////////////////////////////////////////////////////////////////////
+  // STAGE THREE (part two): the solves.
+  //
+  // Both operators are driven from the SAME objects at whatever Nrhs is
+  // asked for -- the matrix elements were built once and survive SetGrid --
+  // so single RHS and multiRHS are the same code path with a different grid.
+  //////////////////////////////////////////////////////////////////////
+  GRID_ASSERT(DenseCC != nullptr);   // the PGCR bottom is not ported yet
+
+  typedef PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> FineSmoother_t;
+
+  ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv);
+  TrivialPrecon<LatticeFermionD> simple_fine;
+  TrivialPrecon<CoarseVector>    simpleC;
+
+  auto RunSolve = [&](int nr)
+  {
+    std::cout << GridLogMessage << "**********************************************" << std::endl;
+    std::cout << GridLogMessage << " THREE-level solve, Nrhs = " << nr << std::endl;
+    // Device-memory budget BEFORE the solve (2026-08-28: NRHS=6 died in hipMalloc at the
+    // first fine-smoother history allocation -- reported as an asynchronous "memory
+    // access fault" unless AMD_SERIALIZE_KERNEL/COPY made it a clean OOM).  The outer
+    // mRHS GCR holds src, sol, r, Az and OuterMmax x (p,q) fine fields PER RHS; the fine
+    // smoother adds FineSmootherMmax x (p,q) once.  The MemoryManager LRU cap
+    // (--device-mem) must be BELOW what is physically left after the non-LRU allocations
+    // (comms buffers, dense slab, stencil buffers), or the device fills before anything
+    // is evicted.  Print the estimate, the LRU state, and the device's own free count.
+    {
+      uint64_t fieldBytes = (uint64_t)FGrid->lSites()*sizeof(typename LatticeFermionD::scalar_object);
+      double outerGB = (double)nr*(4 + 2*OuterMmax)*fieldBytes/1.0e9;
+      double smthGB  = (double)(2*FineSmootherMmax + 4)*fieldBytes/1.0e9;
+      std::cout << GridLogMessage << "Device budget: fine field " << fieldBytes/1.0e6 << " MB; outer GCR history "
+                << nr << " x (4 + 2 x " << OuterMmax << ") fields = " << outerGB << " GB; fine smoother history + temps ~ "
+                << smthGB << " GB; MemoryManager device LRU " << MemoryManager::DeviceCacheBytes()/1.0e9 << " GB now, cap "
+                << MemoryManager::DeviceMaxBytes/1.0e9 << " GB" << std::endl;
+      // Empty the device LRU: every setup-era Lattice copy (coarse null vectors,
+      // coarsening temporaries) goes back to host, so the solve's working set
+      // starts from a clean device and the cap applies to it alone.
+      MemoryManager::EvictAll();
+      // ...and release the allocation caches' held blocks (setup-era deviceVector
+      // scratch that is "free" to the caller but not to hipMalloc).
+      MemoryManager::DropCache();
+      MemoryManager::PrintBytes();
+      acceleratorMem();
+    }
+    std::cout << GridLogMessage << "**********************************************" << std::endl;
+
+    Coordinate cml({nr,1,clatt[0],clatt[1],clatt[2],clatt[3]});
+    Coordinate ccml({nr,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+    GridCartesian *CMrhs  = new GridCartesian(cml, cmsimd,cmmpi);
+    GridCartesian *CCMrhs = new GridCartesian(ccml,cmsimd,cmmpi);
+
+    CoarseOpPV.SetGrid(CMrhs);
+    CoarseOpL2.SetGrid(CCMrhs);
+
+    NonHermitianLinearOperator<CoarseOperator,CoarseVector>            LinOpC (CoarseOpPV);
+    NonHermitianLinearOperator<CoarseCoarseOperator,CoarseCoarseVector> LinOpCC(CoarseOpL2);
+
+    MrhsDenseCCSolve<DenseCC_t,CoarseCoarseVector> ccSolve(*DenseCC,nr);
+
+    ShiftedLinearOperator<CoarseVector> ShiftedC(CoarseSmootherShift, LinOpC);
+    if ( PowerIterations > 0 ) {
+      // Spectral edges the smoother polynomials must respect (see
+      // scripts/gcr_polynomial.py: |R_m|>1 beyond the edge = amplification).
+      PowerIteration<CoarseVector>   ("CoarseSmootherOp(shift="+std::to_string(CoarseSmootherShift)+")", ShiftedC, CMrhs, PowerIterations);
+      PowerIteration<CoarseVector>   ("CoarseOp(unshifted)", LinOpC, CMrhs, PowerIterations);
+      PowerIteration<LatticeFermionD>("FineSmootherOp(shift="+std::to_string(FineSmootherShift)+")", ShiftedPVdagM, Ddwf.FermionGrid(), PowerIterations);
+    }
+    PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
+      CoarseSmootherGCR(0.01,1,ShiftedC,simpleC,CoarseSmootherMmax,CoarseSmootherNstep);
+    CoarseSmootherGCR.Level(2); CoarseSmootherGCR.Name("Csmoother"); CoarseSmootherGCR.SetZeroGuess(1);
+
+    SwitchableSmoother<CoarseVector> CoarseSmootherSlot(CoarseSmootherGCR,"Csmoother GCR");
+    MrhsCoarseThreeLevelPrec<CoarseVector,CoarseCoarseVector>
+      L2to3Precon(LinOpC, CoarseSmootherSlot, MrhsProjectorL2, ccSolve,
+                  Coarse5d, CoarseCoarse5d, CCMrhs, nr);
+
+    PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
+      L2PGCR(CoarseSolverTol, CoarseSolverOrder/16, LinOpC, L2to3Precon, CoarseSolverMmax, 16);
+    L2PGCR.Level(2); L2PGCR.Name("Couter"); L2PGCR.SetZeroGuess(1);
+
+    FineSmoother_t SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,FineSmootherMmax,FineSmootherOrder);
+    SmootherGCR.Level(1); SmootherGCR.Name("Fsmoother"); SmootherGCR.SetZeroGuess(1);
+    SmootherGCR.LogCoefficients(SmootherCoeffLog);
+    CoarseSmootherGCR.LogCoefficients(SmootherCoeffLog);
+
+    SwitchableSmoother<LatticeFermionD> FineSmootherSlot(SmootherGCR,"Fsmoother GCR");
+    MrhsTwoLevelMG<LatticeFermionD,CoarseVector,SwitchableSmoother<LatticeFermionD> >
+      ThreeLevelPrecon(PVdagM, FineSmootherSlot, MrhsProjector, L2PGCR, Coarse5d, CMrhs);
+    ThreeLevelPrecon.SetSloppy   = SetFineSloppy;
+    ThreeLevelPrecon.SloppyComms = FineSloppyComms;
+
+    MrhsPGCRNonHermitian<LatticeFermionD>
+      L1PGCR(OuterTol,1000,PVdagM,ThreeLevelPrecon,OuterMmax,OuterNstep);
+    L1PGCR.Level(1); L1PGCR.Name("Fouter"); L1PGCR.SetZeroGuess(1);
+
+    //////////////////////////////////////////////////////////////////////
+    // Smoother modes.  Objects live for the duration of this RunSolve.
+    //////////////////////////////////////////////////////////////////////
+    std::unique_ptr<ChebyshevNonHermitianSmoother<LatticeFermionD> > FineCheb;
+    std::unique_ptr<ChebyshevNonHermitianSmoother<CoarseVector> >    CoarseCheb;
+    std::unique_ptr<GCRReplaySmoother<LatticeFermionD> > FineReplay;
+    std::unique_ptr<GCRReplaySmoother<CoarseVector> >    CoarseReplay;
+    GCRCoefficients recF, recC;
+    {
+      GCRCoefficients::Select sel = GCRCoefficients::Last;
+      if ( PolyRecordSelect=="first" ) sel = GCRCoefficients::First;
+      if ( PolyRecordSelect=="mean" )  sel = GCRCoefficients::Mean;
+      recF.select = sel; recC.select = sel;
+    }
+    std::unique_ptr<FreeMobius5DInverse<WilsonImplD> > Fms;
+    std::unique_ptr<FreeMobius5DInverse<WilsonImplD> > F1dag;
+    std::unique_ptr<FourierPVdagMSmoother>             FineFourier;
+    std::unique_ptr<FineSmoother_t>                    FineFourierGCR;   // fgcr: Fourier as the preconditioner of a low-order fine GCR ("2-deep" smoother)
+    if ( FineSmootherMode == "fourier" || FineSmootherMode == "fgcr" ) {
+      // Free PVdagM inverse in the Landau frame (config already rotated).
+      // Block build is (4Ls)^2 per local momentum in memory; Ls=24 at large
+      // local volume is heavy -- fine on the laptop harness.
+      std::vector<Complex> pbc = {1,1,1,1};
+      Fms.reset  (new FreeMobius5DInverse<WilsonImplD>(FGrid,Ls,M5,b,c,FourierSmootherMass,pbc));
+      F1dag.reset(new FreeMobius5DInverse<WilsonImplD>(FGrid,Ls,M5,b,c,1.0,pbc));
+      DaggerBlocks(*F1dag);
+      GainCap(*Fms,  SmootherGainCap);
+      GainCap(*F1dag,SmootherGainCap);
+      FineFourier.reset(new FourierPVdagMSmoother(*Fms,*F1dag));
+      if ( FineSmootherMode == "fourier" ) {
+        FineSmootherSlot.Set(*FineFourier,"Fsmoother Fourier");
+      } else {
+        // fgcr: identical to the "gcr" fine smoother but with the Fourier free-inverse
+        // as the preconditioner instead of simple_fine.  A low-order (FineSmootherOrder)
+        // GCR on ShiftedPVdagM, Fourier-preconditioned -- the outer GCR absorbs the
+        // non-normal transient that makes raw Fourier amplify (cf. (iii): S0 alone rho=13
+        // but S0-as-preconditioner converged FGMRES in 181).  Prints its residuals.
+        FineFourierGCR.reset(new FineSmoother_t(0.0,1,ShiftedPVdagM,*FineFourier,FineSmootherMmax,FineSmootherOrder));
+        FineFourierGCR->Level(1); FineFourierGCR->Name("Fsmoother-FGCR"); FineFourierGCR->SetZeroGuess(1);
+        FineSmootherSlot.Set(*FineFourierGCR,"Fsmoother FourierGCR");
+      }
+    }
+    if ( FineSmootherMode == "cheb" ) {
+      FineCheb.reset(new ChebyshevNonHermitianSmoother<LatticeFermionD>(FineChebLo,FineChebHi,FineSmootherOrder,ShiftedPVdagM));
+      FineCheb->Verbose = PolyVerbose; FineCheb->name = "Fsmoother";
+      FineSmootherSlot.Set(*FineCheb,"Fsmoother Chebyshev");
+    }
+    if ( CoarseSmootherMode == "cheb" ) {
+      CoarseCheb.reset(new ChebyshevNonHermitianSmoother<CoarseVector>(CoarseChebLo,CoarseChebHi,CoarseSmootherNstep,ShiftedC));
+      CoarseCheb->Verbose = PolyVerbose; CoarseCheb->name = "Csmoother";
+      CoarseSmootherSlot.Set(*CoarseCheb,"Csmoother Chebyshev");
+    }
+    // Recording window [PolyRecordStart, PolyRecordStart+PolyRecordIters).
+    // M3 (2026-08-26): the GCR polynomial changes fast over the first outer
+    // steps (per-call |r|/|r0| 0.0034 -> 0.017 over steps 1-4) and the MEAN
+    // of those is a poor smoother (replay 0.033-0.049 per call); record a
+    // settled window instead.
+    if ( PolyRecordStart == 0 ) {
+      if ( FineSmootherMode   == "replay" ) SmootherGCR.SetCoefficientRecorder(&recF);
+      if ( CoarseSmootherMode == "replay" ) CoarseSmootherGCR.SetCoefficientRecorder(&recC);
+    }
+    // Record -> replay, with optional periodic re-recording ("re-record, not
+    // fade away": HDCG refreshed its polynomial every 10 steps, tracking the
+    // evolving spectral content of the residual).  Schedule on outer steps:
+    //   [PolyRecordStart, +PolyRecordIters)          adaptive GCR, recording
+    //   then replay of the selected recorded call;
+    //   if PolyRefresh>0: every PolyRefresh steps, ONE adaptive recording
+    //   step, then replay of that call.
+    auto BuildReplays = [&](void){
+      if ( FineSmootherMode == "replay" ) {
+        SmootherGCR.SetCoefficientRecorder(nullptr);
+        recF.Flush(); recF.Report("Fsmoother");
+        FineReplay.reset(new GCRReplaySmoother<LatticeFermionD>(ShiftedPVdagM,recF));
+        FineReplay->Verbose = PolyVerbose; FineReplay->name = "Fsmoother";
+        FineSmootherSlot.Set(*FineReplay,"Fsmoother replay");
+        SmootherGCR.ReleaseHistory();
+      }
+      if ( CoarseSmootherMode == "replay" ) {
+        CoarseSmootherGCR.SetCoefficientRecorder(nullptr);
+        recC.Flush(); recC.Report("Csmoother");
+        CoarseReplay.reset(new GCRReplaySmoother<CoarseVector>(ShiftedC,recC));
+        CoarseReplay->Verbose = PolyVerbose; CoarseReplay->name = "Csmoother";
+        CoarseSmootherSlot.Set(*CoarseReplay,"Csmoother replay");
+        CoarseSmootherGCR.ReleaseHistory();
+      }
+    };
+    auto StartRecording = [&](int step){
+      if ( FineSmootherMode   == "replay" ) { { auto sel=recF.select; recF = GCRCoefficients(); recF.select=sel; } SmootherGCR.SetCoefficientRecorder(&recF);       FineSmootherSlot.Set(SmootherGCR,"Fsmoother GCR (recording)"); }
+      if ( CoarseSmootherMode == "replay" ) { { auto sel=recC.select; recC = GCRCoefficients(); recC.select=sel; } CoarseSmootherGCR.SetCoefficientRecorder(&recC); CoarseSmootherSlot.Set(CoarseSmootherGCR,"Csmoother GCR (recording)"); }
+      std::cout << GridLogMessage << "Smoother coefficient recording starts at outer step " << step << std::endl;
+    };
+    int switchStep = PolyRecordStart + PolyRecordIters;
+    L1PGCR.OnStep = [&](int step){
+      if ( FineSmootherMode != "replay" && CoarseSmootherMode != "replay" ) return;
+      if ( step == PolyRecordStart && PolyRecordStart > 0 ) StartRecording(step);
+      if ( step == switchStep ) { BuildReplays(); return; }
+      if ( PolyRefresh > 0 && step > switchStep ) {
+        int since = step - switchStep;
+        if ( since % PolyRefresh == 0 )     { StartRecording(step); }       // one adaptive, recorded step
+        if ( since % PolyRefresh == 1 )     { BuildReplays(); }             // then replay it
+      }
+    };
+    std::cout << GridLogMessage << "Smoother modes: fine " << FineSmootherMode << "  coarse " << CoarseSmootherMode
+              << (FineSmootherMode=="replay"||CoarseSmootherMode=="replay" ? "  (record outer steps "+std::to_string(PolyRecordStart)+".."+std::to_string(PolyRecordStart+PolyRecordIters)+")" : "")
+              << std::endl;
+
+    std::vector<LatticeFermionD> src(nr,FGrid), sol(nr,FGrid);
+    for(int r=0;r<nr;r++){ gaussian(RNG5,src[r]); sol[r]=Zero(); }
+
+    GridStopWatch w; w.Start();
+    L1PGCR(src,sol);
+    w.Stop();
+    std::cout << GridLogMessage << "3-level solve Nrhs "<<nr<<" total " << w.Elapsed()
+              << "  (per RHS: " << w.useconds()/1.0e6/nr << " s)" << std::endl;
+
+    // The outer operator is exact by policy; assert the state rather than
+    // trust it -- a preconditioner that failed to restore would surface here.
+    SetFineSloppy(0);
+    { LatticeFermionD Ax(FGrid); RealD worst=0.0;
+      for(int r=0;r<nr;r++){ PVdagM.Op(sol[r],Ax); Ax=Ax-src[r];
+        RealD rn=std::sqrt(norm2(Ax)/norm2(src[r]));
+        std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": rhs["<<r<<"] true residual = " << rn << std::endl;
+        worst=std::max(worst,rn); }
+      std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": worst-case residual = " << worst
+                << "   (exact-halo verification)" << std::endl;
+    }
+
+
+    // The operators borrow these grids and build a PaddedCell on them, so
+    // they must let go before the grids are destroyed.
+    CoarseOpPV.ReleaseGrid();
+    CoarseOpL2.ReleaseGrid();
+    delete CMrhs; delete CCMrhs;
+  };
+
+  RunSolve(nrhs);
+  if ( getenv("SOLVE_SRHS")==nullptr || atoi(getenv("SOLVE_SRHS")) ) RunSolve(1);
+
+  std::cout << GridLogMessage << "*** stage three complete: solves done ***" << std::endl;
+
+  Grid_finalize();
+}
diff --git a/systems/WorkArounds.txt b/systems/WorkArounds.txt
index 942e1ec3e..ee678cc6c 100644
--- a/systems/WorkArounds.txt
+++ b/systems/WorkArounds.txt
@@ -206,49 +206,58 @@ Tursa/nsight-sys
 --------------------------------------------------------------------
 
 ============================================================================
-2026-08-28  libfabric memory-registration-cache "memhooks" monitor -- DEFECTIVE
-           (libfabric issue #11451, filed by PB, HPE JIRA opened)
+Piz Daint
 ============================================================================
-STATUS: the ALPS/aarch64 case below is DEMONSTRATED (issue + both fixes verified
-there).  The Frontier case is a HYPOTHESIS as of 2026-08-28: the NO_TRANSLATION
-failure is observed, its attribution to the memhooks MR cache is by mechanism and
-by analogy, and the fix is UNTESTED here.  Test ladder, one knob per run, NRHS=12:
-  (1) FineSloppyComms=0  (2) FI_MR_CACHE_MONITOR=kdreg2  (3) FI_MR_CACHE_MAX_COUNT=0
-  (4) --disable-accelerator-aware-mpi build.
-REVISIT this entry with the result; if (2) does not cure it, remove kdreg2 from the
-job scripts' justification (it stays as OLCF's recommendation regardless).
-  2026-08-28 22:18, job 5372414: step (2) DID NOT cure it -- kdreg2 set (and the
-  environment already carried FI_MR_CACHE_MAX_COUNT=786432), sloppy comms OFF, still
-  NO_TRANSLATION at outer step 60 (Couter, Waitall count=2).
-  CORRECTION 2026-08-28 23:xx (fi_mr(3) man page): FI_MR_CACHE_MONITOR governs SYSTEM
-  memory only; device (HMEM_ROCR) registrations are monitored by
-  FI_MR_ROCR_CACHE_MONITOR_ENABLED=0|1.  Grid's comms window is hipMalloc'd device
-  memory, so the kdreg2 step tested nothing relevant -- the ALPS PLT bug (#11451) and
-  the Frontier NO_TRANSLATION are the same *library* but different monitors.
-  Attributed by gdb on a hung rank (5371826): main thread in MPI_Waitall <- Grid
-  0x5d539c3 (StencilSendToRecvFromComplete/CommsComplete class), i.e. the HALO
-  EXCHANGE; the two aborts' request counts (16 = fine Dhop, 2 = coarse PaddedCell
-  direction) agree.  Next: FI_MR_CACHE_MAX_COUNT=0 (all memory) and
-  FI_MR_ROCR_CACHE_MONITOR_ENABLED=0, one per cell.
 
-Symptom (Frontier, x86, Cray MPICH 8.1.x, ROCm 7.2): device-buffer MPI fails
-with
-   MPI_Waitall ... MPIDI_OFI_handle_cq_error: OFI poll failed
-   (ofi_events.c:MPIDI_OFI_handle_cq_error: Input/output error - NO_TRANSLATION)
-on a LIVE, never-freed hipMalloc buffer (the sloppy-comms compressed halo buffer,
-a static deviceVector) once the process does sustained hipMalloc/hipFree churn
-(MemoryManager eviction at NRHS>=12, EvictAll/DropCache).  Translations cached by
-the provider go stale.
-Symptom (ALPS/CSCS, aarch64 Grace/H200, cray-mpich 8.1.32, libfabric 1.22): the
+https://github.com/ofiwg/libfabric/issues/11451
+
+SEGFAULT in first MPI_Comm_dup on more than one node, after MPI_Init
+ 
+Cause: symptom (ALPS/CSCS, aarch64 Grace/H200, cray-mpich 8.1.32, libfabric 1.22): the
 memhooks monitor intercepts munmap by PATCHING THE PLT at MPI_Init
 (ofi_memhooks_start -> ofi_write_patch with a garbage data_size); the write overruns
-munmap@plt into the NEXT PLT entry (MPI_Comm_dup in Grid's binary), leaving
-`br x15` where an adrp belongs -> segfault on first MPI_Comm_dup.  Diagnosed with
-a hardware watchpoint on the PLT entry during MPI_Init; see the issue for the
+munmap@plt into the NEXT PLT entry (MPI_Comm_dup in Grid's binary)
+leaving
+`br x15` where an adrp belongs ->
+
+Diagnosed with a hardware watchpoint on the PLT entry during MPI_Init; see the issue for the
 gdb transcript.  Reproduce on one node with MPICH_SINGLE_HOST_ENABLED=0.
-Fix (either):   export FI_MR_CACHE_MONITOR=kdreg2      (kernel-driven invalidation; OLCF's
-                                                        own recommendation, NOT the default)
-                export FI_MR_CACHE_MAX_COUNT=0         (no registration cache at all)
+
+Fix (either):
+                export FI_MR_CACHE_MONITOR=kdreg2
+OR                                                 
+                export FI_MR_CACHE_MONITOR=disabled
+                export FI_MR_CACHE_MAX_COUNT=0
+
 On ALPS this is a CORRECTNESS requirement for MPI on Slingshot, not a tuning.
+
+
+
+============================================================================
+Frontier
+============================================================================
+
+a) FI_MR_CACHE_MONITOR=kdreg2 leads to Runtime MPI errors with "NO_TRANSLATION"
+
+Symptom (Frontier, x86, Cray MPICH 8.1.x, ROCm 7.2): device-buffer MPI fails with
+   MPI_Waitall ... MPIDI_OFI_handle_cq_error: OFI poll failed
+   (ofi_events.c:MPIDI_OFI_handle_cq_error: Input/output error - NO_TRANSLATION)
+on a LIVE, never-freed hipMalloc buffer
+
+
+b)              export FI_MR_CACHE_MONITOR=disabled
+                export FI_MR_CACHE_MAX_COUNT=0
+
+Leads to INCORRECT RESULTS
+
+https://github.com/ofiwg/libfabric/issues/12773
+
+https://github.com/ofiwg/libfabric/issues/12775
+
+Fix:
+
+export FI_HMEM_ROCR_USE_DMABUF=0
+
 Every systems/Frontier job and sourceme now sets kdreg2 on OLCF's recommendation;
 whether it resolves the Frontier NO_TRANSLATION is the pending test above.
+
diff --git a/tests/debug/Test_blockcyclic.cc b/tests/debug/Test_blockcyclic.cc
index e5b0d04ff..86fa8708b 100644
--- a/tests/debug/Test_blockcyclic.cc
+++ b/tests/debug/Test_blockcyclic.cc
@@ -20,7 +20,7 @@ Author: Peter Boyle <pboyle@bnl.gov>
 
 //////////////////////////////////////////////////////////////////////////////
 // Regression gate for BlockCyclicLayout -- stage 1 of the 2D distributed
-// dense inverse (documentation/DistributedDenseInverse2D.tex).
+// dense inverse.
 //
 // The layout is pure index arithmetic, so this test is EXHAUSTIVE rather
 // than statistical: every stage sweeps a battery of (N, nb, Pr, Pc)
diff --git a/tests/debug/Test_coarse.cc b/tests/debug/Test_coarse.cc
new file mode 100644
index 000000000..362cbc465
--- /dev/null
+++ b/tests/debug/Test_coarse.cc
@@ -0,0 +1,261 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_coarse.cc
+
+    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.
+
+    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 */
+
+//
+// MultiGeneralCoarsenedOperator against the existing mrhs coarse operator.
+//
+// The reference (the deprecated operator) is constructed on the D+1 grid; the
+// operator under test on the D dimensional grid, with SetGrid() adopting the
+// caller owned D+1 grid and building its padded cell and neighbour table from
+// the D dimensional stencil, with the Nrhs factor multiplied in.
+//
+// Both are given identical matrix elements, so any difference in the apply is
+// the restructured neighbour table. The same geometry object is passed to
+// both: the reference adds one to skip for the rhs direction, the operator
+// under test uses it as is on the D dimensional grid, so both describe the
+// same stencil over the D dimensions.
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const int nbasis = 8;
+
+typedef vSpinColourVector FineObj;
+typedef sTComplexD        CComplexT;   // unvectorised coarse space
+
+typedef DeprecatedMultiGeneralCoarsenedMatrix    <FineObj,CComplexT,nbasis> RefOperator;
+typedef MultiGeneralCoarsenedOperator<FineObj,CComplexT,nbasis> TestOperator;
+
+////////////////////////////////////////////////////////////////////////
+// Identical random matrix elements into both operators
+////////////////////////////////////////////////////////////////////////
+template<class OpA,class OpB>
+void SeedMatrixElements(OpA &A,OpB &B,int npoint,GridSerialRNG &sRNG)
+{
+  typedef typename OpA::calcMatrix calcMatrix;
+
+  deviceVector<calcMatrix> bufA,bufB;
+  for(int p=0;p<npoint;p++){
+
+    // MatrixPointOut/In carry one stencil point in lSite order whichever
+    // internal layout the operator uses.
+    A.MatrixPointOut(p,bufA);
+    B.MatrixPointOut(p,bufB);
+    GRID_ASSERT(bufA.size() == bufB.size());
+
+    int64_t sites = bufA.size();
+    std::vector<calcMatrix> host(sites);
+
+    ComplexD *w = (ComplexD *)&host[0];
+    int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD);
+    for(int64_t i=0;i<words;i++){
+      RealD re,im;
+      random(sRNG,re);
+      random(sRNG,im);
+      w[i] = ComplexD(re-0.5,im-0.5);
+    }
+
+    acceleratorCopyToDevice(&host[0],&bufA[0],sites*sizeof(calcMatrix));
+    A.MatrixPointIn(p,bufA);
+    B.MatrixPointIn(p,bufA);
+  }
+}
+
+template<class Op>
+RealD MatrixChecksum(Op &O,int npoint)
+{
+  typedef typename Op::calcMatrix calcMatrix;
+  RealD sum=0.0;
+  deviceVector<calcMatrix> buf;
+  for(int p=0;p<npoint;p++){
+    O.MatrixPointOut(p,buf);
+    int64_t sites = buf.size();
+    std::vector<calcMatrix> host(sites);
+    acceleratorCopyFromDevice(&buf[0],&host[0],sites*sizeof(calcMatrix));
+    ComplexD *w = (ComplexD *)&host[0];
+    int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD);
+    for(int64_t i=0;i<words;i++) sum += real(w[i])*real(w[i]) + imag(w[i])*imag(w[i]);
+  }
+  return sum;
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int nrhs = 4;
+
+  Coordinate clatt = GridDefaultLatt();
+  Coordinate csimd(Nd,1);                 // the coarse space is unvectorised
+  Coordinate cmpi  = GridDefaultMpi();
+
+  ////////////////////////////////////////////////
+  // D dimensional coarse grid, and D+1 for the reference
+  ////////////////////////////////////////////////
+  GridCartesian *CoarseD = new GridCartesian(clatt,csimd,cmpi);
+
+  std::cout << GridLogMessage << "coarse D   grid "; for(int d=0;d<Nd;d++) std::cout<<clatt[d]<<" ";
+  std::cout << "  Nsimd " << CoarseD->Nsimd() << std::endl;
+
+  ////////////////////////////////////////////////
+  // One geometry object for both, and one D+1 grid
+  // owned here and shared by both operators: fields
+  // conform only across a shared grid object.
+  ////////////////////////////////////////////////
+  NextToNearestStencilGeometry4D geom(CoarseD);
+
+  Coordinate mlatt(1,nrhs), msimd(1,1), mmpi(1,1);
+  for(int d=0;d<Nd;d++){
+    mlatt.push_back(clatt[d]);
+    msimd.push_back(csimd[d]);
+    mmpi .push_back(cmpi[d]);
+  }
+  GridCartesian *CoarseMulti = new GridCartesian(mlatt,msimd,mmpi);
+
+  TestOperator OpTest(geom,CoarseD);
+  RefOperator OpRef(geom,CoarseMulti);
+  OpTest.SetGrid(CoarseMulti);
+
+  std::cout << GridLogMessage << "coarse D+1 grid nrhs " << nrhs
+	    << "  Nsimd " << CoarseMulti->Nsimd() << std::endl;
+
+  std::cout << GridLogMessage << "npoint ref " << OpRef.geom.npoint
+	    << "   npoint test " << OpTest.geom.npoint << std::endl;
+  GRID_ASSERT(OpRef.geom.npoint == OpTest.geom.npoint);
+
+  int npoint = OpRef.geom.npoint;
+
+  ////////////////////////////////////////////////
+  // Identical matrix elements
+  ////////////////////////////////////////////////
+  GridSerialRNG sRNG; sRNG.SeedFixedIntegers(std::vector<int>({7,8,9,10}));
+  SeedMatrixElements(OpRef,OpTest,npoint,sRNG);
+
+  RealD ckRef = MatrixChecksum(OpRef,npoint);
+  RealD ckTest = MatrixChecksum(OpTest,npoint);
+  std::cout << GridLogMessage << "matrix element checksum ref " << ckRef
+	    << "  test " << ckTest << std::endl;
+  GRID_ASSERT( ckRef == ckTest );
+
+  ////////////////////////////////////////////////
+  // Same input, compare the applies
+  ////////////////////////////////////////////////
+  typedef RefOperator::CoarseVector CoarseVector;
+  // RNG on the D dimensional grid fills any D+1 field: the rhs direction is
+  // undistributed and divides cleanly. One RNG serves every Nrhs.
+  GridParallelRNG pRNG(CoarseD); pRNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+
+  CoarseVector in  (CoarseMulti); random(pRNG,in);
+  CoarseVector out1(CoarseMulti);
+  CoarseVector out2(CoarseMulti);
+  CoarseVector err (CoarseMulti);
+
+  OpRef.M(in,out1);
+  OpTest.M(in,out2);
+
+  err = out1 - out2;
+  std::cout << GridLogMessage << "|ref out|^2 = " << norm2(out1)
+	    << "   |test out|^2 = " << norm2(out2) << std::endl;
+  std::cout << GridLogMessage << "|ref - test|^2 = " << norm2(err) << std::endl;
+  GRID_ASSERT( norm2(out1) > 0.0 );
+  GRID_ASSERT( norm2(err) == 0.0 );
+
+  ////////////////////////////////////////////////
+  // SetGrid is idempotent on pointer identity
+  ////////////////////////////////////////////////
+  OpTest.SetGrid(CoarseMulti);
+  GRID_ASSERT( MatrixChecksum(OpTest,npoint) == ckTest );
+  OpTest.M(in,out2);
+  err = out1 - out2;
+  GRID_ASSERT( norm2(err) == 0.0 );
+  std::cout << GridLogMessage << "SetGrid idempotent on identity" << std::endl;
+
+  ////////////////////////////////////////////////
+  // Move to a different Nrhs and back. The matrix
+  // elements are Nrhs independent and must survive
+  // both the release and the rebuild.
+  ////////////////////////////////////////////////
+  // Nrhs 1 is the single RHS case through the multiRHS path, and each slice
+  // of the Nrhs 4 apply must come back unchanged.
+  OpTest.M(in,out2);
+  for(int nr=2;nr>=1;nr--){
+    Coordinate latt2(1,nr), simd2(1,1), mpi2(1,1);
+    for(int d=0;d<Nd;d++){
+      latt2.push_back(clatt[d]);
+      simd2.push_back(csimd[d]);
+      mpi2 .push_back(cmpi[d]);
+    }
+    GridCartesian *CoarseMulti2 = new GridCartesian(latt2,simd2,mpi2);
+
+    OpTest.SetGrid(CoarseMulti2);
+    GRID_ASSERT( OpTest.Nrhs() == nr );
+    GRID_ASSERT( MatrixChecksum(OpTest,npoint) == ckTest );
+
+    CoarseVector in2 (CoarseMulti2);
+    CoarseVector out(CoarseMulti2);
+    for(int r=0;r<nr;r++){
+      CoarseVector slice(CoarseD);
+      ExtractSliceFast(slice,in,r,0);
+      InsertSliceFast(slice,in2,r,0);
+    }
+    OpTest.M(in2,out);
+
+    RealD sdiff=0.0;
+    for(int r=0;r<nr;r++){
+      CoarseVector a(CoarseD),b(CoarseD),e(CoarseD);
+      ExtractSliceFast(a,out ,r,0);
+      ExtractSliceFast(b,out2,r,0);
+      e = a-b;
+      sdiff += norm2(e);
+    }
+    // Not bit exact: a different Nrhs is a different GEMM shape
+    std::cout << GridLogMessage << "Nrhs " << nr << " slices agree with Nrhs "
+	      << nrhs << " : |diff|^2/|out|^2 = " << sdiff/norm2(out) << std::endl;
+    GRID_ASSERT( norm2(out) > 0.0 );
+    GRID_ASSERT( sdiff/norm2(out) < 1.0e-20 );
+
+    OpTest.ReleaseGrid();
+    GRID_ASSERT( MatrixChecksum(OpTest,npoint) == ckTest );  // survives release
+
+    delete CoarseMulti2;
+  }
+
+  OpTest.SetGrid(CoarseMulti);
+  GRID_ASSERT( OpTest.Nrhs() == nrhs );
+
+  OpTest.M(in,out2);
+  err = out1 - out2;
+  std::cout << GridLogMessage << "after Nrhs 4 -> 2 -> 1 -> release -> 4, |ref - test|^2 = "
+	    << norm2(err) << std::endl;
+  GRID_ASSERT( norm2(err) == 0.0 );
+
+  std::cout << GridLogMessage << "Test_coarse: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}
diff --git a/tests/debug/Test_coarse_coarsen.cc b/tests/debug/Test_coarse_coarsen.cc
new file mode 100644
index 000000000..28fa50c9f
--- /dev/null
+++ b/tests/debug/Test_coarse_coarsen.cc
@@ -0,0 +1,279 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_coarse_coarsen.cc
+
+    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.
+
+    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 */
+
+//
+// Coarsen the same fine operator two ways and compare the matrix elements.
+//
+//   ref  : DeprecatedMultiGeneralCoarsenedMatrix::CoarsenOperator, vectorised
+//          coarse space matched to the fine SIMD layout, single RHS fine
+//          applications
+//   test : MultiGeneralCoarsenedOperator::CoarsenOperator on the D+1 grid,
+//          unvectorised sComplexD coarse space, the batch of phased basis
+//          vectors carried in the rhs direction and applied through
+//          MrhsPromotedOperator
+//
+// BLAS_A is written by GridtoBLAS in lSite order, which does not depend on
+// the SIMD layout, so the two are directly comparable.
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const int nbasis = 8;
+const int batch  = 9;
+
+typedef vSpinColourVector FineObj;
+typedef vTComplex         CComplexV;   // vectorised coarse space, reference
+typedef sTComplexD        CComplexS;   // unvectorised coarse space, under test
+
+typedef DeprecatedMultiGeneralCoarsenedMatrix    <FineObj,CComplexV,nbasis> RefOperator;
+typedef MultiGeneralCoarsenedOperator<FineObj,CComplexS,nbasis> TestOperator;
+
+template<class Op>
+void ReadMatrix(Op &O,int npoint,std::vector<std::vector<typename Op::calcMatrix> > &host)
+{
+  typedef typename Op::calcMatrix calcMatrix;
+  host.resize(npoint);
+  // MatrixPointOut carries one stencil point in lSite order whichever internal
+  // layout the operator uses.
+  deviceVector<calcMatrix> buf;
+  for(int p=0;p<npoint;p++){
+    O.MatrixPointOut(p,buf);
+    int64_t sites = buf.size();
+    host[p].resize(sites);
+    acceleratorCopyFromDevice(&buf[0],&host[p][0],sites*sizeof(calcMatrix));
+  }
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  Coordinate flatt = GridDefaultLatt();
+  Coordinate fsimd = GridDefaultSimd(Nd,vComplexD::Nsimd());
+  Coordinate fmpi  = GridDefaultMpi();
+
+  Coordinate block({2,2,2,2});
+  Coordinate clatt(Nd);
+  for(int d=0;d<Nd;d++){
+    GRID_ASSERT(flatt[d]%block[d]==0);
+    clatt[d] = flatt[d]/block[d];
+  }
+
+  GridCartesian         *FineGrid   = new GridCartesian(flatt,fsimd,fmpi);
+  GridRedBlackCartesian *FrbGrid    = SpaceTimeGrid::makeFourDimRedBlackGrid(FineGrid);
+
+  ////////////////////////////////////////////////
+  // Reference coarse space: SIMD layout matched to the fine
+  ////////////////////////////////////////////////
+  Coordinate cvsimd = GridDefaultSimd(Nd,CComplexV::Nsimd());
+  GridCartesian *CoarseV = new GridCartesian(clatt,cvsimd,fmpi);
+
+  // The reference puts all of the SIMD in the rhs direction. CoarsenOperator
+  // does not use the multiRHS grid; it only sizes BLAS_A, so one lane of rhs suffices.
+  int nrhs_ref = CComplexV::Nsimd();
+  Coordinate reflatt(1,nrhs_ref),refsimd(1,CComplexV::Nsimd()),refmpi(1,1);
+  for(int d=0;d<Nd;d++){
+    reflatt.push_back(clatt[d]);
+    refsimd.push_back(1);
+    refmpi .push_back(fmpi[d]);
+  }
+  GridCartesian *CoarseVMulti = new GridCartesian(reflatt,refsimd,refmpi);
+
+  ////////////////////////////////////////////////
+  // Coarse space under test: unvectorised
+  ////////////////////////////////////////////////
+  Coordinate cssimd(Nd,1);
+  GridCartesian *CoarseS = new GridCartesian(clatt,cssimd,fmpi);
+
+  Coordinate cmlatt(1,batch),cmsimd(1,1),cmmpi(1,1);
+  for(int d=0;d<Nd;d++){
+    cmlatt.push_back(clatt[d]);
+    cmsimd.push_back(1);
+    cmmpi .push_back(fmpi[d]);
+  }
+  GridCartesian *CoarseSMulti = new GridCartesian(cmlatt,cmsimd,cmmpi);
+
+  // D+1 fine grid carrying the batch
+  Coordinate fmlatt(1,batch), fmsimd(1,1), fmmpi(1,1);
+  for(int d=0;d<Nd;d++){
+    fmlatt.push_back(flatt[d]);
+    fmsimd.push_back(fsimd[d]);
+    fmmpi .push_back(fmpi[d]);
+  }
+  GridCartesian *FineGridMulti = new GridCartesian(fmlatt,fmsimd,fmmpi);
+
+  std::cout << GridLogMessage << "fine "<<flatt<<"   coarse "<<clatt<<"   batch "<<batch<<std::endl;
+  std::cout << GridLogMessage << "Nsimd   fine "<<FineGrid->Nsimd()
+	    << "   coarse ref "<<CoarseV->Nsimd()
+	    << "   coarse test "<<CoarseS->Nsimd()<<std::endl;
+
+  ////////////////////////////////////////////////
+  // Fine operator
+  ////////////////////////////////////////////////
+  GridParallelRNG pRNG(FineGrid); pRNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+  LatticeGaugeFieldD Umu(FineGrid); SU<Nc>::HotConfiguration(pRNG,Umu);
+
+  RealD mass = 0.1;
+  WilsonFermionD Dw(Umu,*FineGrid,*FrbGrid,mass);
+  MdagMLinearOperator<WilsonFermionD,LatticeFermionD> HermOp(Dw);
+
+  ////////////////////////////////////////////////
+  // One random subspace, two copies: CoarsenOperator orthogonalises in place
+  ////////////////////////////////////////////////
+  Aggregation<FineObj,CComplexV,nbasis> Subspace(CoarseV,FineGrid,0);
+  std::vector<LatticeFermionD> subspace(nbasis,FineGrid);
+  for(int i=0;i<nbasis;i++){
+    random(pRNG,Subspace.subspace[i]);
+    subspace[i] = Subspace.subspace[i];
+  }
+
+  NextToNearestStencilGeometry4D geomV(CoarseV);
+  NextToNearestStencilGeometry4D geomS(CoarseS);
+
+  ////////////////////////////////////////////////
+  // Reference coarsening, matched layouts
+  ////////////////////////////////////////////////
+  RefOperator OpRef(geomV,CoarseVMulti);
+  std::cout << GridLogMessage << "reference CoarsenOperator" << std::endl;
+  OpRef.CoarsenOperator(HermOp,Subspace,CoarseV);
+
+  ////////////////////////////////////////////////
+  // Coarsening under test, D+1 fine applications, unvectorised coarse
+  ////////////////////////////////////////////////
+  TestOperator OpTest(geomS,CoarseS);
+  OpTest.SetGrid(CoarseSMulti);
+
+  MrhsPromotedOperator<LatticeFermionD> MrhsHermOp(HermOp,FineGrid,batch);
+
+  std::cout << GridLogMessage << "test CoarsenOperator (D+1)" << std::endl;
+  OpTest.CoarsenOperator(MrhsHermOp,FineGridMulti,subspace,CoarseS);
+
+  ////////////////////////////////////////////////
+  // Compare matrix elements
+  ////////////////////////////////////////////////
+  int npoint = OpRef.geom.npoint;
+  GRID_ASSERT(npoint == OpTest.geom.npoint);
+  typedef RefOperator::calcMatrix calcMatrix;
+  std::vector<std::vector<calcMatrix> > A1,A2;
+  ReadMatrix(OpRef,npoint,A1);
+  ReadMatrix(OpTest,npoint,A2);
+
+  RealD num=0.0, den=0.0;
+  for(int p=0;p<npoint;p++){
+    GRID_ASSERT(A1[p].size()==A2[p].size());
+    ComplexD *w1 = (ComplexD *)&A1[p][0];
+    ComplexD *w2 = (ComplexD *)&A2[p][0];
+    int64_t words = A1[p].size()*sizeof(calcMatrix)/sizeof(ComplexD);
+    for(int64_t i=0;i<words;i++){
+      ComplexD d = w1[i]-w2[i];
+      num += real(d)*real(d)+imag(d)*imag(d);
+      den += real(w1[i])*real(w1[i])+imag(w1[i])*imag(w1[i]);
+    }
+  }
+  std::cout << GridLogMessage << "|A_ref|^2 = " << den << std::endl;
+  std::cout << GridLogMessage << "|A_ref - A_test|^2 / |A_ref|^2 = " << num/den << std::endl;
+  GRID_ASSERT( den > 0.0 );
+  GRID_ASSERT( num/den < 1.0e-20 );
+
+  ////////////////////////////////////////////////
+  // Same coarsening through the single RHS variant: no mrhs packing, the
+  // batch assembled on the coarse side by the mixed blockProject. Block
+  // Gram-Schmidt is idempotent so the subspace may be reused in place.
+  ////////////////////////////////////////////////
+  TestOperator OpTestSrhs(geomS,CoarseS);
+  OpTestSrhs.SetGrid(CoarseSMulti);
+
+  std::cout << GridLogMessage << "test CoarsenOperator (single RHS fine op)" << std::endl;
+  OpTestSrhs.CoarsenOperator(HermOp,subspace,CoarseS,batch);
+
+  std::vector<std::vector<calcMatrix> > A3;
+  ReadMatrix(OpTestSrhs,npoint,A3);
+
+  RealD nums=0.0;
+  for(int p=0;p<npoint;p++){
+    GRID_ASSERT(A1[p].size()==A3[p].size());
+    ComplexD *w1 = (ComplexD *)&A1[p][0];
+    ComplexD *w3 = (ComplexD *)&A3[p][0];
+    int64_t words = A1[p].size()*sizeof(calcMatrix)/sizeof(ComplexD);
+    for(int64_t i=0;i<words;i++){
+      ComplexD d = w1[i]-w3[i];
+      nums += real(d)*real(d)+imag(d)*imag(d);
+    }
+  }
+  std::cout << GridLogMessage << "|A_ref - A_test_srhs|^2 / |A_ref|^2 = " << nums/den << std::endl;
+  GRID_ASSERT( nums/den < 1.0e-20 );
+
+  ////////////////////////////////////////////////
+  // GetMatrix/SetMatrix round trip through the BLAS layout array. This is
+  // how a bilingual DenseCoarseMatrix will retrieve the elements, and it
+  // needs no SetGrid since the matrix elements are Nrhs independent.
+  ////////////////////////////////////////////////
+  {
+    typedef TestOperator::CoarseMatrix CoarseMatrixS;
+    std::vector<CoarseMatrixS> Aget(npoint,CoarseS);
+    for(int p=0;p<npoint;p++) OpTest.GetMatrix(p,Aget);
+
+    TestOperator OpTestRoundTrip(geomS,CoarseS);
+    for(int p=0;p<npoint;p++) OpTestRoundTrip.SetMatrix(p,Aget);
+
+    std::vector<std::vector<calcMatrix> > A4;
+    ReadMatrix(OpTestRoundTrip,npoint,A4);
+
+    RealD numrt=0.0;
+    for(int p=0;p<npoint;p++){
+      GRID_ASSERT(A2[p].size()==A4[p].size());
+      ComplexD *w2 = (ComplexD *)&A2[p][0];
+      ComplexD *w4 = (ComplexD *)&A4[p][0];
+      int64_t words = A2[p].size()*sizeof(calcMatrix)/sizeof(ComplexD);
+      for(int64_t i=0;i<words;i++){
+        ComplexD d = w2[i]-w4[i];
+        numrt += real(d)*real(d)+imag(d)*imag(d);
+      }
+    }
+    std::cout << GridLogMessage << "GetMatrix/SetMatrix round trip |diff|^2 = " << numrt << std::endl;
+    GRID_ASSERT( numrt == 0.0 );
+  }
+
+  ////////////////////////////////////////////////
+  // and the operator under test applies the matrix it just built
+  ////////////////////////////////////////////////
+  typedef TestOperator::CoarseVector CoarseVectorS;
+  GridParallelRNG cRNG(CoarseS); cRNG.SeedFixedIntegers(std::vector<int>({5,6,7,8}));
+  CoarseVectorS in(CoarseSMulti);  random(cRNG,in);
+  CoarseVectorS out(CoarseSMulti);
+
+  OpTest.M(in,out);
+  std::cout << GridLogMessage << "|in|^2 = " << norm2(in)
+	    << "   |M_test in|^2 = " << norm2(out) << std::endl;
+  GRID_ASSERT( norm2(out) > 0.0 );
+
+  std::cout << GridLogMessage << "Test_coarse_coarsen: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}
diff --git a/tests/debug/Test_mrhs_deflation.cc b/tests/debug/Test_mrhs_deflation.cc
new file mode 100644
index 000000000..e500edaf8
--- /dev/null
+++ b/tests/debug/Test_mrhs_deflation.cc
@@ -0,0 +1,89 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_mrhs_deflation.cc
+
+    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 */
+#include <Grid/Grid.h>
+
+// MultiRHSDeflation: the D+1 multiRHS interface (one permutation pass each
+// way) must reproduce the vector-of-D-fields interface bit for bit up to
+// the summation order of the GEMMs.  Random "eigenvectors" and values: the
+// deflation formula G = E (E^dag R)/lambda does not care that they are not
+// eigenpairs of anything.
+
+using namespace Grid;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int nbasis = 8;
+  const int nev    = 12;
+  const int nrhs   = 5;
+
+  typedef iVector<sTComplexD,nbasis> siteVector;
+  typedef Lattice<siteVector>        CoarseVector;
+
+  // Unvectorised D and D+1 coarse grids, as MGCoarseGrids builds them
+  Coordinate latt({1,4,4,4,4});
+  Coordinate simd({1,1,1,1,1});
+  Coordinate mpi = GridDefaultMpi();
+  Coordinate mpi5({1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  GridCartesian Coarse5d(latt,simd,mpi5);
+
+  Coordinate latt6({nrhs,1,4,4,4,4});
+  Coordinate simd6({1,1,1,1,1,1});
+  Coordinate mpi6({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  GridCartesian CoarseMrhs(latt6,simd6,mpi6);
+
+  GridParallelRNG RNG(&Coarse5d); RNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+
+  std::vector<CoarseVector> evec(nev,&Coarse5d);
+  std::vector<RealD>        eval(nev);
+  for(int e=0;e<nev;e++){ random(RNG,evec[e]); eval[e] = 0.5 + 0.1*e; }
+
+  std::vector<CoarseVector> src(nrhs,&Coarse5d), guess(nrhs,&Coarse5d);
+  for(int r=0;r<nrhs;r++) random(RNG,src[r]);
+
+  MultiRHSDeflation<CoarseVector> Deflator;
+  Deflator.ImportEigenBasis(evec,eval);
+
+  // Reference: the vector interface
+  Deflator.DeflateSources(src,guess);
+
+  // The D+1 interface on the same sources
+  CoarseVector src_mrhs(&CoarseMrhs), guess_mrhs(&CoarseMrhs);
+  for(int r=0;r<nrhs;r++) InsertSliceFast(src[r],src_mrhs,r,0);
+  guess_mrhs = Zero();
+  Deflator.DeflateSources(src_mrhs,guess_mrhs);
+
+  RealD worst = 0.0;
+  CoarseVector g(&Coarse5d), d(&Coarse5d);
+  for(int r=0;r<nrhs;r++){
+    ExtractSliceFast(g,guess_mrhs,r,0);
+    d = g - guess[r];
+    RealD rel = std::sqrt(norm2(d)/norm2(guess[r]));
+    std::cout << GridLogMessage << "rhs " << r << " ||guess_mrhs - guess||/||guess|| = " << rel << std::endl;
+    worst = std::max(worst,rel);
+  }
+  std::cout << GridLogMessage << "MultiRHSDeflation D+1 vs vector interface: worst " << worst << std::endl;
+  GRID_ASSERT( worst < 1.0e-12 );
+  std::cout << GridLogMessage << "PASS" << std::endl;
+
+  Grid_finalize();
+  return 0;
+}
diff --git a/tests/debug/Test_sloppy_dagger.cc b/tests/debug/Test_sloppy_dagger.cc
new file mode 100644
index 000000000..38e3d212e
--- /dev/null
+++ b/tests/debug/Test_sloppy_dagger.cc
@@ -0,0 +1,207 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_sloppy_dagger.cc
+
+    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 */
+
+//
+// Isolate a sloppy-comms dagger defect: on a CPU/GEN multi-rank build the
+// sloppy dagger halo can be wrong by ~5% of norm2, deterministically, which
+// Benchmark_dwf's sloppy pass catches as a failed dagger Cshift check.
+//
+// The sloppy compressed-buffer pool (StencilBuffer::DeviceCommBuf) is a
+// SHARED STATIC across all stencil objects, so scenarios contaminate each
+// other inside one process: this program runs exactly ONE scenario per
+// invocation, selected with --seq <name>:
+//
+//   fresh-dag       sloppy dagger is the operator's first call
+//   nodag-dag       sloppy nodag x3 (one source), then sloppy dagger (same source)
+//   nodag-dag-fresh nodag x3 (fresh source each), then dagger (fresh source)
+//   nodag-dag-newsrc nodag x3 (one source), then dagger (different source)
+//   dag-nodag-dag   dagger, nodag x3, dagger  (dagger-first history)
+//
+// Every call is checked against a separate never-sloppy reference
+// operator (exact Dhop == the Cshift construction at 1e-31, certified by
+// Benchmark_dwf's exact pass).  The last dagger's error is fingerprinted
+// by slice along each MPI-decomposed direction.
+//
+//   OMP_NUM_THREADS=1 mpirun -n 2 ./Test_sloppy_dagger --grid 8.8.8.16 \
+//        --mpi 1.1.1.2 --seq nodag-dag
+//
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::string seq("nodag-dag");
+  if( GridCmdOptionExists(argv,argv+argc,"--seq") )
+    seq = GridCmdOptionPayload(argv,argv+argc,"--seq");
+
+  const int Ls=8;
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+                                                                   GridDefaultSimd(Nd,vComplex::Nsimd()),
+                                                                   GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers({1,2,3,4});
+  GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8});
+
+  LatticeGaugeField Umu(UGrid);
+  SU<Nc>::HotConfiguration(RNG4,Umu);
+
+  RealD mass=0.1, M5=1.8;
+  DomainWallFermionD Dref(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);  // never sloppy
+  DomainWallFermionD Dw  (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);  // under test
+  Dw.SloppyComms(1);
+
+  Coordinate mpi = GridDefaultMpi();
+  LatticeFermionD src(FGrid), exact(FGrid), sloppy(FGrid), err(FGrid);
+
+  // One call: fresh or reused source, dag or not; always checked.
+  LatticeFermionD srcA(FGrid); gaussian(RNG5,srcA);
+  auto call = [&](int dag, int fresh, const char *tag){
+    if ( fresh ) gaussian(RNG5,src); else src = srcA;
+    Dref.Dhop(src,exact,dag);
+    Dw.Dhop(src,sloppy,dag);
+    err = sloppy - exact;
+    std::cout << GridLogMessage << "SEQ[" << seq << "] " << tag
+              << (dag?" DAG   ":" NODAG ") << " norm2(err) " << norm2(err) << std::endl;
+  };
+  // First bad site: coordinate + exact-vs-sloppy hex words (PB: look for
+  // expected/actual mismatching in the high-order half word).
+  auto firstbad = [&](LatticeFermionD &ex, LatticeFermionD &sl){
+    typedef LatticeFermionD::scalar_object sobj;
+    Coordinate gdims = FGrid->GlobalDimensions();
+    int64_t gsites = FGrid->gSites();
+    for(int64_t g=0; g<gsites; g++){
+      Coordinate gcoor(5);
+      Lexicographic::CoorFromIndex(gcoor,g,gdims);
+      sobj se, ss_;
+      peekSite(se,ex,gcoor);
+      peekSite(ss_,sl,gcoor);
+      uint64_t *we=(uint64_t *)&se;
+      uint64_t *ws=(uint64_t *)&ss_;
+      int nw = sizeof(sobj)/8;
+      for(int w=0;w<nw;w++){
+        double de=((double *)we)[w], ds=((double *)ws)[w];
+        if ( fabs(de-ds) > 1.0e-4 ) {
+          std::cout << GridLogMessage << "FIRSTBAD site " << gcoor
+                    << " word " << w << " (spin "<<(w/6)%4<<" col "<<(w/2)%3<<" reim "<<(w%2)<<")"
+                    << "  exact " << std::hex << we[w]
+                    << "  sloppy " << ws[w] << std::dec << std::endl;
+          for(int w2=0;w2<8&&w2<nw;w2++)
+            std::cout << GridLogMessage << "   word["<<w2<<"] exact "<<std::hex<<we[w2]
+                      <<" sloppy "<<ws[w2]<<std::dec<<std::endl;
+          return;
+        }
+      }
+    }
+    std::cout << GridLogMessage << "FIRSTBAD: none above 1e-4" << std::endl;
+  };
+  auto fingerprint = [&](void){
+    for(int mu=0;mu<Nd;mu++){
+      if ( mpi[mu] == 1 ) continue;
+      std::vector<RealD> sn;
+      sliceNorm(sn,err,mu+1);
+      std::cout << GridLogMessage << "SEQ[" << seq << "] err by slice of dim " << mu << ":";
+      for(int t=0;t<(int)sn.size();t++) std::cout << " " << sn[t];
+      std::cout << std::endl;
+    }
+  };
+
+  if        ( seq == "fresh-dag" ) {
+    call(1,0,"call0");
+  } else if ( seq == "nodag-dag" ) {
+    call(0,0,"call0"); call(0,0,"call1"); call(0,0,"call2");
+    call(1,0,"call3");
+  } else if ( seq == "nodag-dag-fresh" ) {
+    call(0,1,"call0"); call(0,1,"call1"); call(0,1,"call2");
+    call(1,1,"call3");
+  } else if ( seq == "nodag-dag-newsrc" ) {
+    call(0,0,"call0"); call(0,0,"call1"); call(0,0,"call2");
+    call(1,1,"call3");
+  } else if ( seq == "dag-nodag-dag" ) {
+    call(1,0,"call0");
+    call(0,0,"call1"); call(0,0,"call2"); call(0,0,"call3");
+    call(1,0,"call4");
+  } else if ( seq == "noleave" ) {
+    // NO exact-operator call between sloppy calls: references precomputed.
+    LatticeFermionD refnodag(FGrid), refdag(FGrid);
+    Dref.Dhop(srcA,refnodag,DaggerNo);
+    Dref.Dhop(srcA,refdag,DaggerYes);
+    for(int i=0;i<2;i++){
+      Dw.Dhop(srcA,sloppy,DaggerNo);
+      err = sloppy - refnodag;
+      std::cout << GridLogMessage << "SEQ[noleave] call" << i << " NODAG  norm2(err) " << norm2(err) << std::endl;
+      if ( i==0 ) firstbad(refnodag,sloppy);
+    }
+    Dw.Dhop(srcA,sloppy,DaggerYes);
+    err = sloppy - refdag;
+    std::cout << GridLogMessage << "SEQ[noleave] call3 DAG    norm2(err) " << norm2(err) << std::endl;
+  } else if ( seq == "twoop" ) {
+    // A SECOND sloppy operator shares the static pool; sequence runs on it
+    // with interleaved exact references (as the honest runs had).
+    DomainWallFermionD Dw2(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+    Dw2.SloppyComms(1);
+    Dw.Dhop(srcA,sloppy,DaggerNo);   // prime the FIRST sloppy op's state
+    for(int i=0;i<3;i++){
+      Dref.Dhop(srcA,exact,DaggerNo);
+      Dw2.Dhop(srcA,sloppy,DaggerNo);
+      err = sloppy - exact;
+      std::cout << GridLogMessage << "SEQ[twoop] call" << i << " NODAG  norm2(err) " << norm2(err) << std::endl;
+    }
+    Dref.Dhop(srcA,exact,DaggerYes);
+    Dw2.Dhop(srcA,sloppy,DaggerYes);
+    err = sloppy - exact;
+    std::cout << GridLogMessage << "SEQ[twoop] call3 DAG    norm2(err) " << norm2(err) << std::endl;
+  } else if ( seq == "self-exact" ) {
+    // Exact call on the SAME object (Dw), opposite sense, then sloppy:
+    // distinguishes per-object from shared state.
+    LatticeFermionD refnodag(FGrid);
+    Dref.Dhop(srcA,refnodag,DaggerNo);              // reference only, early
+    Dw.SloppyComms(0);
+    Dw.Dhop(srcA,sloppy,DaggerYes);                 // exact DAG on Dw itself
+    Dw.SloppyComms(1);
+    Dw.Dhop(srcA,sloppy,DaggerNo);                  // sloppy NODAG: sense mismatch
+    err = sloppy - refnodag;
+    std::cout << GridLogMessage << "SEQ[self-exact] sloppy NODAG after own exact DAG: " << norm2(err) << std::endl;
+  } else if ( seq == "cold" ) {
+    // NO exact Dhop anywhere before the sloppy calls.
+    LatticeFermionD outn(FGrid), outd(FGrid);
+    Dw.Dhop(srcA,outn,DaggerNo);
+    Dw.Dhop(srcA,outd,DaggerYes);
+    LatticeFermionD refnodag(FGrid), refdag(FGrid);
+    Dref.Dhop(srcA,refnodag,DaggerNo);
+    Dref.Dhop(srcA,refdag,DaggerYes);
+    err = outn - refnodag;
+    std::cout << GridLogMessage << "SEQ[cold] first-ever sloppy NODAG: " << norm2(err) << std::endl;
+    err = outd - refdag;
+    std::cout << GridLogMessage << "SEQ[cold] then      sloppy DAG:   " << norm2(err) << std::endl;
+  } else {
+    std::cout << GridLogMessage << "unknown --seq " << seq << std::endl;
+  }
+  fingerprint();
+
+  Grid_finalize();
+}