Submission scripts for distributed schur coarse grid inverse

Grid/algorithms/blas/BatchedBlas.h

@@ -991,6 +991,233 @@ public:
     RealD bytes = 1.0*sizeof(ComplexF)*(m*k+k*n+m*n)*batchCount;
   }
 
+  ///////////////////////////////////////////////////////////////////////////////////
+  // 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
+  // larger row-slab allocations.
+  ///////////////////////////////////////////////////////////////////////////////////
+  void gemmBatched(GridBLASOperation_t OpA,
+                   GridBLASOperation_t OpB,
+                   int m,int n, int k,
+                   ComplexD alpha,
+                   deviceVector<ComplexD*> &Amk, int lda,
+                   deviceVector<ComplexD*> &Bkn, int ldb,
+                   ComplexD beta,
+                   deviceVector<ComplexD*> &Cmn, int ldc)
+  {
+    RealD t2=usecond();
+    int32_t batchCount = Amk.size();
+
+    GRID_ASSERT( lda >= ((OpA==GridBLAS_OP_N) ? m : k) );
+    GRID_ASSERT( ldb >= ((OpB==GridBLAS_OP_N) ? k : n) );
+    GRID_ASSERT( ldc >= m );
+
+    // Cached device constants: copy only on value change (see GridBLASDeviceConstant)
+    static GridBLASDeviceConstant<ComplexD> alpha_c;
+    static GridBLASDeviceConstant<ComplexD> beta_c;
+    ComplexD *alpha_p = alpha_c.put(alpha);
+    ComplexD *beta_p  = beta_c.put(beta);
+    RealD t0=usecond();
+
+    GRID_ASSERT(Bkn.size()==batchCount);
+    GRID_ASSERT(Cmn.size()==batchCount);
+#ifdef GRID_HIP
+    hipblasOperation_t hOpA;
+    hipblasOperation_t hOpB;
+    if ( OpA == GridBLAS_OP_N ) hOpA = HIPBLAS_OP_N;
+    if ( OpA == GridBLAS_OP_T ) hOpA = HIPBLAS_OP_T;
+    if ( OpA == GridBLAS_OP_C ) hOpA = HIPBLAS_OP_C;
+    if ( OpB == GridBLAS_OP_N ) hOpB = HIPBLAS_OP_N;
+    if ( OpB == GridBLAS_OP_T ) hOpB = HIPBLAS_OP_T;
+    if ( OpB == GridBLAS_OP_C ) hOpB = HIPBLAS_OP_C;
+#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
+    auto err = hipblasZgemmBatched(gridblasHandle,
+                                   hOpA,
+                                   hOpB,
+                                   m,n,k,
+                                   (hipDoubleComplex *) &alpha_p[0],
+                                   (hipDoubleComplex **)&Amk[0], lda,
+                                   (hipDoubleComplex **)&Bkn[0], ldb,
+                                   (hipDoubleComplex *) &beta_p[0],
+                                   (hipDoubleComplex **)&Cmn[0], ldc,
+                                   batchCount);
+#else
+    auto err = hipblasZgemmBatched(gridblasHandle,
+                                   hOpA,
+                                   hOpB,
+                                   m,n,k,
+                                   (hipblasDoubleComplex *) &alpha_p[0],
+                                   (hipblasDoubleComplex **)&Amk[0], lda,
+                                   (hipblasDoubleComplex **)&Bkn[0], ldb,
+                                   (hipblasDoubleComplex *) &beta_p[0],
+                                   (hipblasDoubleComplex **)&Cmn[0], ldc,
+                                   batchCount);
+#endif
+    GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS);
+#endif
+#ifdef GRID_CUDA
+    cublasOperation_t hOpA;
+    cublasOperation_t hOpB;
+    if ( OpA == GridBLAS_OP_N ) hOpA = CUBLAS_OP_N;
+    if ( OpA == GridBLAS_OP_T ) hOpA = CUBLAS_OP_T;
+    if ( OpA == GridBLAS_OP_C ) hOpA = CUBLAS_OP_C;
+    if ( OpB == GridBLAS_OP_N ) hOpB = CUBLAS_OP_N;
+    if ( OpB == GridBLAS_OP_T ) hOpB = CUBLAS_OP_T;
+    if ( OpB == GridBLAS_OP_C ) hOpB = CUBLAS_OP_C;
+    auto err = cublasZgemmBatched(gridblasHandle,
+                                  hOpA,
+                                  hOpB,
+                                  m,n,k,
+                                  (cuDoubleComplex *) &alpha_p[0],
+                                  (cuDoubleComplex **)&Amk[0], lda,
+                                  (cuDoubleComplex **)&Bkn[0], ldb,
+                                  (cuDoubleComplex *) &beta_p[0],
+                                  (cuDoubleComplex **)&Cmn[0], ldc,
+                                  batchCount);
+    GRID_ASSERT(err==CUBLAS_STATUS_SUCCESS);
+#endif
+#ifdef GRID_SYCL
+    int64_t m64=m;
+    int64_t n64=n;
+    int64_t k64=k;
+    int64_t lda64=lda;
+    int64_t ldb64=ldb;
+    int64_t ldc64=ldc;
+    int64_t batchCount64=batchCount;
+
+    oneapi::mkl::transpose iOpA;
+    oneapi::mkl::transpose iOpB;
+
+    if ( OpA == GridBLAS_OP_N ) iOpA = oneapi::mkl::transpose::N;
+    if ( OpA == GridBLAS_OP_T ) iOpA = oneapi::mkl::transpose::T;
+    if ( OpA == GridBLAS_OP_C ) iOpA = oneapi::mkl::transpose::C;
+    if ( OpB == GridBLAS_OP_N ) iOpB = oneapi::mkl::transpose::N;
+    if ( OpB == GridBLAS_OP_T ) iOpB = oneapi::mkl::transpose::T;
+    if ( OpB == GridBLAS_OP_C ) iOpB = oneapi::mkl::transpose::C;
+
+    oneapi::mkl::blas::column_major::gemm_batch(*gridblasHandle,
+                                                &iOpA,
+                                                &iOpB,
+                                                &m64,&n64,&k64,
+                                                (ComplexD *) &alpha_p[0],
+                                                (const ComplexD **)&Amk[0], (const int64_t *)&lda64,
+                                                (const ComplexD **)&Bkn[0], (const int64_t *)&ldb64,
+                                                (ComplexD *) &beta_p[0],
+                                                (ComplexD **)&Cmn[0], (const int64_t *)&ldc64,
+                                                (int64_t)1,&batchCount64,std::vector<sycl::event>());
+    synchronise();
+#endif
+#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
+    // Reference implementation: Eigen with explicit outer stride
+    typedef Eigen::Map<Eigen::MatrixXcd,0,Eigen::OuterStride<> > eMat;
+    if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_N) ) {
+      thread_for (p, batchCount, {
+        eMat eAmk(Amk[p],m,k,Eigen::OuterStride<>(lda));
+        eMat eBkn(Bkn[p],k,n,Eigen::OuterStride<>(ldb));
+        eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+        if (std::abs(beta) != 0.0)
+        {
+          eCmn = beta * eCmn + alpha * eAmk * eBkn;
+        }
+        else
+        {
+          eCmn = alpha * eAmk * eBkn;
+        }
+      });
+    } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_N) ) {
+      thread_for (p, batchCount, {
+        eMat eAmk(Amk[p],k,m,Eigen::OuterStride<>(lda));
+        eMat eBkn(Bkn[p],k,n,Eigen::OuterStride<>(ldb));
+        eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+        if (std::abs(beta) != 0.0)
+        {
+          eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn;
+        }
+        else
+        {
+          eCmn = alpha * eAmk.adjoint() * eBkn;
+        }
+      });
+    } else if ( (OpA == GridBLAS_OP_T ) && (OpB == GridBLAS_OP_N) ) {
+      thread_for (p, batchCount, {
+        eMat eAmk(Amk[p],k,m,Eigen::OuterStride<>(lda));
+        eMat eBkn(Bkn[p],k,n,Eigen::OuterStride<>(ldb));
+        eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+        if (std::abs(beta) != 0.0)
+        {
+          eCmn = beta * eCmn + alpha * eAmk.transpose() * eBkn;
+        }
+        else
+        {
+          eCmn = alpha * eAmk.transpose() * eBkn;
+        }
+      });
+    } else if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_C) ) {
+      thread_for (p, batchCount, {
+        eMat eAmk(Amk[p],m,k,Eigen::OuterStride<>(lda));
+        eMat eBkn(Bkn[p],n,k,Eigen::OuterStride<>(ldb));
+        eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+        if (std::abs(beta) != 0.0)
+        {
+          eCmn = beta * eCmn + alpha * eAmk * eBkn.adjoint();
+        }
+        else
+        {
+          eCmn = alpha * eAmk * eBkn.adjoint();
+        }
+      });
+    } else if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_T) ) {
+      thread_for (p, batchCount, {
+        eMat eAmk(Amk[p],m,k,Eigen::OuterStride<>(lda));
+        eMat eBkn(Bkn[p],n,k,Eigen::OuterStride<>(ldb));
+        eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+        if (std::abs(beta) != 0.0)
+        {
+          eCmn = beta * eCmn + alpha * eAmk * eBkn.transpose();
+        }
+        else
+        {
+          eCmn = alpha * eAmk * eBkn.transpose();
+        }
+      });
+    } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_C) ) {
+      thread_for (p, batchCount, {
+        eMat eAmk(Amk[p],k,m,Eigen::OuterStride<>(lda));
+        eMat eBkn(Bkn[p],n,k,Eigen::OuterStride<>(ldb));
+        eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+        if (std::abs(beta) != 0.0)
+        {
+          eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn.adjoint();
+        }
+        else
+        {
+          eCmn = alpha * eAmk.adjoint() * eBkn.adjoint();
+        }
+      });
+    } else if ( (OpA == GridBLAS_OP_T ) && (OpB == GridBLAS_OP_T) ) {
+      thread_for (p, batchCount, {
+        eMat eAmk(Amk[p],k,m,Eigen::OuterStride<>(lda));
+        eMat eBkn(Bkn[p],n,k,Eigen::OuterStride<>(ldb));
+        eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+        if (std::abs(beta) != 0.0)
+        {
+          eCmn = beta * eCmn + alpha * eAmk.transpose() * eBkn.transpose();
+        }
+        else
+        {
+          eCmn = alpha * eAmk.transpose() * eBkn.transpose();
+        }
+      });
+    } else {
+      assert(0);
+    }
+#endif
+    RealD t1=usecond();
+    RealD flops = 8.0*m*n*k*batchCount;
+    RealD bytes = 1.0*sizeof(ComplexD)*(m*k+k*n+m*n)*batchCount;
+  }
+
   ///////////////////////////////////////////////////////////////////////////
   // Single precision real GEMM
   ///////////////////////////////////////////////////////////////////////////

