Update for bug report

Grid/algorithms/blas/BatchedBlas.h

@@ -995,7 +995,7 @@ public:
   ///////////////////////////////////////////////////////////////////////////////////
   // Explicit-leading-dimension complex double GEMM.  Mirror of the ComplexF
   // overload above; motivating use is the fp64 distributed recursive Schur
-  // inversion (RecursiveSchurInverse), whose operands are column windows of
+  // inversion (the distributed Schur recursion), whose operands are column windows of
   // larger row-slab allocations.
   ///////////////////////////////////////////////////////////////////////////////////
   void gemmBatched(GridBLASOperation_t OpA,

Grid/algorithms/iterative/PowerMethod.h

@@ -43,4 +43,69 @@ template<class Field> class PowerMethod
     return evalMaxApprox;
   }
 };
+
+// Non-Hermitian sibling of PowerMethod.  Drives Op() -- the operator the
+// caller actually applies -- rather than HermOp(), so it returns |lambda_max|
+// of a NON-Hermitian operator and reports the spectral-edge diagnostics the
+// Hermitian PowerMethod cannot:
+//   step 0 : |A v|/|v| on the (random) start src -- a one-sample lower bound
+//            on sigma_max(A).
+//   step k : |A 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 => a conjugate pair of equal
+//            modulus at the top).
+//   return : |lambda_max|, with a SUMMARY line flagging non-normality when the
+//            step-0 sigma_max lower bound sits well above the converged
+//            |lambda_max| -- the case where the numerical range extends beyond
+//            the spectrum and a spectrum-based smoother/Chebyshev bound is
+//            unsafe (the field of values must be used instead).
+template<class Field> class PowerMethodNonHermitian
+{
+ public:
+
+  template<typename T> static RealD normalise(T& v)
+  {
+    RealD nn = sqrt(norm2(v));
+    v = v * (1.0/nn);
+    return nn;
+  }
+
+  RealD operator()(LinearOperatorBase<Field> &Op, const Field &src)
+  {
+    GridBase *grid = src.Grid();
+    Field v(grid), Av(grid);
+    v = src;
+
+    RealD ratio = 0.0, ratio0 = 0.0;
+    ComplexD rq(0.0);
+    const int _MAX_ITER_EST_ = 200;
+
+    for (int i=0;i<_MAX_ITER_EST_;i++) {
+      normalise(v);                       // v is now unit
+      Op.Op(v,Av);
+      ratio = sqrt(norm2(Av));            // |A v| = |lambda_max| in the limit
+      rq    = innerProduct(v,Av);         // complex Rayleigh quotient
+
+      if ( i==0 ) {
+        ratio0 = ratio;
+        std::cout << GridLogMessage << "PowerMethodNonHermitian: step 0 (random v): |Av|/|v| = "
+                  << ratio << "   [lower bound on sigma_max]" << std::endl;
+      }
+      if ( (i%10==0) || (i==_MAX_ITER_EST_-1) )
+        std::cout << GridLogMessage << "PowerMethodNonHermitian: step " << i << " |Av|/|v| = " << ratio
+                  << "  Rayleigh (" << real(rq) << "," << imag(rq) << ")" << std::endl;
+
+      v = Av;
+    }
+
+    std::cout << GridLogMessage << "PowerMethodNonHermitian: |lambda_max| ~ " << ratio
+              << "  Rayleigh (" << real(rq) << "," << imag(rq) << ")"
+              << "  phase " << atan2(imag(rq),real(rq)) << " rad"
+              << "  step-0/converged = " << ratio0/ratio
+              << (ratio0/ratio > 1.2 ? "   ** non-normal: sigma_max well above |lambda_max| **"
+                                     : "   (near-normal)")
+              << std::endl;
+    return ratio;
+  }
+};
 }

Grid/algorithms/multigrid/BlockCyclicRedistribute.h

@@ -23,11 +23,47 @@ Author: Peter Boyle <pboyle@bnl.gov>
 
 NAMESPACE_BEGIN(Grid);
 
+///////////////////////////////////////////////////////////////////////////////
+// My rows of a distributed dense matrix in the 1D rank-major row layout:
+// rank r owns contiguous global rows [rowStart[r], rowStart[r+1]) of an
+// N x N matrix, stored rows x cols column major, ld = rows; element (i,j)
+// at data[i + j*ld].  This is the layout the stencil->dense import
+// produces and the apply slab consumes; the 2D inverse slots between them
+// through the redistribution below.
+///////////////////////////////////////////////////////////////////////////////
+class BlockRows
+{
+public:
+  deviceVector<ComplexD> data;
+  int64_t                rows;
+  int64_t                cols;
+  int64_t                ld;
+
+  BlockRows()
+  {
+    rows = 0;
+    cols = 0;
+    ld   = 0;
+  }
+  void Resize(int64_t r, int64_t c)
+  {
+    rows = r;
+    cols = c;
+    ld   = r;
+    data.resize((uint64_t)r*c);
+  }
+  ComplexD *ColumnWindow(int64_t col0)
+  {
+    GRID_ASSERT( col0 >= 0 );
+    GRID_ASSERT( col0 <= cols );
+    return &data[(uint64_t)col0*ld];
+  }
+};
+
 ///////////////////////////////////////////////////////////////////////////////
 // Stage 4 of the 2D distributed dense inverse: redistribution between the
-// 1D rank-major row layout (BlockRows: rank r owns contiguous global rows
-// [rowStart[r], rowStart[r+1]) of an N x N matrix, stored rows x N column
-// major with ld = rows) and the 2D block-cyclic layout.
+// 1D rank-major row layout (BlockRows above) and the 2D block-cyclic
+// layout.
 //
 // This is what lets the EXISTING stencil->dense import, its certificate,
 // the fp32 slab conversion and the apply path all remain byte-for-byte

Grid/algorithms/multigrid/BlockCyclicSchurInverse.h

@@ -28,12 +28,11 @@ NAMESPACE_BEGIN(Grid);
 // Stage 3 of the 2D distributed dense inverse: the recursive Schur
 // complement on a block-cyclic matrix, in place.
 //
-// The nine-step algebra is IDENTICAL to RecursiveSchurInverse (1D); what
-// changes is the decomposition.  The recursion splits the GLOBAL INDEX
-// RANGE at the block boundary nearest the midpoint -- not the rank range --
-// so every rank owns part of every sub-block at every depth, and the
-// ownership gating (inI/inJ, dummy operands, zero-width rank ranges) of the
-// 1D scheme has no analogue here: it is simply gone.
+// The recursion splits the GLOBAL INDEX RANGE at the block boundary
+// nearest the midpoint -- not the rank range as the retired 1D
+// RecursiveSchurInverse did -- so every rank owns part of every sub-block
+// at every depth, and the ownership gating (inI/inJ, dummy operands,
+// zero-width rank ranges) a 1D scheme needs has no analogue here.
 //
 //   I = [c0,m)   J = [m,c1)      (block-aligned, m the mid block boundary)
 //   1. recurse I  :  A11 -> A11inv                        (in place)
@@ -64,8 +63,7 @@ NAMESPACE_BEGIN(Grid);
 // LEAF.  A leaf is a single diagonal block, and block (b,b) of a
 // block-cyclic layout lives ENTIRELY on rank (b%Pr, b%Pc).  The leaf
 // inversion is therefore purely local -- pack the strided block dense,
-// GridBLASInverse, unpack -- with NO communication and no assembly.  The
-// leaf-assembly transport question of the 1D scheme does not arise.
+// GridBLASInverse, unpack -- with NO communication and no assembly.
 // Successive leaves cycle over ranks, so leaf work is naturally spread.
 //
 // COMMUNICATION.  Every transfer in the whole inversion is a
@@ -75,12 +73,9 @@ NAMESPACE_BEGIN(Grid);
 // one optional exception: it performs reductions, and is only ever called
 // explicitly by a caller who wants the numbers.
 //
-// NUMERICS.  No pivoting, exactly as the 1D scheme: every A11 and every
-// Schur complement met on the way down must be non-singular.  The growth
-// telemetry stands in for pivoting; note the recursion splits differently
-// from the 1D rank-range tree, so DIFFERENT sub-blocks are inverted and
-// telemetry values are NOT comparable with the 1D implementation's --
-// re-baseline, do not compare.
+// NUMERICS.  No pivoting: every A11 and every Schur complement met on the
+// way down must be non-singular.  The growth telemetry stands in for
+// pivoting.
 ///////////////////////////////////////////////////////////////////////////////
 
 class BlockCyclicSchurInverse
@@ -92,14 +87,16 @@ public:
   // Telemetry: accumulated LOCALLY, no comms unless ReportTelemetry().
   double                  telLeafMaxInv;
   uint64_t                nLeaf;
-  // BIG LEAVES (SCHUR2D_LEAF_SPAN=s, default 1 = off).  Below span s blocks a
-  // sub-block lives on <= s of the Pr process rows / s of the Pc columns; the
-  // SUMMA rings then run on a few ranks while the rest block in their next
-  // SendToRecvFrom (histogram 2026-08-27: 93% of ring time in the 3.7 MB
-  // single-block panels of exactly these levels).  Instead: gather the
-  // (s*nb)^2 sub-block to one rank, invert locally, scatter back.
-  int                     leafSpan = -1;
-  int                     leafLU   = -1;  // SCHUR2D_LEAF_LU=1: big-leaf inverse via GridBLASInverse::inverseLU (blocked getrf_64 + identity getrs_64) instead of getri_batched
+  // BIG LEAVES.  Below span s blocks a sub-block lives on <= s of the Pr
+  // process rows / s of the Pc columns; when s is small relative to the grid
+  // the SUMMA rings run on a few ranks while the rest block in their next
+  // SendToRecvFrom.  Instead: gather the (s*nb)^2 sub-block to one rank,
+  // invert locally (inverseLU), scatter back.  Fires only while
+  // span < min(Pr,Pc) -- when the whole grid participates the rings are not
+  // degenerate and gathering would only concentrate memory (and at the top
+  // of the tree would gather the whole matrix).  Larger s concentrates more
+  // memory on the root; smaller s loses ring parallelism.
+  int                     leafSpan = 9;
   uint64_t                nBigLeaf = 0;  int64_t maxBigW = 0;
   double                  tBigGather = 0, tBigInv = 0, tBigScatter = 0;
   uint64_t                nNode;
@@ -284,18 +281,9 @@ public:
     }
     tBigGather += usecond();
 
-    // ---- invert on root ----
+    // ---- invert on root: blocked getrf_64 + identity getrs_64 ----
     tBigInv -= usecond();
-    if ( leafLU < 0 ) leafLU = getenv("SCHUR2D_LEAF_LU") ? atoi(getenv("SCHUR2D_LEAF_LU")) : 0;
-    if ( me == root ) {
-      if ( leafLU ) {
-        INV.inverseLU(W, &dense[0]);
-      } else {
-        deviceVector<ComplexD*> bp(1); std::vector<ComplexD*> ptr(1); ptr[0] = &dense[0];
-        acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexD*));
-        INV.inverseBatched(W, bp);
-      }
-    }
+    if ( me == root ) INV.inverseLU(W, &dense[0]);
     tBigInv += usecond();
 
     // ---- scatter ----
@@ -330,8 +318,7 @@ public:
     int64_t span = b1-b0;
     GRID_ASSERT( span >= 1 );
     if ( span == 1 ) { Leaf(A, b0); return; }
-    if ( leafSpan < 0 ) leafSpan = getenv("SCHUR2D_LEAF_SPAN") ? atoi(getenv("SCHUR2D_LEAF_SPAN")) : 1;
-    if ( span <= leafSpan ) { BigLeaf(A, b0, b1); return; }
+    if ( span <= leafSpan && span < std::min(L.Pr,L.Pc) ) { BigLeaf(A, b0, b1); return; }
     GRID_TRACE("SchurNode");
     nNode++;
 
@@ -376,150 +363,11 @@ public:
   // PUBLIC ENTRY.  In-place inverse of the whole matrix.  Scratch (4x the
   // matrix footprint) is allocated here and released on return.
   ///////////////////////////////////////////////////////////////////////////