Grid/algorithms/multigrid/DenseCoarseMatrix.h

@@ -29,6 +29,7 @@ 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 <unordered_map>
 
@@ -68,11 +69,14 @@ NAMESPACE_BEGIN(Grid);
 //
 // 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
 //
-// Eventual internal upgrade (unchanged surface): RecursiveSchur distributed
-// factorisation replacing the single-GCD gather/invert, lifting BOTH the fp32
-// N ~ 90k boss-HBM ceiling AND the CC-grid 256-rank SIMD cap; leaves land on
-// GridBLASInverse::inverseBatched.
+// 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.
 //
 // Tensor-depth agnostic: site scalar objects treated as contiguous ComplexD
 // (iScalar wrappers add no data), so any MG level's coarse operator imports.
@@ -115,6 +119,7 @@ public:
   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(GeneralCoarseOp &Op, GridBase *g)
     : _Op(Op), grid(g)
@@ -126,6 +131,7 @@ 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 << ")"
@@ -276,6 +282,23 @@ public:
       // extract the unpadded field before peeking with unpadded coordinates
       // (exactly as MultiGeneralCoarsenedMatrix::CopyMatrix does).
       CoarseMatrix Aun = _Op.Cell.Extract(_Op._A[p]);
+      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(_Op._A[p])
+                  << " norm2(Extract) " << norm2(Aun)
+                  << " sum|peek|^2 " << gpk << std::endl;
+      }
       autoView(Av, Aun, CpuRead);
       thread_for(ss, lsites, {
         Coordinate ncoor(nd);
@@ -288,16 +311,64 @@ public:
         Msobj m;
         peekLocalSite(m, Av, myLcoor[ss]);
         ComplexD *md = (ComplexD *)&m;
-        for(int a=0; a<nbasis; a++){
-          ComplexF *row = &slab[(uint64_t)(ss*nbasis+a)*N + nsite*nbasis];
-          for(int b=0; b<nbasis; b++)
-            row[b] += ComplexF(md[a*nbasis+b]);   // += : wrapped shifts may collide
+        // 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.
+        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++)
+            row[a] += ComplexF(md[a*nbasis+b]);   // += : wrapped shifts may collide
         }
       });
     }
     t += usecond();
+
+    // Structural diagnostic: identically-zero rows of my slab (a healthy
+    // coarse operator has none; dead rows mean a rank-deficient import
+    // or operator and the inverse will be NaN).
+    int64_t zrows = 0;
+    for(int64_t r=0; r<nrows; r++)
+    {
+      double mx = 0.0;
+      const ComplexF *row = &slab[(uint64_t)r*N];
+      for(int64_t j=0; j<N; j++)
+      {
+        mx = std::max(mx, (double)abs(row[j]));
+      }
+      if ( mx < 1.0e-30 ) zrows++;
+    }
+    RealD gz = (RealD)zrows;
+    grid->GlobalSumVector(&gz, 1);
+
     std::cout << GridLogMessage << "DenseCoarseMatrix: stencil->dense import took "