-  ///////////////////////////////////////////////////////////////////////////
-  // POWER - CLOCK - GROUND.  Before the inverse runs, print the preconditions
-  // that have differed between fast (62 s) and slow (133-141 s) runs of the
-  // SAME inverse, and measure ONE SendToRecvFrom to the ACTUAL ring partners
-  // at three sizes, device and host buffers.  Same code in the harness and in
-  // production, so the two processes are compared on the identical primitive
-  // before any explanation of the SUMMA rings is attempted.
-  // SCHUR2D_PROBE=0 disables (costs ~0.1-0.5 s).
-  ///////////////////////////////////////////////////////////////////////////
-  void Probe(BlockCyclicMatrix &A)
-  {
-    BlockCyclicLayout &L = A.layout;
-    GridBase *grid = A.grid;
-    int me = grid->ThisRank();
-    // --- banner ---
-    int thr = -1;
-#ifdef GRID_COMMS_MPI3
-    MPI_Query_thread(&thr);
-#endif
-    const char *omp = getenv("OMP_NUM_THREADS");
-    MemoryStatus ms = MemoryManager::GetFootprint();
-    std::cout << GridLogMessage << "Schur2D PROBE banner: MPI thread level " << thr
-              << " (0 single,1 funneled,2 serialized,3 multiple)  OMP_NUM_THREADS=" << (omp?omp:"unset")
-              << "  MemoryManager device bytes " << ms.DeviceBytes/1.0e9 << " GB (LRU " << ms.DeviceLRUBytes/1.0e9
-              << " GB, cap " << ms.DeviceMaxBytes/1.0e9 << " GB)"
-              << "  grid " << L.Pr << "x" << L.Pc << " nb " << L.nb << std::endl;
-#ifdef GRID_HIP
-    if ( me==0 ) acceleratorMem();
-#endif
-    // --- ring partners exactly as SUMMA uses them ---
-    int prow=L.prow, pcol=L.pcol, Pr=L.Pr, Pc=L.Pc;
-    struct Ring { const char *name; int dest, src; };
-    Ring rings[2] = { {"ringA(row, q+-1)", prow*Pc + (pcol+1)%Pc,       prow*Pc + (pcol-1+Pc)%Pc},
-                      {"ringB(col, p+-1)", ((prow+1)%Pr)*Pc + pcol,     ((prow-1+Pr)%Pr)*Pc + pcol} };
-    // 2/3/4 MB added 2026-08-27: the SUMMA histogram put 93% of ring time in
-    // [2,4) MB messages at 0.3 GB/s while >=4 MB ran at 11-13 GB/s.
-    const int NSZ = 6;
-    uint64_t sizes[NSZ] = { 64ull*1024, 1024ull*1024, 2048ull*1024, 3072ull*1024, 4096ull*1024, 8ull*1024*1024 };
-    uint64_t maxb = sizes[NSZ-1];
-    deviceVector<char> dsend(maxb), drecv(maxb);
-    std::vector<char>  hsend(maxb), hrecv(maxb);
-    for(int r=0;r<2;r++){
-      if ( (r==0 && Pc==1) || (r==1 && Pr==1) ) continue;
-      int off = grid->IsOffNode(rings[r].dest);
-      for(int si=0;si<NSZ;si++){
-        uint64_t bytes = sizes[si];
-        // warm one, time five, both memory spaces
-        grid->SendToRecvFrom(&dsend[0], rings[r].dest, &drecv[0], rings[r].src, bytes);
-        double t0=usecond();
-        for(int i=0;i<5;i++) grid->SendToRecvFrom(&dsend[0], rings[r].dest, &drecv[0], rings[r].src, bytes);
-        double td=(usecond()-t0)/5.0;
-        grid->SendToRecvFrom(&hsend[0], rings[r].dest, &hrecv[0], rings[r].src, bytes);
-        t0=usecond();
-        for(int i=0;i<5;i++) grid->SendToRecvFrom(&hsend[0], rings[r].dest, &hrecv[0], rings[r].src, bytes);
-        double th=(usecond()-t0)/5.0;
-        // spread over ranks
-        RealD dmax=td, dmin=-td; grid->GlobalMax(dmax); grid->GlobalMax(dmin); dmin=-dmin;
-        std::cout << GridLogMessage << "Schur2D PROBE " << rings[r].name << (off?" OFF-node":" on-node")
-                  << " " << bytes/1024 << " KB: device " << td << " us (" << bytes/td/1.0e3 << " GB/s) [min/max over ranks " << dmin << "/" << dmax << " us]"
-                  << "  host " << th << " us (" << bytes/th/1.0e3 << " GB/s)" << std::endl;
-      }
-    }
-    ///////////////////////////////////////////////////////////////////////
-    // ONE-GCD LEAF MICROBENCHMARK (boss only).  The big-leaf inverse at
-    // W=4320 measured 0.53 s with BOTH getri_batched and getrf_64+getrs_64
-    // (2026-08-27) -- ~0.4 TF/s on a GCD that runs zgemm at ~15.  Time the
-    // three primitives in isolation on a well-conditioned dense matrix so the
-    // leaf's rate can be compared with the machine's, and getrf split from
-    // getrs.  Sizes: the W of span 4 / 9 / 18 leaves on accelerator builds;
-    // tiny on CPU builds (Eigen would take minutes at 4320).
-    ///////////////////////////////////////////////////////////////////////
-    if ( me == 0 ) {
-#if defined(GRID_HIP) || defined(GRID_CUDA) || defined(GRID_SYCL)
-      std::vector<int64_t> Ws({1920, 4320, 8640});
-#else
-      std::vector<int64_t> Ws({240, 480});
-#endif
-      for(int64_t W : Ws){
-        deviceVector<ComplexD> M((uint64_t)W*W), C((uint64_t)W*W);
-        { ComplexD *m = &M[0]; const int64_t WW = W;   // diagonally dominant: (i==j ? W : 0) + cos/sin noise
-          accelerator_for(idx,(uint64_t)W*W,1,{ int64_t j=idx/WW, i=idx-j*WW; double x=0.37*i+0.61*j;
-            m[idx] = ComplexD((i==j)?(double)WW:0.0,0.0) + ComplexD(std::cos(x),std::sin(1.3*x)); });
-          accelerator_barrier(); }
-        double flopLU = 8.0/3.0*(double)W*W*W;          // complex LU  ~ (4 real flops per complex mult-add) * (2/3 n^3)
-        double flopGEMM = 8.0*(double)W*W*W;            // complex GEMM
-        // 1. getri_batched (batch 1)
-        double tb;
-        { deviceVector<ComplexD*> bp(1); std::vector<ComplexD*> ptr(1); ptr[0]=&M[0];
-          acceleratorCopyToDevice(&ptr[0],&bp[0],sizeof(ComplexD*));
-          double t0=usecond(); INV.inverseBatched(W,bp); tb=usecond()-t0; }
-        // 2. inverseLU (getrf_64 + identity getrs_64), timed inside
-        double tl; { double t0=usecond(); INV.inverseLU(W,&M[0]); tl=usecond()-t0; }
-        // 3. one zgemm W x W x W for the machine rate
-        double tg;
-        { deviceVector<ComplexD*> ap(1),bp(1),cp(1); std::vector<ComplexD*> ptr(1);
-          ptr[0]=&M[0]; acceleratorCopyToDevice(&ptr[0],&ap[0],sizeof(ComplexD*)); acceleratorCopyToDevice(&ptr[0],&bp[0],sizeof(ComplexD*));
-          ptr[0]=&C[0]; acceleratorCopyToDevice(&ptr[0],&cp[0],sizeof(ComplexD*));
-          double t0=usecond();
-          SUMMA.BLAS.gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N,(int)W,(int)W,(int)W,ComplexD(1.0,0.0),ap,(int)W,bp,(int)W,ComplexD(0.0,0.0),cp,(int)W);
-          SUMMA.BLAS.synchronise(); tg=usecond()-t0; }
-        std::cout << GridLogMessage << "Schur2D PROBE leaf W=" << W
-                  << ": getri_batched " << tb/1.0e6 << " s (" << flopLU/tb/1.0e6 << " TF/s)"
-                  << "  getrf_64+getrs_64 " << tl/1.0e6 << " s (getrf " << INV.lastGetrfUs/1.0e6 << " getrs " << INV.lastGetrsUs/1.0e6 << ")"
-                  << "  zgemm " << tg/1.0e6 << " s (" << flopGEMM/tg/1.0e6 << " TF/s)" << std::endl;
-      }
-    }
-    ///////////////////////////////////////////////////////////////////////
-    // The SUMMA's conditions, one at a time, at 8 MB on ring B:
-    //  (a) LARGE persistent buffers (the rings use ~0.5 GB Abuf/Bbuf), sending
-    //      from offset 0 and from deep inside the region;
-    //  (b) a pack kernel + accelerator_barrier immediately before each
-    //      message, as the SUMMA does.
-    // Isolated 8 MB messages ran at 11-20 GB/s while the SUMMA averaged 2.1;
-    // whichever variant drops to ~2 GB/s names the condition.
-    ///////////////////////////////////////////////////////////////////////
-    {
-      int r = (Pr>1) ? 1 : 0;
-      uint64_t bytes = sizes[NSZ-1];
-      uint64_t big = 512ull*1024*1024;
-      deviceVector<char> bsend(big), brecv(big);
-      for(int variant=0; variant<3; variant++){
-        uint64_t so = (variant==1) ? big-bytes : 0;               // deep offset in the large region
-        char *sp=&bsend[so], *rp=&brecv[so];
-        grid->SendToRecvFrom(sp, rings[r].dest, rp, rings[r].src, bytes);
-        double t0=usecond();
-        for(int i=0;i<5;i++){
-          if ( variant==2 ) { accelerator_for(k, bytes/8, 1, { ((uint64_t *)sp)[k] = (uint64_t)k; }); accelerator_barrier(); }
-          grid->SendToRecvFrom(sp, rings[r].dest, rp, rings[r].src, bytes);
-        }
-        double t=(usecond()-t0)/5.0;
-        RealD tmax=t, tmin=-t; grid->GlobalMax(tmax); grid->GlobalMax(tmin); tmin=-tmin;
-        const char *vn[3]={"512MB buffer, offset 0","512MB buffer, offset 504MB","pack kernel + barrier before each send"};
-        std::cout << GridLogMessage << "Schur2D PROBE " << rings[r].name << " 8192 KB device, " << vn[variant] << ": "
-                  << t << " us (" << bytes/t/1.0e3 << " GB/s) [min/max over ranks " << tmin << "/" << tmax << " us]" << std::endl;
-      }
-    }
-  }
-
   void Invert(BlockCyclicMatrix &A)
   {
     BlockCyclicLayout &L = A.layout;
     GRID_ASSERT( L.N >= 1 );
     int64_t nblocks = (L.N + L.nb - 1)/L.nb;
-    if ( !(getenv("SCHUR2D_PROBE") && atoi(getenv("SCHUR2D_PROBE"))==0) ) Probe(A);
 
     BlockCyclicMatrix Bt(A.grid, L.N, L.nb, L.Pr, L.Pc);
     BlockCyclicMatrix Ct(A.grid, L.N, L.nb, L.Pr, L.Pc);
@@ -550,8 +398,8 @@ public:
               << std::endl;
     if ( nBigLeaf ) {
       RealD ti = tBigInv/1.0e6; grid->GlobalMax(ti);     // inverse runs on the root of each leaf: report the max over ranks
-      std::cout << GridLogMessage << "BlockCyclicSchurInverse: BIG LEAVES (SCHUR2D_LEAF_SPAN=" << leafSpan
-                << (leafLU>0 ? ", SCHUR2D_LEAF_LU=1: getrf_64+getrs_64" : ", getri_batched") << "): " << nBigLeaf
+      std::cout << GridLogMessage << "BlockCyclicSchurInverse: BIG LEAVES (span " << leafSpan
+                << ", inverseLU): " << nBigLeaf
                 << " leaves, max W " << maxBigW
                 << "  boss secs: gather " << tBigGather/1.0e6 << " scatter " << tBigScatter/1.0e6
                 << "  inverse (max over ranks) " << ti << std::endl;
@@ -576,8 +424,7 @@ public:
     // same sequence of sizes).  Time is wall time inside SendToRecvFrom, so it
     // includes waiting for the partner -- a bucket whose GB/s is far below the
     // probe's for the same size is wait, not wire.
-    std::cout << GridLogMessage << "BlockCyclicSumma ring histogram (boss):  size-bucket  msgs  GB  xfer-secs  GB/s  %time"
-              << (SUMMA.handshake>0 ? "  handshake-secs (partner wait, excluded from xfer)" : "") << std::endl;
+    std::cout << GridLogMessage << "BlockCyclicSumma ring histogram (boss):  size-bucket  msgs  GB  xfer-secs  GB/s  %time" << std::endl;
     std::streamsize oldprec = std::cout.precision();
     for(int b=0;b<SUMMA.NHIST;b++){
       if ( !SUMMA.histN[b] ) continue;
@@ -590,7 +437,6 @@ public:
                 << std::setw(9) << std::setprecision(3) << sec
                 << std::setw(9) << std::setprecision(3) << (sec>0 ? g/sec : 0.0)
                 << std::setw(8) << std::setprecision(3) << (ring>0 ? 100.0*sec/ring : 0.0);
-      if ( SUMMA.handshake>0 ) std::cout << std::setw(12) << std::setprecision(3) << SUMMA.histHsUs[b]/1.0e6;
       std::cout << std::endl;
     }
     std::cout.precision(oldprec);

Grid/algorithms/multigrid/BlockCyclicSumma.h

@@ -146,16 +146,12 @@ public:
   // Per-message-size histogram (bucket = floor(log2 bytes)): count, bytes,
   // microseconds -- decomposes the ring time by packet size so a low average
   // GB/s can be attributed (many small latency-bound messages vs slow large
-  // ones vs partner-wait).  The 8 MB probe runs at 11-20 GB/s; SUMMA averaged 2.
+  // ones vs partner-wait).  Time is wall time inside SendToRecvFrom, so a
+  // bucket far below the wire rate for its size is wait, not wire.
   static const int NHIST=48;
-  uint64_t histN[NHIST]={0}, histBytes[NHIST]={0}; double histUs[NHIST]={0}, histHsUs[NHIST]={0};
-  // SUMMA_HANDSHAKE=1: a 4-byte SendToRecvFrom with the same partner
-  // immediately before each ring message, timed separately (histHsUs).
-  // Handshake time = partner-arrival skew; the remainder = transfer.  Splits
-  // the [2,4) MB bucket's 12 ms/msg (2026-08-27) into wait vs wire.
-  int handshake = -1; int hsTx=0, hsRx=0;
-  void HistAdd(uint64_t bytes, double us, double hs=0.0){ int b=0; while((bytes>>b)>1) b++; histN[b]++; histBytes[b]+=bytes; histUs[b]+=us; histHsUs[b]+=hs; }
-  void ResetTelemetry(void){ tAlloc=tPack=tRingA=tRingB=tGemm=0; bytesRing=nRingMsg=nMultiply=nGemm=0; for(int b=0;b<NHIST;b++){histN[b]=histBytes[b]=0; histUs[b]=histHsUs[b]=0;} }
+  uint64_t histN[NHIST]={0}, histBytes[NHIST]={0}; double histUs[NHIST]={0};
+  void HistAdd(uint64_t bytes, double us){ int b=0; while((bytes>>b)>1) b++; histN[b]++; histBytes[b]+=bytes; histUs[b]+=us; }
+  void ResetTelemetry(void){ tAlloc=tPack=tRingA=tRingB=tGemm=0; bytesRing=nRingMsg=nMultiply=nGemm=0; for(int b=0;b<NHIST;b++){histN[b]=histBytes[b]=0; histUs[b]=0;} }
 
   static int Overlap(int64_t a0,int64_t a1,int64_t b0,int64_t b1)
   { return (a0 < b1) && (b0 < a1); }
@@ -219,7 +215,6 @@ public:
     const uint64_t slotB1 = (uint64_t)nb*nloc_j;         // one panel
     const uint64_t slotB  = (uint64_t)S*slotB1;
     nMultiply++;
-    if ( handshake < 0 ) handshake = getenv("SUMMA_HANDSHAKE") ? atoi(getenv("SUMMA_HANDSHAKE")) : 0;
     tAlloc -= usecond();
     if ( Abuf.size() < std::max<uint64_t>(slotA*Pc,1) ) Abuf.resize( std::max<uint64_t>(slotA*Pc,1) );
     if ( Bbuf.size() < std::max<uint64_t>(slotB*Pr,1) ) Bbuf.resize( std::max<uint64_t>(slotB*Pr,1) );
@@ -286,12 +281,11 @@ public:
         for(int t=1;t<Pc;t++){
           int cs = (pcol - t + 1 + Pc*Pc) % Pc;
           int cr = (pcol - t     + Pc*Pc) % Pc;
-          double ths = 0.0, tm = usecond();
-          if ( handshake ) { grid->SendToRecvFrom((void *)&hsTx, dest, (void *)&hsRx, src, sizeof(int)); ths = usecond()-tm; tm = usecond(); }
+          double tm = usecond();
           grid->SendToRecvFrom((void *)(&Abuf[0]+slotA*cs), dest,
                                (void *)(&Abuf[0]+slotA*cr), src,
                                slotA*sizeof(ComplexD));
-          HistAdd(slotA*sizeof(ComplexD), usecond()-tm, ths);
+          HistAdd(slotA*sizeof(ComplexD), usecond()-tm);
           bytesRing += slotA*sizeof(ComplexD); nRingMsg++;
         }
         tRingA += usecond();
@@ -307,12 +301,11 @@ public:
         for(int t=1;t<Pr;t++){
           int rs = (prow - t + 1 + Pr*Pr) % Pr;
           int rr = (prow - t     + Pr*Pr) % Pr;
-          double ths = 0.0, tm = usecond();
-          if ( handshake ) { grid->SendToRecvFrom((void *)&hsTx, dest, (void *)&hsRx, src, sizeof(int)); ths = usecond()-tm; tm = usecond(); }
+          double tm = usecond();
           grid->SendToRecvFrom((void *)(&Bbuf[0]+slotB*rs), dest,
                                (void *)(&Bbuf[0]+slotB*rr), src,
                                slotB*sizeof(ComplexD));
-          HistAdd(slotB*sizeof(ComplexD), usecond()-tm, ths);
+          HistAdd(slotB*sizeof(ComplexD), usecond()-tm);
           bytesRing += slotB*sizeof(ComplexD); nRingMsg++;
         }
         tRingB += usecond();

Grid/algorithms/multigrid/DenseCoarseMatrix.h

@@ -29,7 +29,6 @@ Author: Peter Boyle <pboyle@bnl.gov>
 
 #include <Grid/algorithms/blas/BatchedBlas.h>
 #include <Grid/algorithms/blas/BatchedInverse.h>
-#include <Grid/algorithms/multigrid/RecursiveSchurInverse.h>
 #include <Grid/algorithms/multigrid/BlockCyclicSchurInverse.h>
 #include <Grid/algorithms/multigrid/BlockCyclicRedistribute.h>
 
@@ -39,46 +38,31 @@ NAMESPACE_BEGIN(Grid);
 
 //////////////////////////////////////////////////////////////////////////////////////
 // DenseCoarseMatrix: a coarsened operator treated as a DENSE matrix -- explicit,
-// row-distributed A^{-1} of a GeneralCoarsenedMatrix.  Library-grade successor of
-// the example-local DistributedDenseInverse (Example_pvdagm_mrhs_3level_dense.cc,
-// FROZEN as the regression baseline).
+// row-distributed A^{-1} of a GeneralCoarsenedMatrix.
 //
-// What is new versus the example class:
 //  - Stencil -> dense DIRECT IMPORT.  The coarse operator IS the dense matrix
 //    unrolled: Dense[(s,a),(s+shift_p,b)] += A[p][s]_{a,b}.  Rows of my sites are
-//    assembled from purely LOCAL _A[p] data: no operator applies, no comms -- the
-//    O(N) probe assembly (93 s at N=69120) is retired.  ACCUMULATE (+=) because on
-//    short axes distinct shifts wrap to the same neighbour.  An IMPORT CERTIFICATE
-//    compares the dense apply against Op.M on a NON-CONSTANT vector (a constant one
-//    cannot see a shift-sign error); DENSE_IMPORT_SIGN=-1 flips the convention
-//    without recompiling.
+//    assembled from purely LOCAL _A[p] data: no operator applies, no comms.
+//    ACCUMULATE (+=) because on short axes distinct shifts wrap to the same
+//    neighbour.  An IMPORT CERTIFICATE compares the dense apply against Op.M on a
+//    NON-CONSTANT vector (a constant one cannot see a shift-sign error).
+//
+//  - Inversion is END-TO-END fp64 through the 2D block-cyclic recursive Schur
+//    complement (BlockCyclicSchurInverse): fp64 rank-major import ->
+//    RowsToCyclic -> in-place recursion (pure point-to-point SUMMA rings and
+//    local leaves; bitwise reproducible) -> CyclicToRows -> ONE terminal
+//    rounding into the fp32 apply slab.  Distributed at every N and P.
+//
 //  - Split-K apply through GridBLAS.gemmBatched with EXPLICIT leading dimensions
 //    (arXiv:2409.03904 fig 11): the tiny-output/huge-K GEMM Y = slab^T X becomes
-//    DENSE_SPLITK chunk-GEMMs by pointer offset into the resident slab (lda = N),
-//    partials reduced in one accelerator_for.  Platform-agnostic: deviceVector +
-//    GridBLAS run the SAME code on HIP/CUDA/SYCL and CPU(Eigen).
-//  - deviceVector everywhere in the apply path; the ONE surviving naked-HIP block
-//    is the boss inversion buffer (quarantined below, documented).
+//    SPLITK chunk-GEMMs by pointer offset into the resident slab (lda = N),
+//    partials reduced in one accelerator_for.  The source vector is assembled by
+//    a cartesian ring ALLGATHER of device buffers (pure P2P; ~8x fewer bytes
+//    than a padded allreduce, and no collective size cliffs).  Platform-agnostic:
+//    deviceVector + GridBLAS run the SAME code on HIP/CUDA/SYCL and CPU(Eigen).
 //
-// Setup: SLAB_FILE=<stem> loads per-rank <stem>.<rank> (header-guarded N/nrows/
-// nbasis -- the interchange format shared with the frozen example; the STEM must
-// encode cfg/mass/blocking/nbasis, only the header is guarded).  Absent: direct
-// import -> import certificate -> chunked zero-fill+GlobalSum gather streamed to
-// the boss GCD -> cgetrf_64 (ILP64) -> rows of A^{-1} via blocked identity
-// cgetrs_64 + broadcast, each rank keeping the rows of its own sites -> save.
-// VERIFY ||A Ainv x - x||/||x|| runs in BOTH paths (and now certifies the DEVICE
-// slab + split-K path, since the single-RHS apply routes through the same core).
-//
-// Env: SLAB_FILE  DENSE_SPLITK (default 32, snapped to a divisor of N)
-//      DENSE_DEVICE_SUM  DENSE_IMPORT_SIGN  DENSE_APPLY_PROFILE  DENSE_CC_CHECK
-//      DENSE_SCHUR (0/absent: single-GCD gather-invert; 1: distributed
-//      recursive Schur; 2: AUDIT -- run BOTH on the same imported A, report
-//      the slab difference, keep the Schur result)  DENSE_PANEL_BYTES
-//
-// The DENSE_SCHUR=1 path is the RecursiveSchurInverse distributed
-// factorisation: it lifts the fp32 N ~ 90k boss-HBM ceiling (the CC-grid
-// 256-rank SIMD cap remains -- separate issue).  Internal only: slab layout,
-// apply path, SLAB_FILE format and VERIFY are identical in every mode.
+// VERIFY ||A Ainv x - x||/||x|| certifies the DEVICE slab + split-K path at the
+// end of Import, since the single-RHS apply routes through the same core.
 //
 // Tensor-depth agnostic: site scalar objects treated as contiguous ComplexD
 // (iScalar wrappers add no data), so any MG level's coarse operator imports.
@@ -111,15 +95,15 @@ public:
   std::vector<int64_t>    myGsite;  // global lex site index of my site ss
   std::vector<ComplexF>   slab;     // nrows x N row-major: A during setup, rows of A^{-1} after
 
-  static const int64_t CHUNKROWS = 1024;  // getrs harvest block (trsm efficiency + fewer broadcasts)
   static const int MRHS_MAX = 32;
+  static const int SPLITK   = 32;   // requested split-K chunk count, snapped DOWN to a divisor of N
 
   // Apply machinery: resident slab + persistent buffers + AOT split-K pointers.
   GridBLAS BLAS;
   deviceVector<ComplexF>  dSlab;
   deviceVector<ComplexF>  dX;       // N x MRHS_MAX
   deviceVector<ComplexF>  dY;       // nrows x MRHS_MAX
-  deviceVector<ComplexF>  dG;       // N x MRHS_MAX lex-major staging for the allgather (devSum==4)
+  deviceVector<ComplexF>  dG;       // N x MRHS_MAX lex-major staging for the allgather
   deviceVector<int>       dLex2Rank;// lex index of a process coordinate -> its rank (allgather block order -> row-block order)
   deviceVector<int64_t>   dRm2G;    // rank-major index (rank*nrows + ss*nbasis + b) -> global column (gsite*nbasis + b) of x / the slab
   int                     myLex;
@@ -130,8 +114,6 @@ public:
   std::vector<ComplexF>   hX;
   std::vector<ComplexF>   hY;
   int NK;                           // split-K chunk count (divides N)
-  int devSum;
-  double schurAuditRel;             // DENSE_SCHUR=2: rel slab diff single-vs-schur (-1 = not run)
 
   DenseCoarseMatrix(GridBase *g)
     : grid(g)
@@ -142,7 +124,6 @@ public:
     N      = grid->gSites() * nbasis;
     lsites = grid->lSites();
     nrows  = (int64_t)lsites * nbasis;
-    schurAuditRel = -1.0;
 
     std::cout << GridLogMessage << "DenseCoarseMatrix: N = " << N
               << " (" << grid->gSites() << " sites x " << nbasis << ")"
@@ -180,47 +161,11 @@ public:
   void Import(CoarseOp &Op)
   {
     double t0 = usecond();
-    ////////////////////////////////////////////////////////////////////
-    // 0. Slab cache: SLAB_FILE=<stem> -> per-rank raw file <stem>.<rank>.
-    //    SAME format as the frozen example (interchange compatible).
-    ////////////////////////////////////////////////////////////////////
-    bool loaded = false;
-    char *sfile = getenv("SLAB_FILE");
-    std::string slabfile;
-    if (sfile) {
-      slabfile = std::string(sfile) + "." + std::to_string(grid->ThisRank());
-      FILE *f = fopen(slabfile.c_str(),"rb");
-      if (f) {
-        int64_t hdr[4] = {0,0,0,0};
-        GRID_ASSERT( fread(hdr,sizeof(int64_t),4,f) == 4 );
-        GRID_ASSERT( hdr[0] == (int64_t)0x44454E5345 );  // magic "DENSE"
-        GRID_ASSERT( hdr[1] == N && hdr[2] == (int64_t)nrows && hdr[3] == (int64_t)nbasis );
-        uint64_t nelem = (uint64_t)nrows * N;
-        GRID_ASSERT( fread(&slab[0], sizeof(ComplexF), nelem, f) == nelem );
-        fclose(f);
-        loaded = true;
-        std::cout << GridLogMessage << "DenseCoarseMatrix: slab loaded from "
-                  << slabfile << " -- skipping import/factor/solve" << std::endl;
-      } else {
-        std::cout << GridLogMessage << "DenseCoarseMatrix: slab cache " << slabfile
-                  << " absent -- full setup, will write it" << std::endl;
-      }
-    }
-    if (!loaded) {
+    {
       ImportDense(Op);        // slab <- my rows of A   (LOCAL, no comms)
       ImportCertificate(Op);  // dense apply == Op.M, before inversion
       InvertDense(Op);        // slab <- my rows of A^{-1}
       double t1 = usecond();
-      if (sfile) {
-        FILE *f = fopen(slabfile.c_str(),"wb");
-        GRID_ASSERT(f != nullptr);
-        int64_t hdr[4] = { (int64_t)0x44454E5345, N, (int64_t)nrows, (int64_t)nbasis };
-        GRID_ASSERT( fwrite(hdr,sizeof(int64_t),4,f) == 4 );
-        uint64_t nelem = (uint64_t)nrows * N;
-        GRID_ASSERT( fwrite(&slab[0], sizeof(ComplexF), nelem, f) == nelem );
-        fclose(f);
-        std::cout << GridLogMessage << "DenseCoarseMatrix: slab written to " << slabfile << std::endl;
-      }
       std::cout << GridLogMessage << "DenseCoarseMatrix: import+invert took "
                 << (t1-t0)/1.0e6 << " s" << std::endl;
     }
@@ -233,11 +178,9 @@ public:
       dSlab.resize((uint64_t)nrows*N);
       acceleratorCopyToDevice(&slab[0],&dSlab[0],sbytes);
 
-      // DENSE_SPLITK: requested chunk count, snapped DOWN to a divisor of N.
-      int req = getenv("DENSE_SPLITK") ? atoi(getenv("DENSE_SPLITK")) : 32;
-      if (req < 1) req = 1;
+      // Split-K chunk count, snapped DOWN to a divisor of N.
       NK = 1;
-      for(int j=1;j<=req;j++) if ( (N % j) == 0 ) NK = j;
+      for(int j=1;j<=SPLITK;j++) if ( (N % j) == 0 ) NK = j;
       int64_t Kc = N / NK;
 
       dX.resize((uint64_t)N*MRHS_MAX);
@@ -255,12 +198,7 @@ public:
       for(int j=0;j<NK;j++) h[j] = &dPartial[0] + (uint64_t)j*nrows*MRHS_MAX;  // compact, ldc=nrows
       acceleratorCopyToDevice(&h[0],&cptrs[0],NK*sizeof(ComplexF*));
 
-      devSum = getenv("DENSE_DEVICE_SUM") ? atoi(getenv("DENSE_DEVICE_SUM")) : 0;
-      const char *sumName[5] = {"host allreduce","DEVICE-buffer allreduce (GPU-aware MPI)",
-                                "DEVICE cartesian ring allreduce (P2P)","DEVICE flat ring allreduce (P2P)",
-                                "DEVICE cartesian ring ALLGATHER (P2P, ~8x fewer bytes than the padded allreduce)"};
-      GRID_ASSERT(devSum>=0 && devSum<=4);
-      if ( devSum==4 ) {
+      {
         dG.resize((uint64_t)N*MRHS_MAX);
         // allgather delivers blocks in lexicographic-coordinate order; the row
         // blocks of x are in RANK order.  Same table as BuildRankMajorMap.
@@ -273,8 +211,7 @@ public:
         GRID_ASSERT( l2r[myLex] == grid->ThisRank() );
         // x and the slab columns are in GLOBAL-SITE order (hX[myGsite*nbasis+b]);
         // the gathered blocks are in RANK-MAJOR order (rank*nrows + ss*nbasis + b).
-        // On one rank the two coincide, which is how the laptop passed while
-        // Frontier VERIFYed 0.9965 (2026-08-27).  Scatter through the inverse map.
+        // The two coincide only on one rank, so scatter through the inverse map.
         std::vector<int64_t> g2rm; BuildRankMajorMap(g2rm);
         std::vector<int64_t> rm2g(N); for(int64_t g=0; g<N; g++) rm2g[g2rm[g]] = g;
         for(int ss=0; ss<lsites; ss++) GRID_ASSERT( rm2g[(int64_t)grid->ThisRank()*nrows + (int64_t)ss*nbasis] == myGsite[ss]*nbasis );
@@ -282,8 +219,8 @@ public:
         acceleratorCopyToDevice(&rm2g[0], &dRm2G[0], N*sizeof(int64_t));
       }
       std::cout << GridLogMessage << "DenseCoarseMatrix: slab resident on device ("
-                << sbytes/1024./1024. << " MB/rank), split-K NK=" << NK << " (Kc=" << Kc << "); "
-                << sumName[devSum] << std::endl;
+                << sbytes/1024./1024. << " MB/rank), split-K NK=" << NK << " (Kc=" << Kc
+                << "); DEVICE cartesian ring ALLGATHER (P2P)" << std::endl;
     }
 
     ////////////////////////////////////////////////////////////////////
@@ -355,8 +292,6 @@ public:
   {
     double t = -usecond();
     Coordinate gdims = grid->GlobalDimensions();
-    int sign = getenv("DENSE_IMPORT_SIGN") ? atoi(getenv("DENSE_IMPORT_SIGN")) : 1;
-    GRID_ASSERT( sign==1 || sign==-1 );
 
     uint64_t nelem = (uint64_t)nrows * N;
     thread_for(i, nelem, { slab[i] = ComplexF(0.0,0.0); });
@@ -367,28 +302,11 @@ public:
       // extract the unpadded field before peeking with unpadded coordinates
       // (exactly as MultiGeneralCoarsenedMatrix::CopyMatrix does).
       CoarseMatrix Aun(grid);  Op.ExtractMatrix(p,Aun);
-      if ( getenv("DENSE_IMPORT_DEBUG") ) {
-        // Peek-path vs field-norm audit: sum |peekLocalSite|^2 must match norm2
-        double pk = 0.0;
-        autoView(Adbg, Aun, CpuRead);
-        for(int ss=0; ss<lsites; ss++){
-          Msobj m;
-          peekLocalSite(m, Adbg, myLcoor[ss]);
-          ComplexD *md = (ComplexD *)&m;
-          for(int i=0; i<nbasis*nbasis; i++) pk += md[i].real()*md[i].real() + md[i].imag()*md[i].imag();
-        }
-        RealD gpk = pk;
-        grid->GlobalSumVector(&gpk, 1);
-        std::cout << GridLogMessage << "DenseCoarseMatrix: DEBUG p=" << p
-                  << " norm2(_A[p]) " << norm2(Aun)
-                  << " norm2(Extract) " << norm2(Aun)
-                  << " sum|peek|^2 " << gpk << std::endl;
-      }
       autoView(Av, Aun, CpuRead);
       thread_for(ss, lsites, {
         Coordinate ncoor(nd);
         for(int d=0; d<nd; d++){
-          int64_t g = grid->_lstart[d] + myLcoor[ss][d] + sign*shift[d];
+          int64_t g = grid->_lstart[d] + myLcoor[ss][d] + shift[d];
           ncoor[d] = (int)((g % gdims[d] + gdims[d]) % gdims[d]);
         }
         int64_t nsite;
@@ -399,11 +317,7 @@ public:
         // The operator contracts out(s,b) = sum_a A[p](s)(a,b) in(nbr,a)
         // (GeneralCoarsenedMatrix.h Mult kernel): the stored site matrix
         // acts TRANSPOSED, so element (a,b) lands at dense row (s,b),
-        // column (nbr,a).  BUG LEDGER 2026-08-14: the original mapping
-        // wrote (s,a),(nbr,b) -- caught by the IMPORT CERTIFICATE on its
-        // FIRST fresh-import exercise (Test_schur_dense_coarse); every
-        // production slab predates this path (probe-import SLAB_FILEs),
-        // so no production output is suspect.
+        // column (nbr,a).
         for(int b=0; b<nbasis; b++){
           ComplexF *row = &slab[(uint64_t)(ss*nbasis+b)*N + nsite*nbasis];
           for(int a=0; a<nbasis; a++)
@@ -499,250 +413,15 @@ public:
               << rel << std::endl;
     if ( rel >= 1.0e-3 ) {
       std::cout << GridLogMessage << "DenseCoarseMatrix: IMPORT CERTIFICATE FAILED. If O(1), the "
-                << "stencil shift-sign convention is opposite: rerun with DENSE_IMPORT_SIGN=-1"
+                << "stencil shift-sign convention of the coarse operator has changed: "
+                << "the import in ImportDense/ImportDenseFP64 must change with it"
                 << std::endl;
     }
     GRID_ASSERT(rel < 1.0e-3);
   }
 
   ////////////////////////////////////////////////////////////////////
-  // 3. Invert dispatcher.  slab holds my rows of A on entry, my rows
-  //    of A^{-1} on exit, in every mode.
-  //      DENSE_SCHUR absent/0 : single-GCD gather-invert (the oracle)
-  //      DENSE_SCHUR=1        : distributed recursive Schur
-  //      DENSE_SCHUR=2        : AUDIT -- both on the same A; report the
-  //                             slab difference; keep the Schur result
-  //                             (so VERIFY certifies the new path).
-  ////////////////////////////////////////////////////////////////////
-  template<class CoarseOp>
-  void InvertDense(CoarseOp &Op)
-  {
-    char *sc  = getenv("DENSE_SCHUR");
-    int  mode = sc ? atoi(sc) : 0;
-
-    if ( mode == 0 )
-    {
-      InvertDenseSingle();
-      return;
-    }
-    if ( mode == 1 )
-    {
-      InvertDenseSchur(Op);
-      return;
-    }
-    GRID_ASSERT( mode == 2 );
-    std::vector<ComplexF> Aimp(slab);       // imported A
-    InvertDenseSingle();
-    std::vector<ComplexF> ref(slab);        // Ainv, single path
-    slab = Aimp;
-    InvertDenseSchur(Op);                   // slab = Ainv, Schur path
-
-    // NaN-PROOF comparison: max() masks NaN, so count non-finite
-    // entries in each result explicitly.
-    double  mx = 0.0;
-    double  mr = 0.0;
-    int64_t badschur  = 0;
-    int64_t badsingle = 0;
-    for(uint64_t i=0; i<(uint64_t)nrows*N; i++)
-    {
-      double as = abs(ComplexD(slab[i]));
-      double ar = abs(ComplexD(ref[i]));
-      if ( !std::isfinite(as) ) badschur++;
-      if ( !std::isfinite(ar) ) badsingle++;
-      if ( std::isfinite(as) && std::isfinite(ar) )
-      {
-        mx = std::max(mx, (double)abs(ComplexD(slab[i]) - ComplexD(ref[i])));
-        mr = std::max(mr, ar);
-      }
-    }
-    RealD gmx = mx;
-    RealD gmr = mr;
-    RealD gbs = (RealD)badschur;
-    RealD gbr = (RealD)badsingle;
-    grid->GlobalMax(gmx);
-    grid->GlobalMax(gmr);
-    grid->GlobalSumVector(&gbs, 1);
-    grid->GlobalSumVector(&gbr, 1);
-    schurAuditRel = gmx/gmr;
-    std::cout << GridLogMessage << "DenseCoarseMatrix: DENSE_SCHUR=2 AUDIT "
-              << "max|Ainv_schur - Ainv_single| = " << gmx
-              << "  relative " << schurAuditRel
-              << "  non-finite: schur " << (int64_t)gbs << " single " << (int64_t)gbr
-              << "  (two fp32 roundings of the same inverse; expect ~ growth * eps32)"
-              << std::endl;
-    GRID_ASSERT( gbs == 0 );
-    GRID_ASSERT( gbr == 0 );
-  }
-
-  ////////////////////////////////////////////////////////////////////
-  // 3a. Single-GCD invert: chunked zero-fill+GlobalSum gather of A
-  //    streamed to the boss GCD, cgetrf_64 (ILP64), rows of A^{-1} by
-  //    blocked identity cgetrs_64 + broadcast; each rank keeps its own
-  //    rows (in `slab`, overwriting A).  Proven path; the SCHUR oracle.
-  ////////////////////////////////////////////////////////////////////
-  void InvertDenseSingle(void)
-  {
-    double t1 = usecond();
-    int boss = grid->IsBoss();
-    std::vector<ComplexF> Afull;
-#ifdef GRID_HIP
-    // QUARANTINED naked HIP: the boss-only N^2 inversion buffer (34GB at
-    // N=65536) must come from raw HBM; EvictAll flushes the device-copy
-    // layer to make the window.  (FreePool of the allocator free-list
-    // awaits the type-dispatched fix.)  Confined to setup; the apply path
-    // is pure Grid primitives.
-    rocblas_float_complex *dA = nullptr;
-    rocblas_float_complex *dB = nullptr;
-    int64_t *dIpiv = nullptr;
-    uint64_t Abytes = (uint64_t)N * N * sizeof(ComplexF);
-    MemoryManager::EvictAll();
-    if (boss) {
-      auto aerr = hipMalloc((void **)&dA, Abytes);
-      if (aerr != hipSuccess) {
-        std::cout << GridLogMessage << "DenseCoarseMatrix: hipMalloc of "
-                  << Abytes/1024./1024./1024. << " GB FAILED -- reduce --device-mem" << std::endl;
-        GRID_ASSERT(aerr == hipSuccess);
-      }
-      std::cout << GridLogMessage << "DenseCoarseMatrix: device inversion buffer allocated ("
-                << Abytes/1024./1024./1024. << " GB)" << std::endl;
-    }
-#else
-    if (boss) Afull.resize((uint64_t)N * N);
-#endif
-    {
-      std::unordered_map<int64_t,int> rowmap;   // global row -> my slab row
-      for(int ss=0; ss<lsites; ss++)
-        for(int a=0; a<nbasis; a++)
-          rowmap[ myGsite[ss]*nbasis + a ] = ss*nbasis + a;
-
-      std::vector<ComplexF> chunk((uint64_t)CHUNKROWS * N);
-      for(int64_t row0=0; row0<N; row0+=CHUNKROWS){
-        int64_t nrow = std::min(CHUNKROWS, N-row0);
-        uint64_t nelem = (uint64_t)nrow * N;
-        for(uint64_t i=0;i<nelem;i++) chunk[i]=ComplexF(0.0,0.0);
-        for(int64_t r=row0; r<row0+nrow; r++){
-          auto it = rowmap.find(r);
-          if (it != rowmap.end()) {
-            uint64_t src = (uint64_t)(it->second) * N;
-            uint64_t dst = (uint64_t)(r-row0) * N;
-            for(int64_t j=0;j<N;j++) chunk[dst+j] = slab[src+j];
-          }
-        }
-        grid->GlobalSumVector(&chunk[0], (int)nelem);
-        if (boss) {
-#ifdef GRID_HIP
-          GRID_ASSERT( hipMemcpy((char *)dA + (uint64_t)row0*N*sizeof(ComplexF),
-                                 &chunk[0], nelem*sizeof(ComplexF),
-                                 hipMemcpyHostToDevice) == hipSuccess );
-#else
-          uint64_t dst = (uint64_t)row0 * N;
-          for(uint64_t i=0;i<nelem;i++) Afull[dst+i] = chunk[i];
-#endif
-        }
-      }
-    }
-    double t2 = usecond();
-    std::cout << GridLogMessage << "DenseCoarseMatrix: gather to boss took "
-              << (t2-t1)/1.0e6 << " s" << std::endl;
-
-    ////////////////////////////////////////////////////////////////////
-    // Factor in place on the boss (fp32, ILP64).  Row-major buffer handed
-    // to column-major LAPACK => LU of A^T.
-    ////////////////////////////////////////////////////////////////////
-    if (boss) {
-#ifdef GRID_HIP
-      std::cout << GridLogMessage << "DenseCoarseMatrix: rocSOLVER cgetrf_64 (ILP64 LU) N=" << N
-                << " in place on resident device buffer" << std::endl;
-      rocblas_handle handle = GridBLASInverse::Handle();
-      int64_t *dInfo;
-      GRID_ASSERT( hipMalloc((void **)&dIpiv, N*sizeof(int64_t)) == hipSuccess );
-      GRID_ASSERT( hipMalloc((void **)&dInfo, sizeof(int64_t))   == hipSuccess );
-      auto st1 = rocsolver_cgetrf_64(handle, (int64_t)N, (int64_t)N, dA, (int64_t)N, dIpiv, dInfo);
-      GRID_ASSERT( hipDeviceSynchronize() == hipSuccess );
-      int64_t info_h = -1;
-      GRID_ASSERT( hipMemcpy(&info_h, dInfo, sizeof(int64_t),
-                             hipMemcpyDeviceToHost) == hipSuccess );
-      std::cout << GridLogMessage << "DenseCoarseMatrix: cgetrf_64 status " << (int)st1
-                << " info = " << (int)info_h << std::endl;
-      GRID_ASSERT(st1 == rocblas_status_success);
-      GRID_ASSERT(info_h == 0);
-      GRID_ASSERT( hipFree(dInfo) == hipSuccess );
-      GRID_ASSERT( hipMalloc((void **)&dB, (uint64_t)CHUNKROWS*N*sizeof(ComplexF)) == hipSuccess );
-      // dA holds the LU of A^T; rows of A^{-1} are produced blockwise below via
-      // cgetrs_64 on identity-column blocks: A^T X = E => X columns = rows of
-      // A^{-1}, in exactly the linear layout the harvest expects.
-#else
-      // Eigen fallback: small local CPU tests only.
-      std::cout << GridLogMessage << "DenseCoarseMatrix: Eigen fallback inversion N=" << N
-                << (N > 10000 ? "  (WARNING: SLOW; use the HIP/rocSOLVER path)" : "")
-                << std::endl;
-      typedef Eigen::Matrix<std::complex<float>,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> MatF;
-      Eigen::Map<MatF> A(reinterpret_cast<std::complex<float>*>(&Afull[0]), N, N);
-      MatF Ainv = A.inverse();
-      A = Ainv;
-#endif
-    }
-    double t3 = usecond();
-    std::cout << GridLogMessage << "DenseCoarseMatrix: factorisation took "
-              << (t3-t2)/1.0e6 << " s" << std::endl;
-
-    ////////////////////////////////////////////////////////////////////
-    // Blocked solve + broadcast: rows of A^{-1} chunk by chunk; each
-    // rank keeps the rows of its own sites (ownership-aligned).
-    ////////////////////////////////////////////////////////////////////
-    {
-      std::unordered_map<int64_t,int> rowmap;
-      for(int ss=0; ss<lsites; ss++)
-        for(int a=0; a<nbasis; a++)
-          rowmap[ myGsite[ss]*nbasis + a ] = ss*nbasis + a;
-
-      std::vector<ComplexF> chunk((uint64_t)CHUNKROWS * N);
-      for(int64_t row0=0; row0<N; row0+=CHUNKROWS){
-        int64_t nrow = std::min(CHUNKROWS, N-row0);
-        uint64_t nelem = (uint64_t)nrow * N;
-        if (boss) {
-#ifdef GRID_HIP
-          // Identity block E: column j = e_{row0+j}; solve A^T X = E so X's
-          // columns are rows [row0,row0+nrow) of A^{-1}.
-          for(uint64_t i=0;i<nelem;i++) chunk[i] = ComplexF(0.0,0.0);
-          for(int64_t j=0;j<nrow;j++) chunk[(uint64_t)j*N + (uint64_t)(row0+j)] = ComplexF(1.0,0.0);
-          GRID_ASSERT( hipMemcpy(dB, &chunk[0], nelem*sizeof(ComplexF), hipMemcpyHostToDevice) == hipSuccess );
-          auto strs = rocsolver_cgetrs_64(GridBLASInverse::Handle(), rocblas_operation_none,
-                                          (int64_t)N, (int64_t)nrow,
-                                          dA, (int64_t)N, dIpiv, dB, (int64_t)N);
-          GRID_ASSERT(strs == rocblas_status_success);
-          GRID_ASSERT( hipDeviceSynchronize() == hipSuccess );
-          GRID_ASSERT( hipMemcpy(&chunk[0], dB, nelem*sizeof(ComplexF), hipMemcpyDeviceToHost) == hipSuccess );
-#else
-          uint64_t src = (uint64_t)row0 * N;
-          for(uint64_t i=0;i<nelem;i++) chunk[i] = Afull[src+i];
-#endif
-        }
-        grid->Broadcast(0, &chunk[0], nelem*sizeof(ComplexF));
-        for(int64_t r=row0; r<row0+nrow; r++){
-          auto it = rowmap.find(r);
-          if (it != rowmap.end()) {
-            uint64_t dst = (uint64_t)(it->second) * N;
-            uint64_t src = (uint64_t)(r-row0) * N;
-            for(int64_t j=0;j<N;j++) slab[dst+j] = chunk[src+j];
-          }
-        }
-      }
-    }
-#ifdef GRID_HIP
-    if (boss) {
-      if (dA)    GRID_ASSERT( hipFree(dA)    == hipSuccess );
-      if (dB)    GRID_ASSERT( hipFree(dB)    == hipSuccess );
-      if (dIpiv) GRID_ASSERT( hipFree(dIpiv) == hipSuccess );
-    }
-#endif
-    double t4 = usecond();
-    std::cout << GridLogMessage << "DenseCoarseMatrix: blocked getrs solve+scatter took "
-              << (t4-t3)/1.0e6 << " s" << std::endl;
-  }
-
-  ////////////////////////////////////////////////////////////////////
-  // 3b. Global column -> rank-major column map, computed LOCALLY.
+  // 3a. Global column -> rank-major column map, computed LOCALLY.
   //    Rank-major ordering: rank q's rows/columns are the contiguous
   //    block [q*nrows, (q+1)*nrows), ordered by q's local site index
   //    (uniform local volumes make ownership arithmetic exact).
@@ -789,7 +468,7 @@ public:
   }
 
   ////////////////////////////////////////////////////////////////////
-  // 3c. Direct stencil -> fp64 rank-major import of MY ROWS of A (the
+  // 3b. Direct stencil -> fp64 rank-major import of MY ROWS of A (the
   //    end-to-end fp64 path: the stencil source IS ComplexD; nothing is
   //    rounded through fp32 on the way into the inversion).  Same
   //    loop/sign/accumulate/transposed-contraction discipline as
@@ -805,8 +484,6 @@ public:
   void ImportDenseFP64(CoarseOp &Op, BlockRows &S, std::vector<int64_t> &g2rm)
   {
     Coordinate gdims = grid->GlobalDimensions();
-    int sign = getenv("DENSE_IMPORT_SIGN") ? atoi(getenv("DENSE_IMPORT_SIGN")) : 1;
-    GRID_ASSERT( sign==1 || sign==-1 );
 
     std::vector<ComplexD> h((uint64_t)nrows*N, ComplexD(0.0,0.0));
     for(int p=0; p<Op.Geometry().npoint; p++)
@@ -818,7 +495,7 @@ public:
         Coordinate ncoor(nd);
         for(int d=0; d<nd; d++)
         {
-          int64_t g = grid->_lstart[d] + myLcoor[ss][d] + sign*shift[d];
+          int64_t g = grid->_lstart[d] + myLcoor[ss][d] + shift[d];
           ncoor[d] = (int)((g % gdims[d] + gdims[d]) % gdims[d]);
         }
         int64_t nsite;
@@ -868,14 +545,13 @@ public:
   }
 
   ////////////////////////////////////////////////////////////////////
-  // 3d. Distributed recursive Schur invert, END-TO-END fp64 (decision
-  //    2026-08-14): stencil (ComplexD) -> fp64 rank-major import ->
-  //    fp64 recursion -> ONE terminal rounding into the fp32 apply
-  //    slab.  Everything downstream (device residency, split-K apply,
-  //    VERIFY, SLAB_FILE) is untouched.
+  // 3c. The inverse: distributed recursive Schur, END-TO-END fp64.
+  //    stencil (ComplexD) -> fp64 rank-major import -> fp64 recursion ->
+  //    ONE terminal rounding into the fp32 apply slab.  Everything
+  //    downstream (device residency, split-K apply, VERIFY) is fp32.
   ////////////////////////////////////////////////////////////////////
   template<class CoarseOp>
-  void InvertDenseSchur(CoarseOp &Op)
+  void InvertDense(CoarseOp &Op)
   {
     double t1 = usecond();
     int P  = grid->ProcessorCount();
@@ -905,26 +581,17 @@ public:
     ImportDenseFP64(Op, S, g2rm);
 
     ////////////////////////////////////////////////////////////////
-    // DENSE_SCHUR2D=1 : invert via the 2D block-cyclic recursion
-    // (BlockCyclicSchurInverse) instead of the 1D rank-range one.
-    // The SAME imported rank-major rows S go in and come back, so the
-    // import certificate above and the slab rounding / VERIFY below are
-    // identical for both paths: a clean A/B on one imported operator.
-    //
-    // Everything in the 2D path -- redistribution, SUMMA rings, leaf --
-    // is point-to-point SendToRecvFrom; no collectives at all.
-    // DENSE_NB overrides the block size (default: rows-per-rank, which
-    // makes the redistribution edges maximally regular).
+    // The 2D block-cyclic recursion (BlockCyclicSchurInverse).
+    // Everything in it -- redistribution, SUMMA rings, leaves -- is
+    // point-to-point SendToRecvFrom; no collectives at all.  Block size
+    // nb = rows-per-rank makes the redistribution edges maximally
+    // regular.
     ////////////////////////////////////////////////////////////////
-    int use2d = getenv("DENSE_SCHUR2D") ? atoi(getenv("DENSE_SCHUR2D")) : 0;
-    int64_t panelBytes = getenv("DENSE_PANEL_BYTES") ? atol(getenv("DENSE_PANEL_BYTES"))
-                                                     : (int64_t)1024*1024*1024;   // 1D path only
     double t2, t3;
-    if ( use2d )
     {
       int Pr,Pc;
       BlockCyclicLayout::ChooseProcessGrid(P, Pr, Pc);
-      int64_t nb = getenv("DENSE_NB") ? atol(getenv("DENSE_NB")) : nrows;
+      int64_t nb = nrows;
       GRID_ASSERT( nb >= 1 );
       std::cout << GridLogMessage << "DenseCoarseMatrix: 2D SCHUR invert, process grid "
                 << Pr << " x " << Pc << "  nb " << nb
@@ -938,14 +605,6 @@ public:
       t3 = usecond();
       RSI2.ReportTelemetry(grid);
     }
-    else
-    {
-      RecursiveSchurInverse RSI(grid, N, rowStart, panelBytes);
-      t2 = usecond();
-      RSI.Invert(S);
-      t3 = usecond();
-      RSI.ReportTelemetry();
-    }
 
     // The single terminal rounding: fp64 inverse -> fp32 apply slab
     // (row-major, global columns)
@@ -962,28 +621,26 @@ public:
     }
     double t4 = usecond();
     std::cout << GridLogMessage << "DenseCoarseMatrix: SCHUR fp64 distributed invert took "
-              << (t4-t1)/1.0e6 << " s (recursion " << (t3-t2)/1.0e6 << " s), panelBytes "
-              << panelBytes << std::endl;
+              << (t4-t1)/1.0e6 << " s (recursion " << (t3-t2)/1.0e6 << " s)" << std::endl;
   }
 
   ////////////////////////////////////////////////////////////////////
   // CORE apply on packed data: hX[N x nr] (zero-filled, my sites only)
-  // -> allreduce -> split-K GEMM against the resident slab -> reduce
-  // partials -> hY[nrows x nr] (column major).  Platform-agnostic:
+  // -> ring allgather -> split-K GEMM against the resident slab ->
+  // reduce partials -> hY[nrows x nr] (column major).  Platform-agnostic:
   // deviceVector + GridBLAS (Eigen fallback on CPU builds).
-  // fp32 allreduce is EXACT: zero-fill assembly gives every element
-  // exactly one contributing rank.
+  // tprof (optional): per-phase microseconds {allgather, H2D, gemm+reduce,
+  // D2H}, printed by the caller on GridLogPerformance.
   ////////////////////////////////////////////////////////////////////
   void SlabApplyPacked(int nr, double *tprof)
   {
     GRID_TRACE("DenseSlabApply");
     GRID_ASSERT(nr <= MRHS_MAX);
-    uint64_t nX = (uint64_t)N * nr;
     uint64_t nY = (uint64_t)nrows * nr;
     int64_t  Kc = N / NK;
     double t1 = usecond();
     double t2, t3;
-    if (devSum==4) {
+    {
       // ALLGATHER: x is not a reduction -- every rank owns the rows of x at
       // global columns myGsite[ss]*nbasis+b (scattered by site coordinate, NOT
       // a contiguous block) and needs all of it.  Only MY rows go host->device
@@ -1011,31 +668,6 @@ public:
         });
       }
       t3 = usecond();
-    } else if (devSum) {
-      { GRID_TRACE("DenseH2D");
-        acceleratorCopyToDevice(&hX[0],&dX[0],nX*sizeof(ComplexF));
-      }
-      t2 = usecond();
-      { GRID_TRACE("DenseAllreduce");
-        // DENSE_DEVICE_SUM=1 : device-buffer MPI_Allreduce (Cray MPICH aborts
-        //                      above ~8 MB: 12 RHS at N=138240 is 13.3 MB)
-        // DENSE_DEVICE_SUM=2 : CartesianRingAllReduce, P2P only, no size cliff
-        // DENSE_DEVICE_SUM=3 : flat RingAllReduce, P2P only
-        // DENSE_DEVICE_SUM=4 : CartesianRingAllGather (branch above)
-        if      (devSum==2) CartesianRingAllReduce(grid,(ComplexF *)&dX[0],nX);
-        else if (devSum==3) RingAllReduce(grid,(ComplexF *)&dX[0],nX);
-        else                grid->GlobalSumVector((ComplexF *)&dX[0], (int)nX);
-      }
-      t3 = usecond();
-    } else {
-      { GRID_TRACE("DenseAllreduce");
-        grid->GlobalSumVector(&hX[0], (int)nX);
-      }
-      t2 = usecond();
-      { GRID_TRACE("DenseH2D");
-        acceleratorCopyToDevice(&hX[0],&dX[0],nX*sizeof(ComplexF));
-      }
-      t3 = usecond();
     }
     // Y = op(slab,T) . X : row-major slab (nrows x N) == col-major A^T
     // (N x nrows, lda=N) => transpose gives the nrows x N operator.
@@ -1065,8 +697,8 @@ public:
     }
     double t5 = usecond();
     if (tprof) {
-      tprof[0] = devSum ? (t3-t2) : (t2-t1);   // allreduce
-      tprof[1] = devSum ? (t2-t1) : (t3-t2);   // H2D
+      tprof[0] = t3-t2;                        // allgather
+      tprof[1] = t2-t1;                        // H2D
       tprof[2] = t4-t3;                        // gemm+reduce
       tprof[3] = t5-t4;                        // D2H
     }
@@ -1208,16 +840,14 @@ public:
     double t6 = usecond();
     std::cout << GridLogMessage << "DenseCoarseMatrix: apply6D " << nr << " rhs took "
               << (t6-t0)/1000.0 << " ms" << std::endl;
-    if ( getenv("DENSE_APPLY_PROFILE") ) {
-      std::cout << GridLogMessage << "DenseCoarseMatrix: apply6D profile:"
-                << " pack "        << (t1-t0)/1000.0
-                << (devSum==4 ? "  allgather " : "  allreduce ")  << tprof[0]/1000.0
-                << "  H2D "        << tprof[1]/1000.0
-                << "  gemm+reduce "<< tprof[2]/1000.0
-                << "  D2H "        << tprof[3]/1000.0
-                << "  unpack "     << (t6-t5)/1000.0
-                << "  ms" << std::endl;
-    }
+    std::cout << GridLogPerformance << "DenseCoarseMatrix: apply6D profile:"
+              << " pack "        << (t1-t0)/1000.0
+              << "  allgather "  << tprof[0]/1000.0
+              << "  H2D "        << tprof[1]/1000.0
+              << "  gemm+reduce "<< tprof[2]/1000.0
+              << "  D2H "        << tprof[3]/1000.0
+              << "  unpack "     << (t6-t5)/1000.0
+              << "  ms" << std::endl;
   }
 };
 

Grid/algorithms/multigrid/MrhsMultiGrid.h

@@ -0,0 +1,256 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/MrhsMultiGrid.h
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+#include <Grid/algorithms/deflation/MultiRHSBlockProject.h>
+
+NAMESPACE_BEGIN(Grid);
+
+//////////////////////////////////////////////////////////////////////
+// mrhs LinearFunction interface: vector-of-fields in, vector out.
+// The outer level carries mrhs as std::vector<Field>; below it mrhs
+// is PACKED into a single D+1 field (rhs = dim 0) and the coarse
+// classes are plain LinearFunctions on that.
+//////////////////////////////////////////////////////////////////////
+template<class Field>
+class MrhsLinearFunction {
+public:
+  virtual void operator()(std::vector<Field> &in, std::vector<Field> &out) = 0;
+};
+
+//////////////////////////////////////////////////////////////////////
+// Single-polynomial mrhs PGCR: one step length and one set of
+// orthogonalisation coefficients SHARED across the right-hand sides
+// (vnorm2/vinnerProduct sum over rhs).  For DWF the dense spectrum
+// makes sharing near-free and the summed reductions amortise; see the
+// mrhs-HDCG paper (arXiv:2409.03904).  Orthogonalisation is classical
+// Gram-Schmidt, batched per rhs (rankInnerProductMulti/axpyMulti) with
+// ONE GlobalSumVector for all coefficients.
+//
+// OnStep is a programmatic per-step hook for algorithmic studies
+// (e.g. driving the GCR coefficient recorder); it is not a consumer
+// parameter.
+//////////////////////////////////////////////////////////////////////
+template<class Field>
+class MrhsPGCRNonHermitian {
+public:
+  RealD Tolerance; Integer MaxIterations; int mmax,nstep,steps,level;
+  int ZeroGuess = 0; int FirstCycle = 0;
+  std::string name = "Level 1";
+  LinearOperatorBase<Field> &Linop;
+  MrhsLinearFunction<Field> &Preconditioner;
+  std::function<void(int)> OnStep;      // called with the outer step count after every step
+  void Level(int lv){ name = "Level " + std::to_string(lv); level=lv; }
+  void Name(std::string n){ name = n; }
+  void SetZeroGuess(int z){ ZeroGuess=z; }
+  MrhsPGCRNonHermitian(RealD tol,Integer maxit,LinearOperatorBase<Field> &_Linop,MrhsLinearFunction<Field> &Prec,int _mmax,int _nstep)
+    : Tolerance(tol),MaxIterations(maxit),Linop(_Linop),Preconditioner(Prec),mmax(_mmax),nstep(_nstep){ level=1; }
+  static RealD vnorm2(std::vector<Field> &x){ RealD s=0; for(auto &f:x) s+=norm2(f); return s; }
+  static ComplexD vinnerProduct(std::vector<Field> &x,std::vector<Field> &y){ ComplexD s(0); for(int r=0;r<(int)x.size();r++) s+=innerProduct(x[r],y[r]); return s; }
+  static void vaxpy(std::vector<Field> &z,ComplexD a,std::vector<Field> &x,std::vector<Field> &y){ for(int r=0;r<(int)z.size();r++) axpy(z[r],a,x[r],y[r]); }
+  void vOp(std::vector<Field> &in,std::vector<Field> &out){ GRID_TRACE("MrhsPGCR::vOp"); for(int r=0;r<(int)in.size();r++) Linop.Op(in[r],out[r]); }
+  void operator()(std::vector<Field> &src,std::vector<Field> &psi){
+    RealD cp,ssq,rsq; int nrhs=src.size(); GridBase *grid=src[0].Grid();
+    ssq=vnorm2(src); rsq=Tolerance*Tolerance*ssq;
+    std::vector<Field> r(nrhs,grid);
+    GridStopWatch T; T.Start(); steps=0; FirstCycle=1;
+    for(int k=0;k<MaxIterations;k++){
+      cp=GCRnStep(src,psi,rsq);
+      std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR("<<mmax<<","<<nstep<<") "<<steps<<" steps cp = "<<cp<<" target "<<rsq<<std::endl;
+      if(cp<rsq){
+        T.Stop(); vOp(psi,r); for(int rr=0;rr<nrhs;rr++) axpy(r[rr],-1.0,src[rr],r[rr]);
+        RealD tr=vnorm2(r);
+        std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR: Converged on iteration "<<steps
+                 <<" computed residual "<<std::sqrt(cp/ssq)<<" true residual "<<std::sqrt(tr/ssq)<<" target "<<Tolerance<<std::endl;
+        std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR Time elapsed: Total "<<T.Elapsed()<<std::endl;
+        return;
+      }
+    }
+    std::cout<<GridLogMessage<<"MrhsPGCR: did not converge"<<std::endl;
+  }
+  RealD GCRnStep(std::vector<Field> &src,std::vector<Field> &psi,RealD rsq){
+    RealD cp; ComplexD a,rq; int nrhs=src.size(); GridBase *grid=src[0].Grid();
+    std::vector<Field> r(nrhs,grid),Az(nrhs,grid);   // Az: restart residual scratch only
+    std::vector< std::vector<Field> > q(mmax,std::vector<Field>(nrhs,grid));
+    std::vector< std::vector<Field> > p(mmax,std::vector<Field>(nrhs,grid));
+    std::vector<RealD> qq(mmax);
+    if (ZeroGuess && FirstCycle) { for(int rr=0;rr<nrhs;rr++){ psi[rr]=Zero(); r[rr]=src[rr]; } }
+    else                         { vOp(psi,Az); for(int rr=0;rr<nrhs;rr++) r[rr]=src[rr]-Az[rr]; }
+    FirstCycle=0;
+    // p[0]=Prec(r), q[0]=A p[0], produced directly in the history slots (no copies)
+    Preconditioner(r,p[0]); vOp(p[0],q[0]); qq[0]=vnorm2(q[0]); cp=vnorm2(r);
+    for(int k=0;k<nstep;k++){
+      steps++; int kp=k+1, peri_k=k%mmax, peri_kp=kp%mmax;
+      if ( OnStep ) OnStep(steps);
+      rq=vinnerProduct(q[peri_k],r); a=rq/qq[peri_k];
+      vaxpy(psi,a,p[peri_k],psi); vaxpy(r,-a,q[peri_k],r); cp=vnorm2(r);
+      std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR step["<<steps<<"]  resid "<<cp<<" target "<<rsq<<std::endl;
+      if((k==nstep-1)||(cp<rsq)) return cp;
+      // New direction straight into its history slot: p=Prec(r), q=A p.
+      Preconditioner(r,p[peri_kp]);
+      vOp(p[peri_kp],q[peri_kp]);
+      int northog=((kp)>(mmax-1))?(mmax-1):(kp);
+      {
+        GRID_TRACE("MrhsPGCR orthog");
+        // Classical Gram-Schmidt: all coefficients against the UN-updated new q
+        // (independent, batchable), then apply.  Complex coefficient: the
+        // operator is non-Hermitian, real(<q_j,Aq>) alone left q's non-orthogonal.
+        // Batched per rhs (one fused kernel + one reduction each), the shared
+        // coefficient summed over rhs on the host, ONE GlobalSumVector.
+        std::vector<ComplexD> bcoef(northog,ComplexD(0.0)), part;
+        for(int rr=0;rr<nrhs;rr++){
+          std::vector<const Field*> qwin(northog);
+          for(int back=0;back<northog;back++){ int peri_back=(k-back)%mmax; GRID_ASSERT((k-back)>=0); qwin[back]=&q[peri_back][rr]; }
+          rankInnerProductMulti(part,qwin,q[peri_kp][rr]);
+          for(int back=0;back<northog;back++) bcoef[back]+=part[back];
+        }
+        if(northog) grid->GlobalSumVector(&bcoef[0],northog);
+        for(int back=0;back<northog;back++){ int peri_back=(k-back)%mmax; bcoef[back]=-bcoef[back]/qq[peri_back]; }
+        for(int rr=0;rr<nrhs;rr++){
+          std::vector<const Field*> qwin(northog), pwin(northog);
+          for(int back=0;back<northog;back++){ int peri_back=(k-back)%mmax; qwin[back]=&q[peri_back][rr]; pwin[back]=&p[peri_back][rr]; }
+          axpyMulti(p[peri_kp][rr],bcoef,pwin);
+          axpyMulti(q[peri_kp][rr],bcoef,qwin);
+        }
+      }
+      qq[peri_kp]=vnorm2(q[peri_kp]);
+    }
+    GRID_ASSERT(0); return cp;
+  }
+};
+
+//////////////////////////////////////////////////////////////////////
+// Dense L3 solve on the packed D+1 coarse-coarse field.  Adapts any
+// class exposing ApplyBatch6D (DenseCoarseMatrix) to a LinearFunction.
+//////////////////////////////////////////////////////////////////////
+template<class DenseType, class CoarseCoarseField>
+class MrhsDenseCCSolve : public LinearFunction<CoarseCoarseField> {
+public:
+  DenseType &_Dense;
+  int _nrhs;
+  MrhsDenseCCSolve(DenseType &D, int nrhs) : _Dense(D), _nrhs(nrhs) {}
+  using LinearFunction<CoarseCoarseField>::operator();
+  virtual void operator()(const CoarseCoarseField &in, CoarseCoarseField &out){
+    _Dense.ApplyBatch6D(in, out, _nrhs);
+  }
+};
+
+//////////////////////////////////////////////////////////////////////
+// L2->L3 half V-cycle on the D+1 coarse field: coarse-coarse correction
+// through the mixed blockProject, then post-smooth.
+//////////////////////////////////////////////////////////////////////
+template<class CoarseField, class CoarseCoarseField>
+class MrhsCoarseThreeLevelPrec : public LinearFunction<CoarseField> {
+public:
+  LinearOperatorBase<CoarseField>          &_CoarseOp;
+  LinearFunction<CoarseField>              &_CoarseSmoother;
+  MultiRHSBlockProject<CoarseField>        &_Projector;
+  LinearFunction<CoarseCoarseField>        &_CoarseCoarseSolve;
+  GridBase *_Coarse5d, *_CoarseCoarse5d, *_CoarseCoarseMrhs;
+  int _nrhs;
+  MrhsCoarseThreeLevelPrec(LinearOperatorBase<CoarseField> &CoarseOp,
+                           LinearFunction<CoarseField> &CoarseSmoother,
+                           MultiRHSBlockProject<CoarseField> &Projector,
+                           LinearFunction<CoarseCoarseField> &CoarseCoarseSolve,
+                           GridBase *Coarse5d, GridBase *CoarseCoarse5d, GridBase *CoarseCoarseMrhs, int nrhs)
+    : _CoarseOp(CoarseOp), _CoarseSmoother(CoarseSmoother), _Projector(Projector),
+      _CoarseCoarseSolve(CoarseCoarseSolve),
+      _Coarse5d(Coarse5d), _CoarseCoarse5d(CoarseCoarse5d), _CoarseCoarseMrhs(CoarseCoarseMrhs), _nrhs(nrhs) {}
+  using LinearFunction<CoarseField>::operator();
+  virtual void operator()(const CoarseField &in, CoarseField &out) {
+    CoarseField vec1(in.Grid());
+    CoarseField vec2(in.Grid());
+    out = in;
+    _CoarseOp.Op(out,vec1);  sub(vec1,in,vec1);
+
+    // restrict, through the mixed blockProject: D+1 coarse in, D+1 cc out
+    CoarseCoarseField CCsrc(_CoarseCoarseMrhs);
+    CoarseCoarseField CCsol(_CoarseCoarseMrhs);
+    _Projector.blockProject(vec1,CCsrc);
+
+    _CoarseCoarseSolve(CCsrc,CCsol);
+
+    _Projector.blockPromote(vec1,CCsol);
+    add(out,out,vec1);
+
+    _CoarseOp.Op(out,vec1);  sub(vec1,in,vec1);
+    _CoarseSmoother(vec1,vec2);
+    add(out,out,vec2);
+  }
+};
+
+//////////////////////////////////////////////////////////////////////
+// L1->L2 mrhs V-cycle.  The whole V-cycle is preconditioner: its fine
+// residuals and the smoother may run with sloppy halos; the caller
+// (the outer Krylov) gets the exact operator back on exit.  SetSloppy
+// is wired by the composer to PVdagMLinearOperator::SloppyComms (a
+// no-op by default), replacing the file-scope global the example used.
+//////////////////////////////////////////////////////////////////////
+template<class FineField, class MrhsCoarseVector, class FineSmoother>
+class MrhsTwoLevelMG : public MrhsLinearFunction<FineField> {
+public:
+  typedef MrhsCoarseVector CoarseVector;
+  LinearOperatorBase<FineField>   &_FineOperator;
+  FineSmoother                    &_PostSmoother;
+  MultiRHSBlockProject<FineField> &_Projector;
+  LinearFunction<CoarseVector>    &_CoarseSolve;
+  GridBase *_CoarseGrid, *_CoarseGridMrhs;
+  std::function<void(int)> SetSloppy = [](int){};
+  int SloppyComms = 0;                 // value passed to SetSloppy on entry
+  MrhsTwoLevelMG(LinearOperatorBase<FineField> &FineOp, FineSmoother &Post,
+                 MultiRHSBlockProject<FineField> &Projector, LinearFunction<CoarseVector> &CoarseSolve,
+                 GridBase *CoarseGrid, GridBase *CoarseGridMrhs)
+    : _FineOperator(FineOp),_PostSmoother(Post),_Projector(Projector),_CoarseSolve(CoarseSolve),
+      _CoarseGrid(CoarseGrid),_CoarseGridMrhs(CoarseGridMrhs){}
+  virtual void operator()(std::vector<FineField> &in, std::vector<FineField> &out){
+    GRID_TRACE("MGVcycle");
+    SetSloppy(SloppyComms);
+    int nrhs=in.size(); GridBase *fgrid=in[0].Grid();
+    std::vector<FineField> vec1(nrhs,fgrid),vec2(nrhs,fgrid);
+    for(int r=0;r<nrhs;r++) out[r]=in[r];
+    { GRID_TRACE("MGFineResidual");
+      for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); }
+    }
+    // fine vector -> D+1 coarse, via the mixed blockProject
+    CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs);
+    { GRID_TRACE("MGProject");
+      _Projector.blockProject(vec1,CsrcMrhs);
+    }
+    CsolMrhs=Zero();
+    { GRID_TRACE("MGCoarseSolve");
+      _CoarseSolve(CsrcMrhs,CsolMrhs);
+    }
+    { GRID_TRACE("MGPromote");
+      _Projector.blockPromote(vec1,CsolMrhs);
+      for(int r=0;r<nrhs;r++) add(out[r],out[r],vec1[r]);
+    }
+    { GRID_TRACE("MGFineResidual2");
+      for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); }
+    }
+    { GRID_TRACE("MGPostSmooth");
+      for(int r=0;r<nrhs;r++){
+        _PostSmoother(vec1[r],vec2[r]); add(out[r],out[r],vec2[r]);
+      }
+    }
+    SetSloppy(0);
+  }
+};
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/MultiGrid.h

@@ -35,3 +35,8 @@ Author: Peter Boyle <pboyle@bnl.gov>
 #include <Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHSV2.h>
 #include <Grid/algorithms/multigrid/MrhsPromotedOperator.h>
 #include <Grid/algorithms/multigrid/Smoothers.h>
+#include <Grid/algorithms/multigrid/PVdagMMultiGridParams.h>
+// PVdagMOperators.h / MrhsMultiGrid.h / PVdagMMultiGrid.h /
+// DenseCoarseMatrix.h are NOT in this umbrella: consumers of the PVdagM
+// chain include PVdagMMultiGrid.h explicitly (it pulls the dense stack
+// and BLAS).

Grid/algorithms/multigrid/PVdagMMultiGrid.h

@@ -0,0 +1,536 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/PVdagMMultiGrid.h
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+#include <Grid/algorithms/multigrid/DenseCoarseMatrix.h>
+#include <Grid/algorithms/multigrid/PVdagMOperators.h>
+#include <Grid/algorithms/multigrid/MrhsMultiGrid.h>
+
+NAMESPACE_BEGIN(Grid);
+
+//////////////////////////////////////////////////////////////////////////////////////
+// The three-level mrhs PVdagM multigrid, as objects:
+//
+//   MGCoarseGrids              the derived coarse grids (owned here, since
+//                              conformability is pointer identity and someone
+//                              must hold them; declare it BEFORE anything that
+//                              borrows from it)
+//   PVdagMMultiGridCoarsening  ALL the coarsening information for one gauge
+//                              configuration: the raw near-null basis, the
+//                              transfer operators it induces, the Galerkin
+//                              coarse operator at each level, and the dense
+//                              bottom inverse.  Everything here is a function
+//                              of the gauge field -- the state HMC must
+//                              rebuild or maintain as U evolves -- and
+//                              everything is solver-independent.
+//   PVdagMMultiGridSolver      the solve chain composed on a borrowed
+//                              coarsening: smoothers, coarse Krylov, V-cycle,
+//                              outer mrhs PGCR.
+//
+// Scope discipline: grids outlive the coarsening outlives the solver.
+//////////////////////////////////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////////////////
+// Subspace I/O: bare-vector scidac records.  Loaded vectors are RAW --
+// deliberately NOT re-orthogonalised: CoarsenOperator block
+// orthonormalises in place, and projecting a block-orthonormal vector
+// onto its own block-orthonormalised aggregation gives e_k with the
+// near-null content silently gone.  GramGuard below catches that.
+//////////////////////////////////////////////////////////////////////
+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
+}
+
+//////////////////////////////////////////////////////////////////////
+// ||<v|v> - I||_F over a set of coarse vectors.  Small means the raw
+// near-null content survived the projection.
+//////////////////////////////////////////////////////////////////////
+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 );
+}
+
+//////////////////////////////////////////////////////////////////////
+// The derived grids of the three-level chain, owned in one place.
+// The coarse space is UNVECTORISED (sComplex scalar, simd {1,..,1});
+// the 5D/6D grids are built directly so the SIMD layout is ours.
+// rhs/batch is dim 0 of the 6D grids, undistributed -- no divisibility
+// constraint on nrhs.
+//////////////////////////////////////////////////////////////////////
+class MGCoarseGrids {
+public:
+  GridCartesian *FGrid;                     // borrowed
+  int Ls;
+  int batch;
+  Coordinate clatt;                         // 4d coarse lattice
+  Coordinate cclatt;                        // 4d coarse-coarse lattice
+  Coordinate c5simd, c5mpi;                 // 5D coarse simd/mpi
+  Coordinate cmsimd, cmmpi;                 // 6D coarse simd/mpi
+  // owned:
+  GridCartesian *Coarse5d;
+  GridCartesian *CoarseBatch;               // 6D at the coarsening batch
+  GridCartesian *CoarseCoarse5d;
+  GridCartesian *CoarseCoarseBatch;
+
+  MGCoarseGrids(GridCartesian *_FGrid, const MGSetupParams &P)
+    : FGrid(_FGrid)
+  {
+    Coordinate fdims = FGrid->FullDimensions();      // {Ls, x,y,z,t}
+    Coordinate fmpi  = FGrid->_processors;
+    GRID_ASSERT( fdims.size() == 5 );
+    Ls    = fdims[0];
+    batch = P.CoarsenBatch;
+
+    clatt.resize(4); cclatt.resize(4);
+    for(int d=0;d<4;d++){
+      GRID_ASSERT( fdims[d+1] % P.Block[d] == 0 );
+      clatt[d] = fdims[d+1] / P.Block[d];
+    }
+    for(int d=0;d<4;d++){
+      GRID_ASSERT( clatt[d] % P.Block2[d] == 0 );
+      cclatt[d] = clatt[d] / P.Block2[d];
+    }
+    std::cout << GridLogMessage << "MGCoarseGrids: Block  " << P.Block  << "  coarse lattice        " << clatt  << std::endl;
+    std::cout << GridLogMessage << "MGCoarseGrids: Block2 " << P.Block2 << "  coarse-coarse lattice " << cclatt << std::endl;
+
+    Coordinate c5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]});
+    c5simd = Coordinate({1,1,1,1,1});
+    c5mpi  = Coordinate({1,fmpi[1],fmpi[2],fmpi[3],fmpi[4]});
+    Coarse5d = new GridCartesian(c5latt,c5simd,c5mpi);
+
+    cmsimd = Coordinate({1,1,1,1,1,1});
+    cmmpi  = Coordinate({1,1,fmpi[1],fmpi[2],fmpi[3],fmpi[4]});
+    Coordinate cblatt({batch,1,clatt[0],clatt[1],clatt[2],clatt[3]});
+    CoarseBatch = new GridCartesian(cblatt,cmsimd,cmmpi);
+
+    Coordinate cc5latt({1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+    CoarseCoarse5d = new GridCartesian(cc5latt,c5simd,c5mpi);
+
+    Coordinate ccblatt({batch,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+    CoarseCoarseBatch = new GridCartesian(ccblatt,cmsimd,cmmpi);
+  }
+  ~MGCoarseGrids()
+  {
+    delete CoarseCoarseBatch;
+    delete CoarseCoarse5d;
+    delete CoarseBatch;
+    delete Coarse5d;
+  }
+};
+
+//////////////////////////////////////////////////////////////////////
+// ALL the coarsening information for one gauge configuration.
+//
+// Owns: the RAW near-null bases (fine and coarse -- retained so the
+// coarse operators can be REBUILT on a changed gauge field with a
+// fixed basis, the cheap HMC maintenance step; DiscardBasis() frees
+// them for valence use), the two coarsened operators, the two block
+// projectors, the dense bottom inverse (raw owning pointer: null in
+// the constructor, allocated by BuildDenseBottom, deleted here), and
+// the transient per-Nrhs solve grids (created by SetNrhs, which owns
+// the ReleaseGrid -> delete -> new -> SetGrid ordering in ONE place).
+//
+// Borrows: the grid bundle.  Declare MGCoarseGrids first.
+//////////////////////////////////////////////////////////////////////
+template<class Fobj,class CComplex,int nbasis>
+class PVdagMMultiGridCoarsening {
+public:
+  typedef Lattice<Fobj>                                                 FineField;
+  typedef MultiGeneralCoarsenedOperatorV2<Fobj,CComplex,nbasis>         CoarseOperator;
+  typedef typename CoarseOperator::CoarseVector                         CoarseVector;
+  typedef typename CoarseVector::vector_object                          CoarseSiteObj;
+  typedef iScalar<CComplex>                                             CComplex2;  // coarsening deepens the nest by one iScalar
+  typedef MultiGeneralCoarsenedOperatorV2<CoarseSiteObj,CComplex2,nbasis> CoarseCoarseOperator;
+  typedef typename CoarseCoarseOperator::CoarseVector                   CoarseCoarseVector;
+  typedef DenseCoarseMatrix<CComplex2,nbasis>                           DenseBottom;
+
+  MGCoarseGrids                     &Grids;      // borrowed
+  MGSetupParams                      Params;
+  NextToNearestStencilGeometry5D     geom;
+  NextToNearestStencilGeometry5D     geom2;
+  CoarseOperator                     CoarseOpPV;
+  CoarseCoarseOperator               CoarseOpL2;
+  MultiRHSBlockProject<FineField>    MrhsProjector;
+  MultiRHSBlockProject<CoarseVector> MrhsProjectorL2;
+  DenseBottom                       *DenseCC;
+  std::vector<FineField>             rawNull;    // RAW fine near-null basis
+  std::vector<CoarseVector>          rawPsi;     // RAW coarse near-null basis
+  GridCartesian                     *CMrhs;      // transient solve grids, owned
+  GridCartesian                     *CCMrhs;
+  int                                nrhs;
+
+  PVdagMMultiGridCoarsening(MGCoarseGrids &_Grids, const MGSetupParams &P)
+    : Grids(_Grids),
+      Params(P),
+      geom (_Grids.Coarse5d),
+      geom2(_Grids.CoarseCoarse5d),
+      CoarseOpPV(geom ,_Grids.Coarse5d),
+      CoarseOpL2(geom2,_Grids.CoarseCoarse5d),
+      DenseCC(nullptr),
+      CMrhs(nullptr),
+      CCMrhs(nullptr),
+      nrhs(-1)
+  {};
+
+  ~PVdagMMultiGridCoarsening()
+  {
+    if ( DenseCC ) delete DenseCC;
+    CoarseOpPV.ReleaseGrid();     // borrowers let go before their grids die
+    CoarseOpL2.ReleaseGrid();
+    if ( CMrhs  ) delete CMrhs;
+    if ( CCMrhs ) delete CCMrhs;
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // The RAW fine basis: load from Params.SubspaceFile if it exists,
+  // else GCR inverse iteration from noise (and save if a file name
+  // was given).  The Aggregation is scaffolding for CreateSubspaceGCR
+  // only -- that runs entirely on the fine grid, so the coarse grid
+  // it holds is never dereferenced.
+  ////////////////////////////////////////////////////////////////////
+  void GetSubspace(GridParallelRNG &RNG, LinearOperatorBase<FineField> &FineOp)
+  {
+    uint64_t file_exists=0;
+    if ( Params.SubspaceFile.length() ) {
+      if ( Grids.FGrid->IsBoss() ){ std::ifstream f(Params.SubspaceFile); file_exists=f.good()?1:0; }
+      Grids.FGrid->GlobalSum(file_exists);
+    }
+    rawNull.clear();
+    rawNull.reserve(nbasis);
+    for(int k=0;k<nbasis;k++) rawNull.push_back(FineField(Grids.FGrid));
+    if ( file_exists ){
+      std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: loading subspace "
+                << Params.SubspaceFile << " (kept RAW)" << std::endl;
+      loadSubspace(rawNull, Params.SubspaceFile);
+    } else {
+      std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: GCR subspace generation" << std::endl;
+      Aggregation<Fobj,CComplex,nbasis> Agg(Grids.Coarse5d,Grids.FGrid,0);
+      Agg.CreateSubspaceGCR(RNG,FineOp,nbasis);
+      for(int k=0;k<nbasis;k++)
+	rawNull[k]=Agg.subspace[k];
+      if ( Params.SubspaceFile.length() )
+	saveSubspace(rawNull, Params.SubspaceFile);
+    }
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // Both coarsenings from the retained RAW basis.  Calling this again
+  // on a CHANGED fine operator is the fixed-basis rebuild: for the
+  // PVdagM coarsening the coarse operator with a fixed basis is exact
+  // in the changed links.  CoarsenOperator block-orthonormalises its
+  // input in place, so working copies are taken and the RAW vectors
+  // survive for the next rebuild.
+  //
+  // The fine operator must be the PVdagM wrapper (it needs
+  // SloppyComms: the coarsening builds the PRECONDITIONER, so its
+  // halos follow the FineSloppyComms policy, restored EXACT on exit).
+  ////////////////////////////////////////////////////////////////////
+  template<class PVdagMOp>
+  void Coarsen(PVdagMOp &FineOp)
+  {
+    GRID_ASSERT( rawNull.size() == nbasis );
+
+    // L1: working copy of the raw fine basis, orthonormalised in place
+    std::vector<FineField> sub(rawNull);
+
+    CoarseOpPV.SetGrid(Grids.CoarseBatch);
+
+    std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: L1 CoarsenOperator, batch "
+              << Grids.batch << std::endl;
+
+    FineOp.SloppyComms(Params.FineSloppyComms);
+    CoarseOpPV.CoarsenOperator(FineOp,sub,Grids.Coarse5d,Grids.batch);
+    FineOp.SloppyComms(0);
+
+    // Transfer operators from the orthonormalised basis; then the
+    // Galerkin images of the RAW basis define the L2 null space.
+    MrhsProjector.Allocate(nbasis,Grids.FGrid,Grids.Coarse5d);
+    MrhsProjector.ImportBasis(sub);
+    sub.clear(); sub.shrink_to_fit();
+
+    std::vector<CoarseVector> psi(nbasis,Grids.Coarse5d);
+    MrhsProjector.blockProject(rawNull,psi);
+    GramGuard("psi_coarse",psi,Grids.Coarse5d);
+
+    rawPsi.clear();
+    rawPsi.reserve(nbasis);
+    for(int k=0;k<nbasis;k++) rawPsi.push_back(psi[k]);
+
+    // L2: the V2 L1 operator is natively multiRHS at the batch grid
+    CoarseOpL2.SetGrid(Grids.CoarseCoarseBatch);
+    NonHermitianLinearOperator<CoarseOperator,CoarseVector> LinOpCoarse(CoarseOpPV);
+    std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: L2 CoarsenOperator, batch "
+              << Grids.batch << std::endl;
+    CoarseOpL2.CoarsenOperator(LinOpCoarse,Grids.CoarseBatch,psi,Grids.CoarseCoarse5d);
+
+    MrhsProjectorL2.Allocate(nbasis,Grids.Coarse5d,Grids.CoarseCoarse5d);
+    MrhsProjectorL2.ImportBasis(psi);              // now block orthonormal
+    {
+      std::vector<CoarseCoarseVector> psi_cc(nbasis,Grids.CoarseCoarse5d);
+      MrhsProjectorL2.blockProject(rawPsi,psi_cc); // RAW vectors in
+      GramGuard("psi_cc",psi_cc,Grids.CoarseCoarse5d);
+    }
+    nrhs = -1;                                     // both operators left at the batch grids
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // Dense bottom: import the L2 matrix elements, invert (2D
+  // block-cyclic, fp64), then the check Import cannot run itself --
+  // the mrhs operator applies only on a D+1 grid, so drive it at
+  // Nrhs 1 through a slice: ||A Ainv x - x||/||x||.
+  ////////////////////////////////////////////////////////////////////
+  void BuildDenseBottom(void)
+  {
+    if ( DenseCC ) delete DenseCC;
+    std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: L3 dense bottom import" << std::endl;
+    DenseCC = new DenseBottom(Grids.CoarseCoarse5d);
+    DenseCC->Import(CoarseOpL2);
+
+    Coordinate cc1latt({1,1,Grids.cclatt[0],Grids.cclatt[1],Grids.cclatt[2],Grids.cclatt[3]});
+    GridCartesian *CoarseCoarseOne = new GridCartesian(cc1latt,Grids.cmsimd,Grids.cmmpi);
+    CoarseOpL2.SetGrid(CoarseCoarseOne);
+
+    CoarseCoarseVector x(Grids.CoarseCoarse5d),y(Grids.CoarseCoarse5d),z(Grids.CoarseCoarse5d);
+    GridParallelRNG dRNG(Grids.CoarseCoarse5d); dRNG.SeedFixedIntegers(std::vector<int>({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 << "PVdagMMultiGridCoarsening: L3 dense ||A Ainv x - x||/||x|| = "
+              << rel << std::endl;
+    GRID_ASSERT( rel < 1.0e-2 );
+
+    CoarseOpL2.ReleaseGrid();           // let go before the grid dies
+    delete CoarseCoarseOne;
+    nrhs = -1;
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // Retarget both operators to a solve Nrhs.  The matrix elements are
+  // Nrhs independent and survive; only the D+1 grids and BLAS buffers
+  // change.  The ordering discipline (borrowers release BEFORE their
+  // grids are destroyed) lives here and nowhere else.
+  ////////////////////////////////////////////////////////////////////
+  void SetNrhs(int nr)
+  {
+    if ( nr == nrhs ) return;
+    CoarseOpPV.ReleaseGrid();
+    CoarseOpL2.ReleaseGrid();
+    if ( CMrhs  ) { delete CMrhs;  CMrhs =nullptr; }
+    if ( CCMrhs ) { delete CCMrhs; CCMrhs=nullptr; }
+    Coordinate cml ({nr,1,Grids.clatt[0], Grids.clatt[1], Grids.clatt[2], Grids.clatt[3]});
+    Coordinate ccml({nr,1,Grids.cclatt[0],Grids.cclatt[1],Grids.cclatt[2],Grids.cclatt[3]});
+    CMrhs  = new GridCartesian(cml, Grids.cmsimd,Grids.cmmpi);
+    CCMrhs = new GridCartesian(ccml,Grids.cmsimd,Grids.cmmpi);
+    CoarseOpPV.SetGrid(CMrhs);
+    CoarseOpL2.SetGrid(CCMrhs);
+    nrhs = nr;
+    std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: operators at Nrhs " << nr << std::endl;
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // Free the retained RAW bases (valence use: coarsening is final for
+  // this configuration and the fine basis is ~GB-scale host memory).
+  ////////////////////////////////////////////////////////////////////
+  void DiscardBasis(void)
+  {
+    rawNull.clear(); rawNull.shrink_to_fit();
+    rawPsi.clear();  rawPsi.shrink_to_fit();
+  }
+};
+
+//////////////////////////////////////////////////////////////////////
+// The solve chain on a borrowed coarsening: adaptive PGCR smoothers
+// at both levels (the one correct route for this non-Hermitian
+// chain), dense bottom, mrhs V-cycle, outer shared-coefficient mrhs
+// PGCR.  Construction retargets the coarsening to this Nrhs.
+//
+// Halo policy: the whole V-cycle is preconditioner and runs with
+// sloppy halos when FineSloppyComms is set; the outer Krylov is EXACT
+// (its applications define what "converged" means), asserted by the
+// FINAL true-residual report in Solve.
+//////////////////////////////////////////////////////////////////////
+template<class Matrix,class Coarsening>
+class PVdagMMultiGridSolver {
+public:
+  typedef typename Coarsening::FineField            FineField;
+  typedef typename Coarsening::CoarseOperator       CoarseOperator;
+  typedef typename Coarsening::CoarseVector         CoarseVector;
+  typedef typename Coarsening::CoarseCoarseOperator CoarseCoarseOperator;
+  typedef typename Coarsening::CoarseCoarseVector   CoarseCoarseVector;
+  typedef typename Coarsening::DenseBottom          DenseBottom;
+  typedef PrecGeneralisedConjugateResidualNonHermitian<FineField>    FineSmoother_t;
+  typedef PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> CoarseKrylov_t;
+
+  Coarsening            &C;
+  PVdagMMultiGridParams  Params;
+  int                    nrhs;
+  int                    _regrid;    // FIRST member-like init: SetNrhs before members capture grids
+
+  PVdagMLinearOperator<Matrix,FineField>        PVdagM;
+  ShiftedPVdagMLinearOperator<Matrix,FineField> ShiftedPVdagM;
+  TrivialPrecon<FineField>                      simple_fine;
+  TrivialPrecon<CoarseVector>                   simpleC;
+  NonHermitianLinearOperator<CoarseOperator,CoarseVector> LinOpC;
+  ShiftedLinearOperator<CoarseVector>           ShiftedC;
+  MrhsDenseCCSolve<DenseBottom,CoarseCoarseVector> ccSolve;
+  CoarseKrylov_t                                CoarseSmootherGCR;
+  MrhsCoarseThreeLevelPrec<CoarseVector,CoarseCoarseVector> L2to3Precon;
+  CoarseKrylov_t                                L2PGCR;
+  FineSmoother_t                                SmootherGCR;
+  MrhsTwoLevelMG<FineField,CoarseVector,FineSmoother_t> ThreeLevelPrecon;
+  MrhsPGCRNonHermitian<FineField>               L1PGCR;
+
+  PVdagMMultiGridSolver(Matrix &Ddwf, Matrix &Dpv,
+                        Coarsening &_C, const PVdagMMultiGridParams &P, int nr)
+    : C(_C), Params(P), nrhs(nr),
+      _regrid((C.SetNrhs(nr),0)),              // members below capture C.CMrhs/C.CCMrhs
+      PVdagM(Ddwf,Dpv),
+      ShiftedPVdagM(P.FineSmoother.Shift,Ddwf,Dpv),
+      LinOpC(C.CoarseOpPV),
+      ShiftedC(P.CoarseSmoother.Shift,LinOpC),
+      ccSolve(*C.DenseCC,nr),
+      // Coarse smoother: one preconditioned restart of nstep GCR steps.
+      CoarseSmootherGCR(0.01,1,ShiftedC,simpleC,P.CoarseSmoother.Mmax,P.CoarseSmoother.Nstep),
+      L2to3Precon(LinOpC,CoarseSmootherGCR,C.MrhsProjectorL2,ccSolve,
+                  C.Grids.Coarse5d,C.Grids.CoarseCoarse5d,C.CCMrhs,nr),
+      // Coarse Krylov: Order/16 restarts of 16 steps.
+      L2PGCR(P.CoarseSolver.Tol,P.CoarseSolver.Order/16,LinOpC,L2to3Precon,P.CoarseSolver.Mmax,16),
+      // Fine smoother: one restart of nstep GCR steps, tolerance 0 = fixed work.
+      SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,P.FineSmoother.Mmax,P.FineSmoother.Nstep),
+      ThreeLevelPrecon(PVdagM,SmootherGCR,C.MrhsProjector,L2PGCR,C.Grids.Coarse5d,C.CMrhs),
+      L1PGCR(P.Outer.Tol,P.Outer.MaxIterations,PVdagM,ThreeLevelPrecon,P.Outer.Mmax,P.Outer.Nstep)
+  {
+    GRID_ASSERT( C.DenseCC != nullptr );
+    
+    CoarseSmootherGCR.Level(2);
+    CoarseSmootherGCR.Name("Csmoother");
+    CoarseSmootherGCR.SetZeroGuess(1);
+
+    L2PGCR.Level(2);
+    L2PGCR.Name("Couter");
+    L2PGCR.SetZeroGuess(1);
+    
+    SmootherGCR.Level(1);
+    SmootherGCR.Name("Fsmoother");
+    SmootherGCR.SetZeroGuess(1);
+
+    L1PGCR.Level(1);
+    L1PGCR.Name("Fouter");
+    L1PGCR.SetZeroGuess(1);
+    
+    // The V-cycle is preconditioner: sloppy halos inside, exact restored on exit.
+    ThreeLevelPrecon.SetSloppy    = [this](int s){ PVdagM.SloppyComms(s); };
+    ThreeLevelPrecon.SloppyComms  = Params.Setup.FineSloppyComms;
+    PVdagM.SloppyComms(0);
+    std::cout << GridLogMessage << "PVdagMMultiGridSolver: Nrhs " << nr
+              << ", fine halo policy: preconditioner+coarsening "
+              << (Params.Setup.FineSloppyComms ? "SLOPPY (fp32 wire)" : "exact")
+              << ", outer Krylov EXACT" << std::endl;
+  }
+
+  void Solve(std::vector<FineField> &src, std::vector<FineField> &sol)
+  {
+    GRID_ASSERT( (int)src.size() == nrhs );
+    // Start the solve from a clean device: evict setup-era Lattice copies
+    // and release the allocation caches' held blocks, so the LRU cap
+    // applies to the solve's own working set (OOM policy, Frontier).
+    MemoryManager::EvictAll();
+    MemoryManager::DropCache();
+
+    GridStopWatch w; w.Start();
+    L1PGCR(src,sol);
+    w.Stop();
+    std::cout << GridLogMessage << "PVdagMMultiGridSolver: Nrhs "<<nrhs<<" total " << w.Elapsed()
+              << "  (per RHS: " << w.useconds()/1.0e6/nrhs << " 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.
+    PVdagM.SloppyComms(0);
+    {
+      FineField Ax(src[0].Grid());
+      RealD worst=0.0;
+      for(int r=0;r<nrhs;r++){
+	PVdagM.Op(sol[r],Ax);
+	Ax=Ax-src[r];
+        RealD rn=std::sqrt(norm2(Ax)/norm2(src[r]));
+        std::cout << GridLogMessage << "FINAL Nrhs "<<nrhs<<": rhs["<<r<<"] true residual = " << rn << std::endl;
+        worst=std::max(worst,rn);
+      }
+      std::cout << GridLogMessage << "FINAL Nrhs "<<nrhs<<": worst-case residual = " << worst
+                << "   (exact-halo verification)" << std::endl;
+    }
+  }
+};
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/PVdagMMultiGridParams.h

@@ -0,0 +1,155 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/PVdagMMultiGridParams.h
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+//////////////////////////////////////////////////////////////////////////////////////
+// Serialisable parameters for the mrhs PVdagM multigrid solver chain.
+//
+// Named for the multigrid FORM it configures, as WilsonMGParams is for the
+// Wilson MG tests; an MrhsHDCGParams will compose the same generic MG*
+// sub-structs for the Hermitian chain.  Read from XML (or JSON -- the
+// serialisation macros give both) via ReadPVdagMMultiGridParams below and
+// printed at startup by the macro's operator<<, so every log identifies its
+// own run.  The constructor defaults ARE a tuned operating point, not
+// arbitrary: an unconfigured run reproduces it, so changing them changes what
+// an unconfigured run does.  Smoother mmax == nstep (full GCR history).
+//
+// There are NO environment-variable controls anywhere in this subsystem:
+// parameters live here, library-internal constants are hard defaults in their
+// classes, verbosity is the --log channels.
+// Research instruments (GCR coefficient recording/replay, Chebyshev and
+// stationary smoother variants, power iteration) are deliberately NOT part
+// of this interface: they remain programmatic, for algorithmic studies,
+// not consumer knobs.
+//////////////////////////////////////////////////////////////////////////////////////
+
+// The smoother is the adaptive shifted PGCR -- the one correct route for
+// this non-Hermitian chain (stationary replay and Chebyshev were explored
+// and did not win here; Chebyshev remains the HERMITIAN chain's smoother).
+struct MGSmootherParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(MGSmootherParams,
+                                  RealD, Shift,
+                                  int,   Nstep,
+                                  int,   Mmax);
+  MGSmootherParams(RealD shift=0.1, int nstep=6, int mmax=6)
+    : Shift(shift), Nstep(nstep), Mmax(mmax) {};
+};
+
+struct MGCoarseSolverParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(MGCoarseSolverParams,
+                                  RealD, Tol,
+                                  int,   Order,
+                                  int,   Mmax);
+  MGCoarseSolverParams() : Tol(0.05), Order(200), Mmax(16) {};
+};
+
+struct MGOuterParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(MGOuterParams,
+                                  RealD, Tol,
+                                  int,   MaxIterations,
+                                  int,   Mmax,
+                                  int,   Nstep);
+  MGOuterParams() : Tol(1.0e-8), MaxIterations(1000), Mmax(4), Nstep(8) {};
+};
+
+struct MGDenseParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(MGDenseParams,
+                                  int, LeafSpan);   // big-leaf span in blocks; see BlockCyclicSchurInverse
+  MGDenseParams() : LeafSpan(9) {};
+};
+
+struct MGSetupParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(MGSetupParams,
+                                  std::vector<int>, Block,           // fine -> coarse blocking
+                                  std::vector<int>, Block2,          // coarse -> coarse-coarse blocking
+                                  int,              CoarsenBatch,
+                                  std::string,      SubspaceFile,    // scidac; empty = create from noise, no I/O
+                                  int,              FineSloppyComms);// fp32 wire INSIDE the preconditioner only
+  MGSetupParams()
+    : Block({2,2,3,3}), Block2({4,4,2,4}), CoarsenBatch(9),
+      SubspaceFile(""), FineSloppyComms(1) {};
+};
+
+struct PVdagMMultiGridParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(PVdagMMultiGridParams,
+                                  MGSetupParams,        Setup,
+                                  MGSmootherParams,     FineSmoother,
+                                  MGSmootherParams,     CoarseSmoother,
+                                  MGCoarseSolverParams, CoarseSolver,
+                                  MGOuterParams,        Outer,
+                                  MGDenseParams,        Dense);
+  PVdagMMultiGridParams()
+    : FineSmoother  (0.1, 6, 6),
+      CoarseSmoother(0.1, 2, 2) {};
+};
+
+inline void CheckValidity(const PVdagMMultiGridParams &P)
+{
+  GRID_ASSERT( P.Setup.Block.size()  == 4 );
+  GRID_ASSERT( P.Setup.Block2.size() == 4 );
+  GRID_ASSERT( P.Setup.CoarsenBatch  >= 1 );
+  GRID_ASSERT( P.FineSmoother.Nstep   > 0 );
+  GRID_ASSERT( P.CoarseSmoother.Nstep > 0 );
+  GRID_ASSERT( P.CoarseSolver.Tol     > 0.0 );
+  GRID_ASSERT( P.Outer.Tol            > 0.0 );
+  GRID_ASSERT( P.Dense.LeafSpan      >= 1 );
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+// Read from an XML file; on a MISSING file the boss writes <file>.templ from
+// the defaults and returns false so the caller can Grid_finalize() and exit
+// -- the produced template is the documentation of every parameter.  An
+// empty filename means "defaults": no I/O, returns true.  Always prints the
+// effective parameters, so the log describes the run.
+//
+// Drivers fetch the filename themselves, e.g.
+//   std::string pfile("");
+//   if( GridCmdOptionExists(argv,argv+argc,"--multigrid-params") )
+//     pfile = GridCmdOptionPayload(argv,argv+argc,"--multigrid-params");
+//   if( !ReadPVdagMMultiGridParams(Params, pfile) ) { Grid_finalize(); return 0; }
+//////////////////////////////////////////////////////////////////////////////////////
+inline bool ReadPVdagMMultiGridParams(PVdagMMultiGridParams &P, const std::string &xmlfile)
+{
+  if ( xmlfile.length() ) {
+    bool good;
+    { std::ifstream f(xmlfile); good = f.good(); }
+    if ( !good ) {
+      std::cout << GridLogMessage << "PVdagMMultiGridParams: " << xmlfile << " does not exist" << std::endl;
+      if ( GlobalSharedMemory::WorldRank == 0 ) {
+        XmlWriter WR(xmlfile+".templ");
+        write(WR, "PVdagMMultiGridParams", P);
+        std::cout << GridLogMessage << "PVdagMMultiGridParams: template written to "
+                  << xmlfile << ".templ" << std::endl;
+      }
+      return false;
+    }
+    XmlReader RD(xmlfile);
+    read(RD, "PVdagMMultiGridParams", P);
+  }
+  CheckValidity(P);
+  std::cout << GridLogMessage << "PVdagMMultiGridParams ("
+            << (xmlfile.length() ? xmlfile : std::string("defaults")) << "):" << std::endl;
+  std::cout << P << std::endl;
+  return true;
+}
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/PVdagMOperators.h

@@ -0,0 +1,88 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/PVdagMOperators.h
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+//////////////////////////////////////////////////////////////////////
+// A = PV^dag M (non-Hermitian).  The operator the PVdagM multigrid
+// coarsens and solves.
+//
+// SloppyComms delegates to BOTH underlying fermion operators, so one
+// call controls every fine halo in the chain -- including any shifted
+// wrapper built on the same two matrices.  Halo-precision POLICY:
+// reduced-precision (fp32 wire) halos belong in the PRECONDITIONER --
+// the smoother, the V-cycle's own residuals, the coarsening -- and
+// NEVER in the outer Krylov, whose applications define what
+// "converged" means.  The operators default to EXACT; the V-cycle
+// turns sloppiness on for its own scope only (MrhsTwoLevelMG).
+//////////////////////////////////////////////////////////////////////
+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 SloppyComms(int sloppy) { _Mat.SloppyComms(sloppy); _PV.SloppyComms(sloppy); }
+  void OpDiag (const Field &in, Field &out) { GRID_ASSERT(0); }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) { GRID_ASSERT(0); }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){ GRID_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); }
+};
+
+//////////////////////////////////////////////////////////////////////
+// A + shift, built directly on the two matrices (fine level: avoids
+// an extra fine-field pass through a generic shifted wrapper).
+//////////////////////////////////////////////////////////////////////
+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) { GRID_ASSERT(0); }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) { GRID_ASSERT(0); }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){ GRID_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(in,tmp); _Mat.Mdag(tmp,out); out = out + shift*in; }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); }
+  void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+//////////////////////////////////////////////////////////////////////
+// Op + shift for any operator (coarse levels).
+//////////////////////////////////////////////////////////////////////
+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) { GRID_ASSERT(0); }
+  void OpDir   (const Field &in, Field &out,int dir,int disp) { GRID_ASSERT(0); }
+  void OpDirAll (const Field &in, std::vector<Field> &out) { GRID_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){ GRID_ASSERT(0); }
+  void HermOp  (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/RecursiveSchurInverse.h

@@ -1,874 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: RecursiveSchurInverse.h
-
-    Copyright (C) 2026
-
-Author: Peter Boyle <pboyle@bnl.gov>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#pragma once
-
-#include <Grid/algorithms/blas/BatchedBlas.h>
-#include <Grid/algorithms/blas/BatchedInverse.h>
-
-#include <algorithm>
-
-NAMESPACE_BEGIN(Grid);
-
-///////////////////////////////////////////////////////////////////////////////
-// RecursiveSchurInverse: distributed dense inversion by recursive Schur
-// complement over a binary rank-range tree.
-//
-// Contract: rank r owns global rows [rowStart[r], rowStart[r+1]) of an
-// N x N matrix in rank-major ordering, and receives its rows of the
-// inverse in the same layout.  Consumes GridBase collectives, GridBLAS
-// and GridBLASInverse only; Eigen reference backends permit CPU unit
-// testing under mpirun (Test_schur_inverse).
-//
-// Arithmetic is fp64 throughout; the caller rounds once into fp32 storage.
-//
-// Execution: SPMD full-tree walk.  Every rank makes the identical call
-// sequence; data participation is ownership-gated; all collectives are
-// world-wide, so no deadlock surface exists.
-//
-// Storage: BlockRows is column-major, ld = rows; element (i,j) at
-// data[i + j*ld]; a column window is the contiguous slice at data[col0*ld].
-///////////////////////////////////////////////////////////////////////////////
-
-///////////////////////////////////////////////////////////////////////////////
-// My rows of a distributed dense matrix: rows x cols, column major, ld = rows.
-///////////////////////////////////////////////////////////////////////////////
-class BlockRows
-{
-public:
-  deviceVector<ComplexD> data;
-  int64_t                rows;
-  int64_t                cols;
-  int64_t                ld;
-
-  BlockRows()
-  {
-    rows = 0;
-    cols = 0;
-    ld   = 0;
-  }
-  void Resize(int64_t r, int64_t c)
-  {
-    rows = r;
-    cols = c;
-    ld   = r;
-    data.resize((uint64_t)r*c);
-  }
-  ComplexD *ColumnWindow(int64_t col0)
-  {
-    GRID_ASSERT( col0 >= 0 );
-    GRID_ASSERT( col0 <= cols );
-    return &data[(uint64_t)col0*ld];
-  }
-};
-
-class RecursiveSchurInverse
-{
-public:
-  GridBase             *grid;       // world collectives only
-  int64_t               N;          // global matrix dimension
-  int                   P;          // ranks
-  int                   me;         // this rank
-  std::vector<int64_t>  rowStart;   // P+1 entries: rank-major row ownership
-  int64_t               myRow0;
-  int64_t               myNrows;
-  int64_t               panelBytes; // gather panel budget (DENSE_PANEL_BYTES)
-
-  GridBLAS              BLAS;
-  GridBLASInverse       INV;
-
-  // Growth telemetry: one entry per merge node, walk order
-  std::vector<double>   telNormB;      // ||B||_F = ||A11inv A12||_F
-  std::vector<double>   telSratio;     // ||S||_F / ||A22||_F
-  double                telLeafMaxInv; // max |(leaf inverse)_ij| over leaves
-
-  // Phase timers/counters, reported by ReportTelemetry
-  double                tMemset;       // device panel zero-fill
-  double                tDeposit;      // owner rows -> panel (device kernel)
-  double                tAllreduce;    // panel collective (see tBarrier)
-  double                tBarrier;      // arrival skew, when DENSE_BARRIER_PROBE
-  double                tRepack;       // rank-major -> panel (gather path)
-  double                tGemm;         // strided gemm + synchronise
-  double                tLeaf;         // leaf inversions
-  double                tARmin;        // fastest single panel collective
-  double                tARmax;        // slowest single panel collective
-  std::vector<double>   tARall;        // every panel collective, for percentiles
-  uint64_t              bytesAllreduce;
-  uint64_t              nAllreduce;    // panel collectives
-  uint64_t              nGatherGemm;   // GatherGemm calls
-  uint64_t              nGather;       // gather-path collectives (debug gate)
-  uint64_t              nGatherV;      // ... of which MPI_Allgatherv
-  uint64_t              nBcast;        // ... of which Bcast-assembled
-
-  // Persistent grow-only device panel; assembly and collectives are
-  // device-resident.  Device builds require GPU-aware MPI.
-  deviceVector<ComplexD> dPanelBuf;
-  // Rank-major receive staging for the AllGatherV path, plus the per-panel-row
-  // owner maps that drive the repack.  Grow-only, same discipline as dPanelBuf:
-  // a fresh device allocation per collective is catastrophic (measured).
-  deviceVector<ComplexD> dRecvBuf;
-  deviceVector<int64_t>  dOwnerOff;    // panel row -> owner's first panel row
-  deviceVector<int64_t>  dOwnerRows;   // panel row -> owner's row count
-
-  // DENSE_GATHER            : transport for the panel assembly.
-  //     0 = zero-fill + GlobalSumVector  (default; moves the payload twice)
-  //     1 = MPI_Allgatherv               (KNOWN BROKEN here, see guard below)
-  //     2 = C sequential MPI_Bcast, one per contributing rank.  Bcast has no
-  //         count vector, so the zero-count shape that breaks Allgatherv
-  //         cannot arise.  Costs C collectives instead of 1 and the roots do
-  //         not transmit concurrently, so the byte cost is ~2*panel per rank
-  //         -- the same as the allreduce.  It wins only if Bcast is faster
-  //         PER BYTE than Allreduce, which is why it must be measured.
-  // DENSE_GATHER_MIN_BYTES  : panels below this stay on the allreduce path.
-  //                           Default 0 (always gather).  The Schur tree puts
-  //                           94% of its bytes in panels with ~3.7 MB per-rank
-  //                           messages, i.e. squarely bandwidth bound; only the
-  //                           two deepest levels (1.2% of bytes, 16-65 KB per
-  //                           rank) are latency bound.  This exists so that
-  //                           tail can be excluded with an env var rather than
-  //                           a rebuild, if it ever proves to matter.
-  // DENSE_BARRIER_PROBE     : Barrier() before each panel collective.
-  //   DEFAULTS ON, and it is an optimisation rather than instrumentation.
-  //   Measured at 288 ranks, N=138240: comms 380.3 s -> 314.7 s and the whole
-  //   invert 394.3 s -> 327.7 s, with the worst single collective falling
-  //   2474.9 ms -> 256.2 ms and effective rate 2.40 -> 4.27 GB/s.  Convoy
-  //   accounting alone predicts NO saving (barrier and collective both
-  //   complete at max-over-ranks), so the mechanism is that the collective
-  //   itself degrades under skewed arrival, not that skew is rebooked.
-  //   Set DENSE_BARRIER_PROBE=0 to recover the old behaviour.  It also
-  //   separates arrival skew (tBarrier) from transfer (tAllreduce).
-  int                    useGather;
-  int64_t                gatherMinBytes;
-  int                    barrierProbe;
-  // DENSE_GATHER_DEBUG=N : dump the descriptor of the first N AllGatherV
-  // calls from rank 0.  The tiling invariant below is checked ALWAYS -- it
-  // costs one O(P) host loop against a collective, and a counts array that
-  // does not tile the panel is exactly the failure we are hunting.
-  int                    gatherDebug;
-
-  ///////////////////////////////////////////////////////////////////////////
-  // Ownership-table validation: a proper partition of [0,N).
-  // Static and communicator-free so synthetic tables unit-test directly.
-  ///////////////////////////////////////////////////////////////////////////
-  static void CheckRowStart(const std::vector<int64_t> &table, int64_t N)
-  {
-    int P = (int)table.size() - 1;
-    GRID_ASSERT( P >= 1 );
-    GRID_ASSERT( table[0] == 0 );
-    GRID_ASSERT( table[P] == N );
-    for(int r=0; r<P; r++)
-    {
-      GRID_ASSERT( table[r+1] >= table[r] );  // zero-row ranks permitted
-    }
-  }
-
-  ///////////////////////////////////////////////////////////////////////////
-  // Build the ownership table from each rank's local row count: zero-fill
-  // allgather (the standing comms idiom) then prefix sum.  Every rank
-  // returns the identical table.
-  ///////////////////////////////////////////////////////////////////////////
-  static std::vector<int64_t> MakeRowStart(GridBase *g, int64_t myNrows)
-  {
-    int P  = g->ProcessorCount();
-    int me = g->ThisRank();
-
-    std::vector<uint64_t> counts(P, 0);
-    counts[me] = (uint64_t)myNrows;
-    g->GlobalSumVector(&counts[0], P);
-
-    std::vector<int64_t> table(P+1);
-    table[0] = 0;
-    for(int r=0; r<P; r++)
-    {
-      table[r+1] = table[r] + (int64_t)counts[r];
-    }
-    CheckRowStart(table, table[P]);
-    return table;
-  }
-
-  RecursiveSchurInverse(GridBase *g,
-                        int64_t N_,
-                        std::vector<int64_t> &rowStart_,
-                        int64_t panelBytes_)
-  {
-    grid       = g;
-    N          = N_;
-    P          = g->ProcessorCount();
-    me         = g->ThisRank();
-    rowStart   = rowStart_;
-    panelBytes = panelBytes_;
-
-    GRID_ASSERT( (int)rowStart.size() == P+1 );
-    CheckRowStart(rowStart, N);
-
-    myRow0  = rowStart[me];
-    myNrows = rowStart[me+1] - rowStart[me];
-
-    telLeafMaxInv = 0.0;
-
-    useGather      = getenv("DENSE_GATHER") ? atoi(getenv("DENSE_GATHER")) : 0;
-    gatherMinBytes = getenv("DENSE_GATHER_MIN_BYTES")
-                   ? atol(getenv("DENSE_GATHER_MIN_BYTES")) : 0;
-    barrierProbe   = getenv("DENSE_BARRIER_PROBE") ? atoi(getenv("DENSE_BARRIER_PROBE")) : 1;
-    gatherDebug    = getenv("DENSE_GATHER_DEBUG") ? atoi(getenv("DENSE_GATHER_DEBUG")) : 0;
-
-    ///////////////////////////////////////////////////////////////////////
-    // DENSE_GATHER=1 is KNOWN BROKEN on Cray MPICH and refuses to start.
-    //
-    // The gather asks MPI_Allgatherv to assemble a panel to which only the
-    // rank sub-range [rB0,rB1) contributes; every other rank passes count 0.
-    // At Schur depth d that is P/2^(d+1) contributors out of P, i.e. 9 of 288
-    // at depth 4 and ~1 of 288 at depth 7.  Measured consequences at 288
-    // ranks (tests/debug/Test_allgather):
-    //   T5  288/288 contributing, 18.4 MB on device : sub-second.
-    //   T6  144/288 contributing,  9.2 MB on device : ~53 s.
-    // and in production, 9/288 contributing at 299 MB hangs and then aborts
-    // inside PMPI_Allgatherv ("req != NULL", mpir_request.h:508).
-    //
-    // The primitive itself is sound -- Test_allgather T1-T6 pass, and T3/T4
-    // verify it BITWISE against the zero-fill+GlobalSum idiom.  What fails is
-    // MPI's handling of a collective in which almost every rank contributes
-    // nothing.  Fixing it needs either P2P (a bisection allgather over
-    // [r0,r1)) or the 2D block-cyclic layout, which removes the shape
-    // altogether by giving every rank part of every sub-block.
-    //
-    // Fail here, in the first second, rather than 100 s and 36 nodes into a
-    // job.  DENSE_GATHER_FORCE=1 proceeds anyway for debugging.
-    ///////////////////////////////////////////////////////////////////////
-    if ( (useGather==1) && !getenv("DENSE_GATHER_FORCE") )
-    {
-      std::cout << GridLogError
-                << "DENSE_GATHER=1 is disabled: MPI_Allgatherv on this MPI is"
-                << " pathological when most ranks contribute count 0 (see the"
-                << " note at RecursiveSchurInverse.h, and Test_allgather T5 vs"
-                << " T6).  Unset DENSE_GATHER, or set DENSE_GATHER_FORCE=1 to"
-                << " proceed anyway." << std::endl;
-      GRID_ASSERT(0 && "DENSE_GATHER=1 known broken: see comment above");
-    }
-  }
-
-  ///////////////////////////////////////////////////////////////////////////
-  // The communication primitive.
-  //
-  //   C(:, colC : colC+widthB) <-   beta  * C(:, colC : colC+widthB)
-  //                               + alpha * A(:, colA : colA+widthA) * Bsub
-  //
-  // Bsub is the widthA x widthB sub-block of a row-distributed operand
-  // owned by ranks [rB0, rB1): owner r contributes its rows of
-  // B(:, colB : colB+widthB) at sub-block row offset
-  // rowStart[r] - rowStart[rB0], gathered in panelBytes row-chunks by
-  // device zero-fill + deposit kernel + GlobalSumVector.
-  //
-  // Every rank calls; non-owners of B add zeros; ranks with A.rows == 0
-  // skip local compute but make every collective call.  Column offsets
-  // are local buffer offsets -- non-participants pass 0.
-  ///////////////////////////////////////////////////////////////////////////
-  void GatherGemm(ComplexD alpha,
-                  BlockRows &A, int64_t colA, int64_t widthA,
-                  int rB0, int rB1,
-                  BlockRows &B, int64_t colB, int64_t widthB,
-                  ComplexD beta,
-                  BlockRows &C, int64_t colC)
-  {
-    GRID_ASSERT( rB0 >= 0 );
-    GRID_ASSERT( rB1 >  rB0 );
-    GRID_ASSERT( rB1 <= P );
-
-    int64_t k = rowStart[rB1] - rowStart[rB0];
-    int64_t m = A.rows;
-    int64_t n = widthB;
-    GRID_ASSERT( widthA == k );
-    GRID_ASSERT( n >= 1 );
-
-    int owner = ( me >= rB0 ) && ( me < rB1 ) && ( B.rows > 0 );
-    int64_t myOff = 0;
-    if ( owner )
-    {
-      myOff = rowStart[me] - rowStart[rB0];
-    }
-
-    if ( m > 0 )
-    {
-      GRID_ASSERT( colA + widthA <= A.cols );
-      GRID_ASSERT( colC + widthB <= C.cols );
-      GRID_ASSERT( C.rows == m );
-    }
-
-    nGatherGemm++;
-    GRID_TRACE("GatherGemm");
-
-    if ( owner )
-    {
-      GRID_ASSERT( colB + widthB <= B.cols );
-    }
-
-    // COLUMN chunking, not row chunking.  Each rank's contribution to a
-    // column chunk is rows_r x nchunk at ld = B.rows -- i.e. exactly the
-    // contiguous window B.ColumnWindow(colB+j0) -- so the AllGatherV send
-    // needs no pack.  Row chunking would slice each column and force one.
-    // It also makes every chunk a full-k product, so beta applies directly
-    // and the cross-chunk accumulation disappears.
-    int64_t bufs = useGather ? 2 : 1;         // panel, plus receive staging
-    int64_t nc = panelBytes / ( bufs * (int64_t)sizeof(ComplexD) * k );
-    if ( nc < 1 ) nc = 1;
-    if ( nc > n ) nc = n;
-    GRID_ASSERT( k*nc < 2147483647L );        // collective counts are int
-
-    deviceVector<ComplexD> &dPanel = dPanelBuf;
-    if ( dPanel.size() < (uint64_t)k*nc ) dPanel.resize((uint64_t)k*nc);
-    deviceVector<ComplexD*> ap(1);
-    deviceVector<ComplexD*> bp(1);
-    deviceVector<ComplexD*> cp(1);
-    std::vector<ComplexD*>  ptr(1);
-
-    // Panel-row -> owner maps for the repack.  Depend only on [rB0,rB1),
-    // so they are built once per call rather than once per chunk.
-    if ( useGather )
-    {
-      if ( dRecvBuf.size()   < (uint64_t)k*nc ) dRecvBuf.resize((uint64_t)k*nc);
-      if ( dOwnerOff.size()  < (uint64_t)k )    dOwnerOff.resize((uint64_t)k);
-      if ( dOwnerRows.size() < (uint64_t)k )    dOwnerRows.resize((uint64_t)k);
-      std::vector<int64_t> hoff(k), hrows(k);
-      for(int r=rB0; r<rB1; r++)
-      {
-        int64_t off_r  = rowStart[r]   - rowStart[rB0];
-        int64_t rows_r = rowStart[r+1] - rowStart[r];
-        for(int64_t i=0;i<rows_r;i++){ hoff[off_r+i]=off_r; hrows[off_r+i]=rows_r; }
-      }
-      acceleratorCopyToDevice(&hoff[0], &dOwnerOff[0], (uint64_t)k*sizeof(int64_t));
-      acceleratorCopyToDevice(&hrows[0],&dOwnerRows[0],(uint64_t)k*sizeof(int64_t));
-      if ( owner ) GRID_ASSERT( B.rows == rowStart[me+1]-rowStart[me] );
-    }
-
-    for(int64_t j0=0; j0<n; j0+=nc)
-    {
-      int64_t nchunk = std::min(nc, n-j0);
-      uint64_t panelWords = (uint64_t)k*nchunk;
-      uint64_t panelBytesThis = panelWords*sizeof(ComplexD);
-
-      // The MODE (0/1/2), not a boolean: `useGather && cond` collapses 2 to 1
-      // and silently dispatched DENSE_GATHER=2 onto the known-broken
-      // AllGatherV path (Frontier hang, 2026-08-24).  Identical on all ranks.
-      int gatherThis = ( (int64_t)panelBytesThis >= gatherMinBytes ) ? useGather : 0;
-
-      // Arrival skew is charged to the barrier, so that tAllreduce measures
-      // transfer alone.  Diagnostic only; off by default.
-      if ( barrierProbe ) { double tb = -usecond(); grid->Barrier(); tBarrier += tb+usecond(); }
-
-      double tar = -usecond();
-      if ( gatherThis )
-      {
-        //////////////////////////////////////////////////////////////////
-        // AllGatherV: move the payload once, no arithmetic, no zero fill.
-        // Receive is rank major -- a concatenation of rows_r x nchunk
-        // column-major blocks -- which is the correct ROW order but the
-        // wrong LAYOUT, so one device repack follows.
-        //////////////////////////////////////////////////////////////////
-        std::vector<int> counts(P,0), displs(P,0);
-        for(int r=rB0; r<rB1; r++)
-        {
-          counts[r] = (int)((rowStart[r+1]-rowStart[r])*nchunk);
-          displs[r] = (int)((rowStart[r]  -rowStart[rB0])*nchunk);
-        }
-        //////////////////////////////////////////////////////////////////
-        // The counts MUST tile the panel exactly: every one of the
-        // k*nchunk words has exactly one owner, or MPI is handed an
-        // inconsistent descriptor.  Always checked; O(P) against a
-        // collective is free.
-        //////////////////////////////////////////////////////////////////
-        int64_t csum = 0; int nz = 0;
-        for(int r=0;r<P;r++){ csum += counts[r]; if ( counts[r] ) nz++; }
-        GRID_ASSERT( csum == (int64_t)panelWords );
-        GRID_ASSERT( displs[rB1-1] + counts[rB1-1] == (int)panelWords );
-        // Gate on the count of GATHERS, not of all collectives.  nAllreduce
-        // is already in the hundreds by the first gather -- the recursion is
-        // depth first and everything below the size threshold takes the
-        // allreduce path -- so gating on it prints nothing at any sane value.
-        nGather++;
-        if ( gatherDebug && ((int)nGather <= gatherDebug) && (me==0) ) {
-          std::cout << GridLogMessage << "GATHER["<<nGather-1<<"]"
-                    << "  ranks ["<<rB0<<","<<rB1<<")"
-                    << "  k "<<k<<"  n "<<n<<"  nchunk "<<nchunk
-                    << "  panelWords "<<panelWords
-                    << " ("<<panelBytesThis/1024/1024<<" MB)"
-                    << "  contributors "<<nz<<"/"<<P
-                    << "  maxcount "<<*std::max_element(counts.begin(),counts.end())
-                    << "  rank0: count "<<counts[me]<<(owner?" owner":" non-owner")
-                    << std::endl;
-        }
-        if ( gatherThis == 1 )
-        {
-          nGatherV++;
-          void *send = owner ? (void *)B.ColumnWindow(colB+j0) : (void *)&dRecvBuf[0];
-          grid->AllGatherV(send, counts[me],
-                           (void *)&dRecvBuf[0], counts, displs, sizeof(ComplexD));
-        }
-        else
-        {
-          ////////////////////////////////////////////////////////////////
-          // C sequential broadcasts into the SAME rank-major staging
-          // buffer, so the repack below is shared with the AllGatherV arm.
-          // The root must already hold its own block: B.ColumnWindow is
-          // contiguous rows_r x nchunk, and so is its slot in dRecvBuf.
-          // Every rank issues all C broadcasts in the same order, so
-          // collective matching stays positional and safe.
-          ////////////////////////////////////////////////////////////////
-          nBcast++;
-          if ( owner )
-          {
-            int64_t off_me = rowStart[me] - rowStart[rB0];
-            acceleratorCopyDeviceToDevice((void *)B.ColumnWindow(colB+j0),
-                                          (void *)&dRecvBuf[off_me*nchunk],
-                                          (uint64_t)B.rows*nchunk*sizeof(ComplexD));
-          }
-          for(int r=rB0;r<rB1;r++)
-          {
-            int64_t off_r  = rowStart[r]   - rowStart[rB0];
-            int64_t rows_r = rowStart[r+1] - rowStart[r];
-            if ( rows_r )
-              grid->Broadcast(r,(void *)&dRecvBuf[off_r*nchunk],
-                              (uint64_t)rows_r*nchunk*sizeof(ComplexD));
-          }
-        }
-        tar += usecond();
-
-        tRepack -= usecond();
-        {
-          ComplexD *pan = &dPanel[0];
-          ComplexD *rcv = &dRecvBuf[0];
-          int64_t  *ooff  = &dOwnerOff[0];
-          int64_t  *orows = &dOwnerRows[0];
-          int64_t   kk = k, ncw = nchunk;
-          accelerator_for(idx, panelWords, 1, {
-            int64_t j = idx / kk;
-            int64_t p = idx - j*kk;
-            int64_t o = ooff[p];
-            int64_t rr= orows[p];
-            pan[idx] = rcv[(uint64_t)(o*ncw + (p-o) + j*rr)];
-          });
-        }
-        tRepack += usecond();
-      }
-      else
-      {
-        //////////////////////////////////////////////////////////////////
-        // Zero fill + deposit + GlobalSumVector.  Exact because the zero
-        // fill leaves exactly one contributing rank per element, but it
-        // moves the payload twice and reduces over zeros.
-        //////////////////////////////////////////////////////////////////
-        tar += usecond();
-        tMemset -= usecond();
-        acceleratorMemSet(&dPanel[0], 0, panelBytesThis);
-        tMemset += usecond();
-        if ( owner )
-        {
-          int64_t   brows = B.rows;
-          int64_t   kk    = k;
-          int64_t   dof   = myOff;
-          ComplexD *src   = B.ColumnWindow(colB+j0);
-          ComplexD *dst   = &dPanel[0];
-          tDeposit -= usecond();
-          accelerator_for(idx, (uint64_t)(brows*nchunk), 1, {
-            int64_t j = idx / brows;
-            int64_t i = idx - j*brows;
-            dst[(uint64_t)(dof + i + j*kk)] = src[(uint64_t)(i + j*brows)];
-          });
-          tDeposit += usecond();
-        }
-        tar = -usecond();
-        grid->GlobalSumVector(&dPanel[0], (int)panelWords);
-        tar += usecond();
-      }
-      tAllreduce += tar;
-      tARmin = std::min(tARmin, tar);
-      tARmax = std::max(tARmax, tar);
-      tARall.push_back(tar);
-      bytesAllreduce += panelBytesThis;
-      nAllreduce++;
-
-      if ( m > 0 )
-      {
-        // Full-k product into this column chunk of C: beta applies directly.
-        ptr[0] = A.ColumnWindow(colA);
-        acceleratorCopyToDevice(&ptr[0], &ap[0], sizeof(ComplexD*));
-        ptr[0] = &dPanel[0];
-        acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexD*));
-        ptr[0] = C.ColumnWindow(colC+j0);
-        acceleratorCopyToDevice(&ptr[0], &cp[0], sizeof(ComplexD*));
-
-        tGemm -= usecond();
-        BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N,
-                         (int)m, (int)nchunk, (int)k,
-                         alpha,    ap, (int)A.ld,
-                                   bp, (int)k,
-                         beta,     cp, (int)C.ld);
-        BLAS.synchronise();
-        tGemm += usecond();
-      }
-    }
-  }
-
-  ///////////////////////////////////////////////////////////////////////////
-  // Local Frobenius norm-squared of a full-height column window.
-  // NO comms; callers GlobalSum the result.  Host staging, setup-scale.
-  ///////////////////////////////////////////////////////////////////////////
-  double FrobNorm2Local(BlockRows &X, int64_t col0, int64_t w)
-  {
-    if ( X.rows == 0 ) return 0.0;
-    GRID_ASSERT( col0 + w <= X.cols );
-    uint64_t len = (uint64_t)X.rows*w;
-    std::vector<ComplexD> h(len);
-    acceleratorCopyFromDevice(X.ColumnWindow(col0), &h[0], len*sizeof(ComplexD));
-    // Member real()/imag(): portable across std::complex (CPU) and
-    // thrust::complex (HIP), where std::norm does not resolve.
-    double s = 0.0;
-    for(uint64_t i=0; i<len; i++)
-    {
-      double re = h[i].real();
-      double im = h[i].imag();
-      s += re*re + im*im;
-    }
-    return s;
-  }
-
-  ///////////////////////////////////////////////////////////////////////////
-  // dst(:, dcol0 : dcol0+w) = - src(:, 0:w).  Both operands have ld == rows
-  // so full-height windows are contiguous: flat elementwise device copy.
-  ///////////////////////////////////////////////////////////////////////////
-  void NegateCopy(BlockRows &src, BlockRows &dst, int64_t dcol0, int64_t w)
-  {
-    GRID_ASSERT( src.rows == dst.rows );
-    GRID_ASSERT( w <= src.cols );
-    GRID_ASSERT( dcol0 + w <= dst.cols );
-    if ( src.rows == 0 ) return;
-    uint64_t len = (uint64_t)src.rows*w;
-    ComplexD *s = &src.data[0];
-    ComplexD *d = dst.ColumnWindow(dcol0);
-    accelerator_for(i, len, 1, {
-      d[i] = -s[i];
-    });
-  }
-
-  ///////////////////////////////////////////////////////////////////////////
-  // Leaf inversion: local, in place on the contiguous diagonal window.
-  // No collectives.
-  ///////////////////////////////////////////////////////////////////////////
-  void LeafInvert(int64_t col0, int64_t width, BlockRows &Arows)
-  {
-    GRID_TRACE("SchurLeaf");
-    GRID_ASSERT( width == Arows.rows );
-    GRID_ASSERT( col0 + width <= Arows.cols );
-    int64_t  w   = width;
-    uint64_t len = (uint64_t)w*w;
-    tLeaf -= usecond();
-
-    deviceVector<ComplexD*> bp(1);
-    std::vector<ComplexD*>  ptr(1);
-    ptr[0] = Arows.ColumnWindow(col0);
-    acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexD*));
-    INV.inverseBatched(w, bp);
-
-    // Telemetry: max |element| of the leaf inverse
-    {
-      std::vector<ComplexD> h(len);
-      acceleratorCopyFromDevice(Arows.ColumnWindow(col0), &h[0], len*sizeof(ComplexD));
-      double mx = 0.0;
-      for(uint64_t i=0; i<len; i++)
-      {
-        double re = h[i].real();
-        double im = h[i].imag();
-        mx = std::max(mx, re*re + im*im);
-      }
-      telLeafMaxInv = std::max(telLeafMaxInv, std::sqrt(mx));
-    }
-    tLeaf += usecond();
-  }
-
-  ///////////////////////////////////////////////////////////////////////////
-  // The recursion (plan 3.5 / 4B.3).  Inverts the diagonal block of the
-  // rank-major matrix spanned by ranks [r0, r1), living in every member
-  // rank's column window [col0, col0+width) -- IN PLACE.
-  //
-  // SPMD: every rank calls with IDENTICAL (r0, r1, width) and its own
-  // local (col0, Arows); ranks outside [r0, r1) participate in the
-  // collectives only (dummy operands, zero contributions).  The collective
-  // sequence -- 5 GatherGemm calls + 3 scalar GlobalSums per merge node --
-  // is identical on every rank by construction.
-  //
-  //   I = [r0, mid)   J = [mid, r1)   widths WI, WJ
-  //   1. recurse I:      A11 -> A11inv
-  //   2. B = A11inv.A12                     (I rows)
-  //   3. C = A21.A11inv                     (J rows)
-  //   4. S = A22 - A21.B    in place        (J rows)   [alpha=-1, beta=1]
-  //   5. recurse J:      S -> Sinv
-  //   6. T = Sinv.C                         (J rows)
-  //   7. U = B.Sinv                         (I rows)
-  //   8. X11 = A11inv + U.C  in place       (I rows)   [beta=1]
-  //   9. X12 = -U, X21 = -T  local negates; X22 = Sinv already in place
-  ///////////////////////////////////////////////////////////////////////////
-  void SchurNode(int r0, int r1, int64_t col0, int64_t width, BlockRows &Arows)
-  {
-    int span = r1 - r0;
-    GRID_ASSERT( span >= 1 );
-    GRID_ASSERT( width == rowStart[r1] - rowStart[r0] );
-
-    if ( span == 1 )
-    {
-      if ( ( me == r0 ) && ( myNrows > 0 ) )
-      {
-        LeafInvert(col0, width, Arows);
-      }
-      return;
-    }
-
-    int     mid = ( r0 + r1 ) / 2;
-    int64_t WI  = rowStart[mid] - rowStart[r0];
-    int64_t WJ  = rowStart[r1]  - rowStart[mid];
-
-    // Zero-width child ranges (all ranks of a half owning no rows) are a
-    // KNOWN LIMITATION: fail loudly rather than divide mysteriously.
-    GRID_ASSERT( WI > 0 );
-    GRID_ASSERT( WJ > 0 );
-
-    int inI = ( me >= r0  ) && ( me < mid );
-    int inJ = ( me >= mid ) && ( me < r1  );
-
-    ComplexD one  ( 1.0,0.0);
-    ComplexD mone (-1.0,0.0);
-    ComplexD zero ( 0.0,0.0);
-
-    BlockRows dummy;
-
-    // 1. A11 -> A11inv
-    SchurNode(r0, mid, col0, WI, Arows);
-
-    // 2. B = A11inv . A12   (I rows; gather A12 from I owners)
-    BlockRows Bbuf;
-    if ( inI ) Bbuf.Resize(myNrows, WJ);
-    {
-      BlockRows &Aop = inI ? Arows : dummy;
-      BlockRows &Cop = inI ? Bbuf  : dummy;
-      int64_t    cA  = inI ? col0  : 0;
-      GatherGemm(one, Aop, cA, WI,
-                 r0, mid,
-                 Arows, col0+WI, WJ,
-                 zero, Cop, 0);
-    }
-    double nB = FrobNorm2Local(Bbuf, 0, inI ? WJ : 0);
-    grid->GlobalSumVector(&nB, 1);
-    telNormB.push_back(std::sqrt(nB));
-
-    // 3. C = A21 . A11inv   (J rows; gather A11inv from I owners)
-    BlockRows Cbuf;
-    if ( inJ ) Cbuf.Resize(myNrows, WI);
-    {
-      BlockRows &Aop = inJ ? Arows : dummy;
-      BlockRows &Cop = inJ ? Cbuf  : dummy;
-      int64_t    cA  = inJ ? col0  : 0;
-      GatherGemm(one, Aop, cA, WI,
-                 r0, mid,
-                 Arows, col0, WI,
-                 zero, Cop, 0);
-    }
-
-    // 4. S = A22 - A21 . B  in place on my A22 window (J rows)
-    double nA22 = FrobNorm2Local( inJ ? Arows : dummy, inJ ? col0+WI : 0, inJ ? WJ : 0 );
-    grid->GlobalSumVector(&nA22, 1);
-    {
-      BlockRows &Aop = inJ ? Arows   : dummy;
-      BlockRows &Cop = inJ ? Arows   : dummy;
-      int64_t    cA  = inJ ? col0    : 0;
-      int64_t    cC  = inJ ? col0+WI : 0;
-      GatherGemm(mone, Aop, cA, WI,
-                 r0, mid,
-                 Bbuf, 0, WJ,
-                 one, Cop, cC);
-    }
-    double nS = FrobNorm2Local( inJ ? Arows : dummy, inJ ? col0+WI : 0, inJ ? WJ : 0 );
-    grid->GlobalSumVector(&nS, 1);
-    telSratio.push_back( std::sqrt(nS) / ( std::sqrt(nA22) + 1.0e-300 ) );
-
-    // 5. S -> Sinv
-    SchurNode(mid, r1, col0+WI, WJ, Arows);
-
-    // 6. T = Sinv . C       (J rows; gather C from J owners)
-    BlockRows Tbuf;
-    if ( inJ ) Tbuf.Resize(myNrows, WI);
-    {
-      BlockRows &Aop = inJ ? Arows   : dummy;
-      BlockRows &Cop = inJ ? Tbuf    : dummy;
-      int64_t    cA  = inJ ? col0+WI : 0;
-      GatherGemm(one, Aop, cA, WJ,
-                 mid, r1,
-                 Cbuf, 0, WI,
-                 zero, Cop, 0);
-    }
-
-    // 7. U = B . Sinv       (I rows; gather Sinv from J owners)
-    BlockRows Ubuf;
-    if ( inI ) Ubuf.Resize(myNrows, WJ);
-    {
-      BlockRows &Aop = inI ? Bbuf : dummy;
-      BlockRows &Cop = inI ? Ubuf : dummy;
-      GatherGemm(one, Aop, 0, WJ,
-                 mid, r1,
-                 Arows, col0+WI, WJ,
-                 zero, Cop, 0);
-    }
-
-    // 8. X11 = A11inv + U . C  in place (I rows; gather C from J owners)
-    {
-      BlockRows &Aop = inI ? Ubuf  : dummy;
-      BlockRows &Cop = inI ? Arows : dummy;
-      int64_t    cC  = inI ? col0  : 0;
-      GatherGemm(one, Aop, 0, WJ,
-                 mid, r1,
-                 Cbuf, 0, WI,
-                 one, Cop, cC);
-    }
-
-    // 9. Off-diagonal signs, local
-    if ( inI ) NegateCopy(Ubuf, Arows, col0+WI, WJ);
-    if ( inJ ) NegateCopy(Tbuf, Arows, col0,    WI);
-  }
-
-  ///////////////////////////////////////////////////////////////////////////
-  // PUBLIC ENTRY.  Arows: my rows of the rank-major N x N matrix (fp64).
-  // On exit Arows holds my rows of the inverse, still fp64; the caller
-  // owns the single terminal rounding into its fp32 apply storage.
-  ///////////////////////////////////////////////////////////////////////////
-  void Invert(BlockRows &Arows)
-  {
-    GRID_ASSERT( Arows.rows == myNrows );
-    GRID_ASSERT( Arows.cols == N );
-
-    telNormB.resize(0);
-    telSratio.resize(0);
-    telLeafMaxInv = 0.0;
-
-    tMemset        = 0.0;
-    tDeposit       = 0.0;
-    tAllreduce     = 0.0;
-    tBarrier       = 0.0;
-    tRepack        = 0.0;
-    tGemm          = 0.0;
-    tLeaf          = 0.0;
-    tARmin         = 1.0e30;
-    tARmax         = 0.0;
-    tARall.clear();
-    bytesAllreduce = 0;
-    nAllreduce     = 0;
-    nGatherGemm    = 0;
-    nGather        = 0;
-    nGatherV       = 0;
-    nBcast         = 0;
-
-    SchurNode(0, P, 0, N, Arows);
-
-    RealD mx = telLeafMaxInv;
-    grid->GlobalMax(mx);
-    telLeafMaxInv = mx;
-  }
-
-  // All telemetry values are globally reduced or boss-local; safe to
-  // stream on every rank (Grid quiesces stdout to the boss unless
-  // --debug-stdout).  NOTE: tAllreduce INCLUDES wait/imbalance -- a rank
-  // arriving early books its wait here; the min/max spread across ranks
-  // separates true wire time (~min) from skew (max-min).
-  void ReportTelemetry(void)
-  {
-    for(uint64_t i=0; i<telNormB.size(); i++)
-    {
-      std::cout << GridLogPerformance
-                << "SchurNode " << i
-                << " ||B||_F "       << telNormB[i]
-                << " ||S||/||A22|| " << telSratio[i]
-                << std::endl;
-    }
-    std::cout << GridLogPerformance
-              << "Schur leaves max|Ainv| " << telLeafMaxInv
-              << std::endl;
-
-    RealD armax =  tAllreduce;
-    RealD armin = -tAllreduce;
-    grid->GlobalMax(armax);
-    grid->GlobalMax(armin);
-    armin = -armin;
-
-    std::cout << GridLogMessage << "Schur phases (boss rank, seconds):"
-              << "  memset "    << tMemset/1.0e6
-              << "  deposit "   << tDeposit/1.0e6
-              << "  allreduce " << tAllreduce/1.0e6
-              << "  barrier "   << tBarrier/1.0e6
-              << "  repack "    << tRepack/1.0e6
-              << "  gemm "      << tGemm/1.0e6
-              << "  leaf "      << tLeaf/1.0e6
-              << std::endl;
-    std::cout << GridLogMessage << "Schur comms:"
-              << "  [transports run: allreduce " << (nAllreduce - nGatherV - nBcast)
-              << " allgatherv " << nGatherV << " bcast " << nBcast << "]"
-              << ( barrierProbe ? " [barrier probe on]" : "" )
-              << "  GatherGemm calls " << nGatherGemm
-              << "  panel allreduces " << nAllreduce
-              << "  allreduce GB "     << bytesAllreduce/1024./1024./1024.
-              << "  allreduce s min/max over ranks " << armin/1.0e6
-              << " / " << armax/1.0e6
-              << std::endl;
-    std::cout << GridLogMessage << "Schur comms per-call (boss):"
-              << "  min " << (nAllreduce ? tARmin/1.0e3 : 0.0) << " ms"
-              << "  avg " << (nAllreduce ? tAllreduce/nAllreduce/1.0e3 : 0.0) << " ms"
-              << "  max " << tARmax/1.0e3 << " ms"
-              << "   effective " << (tAllreduce>0 ? bytesAllreduce/tAllreduce*1.0e6/1.0e9 : 0.0)
-              << " GB/s"
-              << std::endl;
-    //////////////////////////////////////////////////////////////////////
-    // Distribution, not just min/avg/max.  A barrier that merely REBOOKS
-    // skew shifts the median; a barrier that suppresses pathological
-    // collectives shortens the tail.  p50 vs p99 separates the two from a
-    // single run, which min/avg/max cannot.
-    //////////////////////////////////////////////////////////////////////
-    if ( tARall.size() ) {
-      std::vector<double> v = tARall;
-      std::sort(v.begin(),v.end());
-      size_t n = v.size();
-      auto pc = [&](double f){ size_t i=(size_t)(f*(n-1)); return v[i]/1.0e3; };
-      double med = pc(0.50);
-      // how much time is spent in calls that are gross outliers
-      double tail=0.0; size_t ntail=0;
-      for(size_t i=0;i<n;i++) if ( v[i] > 4.0*med*1.0e3 ) { tail+=v[i]; ntail++; }
-      std::cout << GridLogMessage << "Schur comms distribution (boss):"
-                << "  p50 "  << med
-                << "  p90 "  << pc(0.90)
-                << "  p99 "  << pc(0.99)
-                << " ms   calls > 4x median: " << ntail << "/" << n
-                << " carrying " << tail/1.0e6 << " s"
-                << std::endl;
-    }
-  }
-};
-
-NAMESPACE_END(Grid);