-              << t/1.0e6 << " s  (" << _Op.geom.npoint << " points, local, no comms)" << std::endl;
+              << t/1.0e6 << " s  (" << _Op.geom.npoint << " points, local, no comms)"
+              << "  zero rows " << (int64_t)gz << "/" << N << std::endl;
+
+    // Debug: coordinate pattern of live sites (mechanism fingerprint)
+    if ( (int64_t)gz > 0 )
+    {
+      int shown = 0;
+      for(int ss=0; ss<lsites && shown<24; ss++)
+      {
+        double mx = 0.0;
+        const ComplexF *row = &slab[(uint64_t)(ss*nbasis)*N];
+        for(int64_t j=0; j<N; j++)
+        {
+          mx = std::max(mx, (double)abs(row[j]));
+        }
+        if ( mx > 1.0e-30 )
+        {
+          std::cout << GridLogMessage << "DenseCoarseMatrix: LIVE site ss=" << ss
+                    << " lcoor " << myLcoor[ss] << std::endl;
+          shown++;
+        }
+      }
+    }
   }
 
   ////////////////////////////////////////////////////////////////////
@@ -349,12 +420,80 @@ public:
   }
 
   ////////////////////////////////////////////////////////////////////
-  // 3. 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).  Ported proven path from the frozen example.
+  // 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).
   ////////////////////////////////////////////////////////////////////
   void InvertDense(void)
+  {
+    char *sc  = getenv("DENSE_SCHUR");
+    int  mode = sc ? atoi(sc) : 0;
+
+    if ( mode == 0 )
+    {
+      InvertDenseSingle();
+      return;
+    }
+    if ( mode == 1 )
+    {
+      InvertDenseSchur();
+      return;
+    }
+    GRID_ASSERT( mode == 2 );
+    std::vector<ComplexF> Aimp(slab);       // imported A
+    InvertDenseSingle();
+    std::vector<ComplexF> ref(slab);        // Ainv, single path
+    slab = Aimp;
+    InvertDenseSchur();                     // 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();
@@ -514,6 +653,194 @@ public:
               << (t4-t3)/1.0e6 << " s" << std::endl;
   }
 
+  ////////////////////////////////////////////////////////////////////
+  // 3b. 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).
+  //    MPI_Cart_rank is queried ONCE PER RANK (serial, P calls) into a
+  //    lex-processor table; the per-site sweep is then pure arithmetic.
+  ////////////////////////////////////////////////////////////////////
+  void BuildRankMajorMap(std::vector<int64_t> &g2rm)
+  {
+    int P = grid->ProcessorCount();
+    Coordinate pdims = grid->_processors;
+    Coordinate gdims = grid->GlobalDimensions();
+    Coordinate ldims = grid->LocalDimensions();
+
+    std::vector<int> lexp2rank(P);
+    for(int lp=0; lp<P; lp++)
+    {
+      Coordinate pcoor(nd);
+      Lexicographic::CoorFromIndex(pcoor, lp, pdims);
+      lexp2rank[lp] = grid->RankFromProcessorCoor(pcoor);
+    }
+
+    int64_t gsites = grid->gSites();
+    g2rm.resize(N);
+    thread_for(gsite, gsites, {
+      Coordinate gcoor(nd);
+      Coordinate pcoor(nd);
+      Coordinate lcoor(nd);
+      Lexicographic::CoorFromIndex(gcoor, gsite, gdims);
+      for(int d=0; d<nd; d++)
+      {
+        pcoor[d] = gcoor[d]/ldims[d];
+        lcoor[d] = gcoor[d]-pcoor[d]*ldims[d];
+      }
+      int64_t lexp;
+      int64_t lsite;
+      Lexicographic::IndexFromCoor(pcoor, lexp,  pdims);
+      Lexicographic::IndexFromCoor(lcoor, lsite, ldims);
+      int64_t base = (int64_t)lexp2rank[lexp]*nrows + lsite*nbasis;
+      for(int b=0; b<nbasis; b++)
+      {
+        g2rm[(uint64_t)gsite*nbasis + b] = base + b;
+      }
+    });
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // 3c. 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
+  //    ImportDense; output is column-major rows x N with columns in
+  //    rank-major order (g2rm).
+  //    ALWAYS-ON CERTIFICATE: the fp64 import, rounded, must agree with
+  //    the fp32 slab entry at the corresponding global column, over the
+  //    WHOLE of my rows (few ulp: wrapped-shift collisions accumulate in
+  //    different precision order).  NaN-proof: non-finite entries are
+  //    counted explicitly since max() silently masks NaN.
+  ////////////////////////////////////////////////////////////////////
+  void ImportDenseFP64(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.geom.npoint; p++)
+    {
+      Coordinate shift = _Op.geom.shifts[p];
+      CoarseMatrix Aun = _Op.Cell.Extract(_Op._A[p]);
+      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];
+          ncoor[d] = (int)((g % gdims[d] + gdims[d]) % gdims[d]);
+        }
+        int64_t nsite;
+        Lexicographic::IndexFromCoor(ncoor, nsite, gdims);
+        Msobj m;
+        peekLocalSite(m, Av, myLcoor[ss]);
+        ComplexD *md = (ComplexD *)&m;
+        // Transposed contraction as ImportDense: (a,b) lands at
+        // row (s,b), column (nbr,a); column index in rank-major order.
+        for(int a=0; a<nbasis; a++)
+        {
+          int64_t jj = g2rm[ nsite*nbasis + a ];
+          for(int b=0; b<nbasis; b++)
+          {
+            h[(uint64_t)(ss*nbasis+b) + (uint64_t)jj*nrows] += md[a*nbasis+b];
+          }
+        }
+      });
+    }
+
+    // Certificate vs the fp32 slab (slab holds A at this point)
+    double  mx   = 0.0;
+    int64_t nbad = 0;
+    for(int64_t i=0; i<nrows; i++)
+    {
+      for(int64_t gcol=0; gcol<N; gcol++)
+      {
+        ComplexD d64 = h[(uint64_t)(i + g2rm[gcol]*nrows)];
+        ComplexF f32 = slab[(uint64_t)i*N + gcol];
+        double dev = abs(ComplexD(f32) - d64);
+        if ( !std::isfinite(dev) ) nbad++;
+        else mx = std::max(mx, dev);
+      }
+    }
+    RealD gmx  = mx;
+    RealD gbad = (RealD)nbad;
+    grid->GlobalMax(gmx);
+    grid->GlobalSumVector(&gbad, 1);
+    std::cout << GridLogMessage << "DenseCoarseMatrix: fp64 import certificate "
+              << "max|A64 - A32| = " << gmx
+              << "  non-finite entries " << (int64_t)gbad << std::endl;
+    GRID_ASSERT( gbad == 0 );
+    GRID_ASSERT( gmx < 1.0e-5 );
+
+    S.Resize(nrows, N);
+    acceleratorCopyToDevice(&h[0], &S.data[0], (uint64_t)nrows*N*sizeof(ComplexD));
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // 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.
+  ////////////////////////////////////////////////////////////////////
+  void InvertDenseSchur(void)
+  {
+    double t1 = usecond();
+    int P  = grid->ProcessorCount();
+    int me = grid->ThisRank();
+
+    // Uniform local volumes => contiguous uniform ownership
+    std::vector<int64_t> rowStart(P+1);
+    for(int r=0; r<=P; r++)
+    {
+      rowStart[r] = (int64_t)r*nrows;
+    }
+    GRID_ASSERT( rowStart[P] == N );
+
+    std::vector<int64_t> g2rm;
+    BuildRankMajorMap(g2rm);
+
+    // Self-certifying map: my own global rows land at my rank-major slots
+    for(int ss=0; ss<lsites; ss++)
+    {
+      for(int a=0; a<nbasis; a++)
+      {
+        GRID_ASSERT( g2rm[ myGsite[ss]*nbasis + a ] == (int64_t)me*nrows + ss*nbasis + a );
+      }
+    }
+
+    BlockRows S;
+    ImportDenseFP64(S, g2rm);
+
+    int64_t panelBytes = getenv("DENSE_PANEL_BYTES") ? atol(getenv("DENSE_PANEL_BYTES"))
+                                                     : (int64_t)1024*1024*1024;
+    RecursiveSchurInverse RSI(grid, N, rowStart, panelBytes);
+    double t2 = usecond();
+    RSI.Invert(S);
+    double t3 = usecond();
+    RSI.ReportTelemetry();
+
+    // The single terminal rounding: fp64 inverse -> fp32 apply slab
+    // (row-major, global columns)
+    {
+      std::vector<ComplexD> h((uint64_t)nrows*N);
+      acceleratorCopyFromDevice(&S.data[0], &h[0], (uint64_t)nrows*N*sizeof(ComplexD));
+      thread_for(gcol, N, {
+        int64_t jj = g2rm[gcol];
+        for(int64_t i=0; i<nrows; i++)
+        {
+          slab[(uint64_t)i*N + gcol] = ComplexF(h[(uint64_t)(i + jj*nrows)]);
+        }
+      });
+    }
+    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;
+  }
+
   ////////////////////////////////////////////////////////////////////
   // CORE apply on packed data: hX[N x nr] (zero-filled, my sites only)
   // -> allreduce -> split-K GEMM against the resident slab -> reduce

Grid/algorithms/multigrid/GeneralCoarsenedMatrix.h

@@ -453,8 +453,18 @@ public:
 	autoView( A_v  , _A[k], AcceleratorWrite);
 	autoView( FT_v  , FT[k], AcceleratorRead);
 	accelerator_for(sss, osites, nbasis, {
+#ifdef GRID_SIMT
 	    int j = acceleratorSIMTlane(nbasis);
 	    A_v[sss](i,j) = FT_v[sss](j);
+#else
+	    // CPU build: acceleratorSIMTlane()==0 -- an un-looped SIMT tensor
+	    // index writes ONLY j=0 and silently drops the other nbasis-1
+	    // columns (caught by Test_schur_dense_coarse import certificate,
+	    // 2026-08-14).  Loop explicitly.
+	    for(int j=0;j<nbasis;j++){
+	      A_v[sss](i,j) = FT_v[sss](j);
+	    }
+#endif
         });
       }
       tinv+=usecond();
@@ -630,8 +640,18 @@ public:
 	autoView( A_v  , _A[k], AcceleratorWrite);
 	autoView( FT_v  , FT[k], AcceleratorRead);
 	accelerator_for(sss, osites, nbasis, {
+#ifdef GRID_SIMT
 	    int j = acceleratorSIMTlane(nbasis);
 	    A_v[sss](i,j) = FT_v[sss](j);
+#else
+	    // CPU build: acceleratorSIMTlane()==0 -- an un-looped SIMT tensor
+	    // index writes ONLY j=0 and silently drops the other nbasis-1
+	    // columns (caught by Test_schur_dense_coarse import certificate,
+	    // 2026-08-14).  Loop explicitly.
+	    for(int j=0;j<nbasis;j++){
+	      A_v[sss](i,j) = FT_v[sss](j);
+	    }
+#endif
         });
       }
       tinv+=usecond();

Grid/algorithms/multigrid/RecursiveSchurInverse.h