Grid/parallelIO/BinaryIO.h

@@ -266,9 +266,10 @@ class BinaryIO {
   // GRID_BINARYIO_NOAGGREGATE falls back to plain lexicographic I/O.
   static int DefaultControl(void)
   {
-    static int ctrl = getenv("GRID_BINARYIO_NOAGGREGATE")
-                    ? BINARYIO_LEXICOGRAPHIC
-                    : BINARYIO_LEXICOGRAPHIC|BINARYIO_AGGREGATE;
+    //    static int ctrl = getenv("GRID_BINARYIO_NOAGGREGATE")
+    //                    ? BINARYIO_LEXICOGRAPHIC
+    //                    : BINARYIO_LEXICOGRAPHIC|BINARYIO_AGGREGATE;
+    static int ctrl = BINARYIO_LEXICOGRAPHIC|BINARYIO_AGGREGATE;
     return ctrl;
   }
 

Grid/stencil/Stencil.h

@@ -453,9 +453,19 @@ public:
 	double *dbuf =(double *) packet.recv_buf;
 	float  *fbuf =(float  *) packet.compressed_recv_buf;
 
+	// The lane index is a GPU coalescing optimisation: lane = threadIdx.y
+	// puts adjacent threads on adjacent words.  On a CPU build there are no
+	// SIMT lanes and acceleratorSIMTlane() is identically zero, so the plain
+	// form would convert only 1 word of every nsimd; the #else covers all
+	// lanes explicitly.
 	accelerator_forNB(ss,outer,nsimd,{
+#ifdef GRID_SIMT
 	  int lane = acceleratorSIMTlane(nsimd);
 	  dbuf[ss*nsimd+lane] = fbuf[ss*nsimd+lane]; //conversion
+#else
+	  for(int lane=0;lane<nsimd;lane++)
+	    dbuf[ss*nsimd+lane] = fbuf[ss*nsimd+lane]; //conversion
+#endif
 	});
 
       } else if ( sizeof(word)==4){
@@ -467,8 +477,13 @@ public:
 	uint16_t *hbuf =(uint16_t *) packet.compressed_recv_buf;
 
 	accelerator_forNB(ss,outer,nsimd,{
+#ifdef GRID_SIMT
 	  int lane = acceleratorSIMTlane(nsimd);
 	  fbuf[ss*nsimd+lane] = ((uint32_t)hbuf[ss*nsimd+lane])<<16; //copy back and pad each word with zeroes
+#else
+	  for(int lane=0;lane<nsimd;lane++)
+	    fbuf[ss*nsimd+lane] = ((uint32_t)hbuf[ss*nsimd+lane])<<16; //copy back and pad each word with zeroes
+#endif
 	});
 
       } else {
@@ -515,9 +530,15 @@ public:
 	double *dbuf =(double *) packet.send_buf;
 	float  *fbuf =(float  *) packet.compressed_send_buf;
 
+	// CPU lane coverage as in DecompressPacket.
 	accelerator_forNB(ss,outer,nsimd,{
+#ifdef GRID_SIMT
 	  int lane = acceleratorSIMTlane(nsimd);
 	  fbuf[ss*nsimd+lane] = dbuf[ss*nsimd+lane]; // convert fp64 to fp32
+#else
+	  for(int lane=0;lane<nsimd;lane++)
+	    fbuf[ss*nsimd+lane] = dbuf[ss*nsimd+lane]; // convert fp64 to fp32
+#endif
 	});
 
       } else if ( sizeof(word)==4){
@@ -526,8 +547,13 @@ public:
 	uint16_t *hbuf =(uint16_t *) packet.compressed_send_buf;
 	
 	accelerator_forNB(ss,outer,nsimd,{
+#ifdef GRID_SIMT
 	  int lane = acceleratorSIMTlane(nsimd);
 	  hbuf[ss*nsimd+lane] = fbuf[ss*nsimd+lane]>>16; // convert as in Bagel/BFM ; bfloat16 ; s7e8 Intel patent
+#else
+	  for(int lane=0;lane<nsimd;lane++)
+	    hbuf[ss*nsimd+lane] = fbuf[ss*nsimd+lane]>>16; // convert as in Bagel/BFM ; bfloat16 ; s7e8 Intel patent
+#endif
 	});
 
       } else {

examples/Example_pvdagm_multigrid.cc

@@ -0,0 +1,151 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./examples/Example_pvdagm_multigrid.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 */
+
+//
+// The canonical three-level mrhs PVdagM multigrid driver, on the library
+// objects (Grid/algorithms/multigrid/PVdagMMultiGrid.h).
+//
+// Every parameter is in ONE serialisable struct, read from XML:
+//
+//   ./Example_pvdagm_multigrid --grid 48.48.48.96 --mpi ... \
+//                              --pvdagm-params params.xml
+//
+// With no --pvdagm-params the built-in defaults run (hot-start gauge field
+// unless Config is set in the file).  A missing file gets a template
+// written next to it and the program exits: the template documents every
+// parameter.  The effective parameters are always printed, so the log
+// describes its own run.
+//
+// Compile-time: -DNBASIS=8 cuts the basis down for laptop runs.
+//
+#include <Grid/Grid.h>
+#include <Grid/lattice/PaddedCell.h>
+#include <Grid/stencil/GeneralLocalStencil.h>
+#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
+#include <Grid/algorithms/multigrid/PVdagMMultiGrid.h>
+
+using namespace std;
+using namespace Grid;
+
+#ifndef NBASIS
+#define NBASIS 60
+#endif
+
+struct PVdagMDriverParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(PVdagMDriverParams,
+                                  int,         Ls,
+                                  RealD,       Mass,
+                                  RealD,       M5,
+                                  RealD,       MobiusB,
+                                  RealD,       MobiusC,
+                                  std::string, Config,          // empty: hot start
+                                  int,         Nrhs,
+                                  int,         SolveSingleRHS,  // also run Nrhs=1 through the same objects
+                                  PVdagMMultiGridParams, MultiGrid);
+  PVdagMDriverParams()
+    : Ls(24), Mass(0.00078), M5(1.8), MobiusB(1.5), MobiusC(0.5),
+      Config(""), Nrhs(12), SolveSingleRHS(1) {};
+};
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  PVdagMDriverParams P;
+  {
+    std::string pfile("");
+    if( GridCmdOptionExists(argv,argv+argc,"--pvdagm-params") )
+      pfile = GridCmdOptionPayload(argv,argv+argc,"--pvdagm-params");
+    if ( pfile.length() ) {
+      bool good; { std::ifstream f(pfile); good = f.good(); }
+      if ( !good ) {
+        if ( GlobalSharedMemory::WorldRank == 0 ) {
+          XmlWriter WR(pfile+".templ");
+          write(WR, "PVdagMDriver", P);
+          std::cout << GridLogMessage << pfile << " does not exist; template written to "
+                    << pfile << ".templ" << std::endl;
+        }
+        Grid_finalize(); return 0;
+      }
+      XmlReader RD(pfile);
+      read(RD, "PVdagMDriver", P);
+    }
+    CheckValidity(P.MultiGrid);
+    std::cout << GridLogMessage << "PVdagMDriver parameters ("
+              << (pfile.length() ? pfile : std::string("defaults")) << "):" << std::endl;
+    std::cout << P << std::endl;
+  }
+
+  Coordinate latt = GridDefaultLatt();
+  Coordinate mpi  = GridDefaultMpi();
+  Coordinate fsimd= GridDefaultSimd(Nd,vComplex::Nsimd());
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(latt,fsimd,mpi);
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(P.Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(P.Ls,UGrid);
+
+  GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers({1,2,3,4});
+  GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8});
+
+  LatticeGaugeField Umu(UGrid);
+  if ( P.Config.length() ) {
+    std::cout << GridLogMessage << "Reading gauge field " << P.Config << std::endl;
+    FieldMetaData header;
+    NerscIO::readConfiguration(Umu,header,P.Config);
+  } else {
+    std::cout << GridLogMessage << "Hot start gauge field" << std::endl;
+    SU<Nc>::HotConfiguration(RNG4,Umu);
+  }
+
+  MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,P.Mass,P.M5,P.MobiusB,P.MobiusC);
+  MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,   P.M5,P.MobiusB,P.MobiusC);
+
+  typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
+  PVdagM_t PVdagM(Ddwf,Dpv);
+
+  //////////////////////////////////////////////////////////////////////
+  // Grids -> coarsening -> dense bottom.  Scope order is lifetime order.
+  //////////////////////////////////////////////////////////////////////
+  typedef PVdagMMultiGridCoarsening<vSpinColourVector,sTComplexD,NBASIS> Coarsening_t;
+
+  MGCoarseGrids CGrids(FGrid, P.MultiGrid.Setup);
+  Coarsening_t  Coarsening(CGrids, P.MultiGrid.Setup);
+
+  Coarsening.GetSubspace(RNG5, PVdagM);
+  Coarsening.Coarsen(PVdagM);
+  Coarsening.BuildDenseBottom();
+
+  //////////////////////////////////////////////////////////////////////
+  // Solves: mrhs then (optionally) single RHS through the SAME objects.
+  //////////////////////////////////////////////////////////////////////
+  auto RunSolve = [&](int nr)
+  {
+    PVdagMMultiGridSolver<MobiusFermionD,Coarsening_t> Solver(Ddwf,Dpv,Coarsening,P.MultiGrid,nr);
+    std::vector<LatticeFermionD> src(nr,FGrid), sol(nr,FGrid);
+    for(int r=0;r<nr;r++){ gaussian(RNG5,src[r]); sol[r]=Zero(); }
+    Solver.Solve(src,sol);
+  };
+
+  RunSolve(P.Nrhs);
+  if ( P.SolveSingleRHS && P.Nrhs != 1 ) RunSolve(1);
+
+  Grid_finalize();
+}