@@ -0,0 +1,568 @@
+/*************************************************************************************
+
+    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>
+
+NAMESPACE_BEGIN(Grid);
+
+///////////////////////////////////////////////////////////////////////////////
+// RecursiveSchurInverse: distributed dense inversion by recursive Schur
+// complement over a binary rank-range tree.
+//
+// CONTRACT: the caller presents an N x N matrix in RANK-MAJOR row ordering,
+// distributed by rows -- rank r owns global rows [rowStart[r], rowStart[r+1])
+// -- and receives its rows of the INVERSE in the same layout.  This class
+// knows nothing of lattices or coarse operators; it consumes a GridBase for
+// world collectives, GridBLAS for GEMMs and GridBLASInverse for the leaf
+// inversions (all of which have Eigen reference backends, so the whole
+// algorithm unit-tests on a CPU-only laptop build under mpirun).
+//
+// PRECISION (decision 2026-08-14, superseding the fp32-merge design): the
+// ENTIRE inversion runs in fp64 (ComplexD).  The apply-side fp32 gain is
+// taken where it matters -- inside the iterative process -- by rounding the
+// finished inverse ONCE when the caller stores it in the fp32 apply slab.
+// Consequences: merge-growth error accumulates in eps64 and the terminal
+// rounding gives representation-only ~eps32 accuracy independent of growth;
+// the Newton-Schulz refinement and the fp32 escalation ladder are DELETED
+// (resurrectable from git history if a future scale forces reduced-precision
+// merges).  Setup cost: ~2x panel-gather bytes and ~2x transient memory,
+// once per setup; fp64 GEMM runs at fp32 rate on CDNA2/PVC.
+//
+// EXECUTION MODEL: SPMD full-tree walk.  Every rank executes the identical
+// recursion call sequence; participation in DATA is ownership-gated, and
+// every collective is a world-communicator zero-fill GlobalSumVector.  No
+// sub-communicators exist, so no deadlock surface exists.
+//
+// STORAGE CONVENTION (pinned by unit test T1b, Test_schur_inverse.cc):
+// BlockRows is COLUMN-MAJOR with ld = rows, matching the BLAS world:
+// element (i,j) lives at data[ i + j*ld ]; a column window [col0, col0+w)
+// is the contiguous slice starting 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 (diagnostic, not load-bearing at fp64): 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
+
+  ///////////////////////////////////////////////////////////////////////////
+  // 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;
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // THE communication primitive (plan 3.4 / 4B.3).
+  //
+  //   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].  The sub-block is gathered in panelBytes
+  // row-chunks by host zero-fill + world GlobalSumVector.
+  //
+  // SPMD rules: EVERY rank calls (the collectives are world-wide);
+  // non-owners of B add zeros; ranks with A.rows == 0 skip all local
+  // compute but still make every collective call.  Column offsets are
+  // LOCAL buffer offsets -- non-participants pass 0.
+  //
+  // Owners stage their whole B window device->host ONCE (ld == rows makes
+  // the window contiguous); per-chunk deposits are host memcpy runs.
+  ///////////////////////////////////////////////////////////////////////////
+  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 );
+    }
+
+    std::vector<ComplexD> stage;
+    if ( owner )
+    {
+      GRID_ASSERT( colB + widthB <= B.cols );
+      stage.resize((uint64_t)B.rows*n);
+      acceleratorCopyFromDevice(B.ColumnWindow(colB), &stage[0],
+                                (uint64_t)B.rows*n*sizeof(ComplexD));
+    }
+
+    int64_t kc = panelBytes / ( (int64_t)sizeof(ComplexD) * n );
+    if ( kc < 1 ) kc = 1;
+    if ( kc > k ) kc = k;
+    GRID_ASSERT( kc*n < 2147483647L );   // GlobalSumVector count is int
+
+    std::vector<ComplexD>   panel((uint64_t)kc*n);
+    deviceVector<ComplexD>  dPanel((uint64_t)kc*n);
+    deviceVector<ComplexD*> ap(1);
+    deviceVector<ComplexD*> bp(1);
+    deviceVector<ComplexD*> cp(1);
+    std::vector<ComplexD*>  ptr(1);
+
+    for(int64_t k0=0; k0<k; k0+=kc)
+    {
+      int64_t kchunk = std::min(kc, k-k0);
+
+      // PLANNED OPTIMISATION (not yet): single-threaded memset zero-fills
+      // the WHOLE panel; owners then overwrite their segment.  A threaded
+      // zero of only the non-owned rows (thread_for over columns, memset
+      // per column run) halves the host traffic and parallelises it.
+      // Deliberately deferred until the simple version is proven.
+      memset(&panel[0], 0, (uint64_t)kchunk*n*sizeof(ComplexD));
+      if ( owner )
+      {
+        int64_t i0 = std::max(k0,        myOff);
+        int64_t i1 = std::min(k0+kchunk, myOff+B.rows);
+        if ( i1 > i0 )
+        {
+          int64_t len = i1-i0;
+          thread_for(j, n, {
+            memcpy(&panel[(uint64_t)((i0-k0)    + j*kchunk)],
+                   &stage[(uint64_t)((i0-myOff) + j*B.rows)],
+                   len*sizeof(ComplexD));
+          });
+        }
+      }
+      grid->GlobalSumVector(&panel[0], (int)(kchunk*n));
+
+      if ( m > 0 )
+      {
+        acceleratorCopyToDevice(&panel[0], &dPanel[0],
+                                (uint64_t)kchunk*n*sizeof(ComplexD));
+
+        ComplexD beta_use = ( k0==0 ) ? beta : ComplexD(1.0,0.0);
+
+        ptr[0] = A.ColumnWindow(colA + k0);
+        acceleratorCopyToDevice(&ptr[0], &ap[0], sizeof(ComplexD*));
+        ptr[0] = &dPanel[0];
+        acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexD*));
+        ptr[0] = C.ColumnWindow(colC);
+        acceleratorCopyToDevice(&ptr[0], &cp[0], sizeof(ComplexD*));
+
+        BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N,
+                         (int)m, (int)n, (int)kchunk,
+                         alpha,    ap, (int)A.ld,
+                                   bp, (int)kchunk,
+                         beta_use, cp, (int)C.ld);
+        BLAS.synchronise();
+      }
+    }
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // 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.  Purely LOCAL -- the calling rank owns the whole
+  // width x width leaf (width == my row count); no collectives, so the
+  // SPMD walk stays uniform with other ranks doing nothing.  The window
+  // is contiguous (ld == rows == width): invert IN PLACE via
+  // GridBLASInverse.  Everything is already fp64; no promote/demote.
+  ///////////////////////////////////////////////////////////////////////////
+  void LeafInvert(int64_t col0, int64_t width, BlockRows &Arows)
+  {
+    GRID_ASSERT( width == Arows.rows );
+    GRID_ASSERT( col0 + width <= Arows.cols );
+    int64_t  w   = width;
+    uint64_t len = (uint64_t)w*w;
+
+    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));
+    }
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // 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;
+
+    SchurNode(0, P, 0, N, Arows);
+
+    RealD mx = telLeafMaxInv;
+    grid->GlobalMax(mx);
+    telLeafMaxInv = mx;
+  }
+
+  // All telemetry values are globally reduced: safe to stream on every
+  // rank (Grid quiesces stdout to the boss unless --debug-stdout).
+  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;
+  }
+};
+
+NAMESPACE_END(Grid);

systems/Frontier/schur_f1.slurm

@@ -0,0 +1,82 @@
+#!/bin/bash -l
+
+#SBATCH --job-name=schurF1
+#SBATCH --nodes=1
+#SBATCH --ntasks-per-node=8
+#SBATCH --cpus-per-task=7
+#SBATCH --gpus-per-node=8
+#SBATCH --time=00:15:00
+#SBATCH --account=phy157_dwf
+#SBATCH --gpu-bind=none
+#SBATCH --exclusive
+#SBATCH --mem=0
+#SBATCH -q debug
+
+#####################################################################
+# F1 (schur_recursive_inverse_plan.txt 4C): RecursiveSchurInverse on
+# GPU + GPU-aware comms -- no lattice physics, communicator only.
+# Laptop-proven semantics (Eigen oracle, ranks 1-4); this run
+# adjudicates ONLY the device paths: rocSOLVER zgetri leaves, ComplexD
+# strided gemmBatched merges, device panels, thrust-complex NegateCopy
+# kernel, multi-GCD GlobalSumVector panels.
+#
+# PRE-REGISTERED: every stage PASS; T4 certificates ~1e-14..1e-10
+# (eps64 growth model, as laptop); at P=1-equivalent leaves rocSOLVER
+# != Eigen so |X-Xref| = exact 0 is NOT expected, <=1e-12 is healthy.
+#####################################################################
+
+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]}
+unset ROCR_VISIBLE_DEVICES
+export HIP_VISIBLE_DEVICES=\$GPU
+CMD="numactl -m \$NUMA -N \$NUMA  \$*"
+exec \$CMD
+EOF
+chmod +x ./select_gpu
+
+# NUMA_MAP is only valid at 8 ranks/node (see aggregate_io.slurm);
+# below 8 use the GCD-only wrapper and leave placement to Slurm.
+cat << EOF > select_gcd
+#!/bin/bash
+export GPU_MAP=(0 1 2 3 7 6 5 4)
+export GPU=\${GPU_MAP[\$SLURM_LOCALID]}
+unset ROCR_VISIBLE_DEVICES
+export HIP_VISIBLE_DEVICES=\$GPU
+exec \$*
+EOF
+chmod +x ./select_gcd
+SELECT=$SLURM_SUBMIT_DIR/select_gpu
+SELECT_NONUMA=$SLURM_SUBMIT_DIR/select_gcd
+
+root=$HOME/ParallelIO/systems/Frontier/
+source ${root}/sourceme-rocm7.2.sh
+
+export OMP_NUM_THREADS=7
+export MPICH_GPU_SUPPORT_ENABLED=1
+
+BIN=$SLURM_SUBMIT_DIR/tests/debug/Test_schur_inverse
+
+run () {   # run <ranks-per-node> <grid> <mpi> <comment>
+  local rpn=$1 gr=$2 mp=$3 note=$4 wrap=$SELECT
+  if [ $rpn -ne 8 ]; then wrap=$SELECT_NONUMA; fi
+  echo
+  echo "==================================================================="
+  echo "=== F1  n=$rpn  grid=$gr  mpi=$mp   $note"
+  echo "==================================================================="
+  srun -N1 -n$rpn --ntasks-per-node=$rpn $wrap $BIN --grid $gr --mpi $mp \
+       --accelerator-threads 8 --shm 2048 --device-mem 32000
+  echo "=== exit $?"
+}
+
+# Rank-count sweep: even splits, the uneven n=3 case, and the full node.
+run 2 16.16.16.16 1.1.1.2 "even split, 2 GCD"
+run 3 16.16.16.24 1.1.1.3 "UNEVEN split, 3 GCD"
+run 4 16.16.16.16 1.1.1.4 "even split, 4 GCD"
+run 8 16.16.16.16 1.1.1.8 "full node, 8 GCD, depth-3 tree"
+
+echo
+echo "=== F1 done"

systems/Frontier/schur_f2.slurm

@@ -0,0 +1,92 @@
+#!/bin/bash -l
+
+#SBATCH --job-name=schurF2
+#SBATCH --nodes=1
+#SBATCH --ntasks-per-node=8
+#SBATCH --cpus-per-task=7
+#SBATCH --gpus-per-node=8
+#SBATCH --time=00:30:00
+#SBATCH --account=phy157_dwf
+#SBATCH --gpu-bind=none
+#SBATCH --exclusive
+#SBATCH --mem=0
+#SBATCH -q debug
+
+#####################################################################
+# F2 (schur_recursive_inverse_plan.txt 4C): the full DenseCoarseMatrix
+# glue on GPU with a REAL coarsening -- DWF MdagM + 0.5 shift, random
+# aggregation nbasis=4, per-dim blocks {4,4,3,3}, coarse 3.3.4.4/rank
+# (N = 576 * nranks).  The test binary FORCES DENSE_SCHUR=2 audit mode
+# and DENSE_PANEL_BYTES=65536 internally (setenv in main), so both
+# invert paths run and are compared on the same imported A; multi-chunk
+# panel gathers are exercised through the glue.
+#
+# First GPU execution in anger of: GRID_SIMT-guarded CoarsenOperator
+# population kernel; fresh stencil->dense import (fp32 + fp64) with the
+# transposed-contraction fix; fp64 rank-major import; distributed fp64
+# recursion; single terminal fp32 rounding; boss rocSOLVER single-GCD
+# path as the audit oracle.
+#
+# PRE-REGISTERED (laptop values, expect same magnitudes): operator
+# health herm-dev ~1e-15; import certificate ~1e-8; fp64 import
+# certificate ~5e-7 few-ulp; AUDIT rel ~3e-7 (pure fp32 rounding);
+# VERIFY ~1.6e-7; T6 ALL PASS.
+#####################################################################
+
+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]}
+unset ROCR_VISIBLE_DEVICES
+export HIP_VISIBLE_DEVICES=\$GPU
+CMD="numactl -m \$NUMA -N \$NUMA  \$*"
+exec \$CMD
+EOF
+chmod +x ./select_gpu
+
+# NUMA_MAP is only valid at 8 ranks/node (see aggregate_io.slurm);
+# below 8 use the GCD-only wrapper and leave placement to Slurm.
+cat << EOF > select_gcd
+#!/bin/bash
+export GPU_MAP=(0 1 2 3 7 6 5 4)
+export GPU=\${GPU_MAP[\$SLURM_LOCALID]}
+unset ROCR_VISIBLE_DEVICES
+export HIP_VISIBLE_DEVICES=\$GPU
+exec \$*
+EOF
+chmod +x ./select_gcd
+SELECT=$SLURM_SUBMIT_DIR/select_gpu
+SELECT_NONUMA=$SLURM_SUBMIT_DIR/select_gcd
+
+root=$HOME/ParallelIO/systems/Frontier/
+source ${root}/sourceme-rocm7.2.sh
+
+export OMP_NUM_THREADS=7
+export MPICH_GPU_SUPPORT_ENABLED=1
+
+BIN=$SLURM_SUBMIT_DIR/tests/debug/Test_schur_dense_coarse
+
+run () {   # run <ranks-per-node> <grid> <mpi> <comment>
+  local rpn=$1 gr=$2 mp=$3 note=$4 wrap=$SELECT
+  if [ $rpn -ne 8 ]; then wrap=$SELECT_NONUMA; fi
+  echo
+  echo "==================================================================="
+  echo "=== F2  n=$rpn  grid=$gr  mpi=$mp   $note"
+  echo "==================================================================="
+  srun -N1 -n$rpn --ntasks-per-node=$rpn $wrap $BIN --grid $gr --mpi $mp \
+       --accelerator-threads 8 --shm 2048 --device-mem 32000
+  echo "=== exit $?"
+}
+
+# Uniform local fine volume 12.12.12.12 per rank (t = 12 * nranks).
+# n=2 covers the padded-Extract branch; n=8 the full node, N=4608.
+# The n=3 uneven case is laptop-proven; commented to keep the debug
+# slot short -- uncomment for the full sweep.
+run 2 12.12.12.24 1.1.1.2 "2 GCD, padded Extract on t, N=1152"
+#run 3 12.12.12.36 1.1.1.3 "3 GCD, UNEVEN, N=1728"
+run 8 12.12.12.96 1.1.1.8 "full node, 8 GCD, N=4608"
+
+echo
+echo "=== F2 done"

tests/debug/Test_batched_blas.cc

@@ -250,6 +250,84 @@ int main (int argc, char ** argv)
     GRID_ASSERT(rel < 1.0e-4);
   }
 
+  ////////////////////////////////////////////////////////////////
+  // T3 : ComplexD explicit-ld gemmBatched (the RecursiveSchurInverse
+  // merge primitive): same split-K miniature as T2, double precision.
+  // On device builds this is the FIRST exercise of hipblasZ/cublasZ
+  // gemmBatched through the strided overload.
+  ////////////////////////////////////////////////////////////////
+  {
+    const int nrows = 8;
+    const int nrhs  = 4;
+    const int K     = 256;
+    const int NK    = 8;        // split-K chunks
+    const int Kc    = K/NK;
+
+    std::vector<ComplexD> Ahost((uint64_t)K*nrows);
+    std::vector<ComplexD> Xhost((uint64_t)K*nrhs);
+    for(auto &z : Ahost) z = ComplexD(dist(rng),dist(rng));
+    for(auto &z : Xhost) z = ComplexD(dist(rng),dist(rng));
+
+    deviceVector<ComplexD> Adev(Ahost.size());
+    deviceVector<ComplexD> Xdev(Xhost.size());
+    deviceVector<ComplexD> Yref((uint64_t)nrows*nrhs);
+    deviceVector<ComplexD> Ypart((uint64_t)NK*nrows*nrhs);
+    acceleratorCopyToDevice(&Ahost[0],&Adev[0],Ahost.size()*sizeof(ComplexD));
+    acceleratorCopyToDevice(&Xhost[0],&Xdev[0],Xhost.size()*sizeof(ComplexD));
+
+    ComplexD one (1.0,0.0);
+    ComplexD zero(0.0,0.0);
+
+    // (a) reference: one compact batch=1 call (compact lda == K for OP_T)
+    {
+      deviceVector<ComplexD*> Ap(1), Xp(1), Yp(1);
+      std::vector<ComplexD*> h(1);
+      h[0]=&Adev[0]; acceleratorCopyToDevice(&h[0],&Ap[0],sizeof(ComplexD*));
+      h[0]=&Xdev[0]; acceleratorCopyToDevice(&h[0],&Xp[0],sizeof(ComplexD*));
+      h[0]=&Yref[0]; acceleratorCopyToDevice(&h[0],&Yp[0],sizeof(ComplexD*));
+      blas.gemmBatched(GridBLAS_OP_T,GridBLAS_OP_N,
+                       nrows,nrhs,K,
+                       one, Ap, Xp, zero, Yp);
+      blas.synchronise();
+    }
+
+    // (b) split-K: NK chunk-pointers into the SAME allocations, lda/ldb = K
+    {
+      deviceVector<ComplexD*> Ap(NK), Xp(NK), Yp(NK);
+      std::vector<ComplexD*> h(NK);
+      for(int j=0;j<NK;j++) h[j] = &Adev[(uint64_t)j*Kc];        // K-offset slice
+      acceleratorCopyToDevice(&h[0],&Ap[0],NK*sizeof(ComplexD*));
+      for(int j=0;j<NK;j++) h[j] = &Xdev[(uint64_t)j*Kc];
+      acceleratorCopyToDevice(&h[0],&Xp[0],NK*sizeof(ComplexD*));
+      for(int j=0;j<NK;j++) h[j] = &Ypart[(uint64_t)j*nrows*nrhs];
+      acceleratorCopyToDevice(&h[0],&Yp[0],NK*sizeof(ComplexD*));
+      blas.gemmBatched(GridBLAS_OP_T,GridBLAS_OP_N,
+                       nrows,nrhs,Kc,
+                       one, Ap, /*lda*/ K,
+                            Xp, /*ldb*/ K,
+                       zero,Yp, /*ldc*/ nrows);
+      blas.synchronise();
+    }
+
+    std::vector<ComplexD> Yref_h((uint64_t)nrows*nrhs);
+    std::vector<ComplexD> Ypart_h((uint64_t)NK*nrows*nrhs);
+    acceleratorCopyFromDevice(&Yref[0], &Yref_h[0], Yref_h.size()*sizeof(ComplexD));
+    acceleratorCopyFromDevice(&Ypart[0],&Ypart_h[0],Ypart_h.size()*sizeof(ComplexD));
+
+    double maxdev = 0.0;
+    double maxval = 0.0;
+    for(int i=0;i<nrows*nrhs;i++){
+      ComplexD sum(0.0,0.0);
+      for(int j=0;j<NK;j++) sum = sum + Ypart_h[(uint64_t)j*nrows*nrhs + i];
+      maxdev = std::max(maxdev,(double)abs(sum-Yref_h[i]));
+      maxval = std::max(maxval,(double)abs(Yref_h[i]));
+    }
+    double rel = maxdev/maxval;
+    std::cout << GridLogMessage << "T3  split-K strided gemmBatched ComplexD  max rel dev vs compact = "
+              << rel << ( rel < 1.0e-13 ? "   PASS" : "   FAIL" ) << std::endl;
+    GRID_ASSERT(rel < 1.0e-13);
+  }
+
   std::cout << GridLogMessage << "All batched-blas tests PASSED" << std::endl;
 
   Grid_finalize();