systems/Frontier/config-command

@@ -2,9 +2,6 @@ CLIME=`spack find --paths c-lime@2-3-9 | grep c-lime| cut -c 15-`
 
 ../../configure --enable-comms=mpi-auto \
 --with-lime=$CLIME \
---with-slate=$SLATE \
---with-blaspp=$BLASPP \
---with-lapackpp=$LAPACKPP \
 --enable-unified=no \
 --enable-shm=nvlink \
 --enable-tracing=roctx \
@@ -16,7 +13,6 @@ CLIME=`spack find --paths c-lime@2-3-9 | grep c-lime| cut -c 15-`
 --with-openssl=$OPENSSL \
 --with-gmp=$GMP \
 --with-mpfr=$MPFR \
---disable-fermion-reps \
 CXX=hipcc MPICXX=mpicxx \
 CXXFLAGS="-fPIC -I${ROCM_PATH}/include/ -I${MPICH_DIR}/include " \
  LDFLAGS="-L${ROCM_PATH}/lib -L${MPICH_DIR}/lib -lmpi -lmpi_gtl_hsa -lhipblas -lrocsolver -lrocblas -lhipfft -lamdhip64"

systems/Frontier/pvdagm_multigrid.job

@@ -0,0 +1,148 @@
+#!/bin/bash -l
+#SBATCH --job-name=pvdagm-multigrid
+#SBATCH --nodes=36
+#SBATCH --ntasks-per-node=8
+#SBATCH --cpus-per-task=7
+#SBATCH --gpus-per-node=8
+#SBATCH --time=1:30:00
+#SBATCH --account=phy157_dwf
+#SBATCH --gpu-bind=none
+#SBATCH --exclusive
+#SBATCH --mem=0
+#SBATCH -S 0
+
+##############################################################################
+# Production run of the three-level mrhs PVdagM multigrid on the library
+# objects (examples/Example_pvdagm_multigrid.cc), 48^3x96 Ls=24, 288 GCDs.
+#
+# All algorithm/physics parameters are in ONE serialisable struct read from
+# XML (--pvdagm-params); there are NO Grid environment knobs any more.  The
+# library-internal choices the old nrhs_fault_36.job set by env are now the
+# only code path (2D block-cyclic dense inverse, ring-allgather apply,
+# inverseLU big leaves) and are gone from the environment.
+#
+# The binary must be built with -DNBASIS=64 to match the subspace file.
+##############################################################################
+
+root=$HOME/PVdagM/Grid/systems/Frontier
+source $root/sourceme-rocm7.2.sh
+
+# Paths substituted into the XML below.
+LUSTRE=/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle
+SUBSPACE=$LUSTRE/subspace_nb64.scidac
+CONFIG=/ccs/home/paboyle/configs/ckpoint_lat.1000
+RUNDIR=$LUSTRE/runs/${SLURM_JOB_NAME}_${SLURM_JOB_ID}
+mkdir -p $RUNDIR && cd $RUNDIR && echo "RUNDIR $RUNDIR"
+
+
+cat << EOF > select_gpu
+#!/bin/bash
+export GPU_MAP=(0 1 2 3 7 6 5 4)
+export NUMA_MAP=(3 3 1 1 2 2 0 0)
+export GPU=\${GPU_MAP[\$SLURM_LOCALID]}
+export NUMA=\${NUMA_MAP[\$SLURM_LOCALID]}
+export HIP_VISIBLE_DEVICES=\$GPU
+unset ROCR_VISIBLE_DEVICES
+if [ \$SLURM_PROCID = "0" ]; then echo \$*; fi
+exec numactl -m \$NUMA -N \$NUMA \$*
+EOF
+chmod +x ./select_gpu
+
+
+export OMP_NUM_THREADS=7
+ulimit -c 0                             # no 22 GB GPU core dumps
+export FI_MR_CACHE_MONITOR=kdreg2       # site default; device-buffer MPI on Slingshot
+export MPICH_GPU_SUPPORT_ENABLED=1
+export MPICH_SMP_SINGLE_COPY_MODE=CMA
+export MPICH_OFI_NIC_POLICY=GPU
+
+#### This raise the probability of failure ####
+export GRID_ALLOC_NCACHE_LARGE=0       # allocator large-ring depth (Grid allocator, not a solver knob) -- default is 8
+##############################################################################
+
+
+module load libfabric
+
+BIN=$root/examples/Example_pvdagm_multigrid
+OPTS="--accelerator-threads 8 --shm 4096 --shm-mpi 0 --device-mem 32000 --comms-overlap"
+vol=48.48.48.96
+MPI_GEOM=3.6.4.4
+
+
+##############################################################################
+# The parameter file.  Values are the banked adaptive operating point carried
+# over from nrhs_fault_36.job.  Tags mirror PVdagMDriverParams /
+# PVdagMMultiGridParams exactly; a mistyped tag deserialises silently to the
+# struct default, so keep this in step with the headers.
+##############################################################################
+cat << EOF > params.xml
+<?xml version="1.0"?>
+<grid>
+  <PVdagMDriver>
+    <Ls>24</Ls>
+    <Mass>0.00078</Mass>
+    <M5>1.8</M5>
+    <MobiusB>1.5</MobiusB>
+    <MobiusC>0.5</MobiusC>
+    <Config>$CONFIG</Config>
+    <Nrhs>12</Nrhs>
+    <SolveSingleRHS>1</SolveSingleRHS>
+    <MultiGrid>
+      <Setup>
+        <Block><elem>2</elem><elem>2</elem><elem>3</elem><elem>3</elem></Block>
+        <Block2><elem>4</elem><elem>4</elem><elem>2</elem><elem>4</elem></Block2>
+        <CoarsenBatch>9</CoarsenBatch>
+        <SubspaceFile>$SUBSPACE</SubspaceFile>
+        <FineSloppyComms>1</FineSloppyComms>
+      </Setup>
+      <FineSmoother>
+        <Shift>0.1</Shift><Nstep>6</Nstep><Mmax>4</Mmax>
+      </FineSmoother>
+      <CoarseSmoother>
+        <Shift>2.0</Shift><Nstep>2</Nstep><Mmax>2</Mmax>
+      </CoarseSmoother>
+      <CoarseSolver>
+        <Tol>0.05</Tol><Order>200</Order><Mmax>8</Mmax>
+      </CoarseSolver>
+      <Outer>
+        <Tol>1e-8</Tol><MaxIterations>1000</MaxIterations><Mmax>6</Mmax><Nstep>12</Nstep>
+      </Outer>
+      <Dense>
+        <LeafSpan>9</LeafSpan>
+      </Dense>
+    </MultiGrid>
+  </PVdagMDriver>
+</grid>
+EOF
+echo "=== params.xml ==="; cat params.xml; echo "=================="
+
+run_cell () {
+  name=$1; shift
+  echo "----- $name : $* -----"
+  export GRID_STDOUT_ROOT=$RUNDIR/mg_$name
+  env "$@" srun -N36 -n288 --kill-on-bad-exit=1 ./select_gpu $BIN --mpi ${MPI_GEOM} --grid $vol $OPTS \
+       --pvdagm-params params.xml --debug-stdout --log Error,Warning,Message,Performance > log.mg.$name 2>&1
+  echo "  exit $?"; sleep 30
+  f=$(grep -l "Memory access fault\|NO_TRANSLATION" $GRID_STDOUT_ROOT/*/Grid.stderr.* 2>/dev/null | head -1)
+  if [ -n "$f" ]; then
+    echo "  FAULT in $f  (rank $(basename $f | sed 's/Grid.stderr.//'))"; tail -20 $f | cut -c1-160
+  else
+    echo "  no fault; convergence + timing:"
+    grep -h "3-level solve\|FINAL Nrhs\|MrhsPGCR: Converged" $GRID_STDOUT_ROOT/0/Grid.stdout.0 | tail -8 | cut -c1-140
+  fi
+}
+
+# the NRHS>=12 fault was intermittent (~1/3) before the
+# dmabuf avoidance, so run it twice to confirm a clean pair.
+
+# sourceme sets FI_HMEM_ROCR_USE_DMABUF=0 -- the standing avoidance of the CXI
+# NO_TRANSLATION fault at NRHS>=12 (libfabric #12775).
+export FI_HMEM_ROCR_USE_DMABUF=0
+run_cell  dmabuf0_run1 
+run_cell  dmabuf0_run2 
+run_cell  dmabuf0_run3 
+
+export FI_HMEM_ROCR_USE_DMABUF=1
+run_cell  dmabuf1_run1 
+run_cell  dmabuf1_run2 
+run_cell  dmabuf1_run3 

systems/Frontier/sourceme-rocm6.4.sh

@@ -14,3 +14,8 @@ export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH
 export LD_LIBRARY_PATH=/opt/rocm-6.4.0/lib/llvm/lib/:$LD_LIBRARY_PATH
 module load emacs
 export FI_MR_CACHE_MONITOR=kdreg2       # REQUIRED for device-buffer MPI on Slingshot: libfabric memhooks monitor (default) is defective, see systems/WorkArounds.txt (libfabric #11451)
+
+# Standing avoidance of the CXI NO_TRANSLATION fault on persistent device
+# buffers under hipMalloc/hipFree churn (libfabric issue #12775; ruling
+# 2026-09-07). 
+export FI_HMEM_ROCR_USE_DMABUF=0

systems/Frontier/sourceme-rocm7.0.sh

@@ -12,3 +12,8 @@ export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$CRAY_LD_LIBRARY_PATH
 
 module load emacs
 export FI_MR_CACHE_MONITOR=kdreg2       # REQUIRED for device-buffer MPI on Slingshot: libfabric memhooks monitor (default) is defective, see systems/WorkArounds.txt (libfabric #11451)
+
+# Standing avoidance of the CXI NO_TRANSLATION fault on persistent device
+# buffers under hipMalloc/hipFree churn (libfabric issue #12775; ruling
+# 2026-09-07). 
+export FI_HMEM_ROCR_USE_DMABUF=0

systems/Frontier/sourceme-rocm7.2.sh

@@ -3,17 +3,17 @@ module unload slate 2>/dev/null
 echo spack
 . /autofs/nccs-svm1_home1/paboyle/spack-frontier/share/spack/setup-env.sh
 
-export SLATE=`spack find --paths slate   | grep ^slate  | awk '{print $2}' `
-export BLASPP=`spack find --paths blaspp | grep ^blaspp | awk '{print $2}' `
-export LAPACKPP=`spack find --paths lapackpp | grep ^lapackpp | awk '{print $2}' `
+#export SLATE=`spack find --paths slate   | grep ^slate  | awk '{print $2}' `
+#export BLASPP=`spack find --paths blaspp | grep ^blaspp | awk '{print $2}' `
+#export LAPACKPP=`spack find --paths lapackpp | grep ^lapackpp | awk '{print $2}' `
 
-echo SLATE    $SLATE
-echo LAPACKPP $LAPACKPP
-echo BLASPP   $BLASPP
+#echo SLATE    $SLATE
+#echo LAPACKPP $LAPACKPP
+#echo BLASPP   $BLASPP
 
-ls $SLATE/lib64/libslate.so $LAPACKPP/lib64/liblapackpp.so 2>/dev/null || ls $SLATE/lib $LAPACKPP/lib
-ldd $SLATE/lib64/libslate.so | grep -E "blaspp|lapackpp|rocblas|amdhip|mpi_cray|libsci"
-ldd $LAPACKPP/lib64/liblapackpp.so | grep -E "blaspp|libsci"
+#ls $SLATE/lib64/libslate.so $LAPACKPP/lib64/liblapackpp.so 2>/dev/null || ls $SLATE/lib $LAPACKPP/lib
+#ldd $SLATE/lib64/libslate.so | grep -E "blaspp|lapackpp|rocblas|amdhip|mpi_cray|libsci"
+#ldd $LAPACKPP/lib64/liblapackpp.so | grep -E "blaspp|libsci"
 
 export CLIME=`spack find --paths c-lime | grep ^c-lime | awk '{print $2}' `
 export MPFR=`spack find --paths mpfr    | grep ^mpfr  | awk '{print $2}' `
@@ -24,9 +24,14 @@ module load cce/21.0.0
 module load cpe/26.03
 module load rocm/7.2.0
 export LD_LIBRARY_PATH=/opt/rocm-7.2.0/lib/llvm/lib/:$LD_LIBRARY_PATH
-export LD_LIBRARY_PATH=$LAPACKPP/lib64:$LD_LIBRARY_PATH
-export LD_LIBRARY_PATH=$BLASPP/lib64:$LD_LIBRARY_PATH
-export LD_LIBRARY_PATH=$SLATE/lib64:$LD_LIBRARY_PATH
+#export LD_LIBRARY_PATH=$LAPACKPP/lib64:$LD_LIBRARY_PATH
+#export LD_LIBRARY_PATH=$BLASPP/lib64:$LD_LIBRARY_PATH
+#export LD_LIBRARY_PATH=$SLATE/lib64:$LD_LIBRARY_PATH
 export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$CRAY_LD_LIBRARY_PATH
 module load emacs
+
+# Standing avoidance of the CXI NO_TRANSLATION fault on persistent device
+# buffers under hipMalloc/hipFree churn (libfabric issue #12775; ruling
+# 2026-09-07).  
 export FI_MR_CACHE_MONITOR=kdreg2       # REQUIRED for device-buffer MPI on Slingshot: libfabric memhooks monitor (default) is defective, see systems/WorkArounds.txt (libfabric #11451)
+export FI_HMEM_ROCR_USE_DMABUF=0

tests/debug/Test_allgather.cc

@@ -244,7 +244,7 @@ int main(int argc, char **argv)
   // T7 : the SAME shape as T6, but assembled by C sequential MPI_Bcast --
   // one broadcast per contributing rank -- instead of one MPI_Allgatherv.
   //
-  // This is the transport of DENSE_GATHER=2.  Bcast takes no count vector,
+  // This is the transport of the retired chunked-Bcast gather.  Bcast takes no count vector,
   // so the zero-count asymmetry that makes T6 run at ~0.18 MB/s and trip
   // mpir_request.h:508 cannot arise.  It costs C collectives rather than 1
   // and the roots do not transmit concurrently, so the byte cost is about

tests/debug/Test_schur2d.cc

@@ -49,8 +49,8 @@ using namespace Grid;
 static int failures = 0;
 
 // Portable |z|: ComplexD is std::complex on CPU builds and thrust::complex
-// under HIP, where std::abs does not resolve (same trap RecursiveSchurInverse
-// documents at FrobNorm2Local).  Member real()/imag() work on both.
+// under HIP, where std::abs does not resolve.  Member real()/imag() work
+// on both.
 static double Cabs(const ComplexD &z)
 {
   double re = z.real(), im = z.imag();

tests/debug/Test_schur2d_redist.cc

@@ -23,7 +23,7 @@ Author: Peter Boyle <pboyle@bnl.gov>
 //
 //   1D rank-major rows -> block cyclic -> Invert -> back to 1D rows
 //
-// which is exactly what DENSE_SCHUR2D runs inside DenseCoarseMatrix.
+// which is exactly what the dense inverse runs inside DenseCoarseMatrix.
 // CPU build under mpirun at n = 1,2,3,4.
 //
 //   T1 : RowsToCyclic against a direct ImportGlobal of the same matrix --
@@ -31,14 +31,9 @@ Author: Peter Boyle <pboyle@bnl.gov>
 //   T2 : round trip rows -> 2D -> rows -- BITWISE, uniform AND non-uniform
 //        rowStart, layouts with ragged trailing blocks.
 //   T3 : full pipeline inverse against a host Gauss-Jordan reference.
-//   T4 : CROSS-IMPLEMENTATION: the same matrix inverted by the 1D
-//        RecursiveSchurInverse and by the 2D pipeline; results compared
-//        element-wise.  Two independent implementations, two independent
-//        decompositions, one answer.
 //////////////////////////////////////////////////////////////////////////////
 
 #include <Grid/Grid.h>
-#include <Grid/algorithms/multigrid/RecursiveSchurInverse.h>
 #include <Grid/algorithms/multigrid/BlockCyclicSchurInverse.h>
 #include <Grid/algorithms/multigrid/BlockCyclicRedistribute.h>
 
@@ -47,8 +42,8 @@ using namespace Grid;
 static int failures = 0;
 
 // Portable |z|: ComplexD is std::complex on CPU builds and thrust::complex
-// under HIP, where std::abs does not resolve (same trap RecursiveSchurInverse
-// documents at FrobNorm2Local).  Member real()/imag() work on both.
+// under HIP, where std::abs does not resolve.  Member real()/imag() work
+// on both.
 static double Cabs(const ComplexD &z)
 {
   double re = z.real(), im = z.imag();
@@ -190,12 +185,11 @@ int main(int argc, char **argv)
   }
 
   ////////////////////////////////////////////////////////////////////////
-  // T3 + T4 : the DENSE_SCHUR2D pipeline against the host reference and
-  // against the INDEPENDENT 1D RecursiveSchurInverse.
+  // T3 : the 2D pipeline against the host reference.
   ////////////////////////////////////////////////////////////////////////
   {
-    bool ok3 = true, ok4 = true;
-    double worst3 = 0.0, worst4 = 0.0;
+    bool ok3 = true;
+    double worst3 = 0.0;
     BlockCyclicSchurInverse RSI2;
     for(auto &g : grids){
       for(auto &c : cfgs){
@@ -232,27 +226,9 @@ int main(int argc, char **argv)
             worst3 = std::max(worst3,d);
             if ( d > 1.0e-9 ) ok3 = false;
           }
-
-        // ---- 1D RecursiveSchurInverse on the same matrix ----
-        {
-          BlockRows Ar;  Ar.Resize(myrows, N);
-          acceleratorCopyToDevice(&h[0], &Ar.data[0], h.size()*sizeof(ComplexD));
-          std::vector<int64_t> rs = rowStart;
-          RecursiveSchurInverse RSI1(grid, N, rs, 1<<20);
-          RSI1.Invert(Ar);
-          std::vector<ComplexD> h1d(h.size());
-          acceleratorCopyFromDevice(&Ar.data[0], &h1d[0], h1d.size()*sizeof(ComplexD));
-          for(int64_t j=0;j<N;j++)
-            for(int64_t i=0;i<myrows;i++){
-              double d = Cabs(h2d[i+j*myrows]-h1d[i+j*myrows])/mxref;
-              worst4 = std::max(worst4,d);
-              if ( d > 1.0e-9 ) ok4 = false;
-            }
-        }
       }
     }
     Report("T3  2D pipeline vs host reference", ok3, "worst "+std::to_string(worst3));
-    Report("T4  2D pipeline vs 1D RecursiveSchurInverse", ok4, "worst "+std::to_string(worst4));
   }
 
   {

tests/debug/Test_schur2d_scale.cc

@@ -20,7 +20,7 @@ Author: Peter Boyle <pboyle@bnl.gov>
 
 //////////////////////////////////////////////////////////////////////////////
 // SCALE rehearsal for the 2D distributed dense inverse: the full
-// DENSE_SCHUR2D pipeline -- 1D rows -> redistribute -> invert ->
+// 2D block-cyclic pipeline -- 1D rows -> redistribute -> invert ->
 // redistribute back -> certificate -- on a SYNTHETIC matrix of any size,
 // with no multigrid machinery, no configuration and no subspace file.
 //
@@ -31,8 +31,8 @@ Author: Peter Boyle <pboyle@bnl.gov>
 // rank; at N=138240 on 288 ranks it is the production problem shape
 // exactly, in a driver that runs in minutes.
 //
-//   S2D_N   : global dimension          (default 720, laptop friendly)
-//   S2D_NB  : block size                (default N/P rows-per-rank if that
+//   --schur2d-global-dimension <n> : N  (default 720, laptop friendly)
+//   --schur2d-block-size <n>       : nb (default N/P rows-per-rank if that
 //                                        is exact, else 48)
 //
 // The matrix is diagonally dominant (the recursion does not pivot); its
@@ -55,8 +55,8 @@ Author: Peter Boyle <pboyle@bnl.gov>
 using namespace Grid;
 
 // Portable |z|: ComplexD is std::complex on CPU builds and thrust::complex
-// under HIP, where std::abs does not resolve (same trap RecursiveSchurInverse
-// documents at FrobNorm2Local).  Member real()/imag() work on both.
+// under HIP, where std::abs does not resolve.  Member real()/imag() work
+// on both.
 static double Cabs(const ComplexD &z)
 {
   double re = z.real(), im = z.imag();
@@ -76,13 +76,6 @@ static ComplexD Fill(int64_t i, int64_t j, int64_t N)
 
 int main(int argc, char **argv)
 {
-  // Environment walk (2026-08-27, systems/Frontier/schur2d_env.job): knobs to
-  // reproduce the example's environment here, one at a time.  Result: thread
-  // level, OMP_NUM_THREADS, device residency and sustained load all NIL; only
-  // "first job step on fresh nodes" (+3 s) is real.
-  //   S2D_BALLAST_GB=x     x GB of Lattice fields made device-resident before the invert
-  //   S2D_PREHEAT_S=x      x seconds of back-to-back zgemm before the invert
-  //   OMP_NUM_THREADS      set in the job, read by nothing here but the runtime
   Grid_init(&argc, &argv);
 
   GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
@@ -91,9 +84,12 @@ int main(int argc, char **argv)
   const int P  = grid->ProcessorCount();
   const int me = grid->ThisRank();
 
-  int64_t N  = getenv("S2D_N")  ? atol(getenv("S2D_N"))  : 720;
+  int64_t N = 720;
+  if ( GridCmdOptionExists(argv,argv+argc,"--schur2d-global-dimension") )
+    N = atol(GridCmdOptionPayload(argv,argv+argc,"--schur2d-global-dimension").c_str());
   int64_t nb;
-  if      ( getenv("S2D_NB") ) nb = atol(getenv("S2D_NB"));
+  if      ( GridCmdOptionExists(argv,argv+argc,"--schur2d-block-size") )
+    nb = atol(GridCmdOptionPayload(argv,argv+argc,"--schur2d-block-size").c_str());
   else if ( N % P == 0 )       nb = N/P;
   else                         nb = 48;
   GRID_ASSERT( N >= 1 ); GRID_ASSERT( nb >= 1 );
@@ -125,56 +121,13 @@ int main(int argc, char **argv)
   double t1 = usecond();
 
   ////////////////////////////////////////////////////////////////////////
-  // The DENSE_SCHUR2D pipeline, phase-timed.  A0 keeps the original for
+  // The 2D pipeline, phase-timed.  A0 keeps the original for
   // the certificate.
   ////////////////////////////////////////////////////////////////////////
   BlockCyclicMatrix A (grid,N,nb,Pr,Pc);
   BlockCyclicMatrix A0(grid,N,nb,Pr,Pc);
   BlockCyclicSchurInverse RSI2;
 
-  // Pre-heat: drive the GCD with back-to-back zgemm for S2D_PREHEAT_S seconds
-  // before the invert.  The example calls the inverse after ~100 s of full
-  // load on all 288 GCDs and its LOCAL kernels run 20-40% slower than the
-  // idle-start harness (GEMM 3.1 vs 2.5 s, leaf 0.76 vs 0.24 s) with the
-  // wires unchanged; thread level / OMP / residency (E1-E5) did not reproduce
-  // that.  If sustained load does, it is clock/power management, not code.
-  if ( getenv("S2D_PREHEAT_S") ) {
-    double secs = atof(getenv("S2D_PREHEAT_S"));
-    const int64_t W = 4320;
-    deviceVector<ComplexD> M((uint64_t)W*W), C((uint64_t)W*W);
-    { ComplexD *m = &M[0]; accelerator_for(idx,(uint64_t)W*W,1,{ m[idx] = ComplexD(1.0e-3*(idx%97),1.0e-3*(idx%89)); }); accelerator_barrier(); }
-    deviceVector<ComplexD*> ap(1),bp(1),cp(1); std::vector<ComplexD*> ptr(1);
-    ptr[0]=&M[0]; acceleratorCopyToDevice(&ptr[0],&ap[0],sizeof(ComplexD*)); acceleratorCopyToDevice(&ptr[0],&bp[0],sizeof(ComplexD*));
-    ptr[0]=&C[0]; acceleratorCopyToDevice(&ptr[0],&cp[0],sizeof(ComplexD*));
-    double t0=usecond(); int n=0; double tlast=0;
-    while ( (usecond()-t0)/1.0e6 < secs ) {
-      double t1=usecond();
-      RSI2.SUMMA.BLAS.gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N,(int)W,(int)W,(int)W,ComplexD(1.0,0.0),ap,(int)W,bp,(int)W,ComplexD(0.0,0.0),cp,(int)W);
-      RSI2.SUMMA.BLAS.synchronise(); tlast=usecond()-t1; n++;
-    }
-    double tfirst = 0; (void)tfirst;
-    std::cout << GridLogMessage << "Test_schur2d_scale: pre-heat " << (usecond()-t0)/1.0e6 << " s, " << n << " zgemm W=" << W
-              << ", last zgemm " << tlast/1.0e6 << " s (" << 8.0*W*W*W/tlast/1.0e6 << " TF/s; idle-start rate 23.7)" << std::endl;
-  }
-
-  // Device ballast: Lattice fields written on the accelerator so they sit in
-  // the MemoryManager's device LRU exactly as the example's fine-grid state does.
-  typedef Lattice<iVector<iVector<vComplexD,Nc>,Ns> > BallastField;
-  std::vector<BallastField> ballast;
-  if ( getenv("S2D_BALLAST_GB") ) {
-    double gb = atof(getenv("S2D_BALLAST_GB"));
-    uint64_t fbytes = (uint64_t)grid->oSites()*sizeof(BallastField::vector_object);
-    int nf = (int)(gb*1.0e9/(double)fbytes + 0.5);
-    ballast.reserve(nf);
-    for(int i=0;i<nf;i++){
-      ballast.emplace_back(grid);
-      autoView(v, ballast[i], AcceleratorWriteDiscard);
-      accelerator_for(ss, grid->oSites(), 1, { v[ss] = Zero(); });
-    }
-    std::cout << GridLogMessage << "Test_schur2d_scale: device ballast " << nf << " fields x " << fbytes/1.0e6
-              << " MB = " << nf*fbytes/1.0e9 << " GB resident (S2D_BALLAST_GB=" << gb << ")" << std::endl;
-  }
-
   BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,A);
   double t2 = usecond();
   if ( A.data.size() )

tests/debug/Test_schur2d_vs_slate.cc

@@ -64,9 +64,10 @@ Author: Peter Boyle <pboyle@bnl.gov>
 // libblaspp shadows the ROCm one via LD_LIBRARY_PATH and throws
 // "device BLAS not available" from host_malloc_pinned.
 //
-//   S2D_N, S2D_NB as in Test_schur2d_scale (default nb = N/P).
-//   S2D_SKIP_GETRI=1 skips the getri leg (host loop; ~4 min at N=138240).
-//   S2D_NOWARM=1 skips the warm-up.
+//   --schur2d-global-dimension, --schur2d-block-size as in
+//   Test_schur2d_scale (default nb = N/P).
+//   --schur2d-skip-getri skips the getri leg (host loop; ~4 min at N=138240).
+//   --schur2d-nowarm skips the warm-up.
 // A third leg, getrf+getrs(I), is SLATE's device-resident inverse route.
 //////////////////////////////////////////////////////////////////////////////
 
@@ -138,8 +139,12 @@ int main(int argc, char **argv)
   const int P  = grid->ProcessorCount();
   const int me = grid->ThisRank();
 
-  int64_t N  = getenv("S2D_N")  ? atol(getenv("S2D_N"))  : 720;
-  int64_t nb = getenv("S2D_NB") ? atol(getenv("S2D_NB")) : ( (N%P==0) ? N/P : 48 );
+  int64_t N = 720;
+  if ( GridCmdOptionExists(argv,argv+argc,"--schur2d-global-dimension") )
+    N = atol(GridCmdOptionPayload(argv,argv+argc,"--schur2d-global-dimension").c_str());
+  int64_t nb = (N%P==0) ? N/P : 48;
+  if ( GridCmdOptionExists(argv,argv+argc,"--schur2d-block-size") )
+    nb = atol(GridCmdOptionPayload(argv,argv+argc,"--schur2d-block-size").c_str());
   int Pr,Pc; BlockCyclicLayout::ChooseProcessGrid(P,Pr,Pc);
 
   std::vector<int64_t> rowStart(P+1); rowStart[0]=0;
@@ -161,10 +166,10 @@ int main(int argc, char **argv)
   // first.  Run a small throwaway inverse through BOTH paths so the timed
   // legs below measure hot code.  Not reported.
   ////////////////////////////////////////////////////////////////////////
-  // S2D_NOWARM=1 skips it (hang localisation).  Stage markers are flushed so
+  // --schur2d-nowarm skips it (hang localisation).  Stage markers are flushed so
   // a hang shows WHERE even through block-buffered stdout.
   auto Stage = [&](const char *s){ std::cout << GridLogMessage << "stage: " << s << std::endl << std::flush; };
-  if ( !getenv("S2D_NOWARM") ) {
+  if ( !GridCmdOptionExists(argv,argv+argc,"--schur2d-nowarm") ) {
     // Fixed tiny size independent of P: the purpose is handle creation and
     // kernel loading, not work.  (8*P at P=288 was N=2304 -> a 122 s SLATE
     // warm-up dominated by 288-way tile broadcasts.)  Ranks beyond the first
@@ -207,7 +212,7 @@ int main(int argc, char **argv)
 #endif
     std::cout << GridLogMessage << "warm-up done (both paths, N=" << Nw << ")" << std::endl << std::flush;
   } else {
-    std::cout << GridLogMessage << "warm-up SKIPPED (S2D_NOWARM)" << std::endl << std::flush;
+    std::cout << GridLogMessage << "warm-up SKIPPED (--schur2d-nowarm)" << std::endl << std::flush;
   }
 
   ////////////////////////////////////////////////////////////////////////
@@ -240,7 +245,7 @@ int main(int argc, char **argv)
   // LEG 2: SLATE, every layout step timed and charged.
   ////////////////////////////////////////////////////////////////////////
 #ifdef HAVE_SLATE
-  if ( !getenv("S2D_SKIP_GETRI") ) {   // S2D_SKIP_GETRI=1: getri is a host loop, minutes at N=138240
+  if ( !GridCmdOptionExists(argv,argv+argc,"--schur2d-skip-getri") ) {   // getri is a host loop, minutes at N=138240
     typedef std::complex<double> scalar_t;
     acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD));
     BlockCyclicMatrix A(grid,N,nb,Pr,Pc), A0(grid,N,nb,Pr,Pc);

tests/debug/Test_schur_dense_coarse.cc

@@ -27,25 +27,19 @@ Author: Peter Boyle <pboyle@bnl.gov>
 /*  END LEGAL */
 
 //
-// T6 of the RecursiveSchurInverse regression chain
-// (schur_recursive_inverse_plan.txt 4B.5): the DenseCoarseMatrix GLUE,
-// on a real (tiny) lattice coarse operator, CPU laptop build.
+// The DenseCoarseMatrix GLUE test, on a real (tiny) lattice coarse
+// operator, CPU laptop build.
 //
 // Builds a genuine GeneralCoarsenedMatrix (DWF MdagM + 0.5 shift for a
-// guaranteed-invertible Galerkin coarse op, random aggregation basis,
-// nbasis=8, 4^4 x Ls/1 blocking) and constructs DenseCoarseMatrix in
-// DENSE_SCHUR=2 AUDIT mode with small DENSE_PANEL_BYTES (multi-panel
-// gathers exercised through the glue).  The constructor then runs, in
-// order, all the certificates this stage exists to check:
-//   - fresh ImportDense (no SLAB_FILE) + IMPORT CERTIFICATE vs Op.M
-//   - InvertDenseSingle (the oracle)
-//   - InvertDenseSchur: self-certifying rank-major map, fp64 diagonal
-//     import certificate vs the fp32 slab, distributed recursion,
-//     growth telemetry
-//   - AUDIT: max|Ainv_schur - Ainv_single| over the full slab
-//   - VERIFY ||A Ainv x - x||/||x|| through the SCHUR result
-// This program adds asserts on the audit number and a random-vector
-// round trip.
+// guaranteed-invertible Galerkin coarse op, random aggregation basis)
+// and runs the whole Import certificate chain through the glue:
+//   - fresh ImportDense + IMPORT CERTIFICATE vs Op.M
+//   - fp64 rank-major import certificate vs the fp32 slab
+//   - the 2D block-cyclic recursion + growth telemetry
+//   - VERIFY ||A Ainv x - x||/||x|| through the device split-K apply
+// This program adds an INDEPENDENT Eigen fp64 host-inverse oracle at
+// small N (built from applies of M to unit vectors, so it shares no
+// code with the import) and a random-vector round trip.
 //
 // Uniform local volume 12.12.12.12 (fine), per-dim blocks {4,4,3,3},
 // coarse 3.3.4.4/rank, nbasis 4 (N = 576n):
@@ -193,13 +187,18 @@ int main (int argc, char ** argv)
 
   ///////////////////////////////////////////////////////////////////////
   // Full-matrix conditioning probe at small N: dense columns by
-  // applying M to unit vectors, fp64 Eigen SVD.
+  // applying M to unit vectors, fp64 Eigen SVD.  eA is kept: it is the
+  // INDEPENDENT oracle for the inverse below (built from applies of M,
+  // sharing no code with the stencil->dense import).
   ///////////////////////////////////////////////////////////////////////
+  int64_t Nprobe = Coarse5d->gSites() * nbasis;
+  Eigen::MatrixXcd eA;
+  bool haveOracle = false;
   {
-    int64_t Nprobe = Coarse5d->gSites() * nbasis;
     if ( Nprobe <= 700 )
     {
-      Eigen::MatrixXcd eA(Nprobe, Nprobe);
+      eA.resize(Nprobe, Nprobe);
+      haveOracle = true;
       CoarseVector e(Coarse5d);
       CoarseVector Me(Coarse5d);
       for(int64_t j=0; j<Nprobe; j++)
@@ -255,26 +254,15 @@ int main (int argc, char ** argv)
   }
 
   ///////////////////////////////////////////////////////////////////////
-  // T6: AUDIT mode, fresh import, multi-panel gathers.  The constructor
-  // runs every certificate in the chain (see banner).
+  // The glue under test: Import runs the whole certificate chain
+  // (import certificate, fp64 certificate, 2D inverse, VERIFY).
   ///////////////////////////////////////////////////////////////////////
-  setenv("DENSE_SCHUR","2",1);
-  setenv("DENSE_PANEL_BYTES","65536",1);
-  unsetenv("SLAB_FILE");
-
-  // Restructured interface: DenseCoarseMatrix<CComplex,nbasis>, constructed
-  // on the grid and fed by Import (which runs the certificate chain).
   typedef DenseCoarseMatrix<vTComplex,nbasis> DenseCC;
   DenseCC dcm(Coarse5d);
   dcm.Import(LittleDiracOp);
 
-  std::cout << GridLogMessage << "T6  audit relative slab difference (schur vs single) = "
-            << dcm.schurAuditRel << std::endl;
-  GRID_ASSERT( dcm.schurAuditRel >= 0.0 );      // audit actually ran
-  GRID_ASSERT( dcm.schurAuditRel <  1.0e-3 );
-
   ///////////////////////////////////////////////////////////////////////
-  // Random-vector round trip through the SCHUR inverse
+  // Random-vector round trip through the inverse
   ///////////////////////////////////////////////////////////////////////
   CoarseVector x(Coarse5d);
   CoarseVector y(Coarse5d);
@@ -288,6 +276,38 @@ int main (int argc, char ** argv)
             << rel << std::endl;
   GRID_ASSERT( rel < 1.0e-2 );
 
+  ///////////////////////////////////////////////////////////////////////
+  // Independent oracle at small N: y from dcm must match the Eigen fp64
+  // solve of eA (dense columns of M itself) on the same x.  Catches an
+  // inverse that is self-consistent with a WRONG import, which the round
+  // trip above cannot (dense and M would share the error).
+  ///////////////////////////////////////////////////////////////////////
+  if ( haveOracle )
+  {
+    Eigen::VectorXcd xs(Nprobe), ys(Nprobe);
+    typedef typename CoarseVector::vector_object::scalar_object csobj;
+    for(int64_t j=0; j<Nprobe; j++)
+    {
+      int64_t gsite = j / nbasis;
+      int     b     = j % nbasis;
+      Coordinate gcoor(Coarse5d->_ndimension);
+      Lexicographic::CoorFromIndex(gcoor, gsite, Coarse5d->GlobalDimensions());
+      csobj s;
+      peekSite(s, x, gcoor);
+      ComplexD zz = ((ComplexD *)&s)[b];
+      xs(j) = std::complex<double>(zz.real(), zz.imag());
+      peekSite(s, y, gcoor);
+      zz = ((ComplexD *)&s)[b];
+      ys(j) = std::complex<double>(zz.real(), zz.imag());
+    }
+    Eigen::VectorXcd yref = eA.fullPivLu().solve(xs);
+    double dev  = (ys - yref).cwiseAbs().maxCoeff();
+    double ymax = yref.cwiseAbs().maxCoeff();
+    std::cout << GridLogMessage << "T6  oracle max|Ainv x - eigen solve|/max|y| = "
+              << dev/ymax << "   (fp32 slab vs fp64 host solve)" << std::endl;
+    GRID_ASSERT( dev/ymax < 1.0e-4 );
+  }
+
   std::cout << GridLogMessage << "Test_schur_dense_coarse: T6 ALL PASS" << std::endl;
 
   Grid_finalize();

tests/debug/Test_schur_inverse.cc

@@ -1,769 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: Test_schur_inverse.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 */
-
-//
-// Staged regression gate for RecursiveSchurInverse (distributed dense
-// inversion by recursive Schur complement) -- the laptop-side certificate
-// chain of schur_recursive_inverse_plan.txt section 4B.5.  Runs on a
-// CPU-only build (Eigen BLAS backends) under mpirun:
-//
-//   mpirun -n 1 ./Test_schur_inverse --grid 8.8.8.8  --mpi 1.1.1.1
-//   mpirun -n 2 ./Test_schur_inverse --grid 8.8.8.8  --mpi 1.1.1.2
-//   mpirun -n 3 ./Test_schur_inverse --grid 8.8.8.12 --mpi 1.1.1.3
-//   mpirun -n 4 ./Test_schur_inverse --grid 8.8.8.8  --mpi 1.1.1.4
-//
-// (n=3 exercises uneven row splits throughout.)  The lattice exists only
-// to furnish the communicator; no field is ever constructed.
-//
-// PRECISION: the inversion runs ENTIRELY in fp64 (decision 2026-08-14,
-// superseding the fp32-merge design); certificates are eps64-scaled.
-// The single terminal fp32 rounding belongs to the caller (tested at the
-// glue level, Test_schur_dense_coarse).
-//
-// Stages present (cumulative -- earlier tests are never removed):
-//   T1a : ownership tables -- CheckRowStart on synthetic uneven partitions,
-//         MakeRowStart allgather vs closed form on the live communicator.
-//   T1b : STORAGE-CONVENTION PIN -- column-major + ld + window-offset
-//         semantics fixed once via identity multiplies through the
-//         explicit-ld gemmBatched, on INTEGER-VALUED data so all three
-//         cases below are EXACT (values well within the mantissa):
-//           (1) alpha=1,beta=0 read from an input column window
-//           (2) alpha=-1,beta=1 accumulate (the S-formation case)
-//           (3) write INTO an output column window, neighbours untouched
-//         No later failure can be a transposition/convention ambiguity.
-//   T2  : GatherGemm vs naive fp64 oracle (owner sub-ranges, alpha-beta
-//         cases, tiny+huge panels, half-participation call shape).
-//   T3  : LeafInvert in-place residual certificate.
-//   T4  : full recursive Invert vs Eigen fp64 oracle, growth-scaled
-//         certification, adversarial near-singular-A11 family with
-//         telemetry-spike assertion.
-//
-// Hard asserts throughout; thresholds pre-registered in the plan.
-//
-#include <Grid/Grid.h>
-#include <Grid/Grid_Eigen_Dense.h>
-#include <Grid/algorithms/multigrid/RecursiveSchurInverse.h>
-
-using namespace std;
-using namespace Grid;
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  GridCartesian Comm(GridDefaultLatt(),
-                     GridDefaultSimd(Nd,vComplex::Nsimd()),
-                     GridDefaultMpi());
-  GridBase *grid = &Comm;
-
-  ////////////////////////////////////////////////////////////////
-  // T1a : ownership tables
-  ////////////////////////////////////////////////////////////////
-  {
-    // Synthetic partitions of N=97 (prime: every P>1 is uneven)
-    const int64_t N = 97;
-    for(int P=1; P<=4; P++)
-    {
-      std::vector<int64_t> table(P+1);
-      table[0] = 0;
-      for(int r=0; r<P; r++)
-      {
-        int64_t nr = N/P + ( (r < (int)(N%P)) ? 1 : 0 );
-        table[r+1] = table[r] + nr;
-      }
-      RecursiveSchurInverse::CheckRowStart(table, N);
-    }
-
-    // Live allgather: deliberately uneven local counts, closed-form oracle
-    int P  = grid->ProcessorCount();
-    int me = grid->ThisRank();
-
-    int64_t myNrows = 3 + me;
-    std::vector<int64_t> table = RecursiveSchurInverse::MakeRowStart(grid, myNrows);
-
-    std::vector<int64_t> expect(P+1);
-    expect[0] = 0;
-    for(int r=0; r<P; r++)
-    {
-      expect[r+1] = expect[r] + (3 + r);
-    }
-    GRID_ASSERT( (int)table.size() == P+1 );
-    for(int r=0; r<=P; r++)
-    {
-      GRID_ASSERT( table[r] == expect[r] );
-    }
-
-    // Constructor smoke: derived ownership matches
-    RecursiveSchurInverse RSI(grid, table[P], table, 1024*1024);
-    GRID_ASSERT( RSI.P       == P );
-    GRID_ASSERT( RSI.me      == me );
-    GRID_ASSERT( RSI.myRow0  == expect[me] );
-    GRID_ASSERT( RSI.myNrows == myNrows );
-
-    std::cout << GridLogMessage
-      << "T1a ownership tables (synthetic P=1..4, live allgather, ctor)   PASS"
-      << std::endl;
-  }
-
-  ////////////////////////////////////////////////////////////////
-  // T1b : storage-convention pin (every rank, local, exact)
-  ////////////////////////////////////////////////////////////////
-  {
-    const int64_t rows = 5;
-    const int64_t cols = 13;
-    const int64_t col0 = 6;   // input window start
-    const int64_t w    = 4;   // window width
-
-    // f(i,j): integer-valued, unique per element
-    auto f = [](int64_t i, int64_t j) -> ComplexD
-    {
-      return ComplexD( (RealD)(1 + i + 10*j), (RealD)(i - j) );
-    };
-
-    BlockRows A;
-    A.Resize(rows, cols);
-    {
-      std::vector<ComplexD> Ahost((uint64_t)rows*cols);
-      for(int64_t j=0; j<cols; j++)
-      {
-        for(int64_t i=0; i<rows; i++)
-        {
-          Ahost[(uint64_t)(i + j*rows)] = f(i,j);
-        }
-      }
-      acceleratorCopyToDevice(&Ahost[0], &A.data[0], (uint64_t)rows*cols*sizeof(ComplexD));
-    }
-
-    // Identity I_w, column major
-    deviceVector<ComplexD> Idev((uint64_t)w*w);
-    {
-      std::vector<ComplexD> Ihost((uint64_t)w*w, ComplexD(0.0,0.0));
-      for(int64_t d=0; d<w; d++)
-      {
-        Ihost[(uint64_t)(d + d*w)] = ComplexD(1.0,0.0);
-      }
-      acceleratorCopyToDevice(&Ihost[0], &Idev[0], (uint64_t)w*w*sizeof(ComplexD));
-    }
-
-    GridBLAS BLAS;
-    ComplexD one    ( 1.0,0.0);
-    ComplexD minus  (-1.0,0.0);
-    ComplexD zero   ( 0.0,0.0);
-
-    deviceVector<ComplexD*> Ap(1);
-    deviceVector<ComplexD*> Bp(1);
-    deviceVector<ComplexD*> Cp(1);
-    std::vector<ComplexD*>  ptr_h(1);
-
-    auto setptr = [&](deviceVector<ComplexD*> &d, ComplexD *p)
-    {
-      ptr_h[0] = p;
-      acceleratorCopyToDevice(&ptr_h[0], &d[0], sizeof(ComplexD*));
-    };
-
-    ////////////////////////////////////////////////////////////
-    // Case 1: C = A(:, col0:col0+w) . I_w      (alpha=1, beta=0)
-    ////////////////////////////////////////////////////////////
-    {
-      deviceVector<ComplexD> Cdev((uint64_t)rows*w);
-      setptr(Ap, A.ColumnWindow(col0));
-      setptr(Bp, &Idev[0]);
-      setptr(Cp, &Cdev[0]);
-
-      BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N,
-                       (int)rows, (int)w, (int)w,
-                       one,  Ap, (int)A.ld,
-                             Bp, (int)w,
-                       zero, Cp, (int)rows);
-      BLAS.synchronise();
-
-      std::vector<ComplexD> Chost((uint64_t)rows*w);
-      acceleratorCopyFromDevice(&Cdev[0], &Chost[0], (uint64_t)rows*w*sizeof(ComplexD));
-      for(int64_t j=0; j<w; j++)
-      {
-        for(int64_t i=0; i<rows; i++)
-        {
-          GRID_ASSERT( Chost[(uint64_t)(i + j*rows)] == f(i, col0+j) );
-        }
-      }
-    }
-
-    ////////////////////////////////////////////////////////////
-    // Case 2: C = C0 - A(:, col0:col0+w) . I_w (alpha=-1, beta=1)
-    // -- the S-formation accumulate; exact on integer data
-    ////////////////////////////////////////////////////////////
-    {
-      auto g = [](int64_t i, int64_t j) -> ComplexD
-      {
-        return ComplexD( (RealD)(100 + i + j), (RealD)7 );
-      };
-      deviceVector<ComplexD> Cdev((uint64_t)rows*w);
-      {
-        std::vector<ComplexD> Chost((uint64_t)rows*w);
-        for(int64_t j=0; j<w; j++)
-        {
-          for(int64_t i=0; i<rows; i++)
-          {
-            Chost[(uint64_t)(i + j*rows)] = g(i,j);
-          }
-        }
-        acceleratorCopyToDevice(&Chost[0], &Cdev[0], (uint64_t)rows*w*sizeof(ComplexD));
-      }
-      setptr(Ap, A.ColumnWindow(col0));
-      setptr(Bp, &Idev[0]);
-      setptr(Cp, &Cdev[0]);
-
-      BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N,
-                       (int)rows, (int)w, (int)w,
-                       minus, Ap, (int)A.ld,
-                              Bp, (int)w,
-                       one,   Cp, (int)rows);
-      BLAS.synchronise();
-
-      std::vector<ComplexD> Chost((uint64_t)rows*w);
-      acceleratorCopyFromDevice(&Cdev[0], &Chost[0], (uint64_t)rows*w*sizeof(ComplexD));
-      for(int64_t j=0; j<w; j++)
-      {
-        for(int64_t i=0; i<rows; i++)
-        {
-          ComplexD expect = g(i,j) - f(i, col0+j);
-          GRID_ASSERT( Chost[(uint64_t)(i + j*rows)] == expect );
-        }
-      }
-    }
-
-    ////////////////////////////////////////////////////////////
-    // Case 3: write INTO a column window of a wider C;
-    // columns outside the window must be untouched
-    ////////////////////////////////////////////////////////////
-    {
-      const int64_t ccols = 6;
-      const int64_t cw0   = 2;   // output window start
-      BlockRows C;
-      C.Resize(rows, ccols);
-      {
-        std::vector<ComplexD> Chost((uint64_t)rows*ccols, ComplexD(-999.0, 999.0));
-        acceleratorCopyToDevice(&Chost[0], &C.data[0], (uint64_t)rows*ccols*sizeof(ComplexD));
-      }
-      setptr(Ap, A.ColumnWindow(col0));
-      setptr(Bp, &Idev[0]);
-      setptr(Cp, C.ColumnWindow(cw0));
-
-      BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N,
-                       (int)rows, (int)w, (int)w,
-                       one,  Ap, (int)A.ld,
-                             Bp, (int)w,
-                       zero, Cp, (int)C.ld);
-      BLAS.synchronise();
-
-      std::vector<ComplexD> Chost((uint64_t)rows*ccols);
-      acceleratorCopyFromDevice(&C.data[0], &Chost[0], (uint64_t)rows*ccols*sizeof(ComplexD));
-      for(int64_t j=0; j<ccols; j++)
-      {
-        for(int64_t i=0; i<rows; i++)
-        {
-          ComplexD got = Chost[(uint64_t)(i + j*rows)];
-          if ( (j >= cw0) && (j < cw0+w) )
-          {
-            GRID_ASSERT( got == f(i, col0 + (j-cw0)) );
-          }
-          else
-          {
-            GRID_ASSERT( got == ComplexD(-999.0, 999.0) );
-          }
-        }
-      }
-    }
-
-    std::cout << GridLogMessage
-      << "T1b storage-convention pin (window read / S-accumulate / window write, exact)   PASS"
-      << std::endl;
-  }
-
-  ////////////////////////////////////////////////////////////////
-  // T2 : GatherGemm vs naive double-precision oracle.
-  //
-  // Every rank generates the SAME full N x N random fp64 operands
-  // from a fixed seed (no comms needed for the oracle), keeps only
-  // its own rows in BlockRows form, and after each GatherGemm call
-  // checks its output window element-by-element against a plain
-  // triple-loop ComplexD accumulation over the same entries.
-  //
-  // Sweep: N in {8, 96, 97}; owner ranges full/upper-half/single;
-  // (alpha,beta) in {(1,0), (-1,1)}; panelBytes tiny (ragged
-  // many-chunk gathers) and huge (single panel).  Sentinel columns
-  // outside the output window must be untouched.  Finally, a
-  // HALF-PARTICIPATION case rehearses the recursion call pattern:
-  // lower ranks own B but pass EMPTY A/C (collectives only).
-  ////////////////////////////////////////////////////////////////
-  {
-    int P  = grid->ProcessorCount();
-    int me = grid->ThisRank();
-
-    std::mt19937 rng(777);
-    std::uniform_real_distribution<double> dist(-1.0,1.0);
-
-    const int64_t nout = 5;   // output width
-    const int64_t colB = 3;   // B window offset
-    const int64_t colC = 2;   // C window offset
-
-    for(int64_t N : {8L, 96L, 97L})
-    {
-      // Ownership: uneven for any P not dividing N
-      std::vector<int64_t> table(P+1);
-      table[0] = 0;
-      for(int r=0; r<P; r++)
-      {
-        int64_t nr = N/P + ( (r < (int)(N%P)) ? 1 : 0 );
-        table[r+1] = table[r] + nr;
-      }
-      int64_t r0   = table[me];
-      int64_t myNr = table[me+1] - table[me];
-
-      // Identical full operands on every rank
-      std::vector<ComplexD> Aglob((uint64_t)N*N);
-      std::vector<ComplexD> Bglob((uint64_t)N*N);
-      for(uint64_t i=0; i<(uint64_t)N*N; i++) Aglob[i] = ComplexD(dist(rng),dist(rng));
-      for(uint64_t i=0; i<(uint64_t)N*N; i++) Bglob[i] = ComplexD(dist(rng),dist(rng));
-
-      // My rows of a full-matrix operand as a BlockRows
-      auto fillRows = [&](BlockRows &X, std::vector<ComplexD> &glob,
-                          int64_t row0, int64_t nr)
-      {
-        X.Resize(nr, N);
-        if ( nr == 0 ) return;
-        std::vector<ComplexD> h((uint64_t)nr*N);
-        for(int64_t j=0; j<N; j++)
-        {
-          for(int64_t i=0; i<nr; i++)
-          {
-            h[(uint64_t)(i + j*nr)] = glob[(uint64_t)((row0+i) + j*N)];
-          }
-        }
-        acceleratorCopyToDevice(&h[0], &X.data[0], (uint64_t)nr*N*sizeof(ComplexD));
-      };
-
-      // Owner-range cases: full span, upper half, single interior rank
-      std::vector<std::pair<int,int> > ranges;
-      ranges.push_back(std::make_pair(0, P));
-      if ( P > 1 ) ranges.push_back(std::make_pair(P/2, P));
-      if ( P > 1 ) ranges.push_back(std::make_pair(1, 2));
-
-      for(auto range : ranges)
-      {
-        int rB0 = range.first;
-        int rB1 = range.second;
-        int64_t ka0 = table[rB0];       // A-column window start = B row span
-        int64_t k   = table[rB1] - table[rB0];
-
-        for(int acase=0; acase<2; acase++)
-        {
-          ComplexD alpha = ( acase==0 ) ? ComplexD( 1.0,0.0) : ComplexD(-1.0,0.0);
-          ComplexD beta  = ( acase==0 ) ? ComplexD( 0.0,0.0) : ComplexD( 1.0,0.0);
-
-          for(int64_t panelBytes : {64L, 1L<<30})
-          {
-            RecursiveSchurInverse RSI(grid, N, table, panelBytes);
-
-            BlockRows A;
-            BlockRows B;
-            BlockRows C;
-            fillRows(A, Aglob, r0, myNr);
-            fillRows(B, Bglob, r0, myNr);
-
-            // Output: sentinel-filled, window at colC
-            const ComplexD sentinel(-999.0, 999.0);
-            const int64_t  ccols = colC + nout + 2;
-            C.Resize(myNr, ccols);
-            std::vector<ComplexD> C0((uint64_t)myNr*ccols, sentinel);
-            if ( acase == 1 )
-            {
-              // beta=1 needs defined window content: g(i,j), integer-valued
-              for(int64_t j=0; j<nout; j++)
-              {
-                for(int64_t i=0; i<myNr; i++)
-                {
-                  C0[(uint64_t)(i + (colC+j)*myNr)] = ComplexD((RealD)(50+i+j), (RealD)-3);
-                }
-              }
-            }
-            if ( myNr > 0 )
-            {
-              acceleratorCopyToDevice(&C0[0], &C.data[0], (uint64_t)myNr*ccols*sizeof(ComplexD));
-            }
-
-            RSI.GatherGemm(alpha, A, ka0, k,
-                           rB0, rB1,
-                           B, colB, nout,
-                           beta, C, colC);
-
-            std::vector<ComplexD> Chost((uint64_t)myNr*ccols);
-            if ( myNr > 0 )
-            {
-              acceleratorCopyFromDevice(&C.data[0], &Chost[0], (uint64_t)myNr*ccols*sizeof(ComplexD));
-            }
-
-            double tol = 1.0e-14 * (double)k;
-            for(int64_t j=0; j<ccols; j++)
-            {
-              for(int64_t i=0; i<myNr; i++)
-              {
-                ComplexD got = Chost[(uint64_t)(i + j*myNr)];
-                if ( (j >= colC) && (j < colC+nout) )
-                {
-                  int64_t jj = j - colC;
-                  ComplexD acc(0.0,0.0);
-                  if ( acase == 1 )
-                  {
-                    acc = C0[(uint64_t)(i + j*myNr)];
-                  }
-                  for(int64_t t=0; t<k; t++)
-                  {
-                    acc += alpha
-                         * Aglob[(uint64_t)((r0+i)  + (ka0+t)*N)]
-                         * Bglob[(uint64_t)((ka0+t) + (colB+jj)*N)];
-                  }
-                  GRID_ASSERT( abs(got - acc) < tol );
-                }
-                else
-                {
-                  GRID_ASSERT( got == sentinel );
-                }
-              }
-            }
-          }
-        }
-      }
-
-      ////////////////////////////////////////////////////////////
-      // Half-participation: owners = [0,ph) hold B; participants
-      // = [ph,P) hold A/C; owners pass EMPTY A/C and column
-      // offset 0 (collectives only) -- the recursion call shape.
-      ////////////////////////////////////////////////////////////
-      if ( P > 1 )
-      {
-        int ph = ( P+1 ) / 2;
-        int64_t ka0 = table[0];
-        int64_t k   = table[ph] - table[0];
-        int participant = ( me >= ph );
-
-        RecursiveSchurInverse RSI(grid, N, table, 64);
-
-        BlockRows A;
-        BlockRows B;
-        BlockRows C;
-        fillRows(B, Bglob, r0, myNr);
-        if ( participant )
-        {
-          fillRows(A, Aglob, r0, myNr);
-          C.Resize(myNr, nout);
-        }
-
-        ComplexD one (1.0,0.0);
-        ComplexD zero(0.0,0.0);
-        int64_t  cA = participant ? ka0 : 0;
-        RSI.GatherGemm(one, A, cA, k,
-                       0, ph,
-                       B, colB, nout,   // owners deposit from their B window
-                       zero, C, 0);
-
-        if ( participant )
-        {
-          std::vector<ComplexD> Chost((uint64_t)myNr*nout);
-          acceleratorCopyFromDevice(&C.data[0], &Chost[0], (uint64_t)myNr*nout*sizeof(ComplexD));
-          double tol = 1.0e-14 * (double)k;
-          for(int64_t j=0; j<nout; j++)
-          {
-            for(int64_t i=0; i<myNr; i++)
-            {
-              ComplexD acc(0.0,0.0);
-              for(int64_t t=0; t<k; t++)
-              {
-                acc += Aglob[(uint64_t)((r0+i)  + (ka0+t)*N)]
-                     * Bglob[(uint64_t)((ka0+t) + (colB+j)*N)];
-              }
-              GRID_ASSERT( abs(Chost[(uint64_t)(i + j*myNr)] - acc) < tol );
-            }
-          }
-        }
-      }
-    }
-
-    std::cout << GridLogMessage
-      << "T2  GatherGemm vs oracle (N=8/96/97, 3 owner ranges, 2 alpha-beta, tiny+huge panels, half-participation)   PASS"
-      << std::endl;
-  }
-
-  ////////////////////////////////////////////////////////////////
-  // T3 : LeafInvert -- in-place fp64 inversion of the contiguous
-  // leaf window.  Purely local, every rank runs its own
-  // uneven-size leaf; residual certificate in ComplexD.
-  ////////////////////////////////////////////////////////////////
-  {
-    int P  = grid->ProcessorCount();
-    int me = grid->ThisRank();
-
-    int64_t  w   = 17 + 3*me;
-    uint64_t len = (uint64_t)w*w;
-
-    std::vector<int64_t> table = RecursiveSchurInverse::MakeRowStart(grid, w);
-    RecursiveSchurInverse RSI(grid, table[P], table, 1<<20);
-
-    // A = w I + R : well conditioned
-    std::mt19937 rng(31 + me);
-    std::uniform_real_distribution<double> dist(-1.0,1.0);
-    std::vector<ComplexD> Ahost(len);
-    for(uint64_t i=0; i<len; i++) Ahost[i] = ComplexD(dist(rng),dist(rng));
-    for(int64_t d=0; d<w; d++)  Ahost[(uint64_t)(d + d*w)] += ComplexD((RealD)w, 0.0);
-
-    BlockRows Ar;
-    Ar.Resize(w, w);
-    acceleratorCopyToDevice(&Ahost[0], &Ar.data[0], len*sizeof(ComplexD));
-    RSI.LeafInvert(0, w, Ar);
-    std::vector<ComplexD> X(len);
-    acceleratorCopyFromDevice(&Ar.data[0], &X[0], len*sizeof(ComplexD));
-
-    double maxdev = 0.0;
-    for(int64_t j=0; j<w; j++)
-    {
-      for(int64_t i=0; i<w; i++)
-      {
-        ComplexD acc(0.0,0.0);
-        for(int64_t t=0; t<w; t++)
-        {
-          acc += Ahost[(uint64_t)(i + t*w)] * X[(uint64_t)(t + j*w)];
-        }
-        if ( i==j ) acc -= ComplexD(1.0,0.0);
-        maxdev = std::max(maxdev, abs(acc));
-      }
-    }
-    GRID_ASSERT( maxdev < 1.0e-13 );
-
-    std::cout << GridLogMessage
-      << "T3  LeafInvert in-place fp64 (residual " << maxdev << ")   PASS" << std::endl;
-  }
-
-  ////////////////////////////////////////////////////////////////
-  // T4 : full recursive Invert vs Eigen fp64 oracle.
-  //
-  // Every rank builds the SAME N x N fp64 matrix from a fixed seed,
-  // keeps its rows, inverts through the full SPMD recursion, then
-  // the test gathers the complete inverse (zero-fill GlobalSum) and
-  // checks BOTH certificates:
-  //     cert1 = || A X - I ||_max          (ComplexD accumulation)
-  //     cert2 = max|X - Xref| / max|Xref|  (Xref = Eigen fp64 inverse)
-  //
-  // Families (eps64-scaled tolerances; the fp32-era growth data
-  // rescales by eps64/eps32 ~ 1.9e-9):
-  //   kappa-moderate : A = R + 3 sqrt(N) I
-  //   kappa-large    : A = R + 0.3 sqrt(N) I
-  //   adversarial    : leading block (rank 0's whole leaf) REPLACED by
-  //                    1e-2 * (R' + 3 sqrt(b) I) inside a well-conditioned
-  //                    A -- the growth spike must REGISTER in telemetry
-  //                    (asserted > 10 when P > 1); at fp64 the certificate
-  //                    barely notices it: that insensitivity IS the point
-  //                    of the fp64 conversion.
-  //
-  // N=64 runs with panelBytes=128 (ragged many-chunk gathers inside
-  // the recursion); larger N with 1 MB panels.
-  ////////////////////////////////////////////////////////////////
-  {
-    int P  = grid->ProcessorCount();
-    int me = grid->ThisRank();
-
-    std::mt19937 rng(2026);
-    std::uniform_real_distribution<double> dist(-1.0,1.0);
-
-    for(int64_t N : {64L, 200L, 513L})
-    {
-      std::vector<int64_t> table(P+1);
-      table[0] = 0;
-      for(int r=0; r<P; r++)
-      {
-        int64_t nr = N/P + ( (r < (int)(N%P)) ? 1 : 0 );
-        table[r+1] = table[r] + nr;
-      }
-      int64_t r0   = table[me];
-      int64_t myNr = table[me+1] - table[me];
-
-      for(int fam=0; fam<3; fam++)
-      {
-        const char *famname = (fam==0) ? "kappa-moderate" :
-                              (fam==1) ? "kappa-large"    : "adversarial-A11";
-        double shift = (fam==1) ? 0.3*std::sqrt((double)N) : 3.0*std::sqrt((double)N);
-        double tol   = (fam==0) ? 1.0e-12 :
-                       (fam==1) ? 1.0e-11 : 5.0e-11;
-
-        // Identical operand on every rank (all draws rank-independent)
-        std::vector<ComplexD> Aglob((uint64_t)N*N);
-        for(uint64_t i=0; i<(uint64_t)N*N; i++) Aglob[i] = ComplexD(dist(rng),dist(rng));
-        for(int64_t d=0; d<N; d++) Aglob[(uint64_t)(d + d*N)] += ComplexD(shift, 0.0);
-        if ( fam == 2 )
-        {
-          // Leading block = rank 0's whole leaf, scaled down 100x but
-          // internally well conditioned (shift scales as sqrt(b): a
-          // FIXED shift makes A11 itself near-singular at large b).
-          int64_t b = ( P > 1 ) ? table[1] : N/4;
-          for(int64_t j=0; j<b; j++)
-          {
-            for(int64_t i=0; i<b; i++)
-            {
-              Aglob[(uint64_t)(i + j*N)] = ComplexD(0.01,0.0)*ComplexD(dist(rng),dist(rng));
-            }
-          }
-          RealD bshift = (RealD)(0.03*std::sqrt((double)b));
-          for(int64_t d=0; d<b; d++) Aglob[(uint64_t)(d + d*N)] += ComplexD(bshift,0.0);
-        }
-
-        // Eigen fp64 oracle.  Explicit re/im conversion at the boundary:
-        // on HIP builds ComplexD is thrust::complex, which has no
-        // operators against Eigen's std::complex.
-        auto toStd = [](const ComplexD &z) -> std::complex<double>
-        {
-          return std::complex<double>(z.real(), z.imag());
-        };
-        Eigen::MatrixXcd eA(N,N);
-        for(int64_t j=0; j<N; j++)
-        {
-          for(int64_t i=0; i<N; i++)
-          {
-            eA(i,j) = toStd(Aglob[(uint64_t)(i + j*N)]);
-          }
-        }
-        Eigen::MatrixXcd Xref = eA.inverse();
-
-        // Distribute, invert
-        int64_t panelBytes = ( N == 64 ) ? 128 : (1<<20);
-        RecursiveSchurInverse RSI(grid, N, table, panelBytes);
-
-        BlockRows Arows;
-        Arows.Resize(myNr, N);
-        {
-          std::vector<ComplexD> h((uint64_t)myNr*N);
-          for(int64_t j=0; j<N; j++)
-          {
-            for(int64_t i=0; i<myNr; i++)
-            {
-              h[(uint64_t)(i + j*myNr)] = Aglob[(uint64_t)((r0+i) + j*N)];
-            }
-          }
-          acceleratorCopyToDevice(&h[0], &Arows.data[0], (uint64_t)myNr*N*sizeof(ComplexD));
-        }
-
-        RSI.Invert(Arows);
-
-        // Gather the full inverse: zero-fill + GlobalSum
-        std::vector<ComplexD> Xfull((uint64_t)N*N, ComplexD(0.0,0.0));
-        {
-          std::vector<ComplexD> h((uint64_t)myNr*N);
-          acceleratorCopyFromDevice(&Arows.data[0], &h[0], (uint64_t)myNr*N*sizeof(ComplexD));
-          for(int64_t j=0; j<N; j++)
-          {
-            for(int64_t i=0; i<myNr; i++)
-            {
-              Xfull[(uint64_t)((r0+i) + j*N)] = h[(uint64_t)(i + j*myNr)];
-            }
-          }
-        }
-        grid->GlobalSumVector(&Xfull[0], (int)(N*N));
-
-        // cert1 = ||A X - I||_max
-        double cert1 = 0.0;
-        for(int64_t j=0; j<N; j++)
-        {
-          for(int64_t i=0; i<N; i++)
-          {
-            ComplexD acc(0.0,0.0);
-            for(int64_t t=0; t<N; t++)
-            {
-              acc += Aglob[(uint64_t)(i + t*N)] * Xfull[(uint64_t)(t + j*N)];
-            }
-            if ( i==j ) acc -= ComplexD(1.0,0.0);
-            cert1 = std::max(cert1, abs(acc));
-          }
-        }
-
-        // cert2 = max|X - Xref| / max|Xref|
-        double maxref = 0.0;
-        double maxdif = 0.0;
-        for(int64_t j=0; j<N; j++)
-        {
-          for(int64_t i=0; i<N; i++)
-          {
-            maxref = std::max(maxref, std::abs(Xref(i,j)));
-            maxdif = std::max(maxdif, std::abs(toStd(Xfull[(uint64_t)(i + j*N)]) - Xref(i,j)));
-          }
-        }
-        double cert2 = maxdif / maxref;
-
-        double maxNormB = 0.0;
-        for(uint64_t i=0; i<RSI.telNormB.size(); i++)
-        {
-          maxNormB = std::max(maxNormB, RSI.telNormB[i]);
-        }
-
-        // GROWTH-SCALED certification, eps64 (the fp32-era model with
-        // eps swapped: cert2 ~ (10-12) ||B||_F sqrt(N) eps; threshold =
-        // 3x margin, floored at the family tolerance).  ||B||_F capped
-        // per family so growth cannot silently excuse a logic error.
-        // cert1 remains a loose absolute bound (an O(1) logic error
-        // gives cert1 ~ 1e2-1e3; fp64 rounding gives ~1e-10).
-        double eps64    = 2.3e-16;
-        double tolModel = 30.0 * std::max(1.0, maxNormB) * std::sqrt((double)N) * eps64;
-        double tolEff   = std::max(tol, tolModel);
-        double capB     = (fam==0) ? 100.0 : (fam==1) ? 2000.0 : 10000.0;
-
-        std::cout << GridLogMessage
-          << "T4  N=" << N << " " << famname
-          << "  ||AX-I||_max " << cert1
-          << "  |X-Xref|/|Xref| " << cert2
-          << "  max||B||_F " << maxNormB
-          << "  tolEff " << tolEff
-          << ( (cert2 < tolEff) && (cert1 < 1.0e-6) ? "   PASS" : "   FAIL" )
-          << std::endl;
-
-        GRID_ASSERT( cert2 < tolEff );
-        GRID_ASSERT( cert1 < 1.0e-6 );
-        GRID_ASSERT( maxNormB < capB );
-        if ( (fam == 2) && (P > 1) )
-        {
-          GRID_ASSERT( maxNormB > 10.0 );   // the spike must REGISTER
-        }
-      }
-    }
-
-    std::cout << GridLogMessage
-      << "T4  recursive Invert vs Eigen oracle (N=64/200/513, 3 families)   PASS"
-      << std::endl;
-  }
-
-  std::cout << GridLogMessage
-    << "Test_schur_inverse: ALL STAGES PASS" << std::endl;
-
-  Grid_finalize();
-}

tests/debug/Test_summa.cc

@@ -50,8 +50,8 @@ using namespace Grid;
 static int failures = 0;
 
 // Portable |z|: ComplexD is std::complex on CPU builds and thrust::complex
-// under HIP, where std::abs does not resolve (same trap RecursiveSchurInverse
-// documents at FrobNorm2Local).  Member real()/imag() work on both.
+// under HIP, where std::abs does not resolve.  Member real()/imag() work
+// on both.
 static double Cabs(const ComplexD &z)
 {
   double re = z.real(), im = z.imag();

tests/solver/Test_multigrid_common.h

@@ -37,9 +37,9 @@ Zero zero;
 // serialization strategy of Grid?
 
 // clang-format off
-struct MultiGridParams : Serializable {
+struct WilsonMGParams : Serializable {
 public:
-  GRID_SERIALIZABLE_CLASS_MEMBERS(MultiGridParams,
+  GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonMGParams,
                                   int,                           nLevels,
                                   std::vector<std::vector<int>>, blockSizes,           // size == nLevels - 1
                                   std::vector<double>,           smootherTol,          // size == nLevels - 1
@@ -54,7 +54,7 @@ public:
                                   int,                           coarseSolverMaxInnerIter);
 
   // constructor with default values
-  MultiGridParams(int                           _nLevels                  = 2,
+  WilsonMGParams(int                           _nLevels                  = 2,
                   std::vector<std::vector<int>> _blockSizes               = {{4, 4, 4, 4}},
                   std::vector<double>           _smootherTol              = {1e-14},
                   std::vector<int>              _smootherMaxOuterIter     = {4},
@@ -82,7 +82,7 @@ public:
 };
 // clang-format on
 
-void checkParameterValidity(MultiGridParams const &params) {
+void checkParameterValidity(WilsonMGParams const &params) {
 
   auto correctSize = params.nLevels - 1;
 
@@ -101,7 +101,7 @@ public:
   std::vector<GridCartesian *>  Grids;
   std::vector<GridParallelRNG>  PRNGs;
 
-  LevelInfo(GridCartesian *FineGrid, MultiGridParams const &mgParams) {
+  LevelInfo(GridCartesian *FineGrid, WilsonMGParams const &mgParams) {
 
     auto nCoarseLevels = mgParams.blockSizes.size();
 
@@ -176,7 +176,7 @@ public:
   int _CurrentLevel;
   int _NextCoarserLevel;
 
-  MultiGridParams &_MultiGridParams;
+  WilsonMGParams &_WilsonMGParams;
   LevelInfo &      _LevelInfo;
 
   FineDiracMatrix & _FineMatrix;
@@ -201,10 +201,10 @@ public:
   // Member Functions
   /////////////////////////////////////////////
 
-  MultiGridPreconditioner(MultiGridParams &mgParams, LevelInfo &LvlInfo, FineDiracMatrix &FineMat, FineDiracMatrix &SmootherMat)
+  MultiGridPreconditioner(WilsonMGParams &mgParams, LevelInfo &LvlInfo, FineDiracMatrix &FineMat, FineDiracMatrix &SmootherMat)
     : _CurrentLevel(mgParams.nLevels - (nCoarserLevels + 1)) // _Level = 0 corresponds to finest
     , _NextCoarserLevel(_CurrentLevel + 1)                   // incremented for instances on coarser levels
-    , _MultiGridParams(mgParams)
+    , _WilsonMGParams(mgParams)
     , _LevelInfo(LvlInfo)
     , _FineMatrix(FineMat)
     , _SmootherMatrix(SmootherMat)
@@ -212,7 +212,7 @@ public:
     , _CoarseMatrix(*_LevelInfo.Grids[_NextCoarserLevel]) {
 
     _NextPreconditionerLevel
-      = std::unique_ptr<NextPreconditionerLevel>(new NextPreconditionerLevel(_MultiGridParams, _LevelInfo, _CoarseMatrix, _CoarseMatrix));
+      = std::unique_ptr<NextPreconditionerLevel>(new NextPreconditionerLevel(_WilsonMGParams, _LevelInfo, _CoarseMatrix, _CoarseMatrix));
 
     resetTimers();
   }
@@ -261,7 +261,7 @@ public:
     conformable(in, out);
 
     // TODO: implement a W-cycle
-    if(_MultiGridParams.kCycle)
+    if(_WilsonMGParams.kCycle)
       kCycle(in, out);
     else
       vCycle(in, out);
@@ -279,13 +279,13 @@ public:
 
     FineVector fineTmp(in.Grid());
 
-    auto maxSmootherIter = _MultiGridParams.smootherMaxOuterIter[_CurrentLevel] * _MultiGridParams.smootherMaxInnerIter[_CurrentLevel];
+    auto maxSmootherIter = _WilsonMGParams.smootherMaxOuterIter[_CurrentLevel] * _WilsonMGParams.smootherMaxInnerIter[_CurrentLevel];
 
     TrivialPrecon<FineVector>                      fineTrivialPreconditioner;
-    FlexibleGeneralisedMinimalResidual<FineVector> fineFGMRES(_MultiGridParams.smootherTol[_CurrentLevel],
+    FlexibleGeneralisedMinimalResidual<FineVector> fineFGMRES(_WilsonMGParams.smootherTol[_CurrentLevel],
                                                               maxSmootherIter,
                                                               fineTrivialPreconditioner,
-                                                              _MultiGridParams.smootherMaxInnerIter[_CurrentLevel],
+                                                              _WilsonMGParams.smootherMaxInnerIter[_CurrentLevel],
                                                               false);
 
     MdagMLinearOperator<FineDiracMatrix, FineVector> fineMdagMOp(_FineMatrix);
@@ -336,19 +336,19 @@ public:
 
     FineVector fineTmp(in.Grid());
 
-    auto smootherMaxIter = _MultiGridParams.smootherMaxOuterIter[_CurrentLevel] * _MultiGridParams.smootherMaxInnerIter[_CurrentLevel];
-    auto kCycleMaxIter   = _MultiGridParams.kCycleMaxOuterIter[_CurrentLevel] * _MultiGridParams.kCycleMaxInnerIter[_CurrentLevel];
+    auto smootherMaxIter = _WilsonMGParams.smootherMaxOuterIter[_CurrentLevel] * _WilsonMGParams.smootherMaxInnerIter[_CurrentLevel];
+    auto kCycleMaxIter   = _WilsonMGParams.kCycleMaxOuterIter[_CurrentLevel] * _WilsonMGParams.kCycleMaxInnerIter[_CurrentLevel];
 
     TrivialPrecon<FineVector>                        fineTrivialPreconditioner;
-    FlexibleGeneralisedMinimalResidual<FineVector>   fineFGMRES(_MultiGridParams.smootherTol[_CurrentLevel],
+    FlexibleGeneralisedMinimalResidual<FineVector>   fineFGMRES(_WilsonMGParams.smootherTol[_CurrentLevel],
                                                               smootherMaxIter,
                                                               fineTrivialPreconditioner,
-                                                              _MultiGridParams.smootherMaxInnerIter[_CurrentLevel],
+                                                              _WilsonMGParams.smootherMaxInnerIter[_CurrentLevel],
                                                               false);
-    FlexibleGeneralisedMinimalResidual<CoarseVector> coarseFGMRES(_MultiGridParams.kCycleTol[_CurrentLevel],
+    FlexibleGeneralisedMinimalResidual<CoarseVector> coarseFGMRES(_WilsonMGParams.kCycleTol[_CurrentLevel],
                                                                   kCycleMaxIter,
                                                                   *_NextPreconditionerLevel,
-                                                                  _MultiGridParams.kCycleMaxInnerIter[_CurrentLevel],
+                                                                  _WilsonMGParams.kCycleMaxInnerIter[_CurrentLevel],
                                                                   false);
 
     MdagMLinearOperator<FineDiracMatrix, FineVector>     fineMdagMOp(_FineMatrix);
@@ -581,7 +581,7 @@ public:
 
   int _CurrentLevel;
 
-  MultiGridParams &_MultiGridParams;
+  WilsonMGParams &_WilsonMGParams;
   LevelInfo &      _LevelInfo;
 
   FineDiracMatrix &_FineMatrix;
@@ -594,9 +594,9 @@ public:
   // Member Functions
   /////////////////////////////////////////////
 
-  MultiGridPreconditioner(MultiGridParams &mgParams, LevelInfo &LvlInfo, FineDiracMatrix &FineMat, FineDiracMatrix &SmootherMat)
+  MultiGridPreconditioner(WilsonMGParams &mgParams, LevelInfo &LvlInfo, FineDiracMatrix &FineMat, FineDiracMatrix &SmootherMat)
     : _CurrentLevel(mgParams.nLevels - (0 + 1))
-    , _MultiGridParams(mgParams)
+    , _WilsonMGParams(mgParams)
     , _LevelInfo(LvlInfo)
     , _FineMatrix(FineMat)
     , _SmootherMatrix(SmootherMat) {
@@ -613,12 +613,12 @@ public:
     conformable(_LevelInfo.Grids[_CurrentLevel], in.Grid());
     conformable(in, out);
 
-    auto coarseSolverMaxIter = _MultiGridParams.coarseSolverMaxOuterIter * _MultiGridParams.coarseSolverMaxInnerIter;
+    auto coarseSolverMaxIter = _WilsonMGParams.coarseSolverMaxOuterIter * _WilsonMGParams.coarseSolverMaxInnerIter;
 
     // On the coarsest level we only have what I above call the fine level, no coarse one
     TrivialPrecon<FineVector>                      fineTrivialPreconditioner;
     FlexibleGeneralisedMinimalResidual<FineVector> fineFGMRES(
-      _MultiGridParams.coarseSolverTol, coarseSolverMaxIter, fineTrivialPreconditioner, _MultiGridParams.coarseSolverMaxInnerIter, false);
+      _WilsonMGParams.coarseSolverTol, coarseSolverMaxIter, fineTrivialPreconditioner, _WilsonMGParams.coarseSolverMaxInnerIter, false);
 
     MdagMLinearOperator<FineDiracMatrix, FineVector> fineMdagMOp(_FineMatrix);
 
@@ -651,7 +651,7 @@ using NLevelMGPreconditioner = MultiGridPreconditioner<Fobj, CComplex, nBasis, n
 
 template<class Fobj, class CComplex, int nBasis, class Matrix>
 std::unique_ptr<MultiGridPreconditionerBase<Lattice<Fobj>>>
-createMGInstance(MultiGridParams &mgParams, LevelInfo &levelInfo, Matrix &FineMat, Matrix &SmootherMat) {
+createMGInstance(WilsonMGParams &mgParams, LevelInfo &levelInfo, Matrix &FineMat, Matrix &SmootherMat) {
 
 #define CASE_FOR_N_LEVELS(nLevels)                                                                                     \
   case nLevels:                                                                                                        \

tests/solver/Test_wilson_mg.cc

@@ -51,7 +51,7 @@ int main(int argc, char **argv) {
 
   RealD mass = -0.25;
 
-  MultiGridParams mgParams;
+  WilsonMGParams mgParams;
   std::string     inputXml{"./mg_params.xml"};
 
   if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) {

tests/solver/Test_wilson_mg_mp.cc

@@ -58,7 +58,7 @@ int main(int argc, char **argv) {
 
   RealD mass = -0.25;
 
-  MultiGridParams mgParams;
+  WilsonMGParams mgParams;
   std::string     inputXml{"./mg_params.xml"};
 
   if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) {

tests/solver/Test_wilsonclover_mg.cc

@@ -54,7 +54,7 @@ int main(int argc, char **argv) {
   RealD csw_r = 1.0;
   RealD csw_t = 1.0;
 
-  MultiGridParams mgParams;
+  WilsonMGParams mgParams;
   std::string     inputXml{"./mg_params.xml"};
 
   if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) {

tests/solver/Test_wilsonclover_mg_lime.cc

@@ -84,7 +84,7 @@ int main(int argc, char **argv) {
   RealD csw_r = 1.0;
   RealD csw_t = 1.0;
 
-  MultiGridParams mgParams;
+  WilsonMGParams mgParams;
   std::string     inputXml{"./mg_params.xml"};
 
   if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) {

tests/solver/Test_wilsonclover_mg_mp.cc

@@ -60,7 +60,7 @@ int main(int argc, char **argv) {
   RealD csw_r = 1.0;
   RealD csw_t = 1.0;
 
-  MultiGridParams mgParams;
+  WilsonMGParams mgParams;
   std::string     inputXml{"./mg_params.xml"};
 
   if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) {