tests/debug/Test_schur_dense_coarse.cc

@@ -0,0 +1,291 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: Test_schur_dense_coarse.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//
+// 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.
+//
+// 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.
+//
+// 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):
+//   mpirun -n 1 ./Test_schur_dense_coarse --grid 12.12.12.12 --mpi 1.1.1.1
+//   mpirun -n 2 ./Test_schur_dense_coarse --grid 12.12.12.24 --mpi 1.1.1.2
+//   mpirun -n 3 ./Test_schur_dense_coarse --grid 12.12.12.36 --mpi 1.1.1.3
+//   mpirun -n 4 ./Test_schur_dense_coarse --grid 12.12.12.48 --mpi 1.1.1.4
+//
+#include <Grid/Grid.h>
+#include <Grid/lattice/PaddedCell.h>
+#include <Grid/stencil/GeneralLocalStencil.h>
+#include <Grid/algorithms/multigrid/DenseCoarseMatrix.h>
+
+using namespace std;
+using namespace Grid;
+
+///////////////////////////////////////////////////////////////////////
+// MdagM + shift: Galerkin projection of a PD operator plus sigma I is
+// safely invertible whatever the (random) subspace quality.
+///////////////////////////////////////////////////////////////////////
+template<class Field>
+class ShiftedHermOpAdaptor : public LinearOperatorBase<Field>
+{
+  LinearOperatorBase<Field> &wrapped;
+  RealD                      shift;
+public:
+  ShiftedHermOpAdaptor(LinearOperatorBase<Field> &wrapme, RealD s)
+    : wrapped(wrapme), shift(s) {};
+
+  void Op(const Field &in, Field &out)
+  {
+    wrapped.HermOp(in, out);
+    out = out + shift*in;
+  }
+  void AdjOp(const Field &in, Field &out)
+  {
+    Op(in, out);
+  }
+  void HermOp(const Field &in, Field &out)
+  {
+    Op(in, out);
+  }
+  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 HermOpAndNorm(const Field &in, Field &out, RealD &n1, RealD &n2) { GRID_ASSERT(0); }
+};
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int Ls     = 4;
+  const int nbasis = 4;
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+                                                                   GridDefaultSimd(Nd,vComplex::Nsimd()),
+                                                                   GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  // Per-dimension blocking {4,4,3,3}: fine 12.12.12.12 -> coarse
+  // 3.3.4.4.  Two constraints meet here (both MEASURED today):
+  //  - coarse dims of 2 hit the probing pathology (health probe below)
+  //  - GEN-simd lanes {1,1,2,2} must land on even coarse dims, so the
+  //    odd production-like 3s go on the lane-free x,y axes (exactly the
+  //    production [3,6,8,8] trick).
+  Coordinate blocks({4,4,3,3});
+  Coordinate clatt = GridDefaultLatt();
+  for(int d=0; d<clatt.size(); d++)
+  {
+    GRID_ASSERT( (clatt[d] % blocks[d]) == 0 );
+    clatt[d] = clatt[d]/blocks[d];
+  }
+  GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt,
+                                                           GridDefaultSimd(Nd,vComplex::Nsimd()),
+                                                           GridDefaultMpi());
+  GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
+
+  std::vector<int> seeds4({1,2,3,4});
+  std::vector<int> seeds5({5,6,7,8});
+  std::vector<int> cseeds({9,10,11,12});
+  GridParallelRNG RNG4(UGrid);    RNG4.SeedFixedIntegers(seeds4);
+  GridParallelRNG RNG5(FGrid);    RNG5.SeedFixedIntegers(seeds5);
+  GridParallelRNG CRNG(Coarse5d); CRNG.SeedFixedIntegers(cseeds);
+
+  LatticeGaugeField Umu(UGrid);
+  SU<Nc>::HotConfiguration(RNG4,Umu);
+
+  RealD mass = 0.1;
+  RealD M5   = 1.8;
+  DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+
+  MdagMLinearOperator<DomainWallFermionD,LatticeFermion> HermDefOp(Ddwf);
+  ShiftedHermOpAdaptor<LatticeFermionD> HOA(HermDefOp, 0.5);
+
+  std::cout << GridLogMessage << "Building random aggregation space, nbasis " << nbasis << std::endl;
+  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
+  Subspace Aggregates(Coarse5d,FGrid,0);
+  Aggregates.CreateSubspaceRandom(RNG5);
+
+  std::cout << GridLogMessage << "Coarsening shifted MdagM" << std::endl;
+  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
+  typedef LittleDiracOperator::CoarseVector CoarseVector;
+  NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d);
+  LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse5d);
+  LittleDiracOp.CoarsenOperator(HOA,Aggregates);
+
+  ///////////////////////////////////////////////////////////////////////
+  // Operator health probes (independent of DenseCoarseMatrix).
+  //
+  // MEASURED PATHOLOGY, banked 2026-08-14: on coarse dims of 2 (fine
+  // 8.8.8.8, block 4 -> coarse 2.2.2.2) the coarsened operator is
+  // rank 16/128 with 112 zero ROWS (output support = 2 of 16 sites)
+  // and Hermiticity violation 0.17 -- the probing construction breaks
+  // on the size-2 torus.  The import certificate cannot see this
+  // (dense and M share _A).  Out of scope here; coarse dims >= 3.
+  //
+  // Cheap any-size probes: output support + Hermiticity via inner
+  // products on random vectors.
+  ///////////////////////////////////////////////////////////////////////
+  {
+    CoarseVector px(Coarse5d);
+    CoarseVector py(Coarse5d);
+    CoarseVector Mx(Coarse5d);
+    CoarseVector My(Coarse5d);
+    random(CRNG, px);
+    random(CRNG, py);
+    LittleDiracOp.M(px, Mx);
+    LittleDiracOp.M(py, My);
+    ComplexD ip1 = innerProduct(py, Mx);        // <y, M x>
+    ComplexD ip3 = innerProduct(px, My);        // <x, M y>
+    RealD hermdev = abs(ip1 - conj(ip3)) / std::sqrt(norm2(Mx)*norm2(py));
+    RealD support = norm2(Mx) / norm2(px);
+    std::cout << GridLogMessage << "Operator health: ||Mx||^2/||x||^2 = " << support
+              << "  herm-dev " << hermdev << std::endl;
+    // Hermitian fine op => exactly Hermitian Galerkin coarse op.
+    // (A measured herm-dev of 1.3e-4 here was the CPU SIMT-lane
+    // CoarsenOperator bug -- fixed 2026-08-14, now 4e-15.  A loud
+    // failure here means _A population is broken again.)
+    GRID_ASSERT( support > 1.0e-3 );
+    GRID_ASSERT( hermdev < 1.0e-10 );
+  }
+
+  ///////////////////////////////////////////////////////////////////////
+  // Full-matrix conditioning probe at small N: dense columns by
+  // applying M to unit vectors, fp64 Eigen SVD.
+  ///////////////////////////////////////////////////////////////////////
+  {
+    int64_t Nprobe = Coarse5d->gSites() * nbasis;
+    if ( Nprobe <= 700 )
+    {
+      Eigen::MatrixXcd eA(Nprobe, Nprobe);
+      CoarseVector e(Coarse5d);
+      CoarseVector Me(Coarse5d);
+      for(int64_t j=0; j<Nprobe; j++)
+      {
+        int64_t gsite = j / nbasis;
+        int     b     = j % nbasis;
+        e = Zero();
+        Coordinate gcoor(Coarse5d->_ndimension);
+        Lexicographic::CoorFromIndex(gcoor, gsite, Coarse5d->GlobalDimensions());
+        typedef typename CoarseVector::vector_object::scalar_object csobj;
+        csobj s;
+        s = Zero();
+        ((ComplexD *)&s)[b] = ComplexD(1.0,0.0);
+        pokeSite(s, e, gcoor);
+        LittleDiracOp.M(e, Me);
+        for(int64_t i=0; i<Nprobe; i++)
+        {
+          int64_t gsi = i / nbasis;
+          int     bi  = i % nbasis;
+          Coordinate gci(Coarse5d->_ndimension);
+          Lexicographic::CoorFromIndex(gci, gsi, Coarse5d->GlobalDimensions());
+          csobj si;
+          peekSite(si, Me, gci);
+          // Explicit re/im at the thrust/std boundary (HIP builds)
+          ComplexD zz = ((ComplexD *)&si)[bi];
+          eA(i,j) = std::complex<double>(zz.real(), zz.imag());
+        }
+      }
+      Eigen::JacobiSVD<Eigen::MatrixXcd> svd(eA);
+      double smax = svd.singularValues()(0);
+      double smin = svd.singularValues()(Nprobe-1);
+      int64_t rank = 0;
+      for(int64_t i=0; i<Nprobe; i++)
+      {
+        if ( svd.singularValues()(i) > 1.0e-10*smax ) rank++;
+      }
+      double herm = (eA - eA.adjoint()).cwiseAbs().maxCoeff();
+      int64_t zrows = 0;
+      int64_t zcols = 0;
+      for(int64_t i=0; i<Nprobe; i++)
+      {
+        if ( eA.row(i).cwiseAbs().maxCoeff() < 1.0e-12 ) zrows++;
+        if ( eA.col(i).cwiseAbs().maxCoeff() < 1.0e-12 ) zcols++;
+      }
+      std::cout << GridLogMessage << "Operator probe: N=" << Nprobe
+                << "  sigma_max " << smax
+                << "  sigma_min " << smin
+                << "  rank " << rank << "/" << Nprobe
+                << "  herm-dev " << herm
+                << "  zero rows/cols " << zrows << "/" << zcols
+                << std::endl;
+    }
+  }
+
+  ///////////////////////////////////////////////////////////////////////
+  // T6: AUDIT mode, fresh import, multi-panel gathers.  The constructor
+  // runs every certificate in the chain (see banner).
+  ///////////////////////////////////////////////////////////////////////
+  setenv("DENSE_SCHUR","2",1);
+  setenv("DENSE_PANEL_BYTES","65536",1);
+  unsetenv("SLAB_FILE");
+
+  typedef DenseCoarseMatrix<vSpinColourVector,vTComplex,nbasis> DenseCC;
+  DenseCC dcm(LittleDiracOp, Coarse5d);
+
+  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
+  ///////////////////////////////////////////////////////////////////////
+  CoarseVector x(Coarse5d);
+  CoarseVector y(Coarse5d);
+  CoarseVector z(Coarse5d);
+  random(CRNG, x);
+  dcm(x, y);
+  LittleDiracOp.M(y, z);
+  z = z - x;
+  RealD rel = std::sqrt(norm2(z)/norm2(x));
+  std::cout << GridLogMessage << "T6  round trip ||A Ainv x - x||/||x|| (random x) = "
+            << rel << std::endl;
+  GRID_ASSERT( rel < 1.0e-2 );
+
+  std::cout << GridLogMessage << "Test_schur_dense_coarse: T6 ALL PASS" << std::endl;
+
+  Grid_finalize();
+}

tests/debug/Test_schur_inverse.cc

@@ -0,0 +1,769 @@
+/*************************************************************************************
+
+    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();
+}