Claude build fixeds

Grid/algorithms/multigrid/BlockCyclic.h

@@ -0,0 +1,215 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/BlockCyclic.h
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+///////////////////////////////////////////////////////////////////////////////
+// BlockCyclicLayout: the index arithmetic of a 2D block-cyclic distribution
+// of an N x N matrix over a Pr x Pc logical process grid with block size nb.
+//
+// This is stage 1 of the 2D distributed dense inverse
+// (documentation/DistributedDenseInverse2D.tex).  It is deliberately
+// COMMUNICATOR-FREE: every mapping is a static pure function of
+// (N, nb, Pr, Pc), so the whole layout is exhaustively unit-testable on one
+// rank with no MPI in the loop (Test_blockcyclic).  A thin instance layer
+// binds a world rank to a grid coordinate and caches local extents.
+//
+// Conventions (fixed here, relied on by every later stage):
+//
+//  * Global block b of a dimension with Pg processes is owned by process
+//    coordinate  b % Pg  (ScaLAPACK csrc=0), and is that process's local
+//    block  b / Pg.
+//  * Rank <-> grid coordinate is ROW MAJOR over the process grid:
+//        rank = p*Pc + q ,  p = rank/Pc ,  q = rank%Pc .
+//    The eventual ring transport must construct its neighbour tables with
+//    the same convention.
+//  * Local storage is COLUMN MAJOR with ld = mloc, matching BlockRows:
+//    local element (i,j) lives at  data[i + j*mloc].
+//  * The trailing partial block (N % nb != 0) belongs to the owner of the
+//    last full-size block position; only that one block is short.
+//
+// Element (gi,gj) therefore lives on grid coordinate
+//    ( (gi/nb) % Pr , (gj/nb) % Pc )
+// at local coordinate
+//    ( ((gi/nb)/Pr)*nb + gi%nb , ((gj/nb)/Pc)*nb + gj%nb ).
+//
+// Everything here is host-side integer arithmetic; nothing allocates.
+///////////////////////////////////////////////////////////////////////////////
+
+class BlockCyclicLayout
+{
+public:
+  ///////////////////////////////////////////////////////////////////////////
+  // Closest-to-square factorisation Pr*Pc == P with Pr <= Pc.
+  // For fixed P the per-rank SUMMA volume  N^2 (1/Pr + 1/Pc)  is minimised
+  // at the most square grid.  P=288 -> 16 x 18.
+  ///////////////////////////////////////////////////////////////////////////
+  static void ChooseProcessGrid(int P, int &Pr, int &Pc)
+  {
+    GRID_ASSERT(P >= 1);
+    Pr = 1;
+    for(int r=1; (int64_t)r*r <= (int64_t)P; r++)
+      if ( P % r == 0 ) Pr = r;
+    Pc = P / Pr;
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Number of rows (or columns) of a dimension of global extent N, block nb,
+  // owned by process coordinate p of Pg.  ScaLAPACK "numroc", csrc=0.
+  ///////////////////////////////////////////////////////////////////////////
+  static int64_t NumLocal(int64_t N, int64_t nb, int p, int Pg)
+  {
+    GRID_ASSERT(N  >= 0);
+    GRID_ASSERT(nb >= 1);
+    GRID_ASSERT(p  >= 0);
+    GRID_ASSERT(p  <  Pg);
+    int64_t nblocks = N / nb;            // full blocks
+    int64_t extra   = N % nb;            // trailing partial block
+    int64_t full    = nblocks / Pg;      // full blocks everyone owns
+    int64_t rem     = nblocks % Pg;      // coords [0,rem) own one more
+    int64_t n = full*nb;
+    if ( p <  (int)rem ) n += nb;        // an extra full block
+    if ( p == (int)rem ) n += extra;     // the partial block, if any
+    return n;
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Global index -> (owner coordinate, local index) in one dimension.
+  ///////////////////////////////////////////////////////////////////////////
+  static void GlobalToLocal(int64_t g, int64_t nb, int Pg,
+                            int &owner, int64_t &loc)
+  {
+    GRID_ASSERT(g >= 0);
+    int64_t b = g / nb;                  // global block
+    owner = (int)(b % Pg);
+    loc   = (b / Pg)*nb + (g % nb);
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // (process coordinate, local index) -> global index in one dimension.
+  // Inverse of GlobalToLocal on the owned set.
+  ///////////////////////////////////////////////////////////////////////////
+  static int64_t LocalToGlobal(int64_t l, int64_t nb, int p, int Pg)
+  {
+    GRID_ASSERT(l >= 0);
+    int64_t lb = l / nb;                 // local block
+    int64_t b  = lb*Pg + p;              // global block
+    return b*nb + (l % nb);
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Instance layer: bind a rank of a Pr x Pc grid.
+  ///////////////////////////////////////////////////////////////////////////
+  int64_t N;         // global matrix dimension (square)
+  int64_t nb;        // block size
+  int     Pr, Pc;    // process grid
+  int     me;        // world rank within the grid, row major
+  int     prow, pcol;// my grid coordinate
+  int64_t mloc,nloc; // my local extents; storage column major, ld = mloc
+
+  BlockCyclicLayout(int64_t N_, int64_t nb_, int Pr_, int Pc_, int me_)
+  {
+    N  = N_;
+    nb = nb_;
+    Pr = Pr_;
+    Pc = Pc_;
+    me = me_;
+    GRID_ASSERT( N  >= 0 );
+    GRID_ASSERT( nb >= 1 );
+    GRID_ASSERT( Pr >= 1 );
+    GRID_ASSERT( Pc >= 1 );
+    GRID_ASSERT( me >= 0 );
+    GRID_ASSERT( me <  Pr*Pc );
+    prow = me / Pc;                      // ROW MAJOR rank convention
+    pcol = me % Pc;
+    mloc = NumLocal(N, nb, prow, Pr);
+    nloc = NumLocal(N, nb, pcol, Pc);
+  }
+
+  // Owning rank of global element (gi,gj), row-major rank convention.
+  int OwnerRank(int64_t gi, int64_t gj) const
+  {
+    int pr,pc; int64_t li,lj;
+    GlobalToLocal(gi, nb, Pr, pr, li);
+    GlobalToLocal(gj, nb, Pc, pc, lj);
+    return pr*Pc + pc;
+  }
+
+  // My local storage offset of global element (gi,gj).
+  // The caller must know I own it; asserted, not assumed.
+  int64_t LocalOffset(int64_t gi, int64_t gj) const
+  {
+    int pr,pc; int64_t li,lj;
+    GlobalToLocal(gi, nb, Pr, pr, li);
+    GlobalToLocal(gj, nb, Pc, pc, lj);
+    GRID_ASSERT( pr == prow );
+    GRID_ASSERT( pc == pcol );
+    return li + lj*mloc;                 // column major, ld = mloc
+  }
+
+  // Do I own global element (gi,gj)?
+  int Owns(int64_t gi, int64_t gj) const
+  {
+    return OwnerRank(gi,gj) == me;
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Block-aligned global range [g0,g1) -> my contiguous local range [l0,l1).
+  //
+  // For fixed owner p the local index is monotone in the global index, so a
+  // coordinate's owned elements of ANY global range are contiguous in local
+  // storage; and for a BLOCK-ALIGNED range the bounds are exactly
+  // NumLocal(g0) and NumLocal(g1), because NumLocal(g,...) counts the owned
+  // elements below g.  This is what lets a windowed product view the local
+  // sub-matrix of a global window as &data[l0 + c0*mloc] with the SAME ld --
+  // no gather, no copy.  Verified exhaustively in Test_blockcyclic T7.
+  //
+  // g0 must be a block multiple; g1 a block multiple or N itself.
+  ///////////////////////////////////////////////////////////////////////////
+  static void RangeToLocal(int64_t g0, int64_t g1,
+                           int64_t N, int64_t nb, int p, int Pg,
+                           int64_t &l0, int64_t &l1)
+  {
+    GRID_ASSERT( 0 <= g0 );
+    GRID_ASSERT( g0 <= g1 );
+    GRID_ASSERT( g1 <= N );
+    GRID_ASSERT( g0 % nb == 0 );
+    GRID_ASSERT( (g1 % nb == 0) || (g1 == N) );
+    l0 = NumLocal(g0, nb, p, Pg);
+    l1 = NumLocal(g1, nb, p, Pg);
+  }
+
+  // Instance forms, rows and columns of my own coordinate.
+  void RowRange(int64_t g0, int64_t g1, int64_t &l0, int64_t &l1) const
+  { RangeToLocal(g0,g1,N,nb,prow,Pr,l0,l1); }
+  void ColRange(int64_t g0, int64_t g1, int64_t &l0, int64_t &l1) const
+  { RangeToLocal(g0,g1,N,nb,pcol,Pc,l0,l1); }
+
+  // Size of global block b (the trailing block may be short).
+  int64_t BlockSize(int64_t b) const
+  {
+    int64_t lo = b*nb;
+    GRID_ASSERT( lo < N );
+    return std::min(N, lo+nb) - lo;
+  }
+};
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/BlockCyclicRedistribute.h

@@ -0,0 +1,247 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/BlockCyclicRedistribute.h
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+#include <Grid/algorithms/multigrid/BlockCyclicSumma.h>
+
+NAMESPACE_BEGIN(Grid);
+
+///////////////////////////////////////////////////////////////////////////////
+// Stage 4 of the 2D distributed dense inverse: redistribution between the
+// 1D rank-major row layout (BlockRows: rank r owns contiguous global rows
+// [rowStart[r], rowStart[r+1]) of an N x N matrix, stored rows x N column
+// major with ld = rows) and the 2D block-cyclic layout.
+//
+// This is what lets the EXISTING stencil->dense import, its certificate,
+// the fp32 slab conversion and the apply path all remain byte-for-byte
+// untouched: the 2D inverse slots between them as
+//
+//    RowsToCyclic -> BlockCyclicSchurInverse::Invert -> CyclicToRows
+//
+// Volume is one matrix pass each way -- N^2/P elements per rank (~1 GB at
+// production), trivial against the inversion itself.
+//
+// Transport: PURE POINT-TO-POINT, like everything else in this stack.
+// Ranks exchange in a round-robin TOURNAMENT (the circle method, on an odd
+// modulus M so it covers every pair exactly once for any P, with byes):
+// at round r, ranks x and y are partners iff x+y == r (mod M).  Each
+// meeting handles both directed edges of the pair in ONE SendToRecvFrom,
+// padded to the larger of the two edge sizes -- SendToRecvFrom carries a
+// single byte count for both directions, and both endpoints compute the
+// same max from the shared descriptors, so there is no asymmetric-size
+// case and no zero-count shape.  Pairs with nothing to exchange skip the
+// round, decided identically at both ends.
+//
+// Element enumeration within an edge is canonical -- ascending global
+// column outer, ascending global row inner -- and each endpoint builds its
+// OWN local offset tables from the shared descriptors, so no index data is
+// ever transmitted.  The round trip is BITWISE exact (pure data movement,
+// no arithmetic): Test_schur2d_redist proves it.
+///////////////////////////////////////////////////////////////////////////////
+
+class BlockCyclicRedistribute
+{
+public:
+  /////////////////////////////////////////////////////////////////////////
+  // The directed edge (1D rank r1, 2D rank r2): global rows of r1's range
+  // whose row-block coordinate is r2's prow; ALL global columns whose
+  // column-block coordinate is r2's pcol.  Every rank can enumerate any
+  // edge from (rowStart, layout) alone.
+  /////////////////////////////////////////////////////////////////////////
+  static void EdgeRows(const std::vector<int64_t> &rowStart, int r1,
+                       const BlockCyclicLayout &L, int r2,
+                       std::vector<int64_t> &rows)
+  {
+    rows.clear();
+    int p = r2 / L.Pc;                       // row-major rank convention
+    for(int64_t i=rowStart[r1]; i<rowStart[r1+1]; i++)
+      if ( (int)((i/L.nb) % L.Pr) == p ) rows.push_back(i);
+  }
+  static void EdgeCols(const BlockCyclicLayout &L, int r2,
+                       std::vector<int64_t> &cols)
+  {
+    cols.clear();
+    int q = r2 % L.Pc;
+    for(int64_t b=0; b*L.nb<L.N; b++){
+      if ( (int)(b % L.Pc) != q ) continue;
+      int64_t g0=b*L.nb, g1=std::min(L.N,(b+1)*L.nb);
+      for(int64_t j=g0;j<g1;j++) cols.push_back(j);
+    }
+  }
+
+  /////////////////////////////////////////////////////////////////////////
+  // Gather/scatter one edge between a matrix (device, column major, ld)
+  // and a dense edge buffer, through device offset tables.
+  //   buffer(a,b) = elem(rows[a], cols[b]),  a fastest.
+  /////////////////////////////////////////////////////////////////////////
+  static void MoveEdge(int toBuffer,
+                       ComplexD *mat, int64_t ld,
+                       const std::vector<int64_t> &roff,  // per-row offset in mat
+                       const std::vector<int64_t> &coff,  // per-col offset in mat
+                       ComplexD *buf)
+  {
+    int64_t nr = roff.size();
+    int64_t nc = coff.size();
+    if ( !(nr && nc) ) return;
+    deviceVector<int64_t> dro(nr), dco(nc);
+    acceleratorCopyToDevice((void *)&roff[0], (void *)&dro[0], nr*sizeof(int64_t));
+    acceleratorCopyToDevice((void *)&coff[0], (void *)&dco[0], nc*sizeof(int64_t));
+    int64_t *ro = &dro[0];
+    int64_t *co = &dco[0];
+    if ( toBuffer ) {
+      accelerator_for(idx, (uint64_t)(nr*nc), 1, {
+        int64_t b = idx / nr;
+        int64_t a = idx - b*nr;
+        buf[a + b*nr] = mat[ ro[a] + co[b]*ld ];
+      });
+    } else {
+      accelerator_for(idx, (uint64_t)(nr*nc), 1, {
+        int64_t b = idx / nr;
+        int64_t a = idx - b*nr;
+        mat[ ro[a] + co[b]*ld ] = buf[a + b*nr];
+      });
+    }
+  }
+
+  /////////////////////////////////////////////////////////////////////////
+  // My offset tables for an edge, on whichever side I am.
+  /////////////////////////////////////////////////////////////////////////
+  static void Offsets1D(const std::vector<int64_t> &rows,
+                        const std::vector<int64_t> &cols,
+                        int64_t row0,
+                        std::vector<int64_t> &roff, std::vector<int64_t> &coff)
+  {
+    roff.resize(rows.size());  coff.resize(cols.size());
+    for(uint64_t a=0;a<rows.size();a++) roff[a] = rows[a]-row0;   // local row
+    for(uint64_t b=0;b<cols.size();b++) coff[b] = cols[b];        // global col
+  }
+  static void Offsets2D(const BlockCyclicLayout &L,
+                        const std::vector<int64_t> &rows,
+                        const std::vector<int64_t> &cols,
+                        std::vector<int64_t> &roff, std::vector<int64_t> &coff)
+  {
+    roff.resize(rows.size());  coff.resize(cols.size());
+    for(uint64_t a=0;a<rows.size();a++){
+      int p; int64_t l;
+      BlockCyclicLayout::GlobalToLocal(rows[a], L.nb, L.Pr, p, l);
+      GRID_ASSERT( p == L.prow );
+      roff[a] = l;
+    }
+    for(uint64_t b=0;b<cols.size();b++){
+      int q; int64_t l;
+      BlockCyclicLayout::GlobalToLocal(cols[b], L.nb, L.Pc, q, l);
+      GRID_ASSERT( q == L.pcol );
+      coff[b] = l;
+    }
+  }
+
+  /////////////////////////////////////////////////////////////////////////
+  // The worker.  dir=+1 : 1D rows -> block cyclic ;  dir=-1 : back.
+  /////////////////////////////////////////////////////////////////////////
+  static void Redistribute(int dir, GridBase *grid,
+                           const std::vector<int64_t> &rowStart,
+                           ComplexD *rows1d, int64_t myrows,
+                           BlockCyclicMatrix &A)
+  {
+    BlockCyclicLayout &L = A.layout;
+    int P  = grid->ProcessorCount();
+    int me = grid->ThisRank();
+    GRID_ASSERT( (int)rowStart.size() == P+1 );
+    GRID_ASSERT( rowStart[P] == L.N );
+    GRID_ASSERT( rowStart[me+1]-rowStart[me] == myrows );
+    int64_t row0 = rowStart[me];
+    int64_t ld1  = myrows ? myrows : 1;
+
+    std::vector<int64_t> rows, cols, roff, coff;
+    deviceVector<ComplexD> sbuf(1), rbuf(1);
+
+    ///////////////////////////////////////////////////////////////////////
+    // Self edge: purely local, via a bounce buffer (shares all the code).
+    ///////////////////////////////////////////////////////////////////////
+    EdgeRows(rowStart, me, L, me, rows);
+    EdgeCols(L, me, cols);
+    if ( rows.size() && cols.size() ){
+      uint64_t ne = rows.size()*cols.size();
+      if ( sbuf.size() < ne ) sbuf.resize(ne);
+      std::vector<int64_t> roff2, coff2;
+      Offsets1D(rows, cols, row0, roff, coff);
+      Offsets2D(L,  rows, cols, roff2, coff2);
+      if ( dir > 0 ) {
+        MoveEdge(1, rows1d,     ld1,     roff,  coff,  &sbuf[0]);
+        MoveEdge(0, &A.data[0], L.mloc,  roff2, coff2, &sbuf[0]);
+      } else {
+        MoveEdge(1, &A.data[0], L.mloc,  roff2, coff2, &sbuf[0]);
+        MoveEdge(0, rows1d,     ld1,     roff,  coff,  &sbuf[0]);
+      }
+    }
+
+    ///////////////////////////////////////////////////////////////////////
+    // Tournament over all pairs: odd modulus M, partner = (r - me) mod M.
+    // Every unordered pair meets exactly once; partner==me or >=P is a bye.
+    ///////////////////////////////////////////////////////////////////////
+    int M = (P%2) ? P : P+1;
+    for(int r=0;r<M;r++){
+      int partner = (int)(((int64_t)r - me + 2L*M) % M);
+      if ( partner == me || partner >= P ) continue;
+
+      // outbound edge: my (dir>0 ? 1D rows : 2D data) -> partner
+      // inbound  edge: partner -> my (dir>0 ? 2D data : 1D rows)
+      std::vector<int64_t> orow, ocol, irow, icol;
+      if ( dir > 0 ) { EdgeRows(rowStart, me,      L, partner, orow); EdgeCols(L, partner, ocol);
+                       EdgeRows(rowStart, partner, L, me,      irow); EdgeCols(L, me,      icol); }
+      else           { EdgeRows(rowStart, partner, L, me,      orow); EdgeCols(L, me,      ocol);
+                       EdgeRows(rowStart, me,      L, partner, irow); EdgeCols(L, partner, icol); }
+
+      uint64_t nout = orow.size()*ocol.size();
+      uint64_t nin  = irow.size()*icol.size();
+      if ( !(nout || nin) ) continue;              // both ends compute this identically
+
+      uint64_t nmax = std::max(nout,nin);          // symmetric padded transfer
+      if ( sbuf.size() < nmax ) sbuf.resize(nmax);
+      if ( rbuf.size() < nmax ) rbuf.resize(nmax);
+
+      if ( nout ){
+        if ( dir > 0 ) { Offsets1D(orow, ocol, row0, roff, coff);
+                         MoveEdge(1, rows1d, ld1, roff, coff, &sbuf[0]); }
+        else           { Offsets2D(L, orow, ocol, roff, coff);
+                         MoveEdge(1, &A.data[0], L.mloc, roff, coff, &sbuf[0]); }
+      }
+      grid->SendToRecvFrom((void *)&sbuf[0], partner,
+                           (void *)&rbuf[0], partner,
+                           nmax*sizeof(ComplexD));
+      if ( nin ){
+        if ( dir > 0 ) { Offsets2D(L, irow, icol, roff, coff);
+                         MoveEdge(0, &A.data[0], L.mloc, roff, coff, &rbuf[0]); }
+        else           { Offsets1D(irow, icol, row0, roff, coff);
+                         MoveEdge(0, rows1d, ld1, roff, coff, &rbuf[0]); }
+      }
+    }
+  }
+
+  static void RowsToCyclic(GridBase *grid, const std::vector<int64_t> &rowStart,
+                           ComplexD *rows1d, int64_t myrows, BlockCyclicMatrix &A)
+  { Redistribute(+1, grid, rowStart, rows1d, myrows, A); }
+
+  static void CyclicToRows(GridBase *grid, const std::vector<int64_t> &rowStart,
+                           BlockCyclicMatrix &A, ComplexD *rows1d, int64_t myrows)
+  { Redistribute(-1, grid, rowStart, rows1d, myrows, A); }
+};
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/BlockCyclicSchurInverse.h

@@ -0,0 +1,291 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/BlockCyclicSchurInverse.h
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+#include <Grid/algorithms/blas/BatchedInverse.h>
+#include <Grid/algorithms/multigrid/BlockCyclicSumma.h>
+
+NAMESPACE_BEGIN(Grid);
+
+///////////////////////////////////////////////////////////////////////////////
+// Stage 3 of the 2D distributed dense inverse: the recursive Schur
+// complement on a block-cyclic matrix, in place.
+//
+// The nine-step algebra is IDENTICAL to RecursiveSchurInverse (1D); what
+// changes is the decomposition.  The recursion splits the GLOBAL INDEX
+// RANGE at the block boundary nearest the midpoint -- not the rank range --
+// so every rank owns part of every sub-block at every depth, and the
+// ownership gating (inI/inJ, dummy operands, zero-width rank ranges) of the
+// 1D scheme has no analogue here: it is simply gone.
+//
+//   I = [c0,m)   J = [m,c1)      (block-aligned, m the mid block boundary)
+//   1. recurse I  :  A11 -> A11inv                        (in place)
+//   2. Bt = A11inv . A12                                  (scratch, I x J)
+//   3. Ct = A21 . A11inv                                  (scratch, J x I)
+//   4. A22 -= A21 . Bt          == S                      (in place)
+//   5. recurse J  :  S -> Sinv                            (in place)
+//   6. Tt = Sinv . Ct                                     (scratch, J x I)
+//   7. Ut = Bt . Sinv                                     (scratch, I x J)
+//   8. A11 += Ut . Ct           == X11                    (in place)
+//   9. A12 = -Ut ,  A21 = -Tt                             (window copies)
+//
+// SCRATCH SHARING.  Four full-size block-cyclic scratch matrices (Bt, Ct,
+// Tt, Ut) serve the ENTIRE tree, used through windows.  This is safe at
+// every depth because of a window-disjointness invariant:
+//
+//   * every temporary of a node has its row range in one half of the
+//     node's window and its column range in the other (I x J or J x I);
+//   * everything any DESCENDANT touches -- its A windows and its own
+//     temporaries -- has BOTH ranges inside a single half (I x I during
+//     step 1, J x J during step 5).
+//
+// Hence a descendant window and a live ancestor temporary always differ in
+// at least one dimension by disjoint ranges.  Only Bt and Ct are live
+// across the step-5 recursion (Tt, Ut are written after it), and both are
+// covered by the invariant.
+//
+// LEAF.  A leaf is a single diagonal block, and block (b,b) of a
+// block-cyclic layout lives ENTIRELY on rank (b%Pr, b%Pc).  The leaf
+// inversion is therefore purely local -- pack the strided block dense,
+// GridBLASInverse, unpack -- with NO communication and no assembly.  The
+// leaf-assembly transport question of the 1D scheme does not arise.
+// Successive leaves cycle over ranks, so leaf work is naturally spread.
+//
+// COMMUNICATION.  Every transfer in the whole inversion is a
+// SendToRecvFrom inside BlockCyclicSumma's rings: pure point-to-point, no
+// collectives on the critical path, deterministic summation order (so
+// repeated inversions are bitwise identical).  ReportTelemetry() is the
+// one optional exception: it performs reductions, and is only ever called
+// explicitly by a caller who wants the numbers.
+//
+// NUMERICS.  No pivoting, exactly as the 1D scheme: every A11 and every
+// Schur complement met on the way down must be non-singular.  The growth
+// telemetry stands in for pivoting; note the recursion splits differently
+// from the 1D rank-range tree, so DIFFERENT sub-blocks are inverted and
+// telemetry values are NOT comparable with the 1D implementation's --
+// re-baseline, do not compare.
+///////////////////////////////////////////////////////////////////////////////
+
+class BlockCyclicSchurInverse
+{
+public:
+  BlockCyclicSumma        SUMMA;
+  GridBLASInverse         INV;
+
+  // Telemetry: accumulated LOCALLY, no comms unless ReportTelemetry().
+  double                  telLeafMaxInv;
+  uint64_t                nLeaf;
+  uint64_t                nNode;
+  double                  tLeaf;
+  double                  tGemm;      // wall in Multiply calls (comms+gemm)
+  double                  tCopy;
+
+  BlockCyclicSchurInverse()
+  {
+    telLeafMaxInv = 0.0;
+    nLeaf = nNode = 0;
+    tLeaf = tGemm = tCopy = 0.0;
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Window copy-scale:  Dst[i0:i1, j0:j1] = alpha * Src[same window].
+  // Both share one layout, so the local bands coincide; pure local kernel.
+  ///////////////////////////////////////////////////////////////////////////
+  void WindowCopyScale(ComplexD alpha,
+                       BlockCyclicMatrix &Src, BlockCyclicMatrix &Dst,
+                       int64_t i0, int64_t i1, int64_t j0, int64_t j1)
+  {
+    BlockCyclicLayout &L = Dst.layout;
+    GRID_ASSERT( Src.layout.N==L.N && Src.layout.nb==L.nb );
+    GRID_ASSERT( Src.layout.Pr==L.Pr && Src.layout.Pc==L.Pc );
+    int64_t li0,li1, lj0,lj1;
+    L.RowRange(i0,i1, li0,li1);
+    L.ColRange(j0,j1, lj0,lj1);
+    int64_t m = li1-li0, n = lj1-lj0;
+    if ( !(m && n) ) return;
+    ComplexD *src = Src.LocalWindow(li0,lj0);
+    ComplexD *dst = Dst.LocalWindow(li0,lj0);
+    int64_t ldS = Src.layout.mloc;
+    int64_t ldD = L.mloc;
+    tCopy -= usecond();
+    accelerator_for(idx, (uint64_t)(m*n), 1, {
+      int64_t jj = idx / m;
+      int64_t ii = idx - jj*m;
+      dst[ii + jj*ldD] = alpha*src[ii + jj*ldS];
+    });
+    tCopy += usecond();
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Leaf: single diagonal block (b,b), entirely on rank (b%Pr, b%Pc).
+  // Local pack -> dense inverse -> unpack; every other rank does nothing
+  // and needs no synchronisation: the next SUMMA's rings pair them up.
+  ///////////////////////////////////////////////////////////////////////////
+  void Leaf(BlockCyclicMatrix &A, int64_t b)
+  {
+    BlockCyclicLayout &L = A.layout;
+    nLeaf++;
+    if ( (int)(b % L.Pr) != L.prow ) return;
+    if ( (int)(b % L.Pc) != L.pcol ) return;
+
+    tLeaf -= usecond();
+    int64_t g0 = b*L.nb;
+    int64_t g1 = std::min(L.N, g0+L.nb);
+    int64_t w  = g1-g0;
+    int64_t lr0,lr1, lc0,lc1;
+    L.RowRange(g0,g1, lr0,lr1);
+    L.ColRange(g0,g1, lc0,lc1);
+    GRID_ASSERT( lr1-lr0 == w );
+    GRID_ASSERT( lc1-lc0 == w );
+
+    // Pack the strided block dense (inverseBatched assumes lda == w).
+    deviceVector<ComplexD> dense((uint64_t)w*w);
+    {
+      ComplexD *src = A.LocalWindow(lr0,lc0);
+      ComplexD *dst = &dense[0];
+      int64_t   ld  = L.mloc;
+      accelerator_for(idx, (uint64_t)(w*w), 1, {
+        int64_t jj = idx / w;
+        int64_t ii = idx - jj*w;
+        dst[ii + jj*w] = src[ii + jj*ld];
+      });
+    }
+    {
+      deviceVector<ComplexD*> bp(1);
+      std::vector<ComplexD*>  ptr(1);
+      ptr[0] = &dense[0];
+      acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexD*));
+      INV.inverseBatched(w, bp);
+    }
+    {
+      ComplexD *src = &dense[0];
+      ComplexD *dst = A.LocalWindow(lr0,lc0);
+      int64_t   ld  = L.mloc;
+      accelerator_for(idx, (uint64_t)(w*w), 1, {
+        int64_t jj = idx / w;
+        int64_t ii = idx - jj*w;
+        dst[ii + jj*ld] = src[ii + jj*w];
+      });
+    }
+    // Growth telemetry, local only.
+    {
+      std::vector<ComplexD> h((uint64_t)w*w);
+      acceleratorCopyFromDevice(&dense[0], &h[0], h.size()*sizeof(ComplexD));
+      double mx = 0.0;
+      for(auto &z : h){
+        double re=z.real(), im=z.imag();
+        mx = std::max(mx, re*re+im*im);
+      }
+      telLeafMaxInv = std::max(telLeafMaxInv, std::sqrt(mx));
+    }
+    tLeaf += usecond();
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // The recursion, on global BLOCK range [b0,b1).  SPMD: every rank calls
+  // with identical arguments; there is no ownership gating to get wrong.
+  ///////////////////////////////////////////////////////////////////////////
+  void SchurNode(BlockCyclicMatrix &A,
+                 BlockCyclicMatrix &Bt, BlockCyclicMatrix &Ct,
+                 BlockCyclicMatrix &Tt, BlockCyclicMatrix &Ut,
+                 int64_t b0, int64_t b1)
+  {
+    BlockCyclicLayout &L = A.layout;
+    int64_t span = b1-b0;
+    GRID_ASSERT( span >= 1 );
+    if ( span == 1 ) { Leaf(A, b0); return; }
+    nNode++;
+
+    int64_t bm = b0 + span/2;
+    int64_t c0 = b0*L.nb;
+    int64_t m  = bm*L.nb;
+    int64_t c1 = std::min(L.N, b1*L.nb);
+    ComplexD one (1.0,0.0), mone(-1.0,0.0), zero(0.0,0.0);
+
+    // 1. A11 -> A11inv
+    SchurNode(A,Bt,Ct,Tt,Ut, b0,bm);
+
+    tGemm -= usecond();
+    // 2. Bt[I,J] = A11inv . A12
+    SUMMA.Multiply(one,  A,  A,  zero, Bt, c0,m,  m,c1,  c0,m );
+    // 3. Ct[J,I] = A21 . A11inv
+    SUMMA.Multiply(one,  A,  A,  zero, Ct, m,c1,  c0,m,  c0,m );
+    // 4. A22 -= A21 . Bt   (the Schur complement, in place)
+    SUMMA.Multiply(mone, A,  Bt, one,  A,  m,c1,  m,c1,  c0,m );
+    tGemm += usecond();
+
+    // 5. S -> Sinv          (Bt, Ct live across this call: see invariant)
+    SchurNode(A,Bt,Ct,Tt,Ut, bm,b1);
+
+    tGemm -= usecond();
+    // 6. Tt[J,I] = Sinv . Ct
+    SUMMA.Multiply(one,  A,  Ct, zero, Tt, m,c1,  c0,m,  m,c1 );
+    // 7. Ut[I,J] = Bt . Sinv
+    SUMMA.Multiply(one,  Bt, A,  zero, Ut, c0,m,  m,c1,  m,c1 );
+    // 8. A11 += Ut . Ct
+    SUMMA.Multiply(one,  Ut, Ct, one,  A,  c0,m,  c0,m,  m,c1 );
+    tGemm += usecond();
+
+    // 9. Off-diagonal signs
+    WindowCopyScale(mone, Ut, A, c0,m,  m,c1);
+    WindowCopyScale(mone, Tt, A, m,c1,  c0,m);
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // PUBLIC ENTRY.  In-place inverse of the whole matrix.  Scratch (4x the
+  // matrix footprint) is allocated here and released on return.
+  ///////////////////////////////////////////////////////////////////////////
+  void Invert(BlockCyclicMatrix &A)
+  {
+    BlockCyclicLayout &L = A.layout;
+    GRID_ASSERT( L.N >= 1 );
+    int64_t nblocks = (L.N + L.nb - 1)/L.nb;
+
+    BlockCyclicMatrix Bt(A.grid, L.N, L.nb, L.Pr, L.Pc);
+    BlockCyclicMatrix Ct(A.grid, L.N, L.nb, L.Pr, L.Pc);
+    BlockCyclicMatrix Tt(A.grid, L.N, L.nb, L.Pr, L.Pc);
+    BlockCyclicMatrix Ut(A.grid, L.N, L.nb, L.Pr, L.Pc);
+
+    telLeafMaxInv = 0.0;
+    nLeaf = nNode = 0;
+    tLeaf = tGemm = tCopy = 0.0;
+
+    SchurNode(A, Bt,Ct,Tt,Ut, 0, nblocks);
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Optional, and the ONLY place any reduction happens: call it if you
+  // want the numbers, never from Invert.
+  ///////////////////////////////////////////////////////////////////////////
+  void ReportTelemetry(GridBase *grid)
+  {
+    RealD mx = telLeafMaxInv;
+    grid->GlobalMax(mx);
+    std::cout << GridLogMessage << "BlockCyclicSchurInverse:"
+              << "  nodes " << nNode << "  leaves " << nLeaf
+              << "  max|leafinv| " << mx
+              << "  (boss secs: gemm+comms " << tGemm/1.0e6
+              << "  leaf " << tLeaf/1.0e6
+              << "  copy " << tCopy/1.0e6 << ")"
+              << std::endl;
+  }
+};
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/BlockCyclicSumma.h

@@ -0,0 +1,299 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./Grid/algorithms/multigrid/BlockCyclicSumma.h
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+#include <Grid/algorithms/multigrid/BlockCyclic.h>
+
+NAMESPACE_BEGIN(Grid);
+
+///////////////////////////////////////////////////////////////////////////////
+// Stage 2 of the 2D distributed dense inverse: the windowed SUMMA product
+//
+//   C[i0:i1, j0:j1]  <-  beta C[i0:i1, j0:j1]
+//                      + alpha A[i0:i1, k0:k1] . B[k0:k1, j0:j1]
+//
+// on block-cyclic matrices sharing one BlockCyclicLayout.  All ranges are
+// BLOCK ALIGNED (multiples of nb, or N itself at the top end): stage 3's
+// recursion splits on block boundaries, so nothing else is ever needed, and
+// alignment makes every local window a contiguous band of local storage
+// (BlockCyclic.h RangeToLocal, Test_blockcyclic T7).
+//
+// Transport is PURE POINT-TO-POINT: SendToRecvFrom on explicit world ranks
+// computed from the row-major rank convention.  No collectives of any kind
+// -- no sub-communicators, no broadcast, no allgather -- by design: the
+// collective pathologies measured on this machine (MPI_Allgatherv at
+// ~0.18 MB/s with skewed counts, mpir_request.h:508 aborts, allreduce at
+// 39% of the P2P rate) motivated this implementation.  SendToRecvFrom is
+// the most exercised device-buffer path in Grid and the only one never
+// implicated.
+//
+// Algorithm: round-based ring allgather SUMMA.  The k range is processed in
+// rounds of Pc consecutive global blocks.  Within a round
+//
+//   * process column c owns at most one A panel  (blocks s with s%Pc == c);
+//     the Pc panels circulate around each process-ROW ring in Pc-1 steps;
+//   * process row r owns up to ceil(Pc/Pr) B panels (blocks s%Pr == r);
+//     they circulate around each process-COLUMN ring in Pr-1 steps;
+//   * every rank then accumulates  Cloc += alpha * Apanel_s . Bpanel_s
+//     for each block s of the round, in ascending s: a fixed summation
+//     order, so REPEATED RUNS ARE BITWISE IDENTICAL (no reduction, no
+//     order ambiguity -- the property the P2P doctrine buys).
+//
+// Ring chunks are PADDED to a fixed size (full nb panels, fixed
+// panels-per-origin): SendToRecvFrom carries one byte count for both
+// directions, so symmetric transfers eliminate every variable-size edge
+// case at a worst-case ~1/Pc extra volume on ragged rounds.  Padding is
+// never read: GEMMs address only the leading nb_s x width of each slot.
+//
+// Per-rank received volume: (k-extent) * (mloc_i + nloc_j) elements --
+// the N^2 (1/Pr + 1/Pc) SUMMA optimum, ~sqrt(P)/2 below the 1D scheme.
+///////////////////////////////////////////////////////////////////////////////
+
+class BlockCyclicMatrix
+{
+public:
+  GridBase              *grid;      // borrowed, never owned
+  BlockCyclicLayout      layout;
+  deviceVector<ComplexD> data;      // column major, ld = layout.mloc
+
+  BlockCyclicMatrix(GridBase *g, int64_t N, int64_t nb, int Pr, int Pc)
+    : grid(g),
+      layout(N, nb, Pr, Pc, g->ThisRank())
+  {
+    GRID_ASSERT( Pr*Pc == g->ProcessorCount() );
+    uint64_t sz = (uint64_t)layout.mloc*layout.nloc;
+    data.resize( sz ? sz : 1 );
+  }
+
+  ComplexD *LocalWindow(int64_t li, int64_t lj)
+  {
+    return &data[0] + li + lj*layout.mloc;
+  }
+
+  /////////////////////////////////////////////////////////////////////////
+  // TEST-SCALE import/export of a replicated global matrix (host, O(N^2)
+  // loops, one collective in Export).  For unit tests and the stage-3
+  // oracle only; production data enters through the direct block-cyclic
+  // import, never through these.
+  /////////////////////////////////////////////////////////////////////////
+  void ImportGlobal(const std::vector<ComplexD> &G)
+  {
+    int64_t N = layout.N;
+    GRID_ASSERT( (int64_t)G.size() == N*N );
+    std::vector<ComplexD> h((uint64_t)layout.mloc*layout.nloc, ComplexD(0.0,0.0));
+    for(int64_t j=0;j<N;j++){
+      for(int64_t i=0;i<N;i++){
+        if ( layout.Owns(i,j) ) h[layout.LocalOffset(i,j)] = G[i + j*N];
+      }
+    }
+    if ( h.size() )
+      acceleratorCopyToDevice(&h[0], &data[0], h.size()*sizeof(ComplexD));
+  }
+  void ExportGlobal(std::vector<ComplexD> &G)
+  {
+    int64_t N = layout.N;
+    G.assign((uint64_t)N*N, ComplexD(0.0,0.0));
+    std::vector<ComplexD> h((uint64_t)layout.mloc*layout.nloc);
+    if ( h.size() )
+      acceleratorCopyFromDevice(&data[0], &h[0], h.size()*sizeof(ComplexD));
+    for(int64_t j=0;j<N;j++){
+      for(int64_t i=0;i<N;i++){
+        if ( layout.Owns(i,j) ) G[i + j*N] = h[layout.LocalOffset(i,j)];
+      }
+    }
+    if ( N ) grid->GlobalSumVector((ComplexD *)&G[0], (int)(N*N)); // zero-fill: exact
+  }
+};
+
+class BlockCyclicSumma
+{
+public:
+  GridBLAS BLAS;
+
+  static int Overlap(int64_t a0,int64_t a1,int64_t b0,int64_t b1)
+  { return (a0 < b1) && (b0 < a1); }
+
+  void Multiply(ComplexD alpha,
+                BlockCyclicMatrix &A,
+                BlockCyclicMatrix &B,
+                ComplexD beta,
+                BlockCyclicMatrix &C,
+                int64_t i0, int64_t i1,
+                int64_t j0, int64_t j1,
+                int64_t k0, int64_t k1)
+  {
+    BlockCyclicLayout &L = C.layout;
+    GridBase *grid = C.grid;
+    const int64_t N  = L.N;
+    const int64_t nb = L.nb;
+    const int     Pr = L.Pr,   Pc   = L.Pc;
+    const int     prow = L.prow, pcol = L.pcol;
+
+    ///////////////////////////////////////////////////////////////////////
+    // Conformability: one layout, one communicator, aligned ranges.
+    ///////////////////////////////////////////////////////////////////////
+    auto same = [&](BlockCyclicLayout &X){
+      GRID_ASSERT( X.N==N ); GRID_ASSERT( X.nb==nb );
+      GRID_ASSERT( X.Pr==Pr ); GRID_ASSERT( X.Pc==Pc );
+      GRID_ASSERT( X.me==L.me );
+    };
+    same(A.layout); same(B.layout);
+    GRID_ASSERT( A.grid==grid ); GRID_ASSERT( B.grid==grid );
+    auto aligned = [&](int64_t g0,int64_t g1){
+      GRID_ASSERT( 0<=g0 ); GRID_ASSERT( g0<=g1 ); GRID_ASSERT( g1<=N );
+      GRID_ASSERT( g0%nb==0 ); GRID_ASSERT( (g1%nb==0)||(g1==N) );
+    };
+    aligned(i0,i1); aligned(j0,j1); aligned(k0,k1);
+    GRID_ASSERT( k1 > k0 );                    // pure scaling not supported here
+    // In-place windows are legal only if the written window is disjoint
+    // from anything read (stage 3 uses this; make violation loud).
+    if ( &C==&A ) GRID_ASSERT( !Overlap(j0,j1,k0,k1) );
+    if ( &C==&B ) GRID_ASSERT( !Overlap(i0,i1,k0,k1) );
+
+    ///////////////////////////////////////////////////////////////////////
+    // My local bands of the three windows (contiguous: T7).
+    ///////////////////////////////////////////////////////////////////////
+    int64_t li0,li1, lj0,lj1;
+    L.RowRange(i0,i1, li0,li1);
+    L.ColRange(j0,j1, lj0,lj1);
+    const int64_t mloc_i = li1-li0;            // my rows of the i window
+    const int64_t nloc_j = lj1-lj0;            // my cols of the j window
+
+    const int64_t kb0 = k0/nb;
+    const int64_t kb1 = (k1+nb-1)/nb;          // block-aligned or ==N: exact
+    const int64_t S   = (Pc + Pr - 1)/Pr;      // max B panels per origin row
+
+    ///////////////////////////////////////////////////////////////////////
+    // Round buffers, padded to fixed slot sizes (see header comment).
+    //   A: Pc slots of mloc_i x nb   (slot c = panel of the block owned by c)
+    //   B: Pr slots of S x (nb x nloc_j)
+    ///////////////////////////////////////////////////////////////////////
+    const uint64_t slotA  = (uint64_t)mloc_i*nb;
+    const uint64_t slotB1 = (uint64_t)nb*nloc_j;         // one panel
+    const uint64_t slotB  = (uint64_t)S*slotB1;
+    deviceVector<ComplexD> Abuf( slotA*Pc ? slotA*Pc : 1 );
+    deviceVector<ComplexD> Bbuf( slotB*Pr ? slotB*Pr : 1 );
+
+    deviceVector<ComplexD *> ap(1), bp(1), cp(1);
+    std::vector<ComplexD *>  ptr(1);
+
+    int firstblock = 1;
+    for(int64_t r0=kb0; r0<kb1; r0+=Pc){       // rounds of Pc k-blocks
+      int64_t r1 = std::min(kb1, r0+Pc);
+
+      /////////////////////////////////////////////////////////////////////
+      // Pack MY panels of this round into my origin slots.
+      /////////////////////////////////////////////////////////////////////
+      for(int64_t s=r0; s<r1; s++){
+        int64_t nb_s = L.BlockSize(s);
+        if ( (int)(s%Pc) == pcol && mloc_i ){  // my A panel: block-column s
+          int64_t lc0, lc1;
+          L.ColRange(s*nb, std::min(N,(s+1)*nb), lc0, lc1);
+          GRID_ASSERT( lc1-lc0 == nb_s );
+          ComplexD *src = A.LocalWindow(li0, lc0);
+          ComplexD *dst = &Abuf[0] + slotA*pcol;
+          int64_t   ld  = A.layout.mloc;
+          int64_t   m   = mloc_i;
+          accelerator_for(idx, (uint64_t)(m*nb_s), 1, {
+            int64_t jj = idx / m;
+            int64_t ii = idx - jj*m;
+            dst[ii + jj*m] = src[ii + jj*ld];
+          });
+        }
+        if ( (int)(s%Pr) == prow && nloc_j ){  // my B panel: block-row s
+          int64_t lr0, lr1;
+          L.RowRange(s*nb, std::min(N,(s+1)*nb), lr0, lr1);
+          GRID_ASSERT( lr1-lr0 == nb_s );
+          int64_t idxs = (s - r0 - ((prow - r0%Pr + Pr) % Pr)) / Pr; // my panel # in round
+          GRID_ASSERT( idxs >= 0 ); GRID_ASSERT( idxs < S );
+          ComplexD *src = B.LocalWindow(lr0, lj0);
+          ComplexD *dst = &Bbuf[0] + slotB*prow + slotB1*idxs;
+          int64_t   ld  = B.layout.mloc;
+          int64_t   nn  = nloc_j;
+          accelerator_for(idx, (uint64_t)(nb_s*nn), 1, {
+            int64_t jj = idx / nb_s;
+            int64_t ii = idx - jj*nb_s;
+            dst[ii + jj*nb_s] = src[ii + jj*ld];
+          });
+        }
+      }
+      accelerator_barrier();
+
+      /////////////////////////////////////////////////////////////////////
+      // Ring allgather along my process ROW: Pc-1 symmetric steps.  At
+      // step t send the slot of origin (pcol-t+1), receive origin (pcol-t).
+      /////////////////////////////////////////////////////////////////////
+      if ( Pc > 1 && slotA ){
+        int dest = prow*Pc + (pcol+1)%Pc;
+        int src  = prow*Pc + (pcol-1+Pc)%Pc;
+        for(int t=1;t<Pc;t++){
+          int cs = (pcol - t + 1 + Pc*Pc) % Pc;
+          int cr = (pcol - t     + Pc*Pc) % Pc;
+          grid->SendToRecvFrom((void *)(&Abuf[0]+slotA*cs), dest,
+                               (void *)(&Abuf[0]+slotA*cr), src,
+                               slotA*sizeof(ComplexD));
+        }
+      }
+      /////////////////////////////////////////////////////////////////////
+      // Ring allgather along my process COLUMN: Pr-1 symmetric steps.
+      /////////////////////////////////////////////////////////////////////
+      if ( Pr > 1 && slotB ){
+        int dest = ((prow+1)%Pr)*Pc + pcol;
+        int src  = ((prow-1+Pr)%Pr)*Pc + pcol;
+        for(int t=1;t<Pr;t++){
+          int rs = (prow - t + 1 + Pr*Pr) % Pr;
+          int rr = (prow - t     + Pr*Pr) % Pr;
+          grid->SendToRecvFrom((void *)(&Bbuf[0]+slotB*rs), dest,
+                               (void *)(&Bbuf[0]+slotB*rr), src,
+                               slotB*sizeof(ComplexD));
+        }
+      }
+
+      /////////////////////////////////////////////////////////////////////
+      // Local update, ascending s: fixed order, bitwise-reproducible.
+      /////////////////////////////////////////////////////////////////////
+      for(int64_t s=r0; s<r1; s++){
+        int64_t nb_s = L.BlockSize(s);
+        if ( !(mloc_i && nloc_j && nb_s) ) { firstblock = 0; continue; }
+        int cA = (int)(s%Pc);
+        int rB = (int)(s%Pr);
+        int64_t idxs = (s - r0 - ((rB - r0%Pr + Pr) % Pr)) / Pr;
+        ComplexD beta_use = firstblock ? beta : ComplexD(1.0,0.0);
+        firstblock = 0;
+
+        ptr[0] = &Abuf[0] + slotA*cA;
+        acceleratorCopyToDevice(&ptr[0], &ap[0], sizeof(ComplexD *));
+        ptr[0] = &Bbuf[0] + slotB*rB + slotB1*idxs;
+        acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexD *));
+        ptr[0] = C.LocalWindow(li0, lj0);
+        acceleratorCopyToDevice(&ptr[0], &cp[0], sizeof(ComplexD *));
+
+        BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N,
+                         (int)mloc_i, (int)nloc_j, (int)nb_s,
+                         alpha,    ap, (int)mloc_i,
+                                   bp, (int)nb_s,
+                         beta_use, cp, (int)C.layout.mloc);
+        BLAS.synchronise();
+      }
+    }
+  }
+};
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/DenseCoarseMatrix.h

@@ -30,6 +30,8 @@ Author: Peter Boyle <pboyle@bnl.gov>
 #include <Grid/algorithms/blas/BatchedBlas.h>
 #include <Grid/algorithms/blas/BatchedInverse.h>
 #include <Grid/algorithms/multigrid/RecursiveSchurInverse.h>
+#include <Grid/algorithms/multigrid/BlockCyclicSchurInverse.h>
+#include <Grid/algorithms/multigrid/BlockCyclicRedistribute.h>
 
 #include <unordered_map>
 
@@ -874,13 +876,48 @@ public:
     BlockRows S;
     ImportDenseFP64(Op, S, g2rm);
 
+    ////////////////////////////////////////////////////////////////
+    // DENSE_SCHUR2D=1 : invert via the 2D block-cyclic recursion
+    // (BlockCyclicSchurInverse) instead of the 1D rank-range one.
+    // The SAME imported rank-major rows S go in and come back, so the
+    // import certificate above and the slab rounding / VERIFY below are
+    // identical for both paths: a clean A/B on one imported operator.
+    //
+    // Everything in the 2D path -- redistribution, SUMMA rings, leaf --
+    // is point-to-point SendToRecvFrom; no collectives at all.
+    // DENSE_NB overrides the block size (default: rows-per-rank, which
+    // makes the redistribution edges maximally regular).
+    ////////////////////////////////////////////////////////////////
+    int use2d = getenv("DENSE_SCHUR2D") ? atoi(getenv("DENSE_SCHUR2D")) : 0;
     int64_t panelBytes = getenv("DENSE_PANEL_BYTES") ? atol(getenv("DENSE_PANEL_BYTES"))
-                                                     : (int64_t)1024*1024*1024;
-    RecursiveSchurInverse RSI(grid, N, rowStart, panelBytes);
-    double t2 = usecond();
-    RSI.Invert(S);
-    double t3 = usecond();
-    RSI.ReportTelemetry();
+                                                     : (int64_t)1024*1024*1024;   // 1D path only
+    double t2, t3;
+    if ( use2d )
+    {
+      int Pr,Pc;
+      BlockCyclicLayout::ChooseProcessGrid(P, Pr, Pc);
+      int64_t nb = getenv("DENSE_NB") ? atol(getenv("DENSE_NB")) : nrows;
+      GRID_ASSERT( nb >= 1 );
+      std::cout << GridLogMessage << "DenseCoarseMatrix: 2D SCHUR invert, process grid "
+                << Pr << " x " << Pc << "  nb " << nb
+                << "  (pure P2P: redistribute + SUMMA rings + local leaves)" << std::endl;
+      BlockCyclicMatrix A2(grid, N, nb, Pr, Pc);
+      BlockCyclicSchurInverse RSI2;
+      t2 = usecond();
+      BlockCyclicRedistribute::RowsToCyclic(grid, rowStart, &S.data[0], nrows, A2);
+      RSI2.Invert(A2);
+      BlockCyclicRedistribute::CyclicToRows(grid, rowStart, A2, &S.data[0], nrows);
+      t3 = usecond();
+      RSI2.ReportTelemetry(grid);
+    }
+    else
+    {
+      RecursiveSchurInverse RSI(grid, N, rowStart, panelBytes);
+      t2 = usecond();
+      RSI.Invert(S);
+      t3 = usecond();
+      RSI.ReportTelemetry();
+    }
 
     // The single terminal rounding: fp64 inverse -> fp32 apply slab
     // (row-major, global columns)

tests/debug/Test_blockcyclic.cc

@@ -0,0 +1,280 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_blockcyclic.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//////////////////////////////////////////////////////////////////////////////
+// Regression gate for BlockCyclicLayout -- stage 1 of the 2D distributed
+// dense inverse (documentation/DistributedDenseInverse2D.tex).
+//
+// The layout is pure index arithmetic, so this test is EXHAUSTIVE rather
+// than statistical: every stage sweeps a battery of (N, nb, Pr, Pc)
+// configurations chosen for their edge cases -- nb=1, nb>N, N%nb!=0,
+// more processes than blocks, prime N, the production shape -- and checks
+// every global index (T2/T3) or every element of an outer-product grid (T4)
+// against a brute-force reference built by walking blocks.
+//
+// No communication: correct at mpirun -n 1, and identical at any rank count.
+//
+//   T1 : NumLocal partitions N (sums over coords; against brute force).
+//   T2 : GlobalToLocal / LocalToGlobal round trip, every index.
+//   T3 : local indices are dense [0,mloc): a bijection, not just a cover.
+//   T4 : 2D ownership: every element has exactly one owner; per-rank counts
+//        equal mloc*nloc; LocalOffset is a bijection onto [0,mloc*nloc).
+//   T5 : block contiguity: within any owned global block, consecutive
+//        global rows are consecutive local rows (what SUMMA panels rely on).
+//   T6 : ChooseProcessGrid: exact factorisation, Pr<=Pc, most-square.
+//////////////////////////////////////////////////////////////////////////////
+
+#include <Grid/Grid.h>
+#include <Grid/algorithms/multigrid/BlockCyclic.h>
+
+using namespace Grid;
+
+static int failures = 0;
+
+static void Report(const std::string &name, bool pass, const std::string &detail="")
+{
+  std::cout << GridLogMessage << "  " << name << (pass ? "  PASS" : "  ** FAIL **");
+  if ( detail.size() ) std::cout << "   " << detail;
+  std::cout << std::endl;
+  if ( !pass ) failures++;
+}
+
+// One-dimensional configurations: {N, nb, Pg}
+struct Cfg1d { int64_t N; int64_t nb; int Pg; };
+static const std::vector<Cfg1d> configs1d = {
+  {  0,  1, 1},          // empty matrix
+  {  1,  1, 1},          // minimal
+  { 16,  4, 4},          // exact tiling
+  { 17,  4, 4},          // trailing partial block
+  { 16,  1, 4},          // nb=1: pure cyclic
+  { 16, 32, 4},          // nb>N: one short block, most coords empty
+  { 16,  4, 7},          // more processes than blocks
+  { 97, 13, 5},          // prime N, awkward everything
+  {138240, 480, 16},     // production N, row dimension of 16x18
+  {138240, 480, 18},     // production N, column dimension of 16x18
+  {138240, 512, 16},     // production N, non-dividing block
+};
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  std::cout << GridLogMessage << "BlockCyclicLayout regression (communicator-free index arithmetic)" << std::endl;
+
+  ////////////////////////////////////////////////////////////////////////
+  // T1 : NumLocal partitions N, and agrees with brute force block-walking.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true;
+    for(auto &c : configs1d){
+      // brute force: walk global blocks, count elements per coordinate
+      std::vector<int64_t> counts(c.Pg, 0);
+      for(int64_t b=0; b*c.nb < c.N || (c.N==0 && b<0); b++){
+        if ( b*c.nb >= c.N ) break;
+        int64_t lo = b*c.nb;
+        int64_t hi = std::min(c.N, lo+c.nb);
+        counts[b % c.Pg] += hi-lo;
+      }
+      int64_t sum = 0;
+      for(int p=0;p<c.Pg;p++){
+        int64_t n = BlockCyclicLayout::NumLocal(c.N, c.nb, p, c.Pg);
+        if ( n != counts[p] ) ok = false;
+        sum += n;
+      }
+      if ( sum != c.N ) ok = false;
+    }
+    Report("T1  NumLocal partitions N, == brute force", ok);
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T2 : round trip, every global index of every configuration.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true;
+    for(auto &c : configs1d){
+      int64_t sweep = std::min<int64_t>(c.N, 200000); // full N except production
+      for(int64_t g=0; g<sweep; g++){
+        int p; int64_t l;
+        BlockCyclicLayout::GlobalToLocal(g, c.nb, c.Pg, p, l);
+        if ( BlockCyclicLayout::LocalToGlobal(l, c.nb, p, c.Pg) != g ) ok = false;
+      }
+      // and the production tail, where the arithmetic could overflow or drift
+      for(int64_t g=std::max<int64_t>(0,c.N-1000); g<c.N; g++){
+        int p; int64_t l;
+        BlockCyclicLayout::GlobalToLocal(g, c.nb, c.Pg, p, l);
+        if ( BlockCyclicLayout::LocalToGlobal(l, c.nb, p, c.Pg) != g ) ok = false;
+      }
+    }
+    Report("T2  GlobalToLocal <-> LocalToGlobal round trip", ok);
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T3 : per coordinate, the local indices hit [0, NumLocal) exactly once.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true;
+    for(auto &c : configs1d){
+      if ( c.N > 4096 ) continue;        // dense bitmap check: small cfgs only
+      for(int p=0;p<c.Pg;p++){
+        int64_t n = BlockCyclicLayout::NumLocal(c.N, c.nb, p, c.Pg);
+        std::vector<int> hit(n, 0);
+        for(int64_t g=0; g<c.N; g++){
+          int pp; int64_t l;
+          BlockCyclicLayout::GlobalToLocal(g, c.nb, c.Pg, pp, l);
+          if ( pp != p ) continue;
+          if ( l < 0 || l >= n ) { ok = false; continue; }
+          hit[l]++;
+        }
+        for(int64_t l=0;l<n;l++) if ( hit[l] != 1 ) ok = false;
+      }
+    }
+    Report("T3  local indices dense and unique on [0,NumLocal)", ok);
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T4 : 2D ownership.  Small grids, every element: exactly one owner,
+  // owner counts equal mloc*nloc, LocalOffset bijective onto storage.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true;
+    struct Cfg2d { int64_t N; int64_t nb; int Pr; int Pc; };
+    std::vector<Cfg2d> cfgs = {
+      { 16, 4, 2, 2},
+      { 17, 4, 2, 3},                    // partial block, non-square grid
+      { 23, 5, 3, 2},
+      { 12, 2, 3, 4},                    // 12 ranks
+      { 30, 7, 4, 2},
+    };
+    for(auto &c : cfgs){
+      int P = c.Pr*c.Pc;
+      std::vector<BlockCyclicLayout> L;
+      for(int r=0;r<P;r++) L.push_back(BlockCyclicLayout(c.N,c.nb,c.Pr,c.Pc,r));
+      // ownership count per rank, and per-rank storage bitmap
+      std::vector<int64_t> owned(P,0);
+      std::vector<std::vector<int>> slot(P);
+      for(int r=0;r<P;r++) slot[r].assign(L[r].mloc*L[r].nloc, 0);
+      for(int64_t i=0;i<c.N;i++){
+        for(int64_t j=0;j<c.N;j++){
+          int r = L[0].OwnerRank(i,j);
+          if ( r < 0 || r >= P ) { ok=false; continue; }
+          // every layout instance must agree on the owner
+          if ( L[r].OwnerRank(i,j) != r ) ok = false;
+          if ( !L[r].Owns(i,j) )          ok = false;
+          owned[r]++;
+          int64_t off = L[r].LocalOffset(i,j);
+          if ( off < 0 || off >= (int64_t)slot[r].size() ) { ok=false; continue; }
+          slot[r][off]++;
+        }
+      }
+      int64_t tot=0;
+      for(int r=0;r<P;r++){
+        if ( owned[r] != L[r].mloc*L[r].nloc ) ok = false;
+        for(auto h : slot[r]) if ( h != 1 ) ok = false;
+        tot += owned[r];
+      }
+      if ( tot != c.N*c.N ) ok = false;
+    }
+    Report("T4  2D ownership: unique owner, counts == mloc*nloc, offsets bijective", ok);
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T5 : block contiguity.  For every owned global block, consecutive
+  // global rows are consecutive local rows: g and g+1 in the same block
+  // must give l and l+1 on the same coordinate.  SUMMA panel extraction
+  // (stage 2) sends whole local blocks as contiguous strides through the
+  // column-major store; that only works if this holds.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true;
+    for(auto &c : configs1d){
+      if ( c.N == 0 ) continue;
+      int64_t sweep = std::min<int64_t>(c.N-1, 100000);
+      for(int64_t g=0; g<sweep; g++){
+        if ( (g+1) % c.nb == 0 ) continue;   // block boundary: owner may change
+        int p0,p1; int64_t l0,l1;
+        BlockCyclicLayout::GlobalToLocal(g,   c.nb, c.Pg, p0, l0);
+        BlockCyclicLayout::GlobalToLocal(g+1, c.nb, c.Pg, p1, l1);
+        if ( p1 != p0     ) ok = false;
+        if ( l1 != l0 + 1 ) ok = false;
+      }
+    }
+    Report("T5  intra-block contiguity (SUMMA panel precondition)", ok);
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T6 : ChooseProcessGrid.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true;
+    for(int P : {1,2,3,4,6,8,12,16,17,64,96,144,256,288,512}){
+      int Pr,Pc;
+      BlockCyclicLayout::ChooseProcessGrid(P,Pr,Pc);
+      if ( Pr*Pc != P ) ok = false;
+      if ( Pr > Pc )    ok = false;
+      // most-square: no divisor r with Pr < r <= sqrt(P)
+      for(int r=Pr+1; (int64_t)r*r <= (int64_t)P; r++)
+        if ( P % r == 0 ) ok = false;
+    }
+    int Pr,Pc;
+    BlockCyclicLayout::ChooseProcessGrid(288,Pr,Pc);
+    if ( !(Pr==16 && Pc==18) ) ok = false;
+    Report("T6  ChooseProcessGrid exact, Pr<=Pc, most-square (288 -> 16x18)", ok);
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T7 : RangeToLocal.  For every block-aligned range of every small
+  // configuration: the owned global indices of [g0,g1) map exactly onto
+  // local [l0,l1), contiguously and in order (brute force).
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true;
+    for(auto &c : configs1d){
+      if ( c.N == 0 || c.N > 4096 ) continue;
+      int64_t nblocks = (c.N + c.nb - 1)/c.nb;
+      for(int p=0;p<c.Pg;p++){
+        for(int64_t b0=0;b0<=nblocks;b0++){
+          for(int64_t b1=b0;b1<=nblocks;b1++){
+            int64_t g0 = b0*c.nb;
+            int64_t g1 = std::min(c.N, b1*c.nb);
+            if ( g0 > c.N ) continue;
+            int64_t l0,l1;
+            BlockCyclicLayout::RangeToLocal(g0,g1,c.N,c.nb,p,c.Pg,l0,l1);
+            // brute force: owned globals in [g0,g1) in ascending order
+            int64_t expect = l0;
+            for(int64_t g=g0; g<g1; g++){
+              int pp; int64_t l;
+              BlockCyclicLayout::GlobalToLocal(g, c.nb, c.Pg, pp, l);
+              if ( pp != p ) continue;
+              if ( l != expect ) ok = false;   // contiguous, in order
+              expect++;
+            }
+            if ( expect != l1 ) ok = false;    // count matches the bounds
+          }
+        }
+      }
+    }
+    Report("T7  RangeToLocal contiguous, ordered, exact bounds", ok);
+  }
+
+  std::cout << GridLogMessage << (failures ? "Test_blockcyclic: FAILURES"
+                                           : "Test_blockcyclic: ALL PASS") << std::endl;
+  Grid_finalize();
+  return failures ? 1 : 0;
+}

tests/debug/Test_general_coarse.cc

@@ -36,8 +36,9 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 using namespace std;
 using namespace Grid;
 
-gridblasHandle_t GridBLAS::gridblasHandle;
-int            GridBLAS::gridblasInit;
+// GridBLAS statics are defined ONCE, in Grid/algorithms/blas/BatchedBlas.cc;
+// a second definition here is a duplicate-symbol link error the moment the
+// archive member is pulled in.
 
 ///////////////////////
 // Tells little dirac op to use MdagM as the .Op()

tests/debug/Test_schur2d.cc

@@ -0,0 +1,204 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_schur2d.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//////////////////////////////////////////////////////////////////////////////
+// Regression gate for BlockCyclicSchurInverse -- stage 3 of the 2D
+// distributed dense inverse.  CPU build under mpirun:
+//
+//   mpirun -n 1 ./Test_schur2d --grid 8.8.8.8  --mpi 1.1.1.1
+//   mpirun -n 2 ./Test_schur2d --grid 8.8.8.8  --mpi 1.1.1.2
+//   mpirun -n 3 ./Test_schur2d --grid 8.8.8.12 --mpi 1.1.1.3
+//   mpirun -n 4 ./Test_schur2d --grid 8.8.8.8  --mpi 1.1.1.4
+//
+// Sweeps all process-grid factorisations of P and a battery of (N,nb)
+// including ragged trailing blocks, a single-leaf matrix (nblocks==1),
+// and nb=3 with many blocks.  The matrices are diagonally dominant --
+// the recursion does not pivot, exactly like the 1D implementation, and
+// the test respects that contract.
+//
+//   T1 : certificate  max|A . Ainv - I|  with the product computed by the
+//        (independently validated) distributed SUMMA on an untouched copy.
+//   T2 : element-wise against a host Gauss-Jordan reference inverse
+//        (partial pivoting, fp64).
+//   T3 : repeated inversion bitwise identical -- the P2P determinism
+//        property, which a collective-reduce implementation cannot offer.
+//////////////////////////////////////////////////////////////////////////////
+
+#include <Grid/Grid.h>
+#include <Grid/algorithms/multigrid/BlockCyclicSchurInverse.h>
+
+using namespace Grid;
+
+static int failures = 0;
+
+static void Report(const std::string &name, bool pass, const std::string &detail="")
+{
+  std::cout << GridLogMessage << "  " << name << (pass ? "  PASS" : "  ** FAIL **");
+  if ( detail.size() ) std::cout << "   " << detail;
+  std::cout << std::endl;
+  if ( !pass ) failures++;
+}
+
+static ComplexD Fill(int64_t i, int64_t j, int salt)
+{
+  double x = std::sin(0.7*i + 1.3*j + 0.31*salt);
+  double y = std::cos(1.9*i - 0.4*j + 0.77*salt);
+  return ComplexD(x,y);
+}
+
+// Diagonally dominant test matrix: no pivoting required at any depth.
+static void MakeMatrix(std::vector<ComplexD> &A, int64_t N, int salt)
+{
+  A.resize((uint64_t)N*N);
+  for(int64_t j=0;j<N;j++)
+    for(int64_t i=0;i<N;i++)
+      A[i+j*N] = Fill(i,j,salt) + ((i==j) ? ComplexD(3.0*N, 0.5) : ComplexD(0.0,0.0));
+}
+
+// Host reference inverse: Gauss-Jordan with partial pivoting, fp64.
+static void HostInverse(std::vector<ComplexD> A, std::vector<ComplexD> &X, int64_t N)
+{
+  X.assign((uint64_t)N*N, ComplexD(0.0,0.0));
+  for(int64_t i=0;i<N;i++) X[i+i*N] = ComplexD(1.0,0.0);
+  for(int64_t c=0;c<N;c++){
+    int64_t piv=c; double mx = std::abs(A[c+c*N]);
+    for(int64_t r=c+1;r<N;r++)
+      if ( std::abs(A[r+c*N]) > mx ){ mx=std::abs(A[r+c*N]); piv=r; }
+    GRID_ASSERT( mx > 0.0 );
+    if ( piv != c )
+      for(int64_t j=0;j<N;j++){
+        std::swap(A[c+j*N],A[piv+j*N]);
+        std::swap(X[c+j*N],X[piv+j*N]);
+      }
+    ComplexD d = ComplexD(1.0,0.0)/A[c+c*N];
+    for(int64_t j=0;j<N;j++){ A[c+j*N]*=d; X[c+j*N]*=d; }
+    for(int64_t r=0;r<N;r++){
+      if ( r==c ) continue;
+      ComplexD f = A[r+c*N];
+      if ( f == ComplexD(0.0,0.0) ) continue;
+      for(int64_t j=0;j<N;j++){
+        A[r+j*N] -= f*A[c+j*N];
+        X[r+j*N] -= f*X[c+j*N];
+      }
+    }
+  }
+}
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+                                                       GridDefaultSimd(Nd, vComplexD::Nsimd()),
+                                                       GridDefaultMpi());
+  const int P = grid->ProcessorCount();
+
+  std::vector<std::pair<int,int>> grids;
+  for(int r=1;r<=P;r++) if ( P%r==0 ) grids.push_back({r,P/r});
+
+  struct Cfg { int64_t N; int64_t nb; };
+  std::vector<Cfg> cfgs = { {24,4}, {26,4}, {17,5}, {30,7}, {8,8}, {33,3}, {40,5} };
+
+  BlockCyclicSchurInverse RSI;
+  BlockCyclicSumma        SUMMA;
+
+  std::cout << GridLogMessage << "BlockCyclicSchurInverse regression: P=" << P
+            << ", " << grids.size() << " process grids, " << cfgs.size()
+            << " layouts" << std::endl;
+
+  ////////////////////////////////////////////////////////////////////////
+  // T1 + T2, one pass: invert, certify with distributed SUMMA on an
+  // untouched copy, and compare against the host reference inverse.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool okC = true, okR = true;
+    double worstC = 0.0, worstR = 0.0;
+    for(auto &g : grids){
+      for(auto &c : cfgs){
+        int64_t N = c.N;
+        std::vector<ComplexD> Ag, Ref, Ainv, Cert;
+        MakeMatrix(Ag, N, 12);
+        HostInverse(Ag, Ref, N);
+
+        BlockCyclicMatrix A (grid,N,c.nb,g.first,g.second);
+        BlockCyclicMatrix A0(grid,N,c.nb,g.first,g.second);
+        BlockCyclicMatrix Ce(grid,N,c.nb,g.first,g.second);
+        A.ImportGlobal(Ag);
+        A0.ImportGlobal(Ag);
+
+        RSI.Invert(A);                       // in place: A now holds Ainv
+
+        // certificate: Ce = A0 . Ainv, distributed
+        SUMMA.Multiply(ComplexD(1.0,0.0),A0,A,ComplexD(0.0,0.0),Ce, 0,N,0,N,0,N);
+        Ce.ExportGlobal(Cert);
+        double dc = 0.0;
+        for(int64_t j=0;j<N;j++)
+          for(int64_t i=0;i<N;i++){
+            ComplexD id = (i==j) ? ComplexD(1.0,0.0) : ComplexD(0.0,0.0);
+            dc = std::max(dc, std::abs(Cert[i+j*N]-id));
+          }
+        worstC = std::max(worstC,dc);
+        if ( dc > 1.0e-10 ) okC = false;
+
+        // reference: element-wise, scaled by the largest inverse entry
+        A.ExportGlobal(Ainv);
+        double mxref = 0.0, dr = 0.0;
+        for(uint64_t i=0;i<Ref.size();i++) mxref = std::max(mxref, std::abs(Ref[i]));
+        for(uint64_t i=0;i<Ref.size();i++) dr = std::max(dr, std::abs(Ainv[i]-Ref[i]));
+        dr /= mxref;
+        worstR = std::max(worstR,dr);
+        if ( dr > 1.0e-9 ) okR = false;
+      }
+    }
+    Report("T1  certificate max|A.Ainv - I|, all grids x layouts", okC,
+           "worst "+std::to_string(worstC));
+    Report("T2  vs host Gauss-Jordan reference (relative)", okR,
+           "worst "+std::to_string(worstR));
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T3 : determinism.  Same import, two inversions, bitwise comparison.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true;
+    for(auto &g : grids){
+      int64_t N = 30, nb = 7;
+      std::vector<ComplexD> Ag, X1, X2;
+      MakeMatrix(Ag, N, 13);
+      BlockCyclicMatrix A(grid,N,nb,g.first,g.second);
+      A.ImportGlobal(Ag); RSI.Invert(A); A.ExportGlobal(X1);
+      A.ImportGlobal(Ag); RSI.Invert(A); A.ExportGlobal(X2);
+      for(uint64_t i=0;i<X1.size();i++)
+        if ( !(X1[i]==X2[i]) ) ok = false;   // BITWISE
+    }
+    Report("T3  repeated inversion bitwise identical", ok);
+  }
+
+  {
+    uint64_t f = failures;  grid->GlobalSum(f);
+    if ( f && !failures )
+      std::cout << GridLogMessage << "  ** failures on OTHER ranks: " << f << " **" << std::endl;
+    failures = (int)f;
+  }
+  std::cout << GridLogMessage << (failures ? "Test_schur2d: FAILURES"
+                                           : "Test_schur2d: ALL PASS") << std::endl;
+  Grid_finalize();
+  return failures ? 1 : 0;
+}

tests/debug/Test_schur2d_redist.cc

@@ -0,0 +1,259 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_schur2d_redist.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//////////////////////////////////////////////////////////////////////////////
+// Regression gate for BlockCyclicRedistribute and the full stage-4 path
+//
+//   1D rank-major rows -> block cyclic -> Invert -> back to 1D rows
+//
+// which is exactly what DENSE_SCHUR2D runs inside DenseCoarseMatrix.
+// CPU build under mpirun at n = 1,2,3,4.
+//
+//   T1 : RowsToCyclic against a direct ImportGlobal of the same matrix --
+//        BITWISE (pure data movement, no arithmetic).
+//   T2 : round trip rows -> 2D -> rows -- BITWISE, uniform AND non-uniform
+//        rowStart, layouts with ragged trailing blocks.
+//   T3 : full pipeline inverse against a host Gauss-Jordan reference.
+//   T4 : CROSS-IMPLEMENTATION: the same matrix inverted by the 1D
+//        RecursiveSchurInverse and by the 2D pipeline; results compared
+//        element-wise.  Two independent implementations, two independent
+//        decompositions, one answer.
+//////////////////////////////////////////////////////////////////////////////
+
+#include <Grid/Grid.h>
+#include <Grid/algorithms/multigrid/RecursiveSchurInverse.h>
+#include <Grid/algorithms/multigrid/BlockCyclicSchurInverse.h>
+#include <Grid/algorithms/multigrid/BlockCyclicRedistribute.h>
+
+using namespace Grid;
+
+static int failures = 0;
+
+static void Report(const std::string &name, bool pass, const std::string &detail="")
+{
+  std::cout << GridLogMessage << "  " << name << (pass ? "  PASS" : "  ** FAIL **");
+  if ( detail.size() ) std::cout << "   " << detail;
+  std::cout << std::endl;
+  if ( !pass ) failures++;
+}
+
+static ComplexD Fill(int64_t i, int64_t j, int salt)
+{
+  double x = std::sin(0.7*i + 1.3*j + 0.31*salt);
+  double y = std::cos(1.9*i - 0.4*j + 0.77*salt);
+  return ComplexD(x,y);
+}
+static void MakeMatrix(std::vector<ComplexD> &A, int64_t N, int salt)
+{
+  A.resize((uint64_t)N*N);
+  for(int64_t j=0;j<N;j++)
+    for(int64_t i=0;i<N;i++)
+      A[i+j*N] = Fill(i,j,salt) + ((i==j) ? ComplexD(3.0*N,0.5) : ComplexD(0.0,0.0));
+}
+static void HostInverse(std::vector<ComplexD> A, std::vector<ComplexD> &X, int64_t N)
+{
+  X.assign((uint64_t)N*N, ComplexD(0.0,0.0));
+  for(int64_t i=0;i<N;i++) X[i+i*N] = ComplexD(1.0,0.0);
+  for(int64_t c=0;c<N;c++){
+    int64_t piv=c; double mx = std::abs(A[c+c*N]);
+    for(int64_t r=c+1;r<N;r++)
+      if ( std::abs(A[r+c*N]) > mx ){ mx=std::abs(A[r+c*N]); piv=r; }
+    GRID_ASSERT( mx > 0.0 );
+    if ( piv != c )
+      for(int64_t j=0;j<N;j++){ std::swap(A[c+j*N],A[piv+j*N]); std::swap(X[c+j*N],X[piv+j*N]); }
+    ComplexD d = ComplexD(1.0,0.0)/A[c+c*N];
+    for(int64_t j=0;j<N;j++){ A[c+j*N]*=d; X[c+j*N]*=d; }
+    for(int64_t r=0;r<N;r++){
+      if ( r==c ) continue;
+      ComplexD f = A[r+c*N];
+      if ( f == ComplexD(0.0,0.0) ) continue;
+      for(int64_t j=0;j<N;j++){ A[r+j*N]-=f*A[c+j*N]; X[r+j*N]-=f*X[c+j*N]; }
+    }
+  }
+}
+
+// Ownership tables over the ranks: uniform-ish, and deliberately lopsided.
+static std::vector<int64_t> MakeRowStart(int64_t N, int P, int lopsided)
+{
+  std::vector<int64_t> t(P+1); t[0]=0;
+  for(int r=0;r<P;r++){
+    int64_t base = N/P;
+    int64_t extra = ( r < (int)(N%P) ) ? 1 : 0;
+    int64_t n = base + extra;
+    if ( lopsided ){                       // shift work toward low ranks
+      if ( r==0 && P>1 ) n = std::min(N, base+extra+3);
+      if ( r==P-1 )      n = N - t[r];     // remainder
+    }
+    t[r+1] = std::min(N, t[r]+n);
+  }
+  t[P]=N;
+  return t;
+}
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+                                                       GridDefaultSimd(Nd, vComplexD::Nsimd()),
+                                                       GridDefaultMpi());
+  const int P  = grid->ProcessorCount();
+  const int me = grid->ThisRank();
+
+  std::vector<std::pair<int,int>> grids;
+  for(int r=1;r<=P;r++) if ( P%r==0 ) grids.push_back({r,P/r});
+
+  struct Cfg { int64_t N; int64_t nb; };
+  std::vector<Cfg> cfgs = { {24,4}, {26,4}, {17,5}, {30,7}, {33,3} };
+
+  std::cout << GridLogMessage << "BlockCyclicRedistribute regression: P=" << P << std::endl;
+
+  ////////////////////////////////////////////////////////////////////////
+  // T1 + T2 : distribution correctness and round trip, bitwise.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok1 = true, ok2 = true;
+    for(int lop=0; lop<2; lop++){
+      for(auto &g : grids){
+        for(auto &c : cfgs){
+          int64_t N = c.N;
+          std::vector<int64_t> rowStart = MakeRowStart(N,P,lop);
+          int64_t myrows = rowStart[me+1]-rowStart[me];
+
+          std::vector<ComplexD> Ag;  MakeMatrix(Ag, N, 21);
+
+          // my 1D rows, column major ld = myrows
+          std::vector<ComplexD> h((uint64_t)std::max<int64_t>(myrows,1)*N);
+          for(int64_t j=0;j<N;j++)
+            for(int64_t i=0;i<myrows;i++)
+              h[i + j*myrows] = Ag[(rowStart[me]+i) + j*N];
+          deviceVector<ComplexD> rows1d(h.size());
+          acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD));
+
+          BlockCyclicMatrix A(grid,N,c.nb,g.first,g.second);
+          BlockCyclicMatrix R(grid,N,c.nb,g.first,g.second);
+          BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,A);
+          R.ImportGlobal(Ag);
+
+          // T1: bitwise against direct import
+          {
+            std::vector<ComplexD> x((uint64_t)A.layout.mloc*A.layout.nloc);
+            std::vector<ComplexD> y(x.size());
+            if ( x.size() ){
+              acceleratorCopyFromDevice(&A.data[0], &x[0], x.size()*sizeof(ComplexD));
+              acceleratorCopyFromDevice(&R.data[0], &y[0], y.size()*sizeof(ComplexD));
+            }
+            for(uint64_t i=0;i<x.size();i++) if ( !(x[i]==y[i]) ) ok1 = false;
+          }
+          // T2: round trip, bitwise
+          {
+            deviceVector<ComplexD> back(h.size());
+            std::vector<ComplexD>  hb(h.size(), ComplexD(0.0,0.0));
+            acceleratorCopyToDevice(&hb[0], &back[0], hb.size()*sizeof(ComplexD));
+            BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&back[0],myrows);
+            acceleratorCopyFromDevice(&back[0], &hb[0], hb.size()*sizeof(ComplexD));
+            for(int64_t j=0;j<N;j++)
+              for(int64_t i=0;i<myrows;i++)
+                if ( !(hb[i+j*myrows]==h[i+j*myrows]) ) ok2 = false;
+          }
+        }
+      }
+    }
+    Report("T1  RowsToCyclic == ImportGlobal, bitwise", ok1);
+    Report("T2  round trip rows->2D->rows, bitwise, incl. lopsided rowStart", ok2);
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T3 + T4 : the DENSE_SCHUR2D pipeline against the host reference and
+  // against the INDEPENDENT 1D RecursiveSchurInverse.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok3 = true, ok4 = true;
+    double worst3 = 0.0, worst4 = 0.0;
+    BlockCyclicSchurInverse RSI2;
+    for(auto &g : grids){
+      for(auto &c : cfgs){
+        int64_t N = c.N;
+        std::vector<int64_t> rowStart = MakeRowStart(N,P,0);
+        int64_t myrows = rowStart[me+1]-rowStart[me];
+
+        std::vector<ComplexD> Ag, Ref;
+        MakeMatrix(Ag, N, 22);
+        HostInverse(Ag, Ref, N);
+
+        std::vector<ComplexD> h((uint64_t)std::max<int64_t>(myrows,1)*N);
+        for(int64_t j=0;j<N;j++)
+          for(int64_t i=0;i<myrows;i++)
+            h[i + j*myrows] = Ag[(rowStart[me]+i) + j*N];
+
+        double mxref = 0.0;
+        for(auto &z : Ref) mxref = std::max(mxref, std::abs(z));
+
+        // ---- 2D pipeline: rows -> cyclic -> invert -> rows ----
+        std::vector<ComplexD> h2d(h.size());
+        {
+          deviceVector<ComplexD> rows1d(h.size());
+          acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD));
+          BlockCyclicMatrix A(grid,N,c.nb,g.first,g.second);
+          BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,A);
+          RSI2.Invert(A);
+          BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&rows1d[0],myrows);
+          acceleratorCopyFromDevice(&rows1d[0], &h2d[0], h2d.size()*sizeof(ComplexD));
+        }
+        for(int64_t j=0;j<N;j++)
+          for(int64_t i=0;i<myrows;i++){
+            double d = std::abs(h2d[i+j*myrows]-Ref[(rowStart[me]+i)+j*N])/mxref;
+            worst3 = std::max(worst3,d);
+            if ( d > 1.0e-9 ) ok3 = false;
+          }
+
+        // ---- 1D RecursiveSchurInverse on the same matrix ----
+        {
+          BlockRows Ar;  Ar.Resize(myrows, N);
+          acceleratorCopyToDevice(&h[0], &Ar.data[0], h.size()*sizeof(ComplexD));
+          std::vector<int64_t> rs = rowStart;
+          RecursiveSchurInverse RSI1(grid, N, rs, 1<<20);
+          RSI1.Invert(Ar);
+          std::vector<ComplexD> h1d(h.size());
+          acceleratorCopyFromDevice(&Ar.data[0], &h1d[0], h1d.size()*sizeof(ComplexD));
+          for(int64_t j=0;j<N;j++)
+            for(int64_t i=0;i<myrows;i++){
+              double d = std::abs(h2d[i+j*myrows]-h1d[i+j*myrows])/mxref;
+              worst4 = std::max(worst4,d);
+              if ( d > 1.0e-9 ) ok4 = false;
+            }
+        }
+      }
+    }
+    Report("T3  2D pipeline vs host reference", ok3, "worst "+std::to_string(worst3));
+    Report("T4  2D pipeline vs 1D RecursiveSchurInverse", ok4, "worst "+std::to_string(worst4));
+  }
+
+  {
+    uint64_t f = failures;  grid->GlobalSum(f);
+    if ( f && !failures )
+      std::cout << GridLogMessage << "  ** failures on OTHER ranks: " << f << " **" << std::endl;
+    failures = (int)f;
+  }
+  std::cout << GridLogMessage << (failures ? "Test_schur2d_redist: FAILURES"
+                                           : "Test_schur2d_redist: ALL PASS") << std::endl;
+  Grid_finalize();
+  return failures ? 1 : 0;
+}

tests/debug/Test_schur2d_scale.cc

@@ -0,0 +1,198 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_schur2d_scale.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//////////////////////////////////////////////////////////////////////////////
+// SCALE rehearsal for the 2D distributed dense inverse: the full
+// DENSE_SCHUR2D pipeline -- 1D rows -> redistribute -> invert ->
+// redistribute back -> certificate -- on a SYNTHETIC matrix of any size,
+// with no multigrid machinery, no configuration and no subspace file.
+//
+// This is the missing rung between the small-N oracle tests (which build
+// the whole matrix on every host, impossible at production N) and the
+// production example (which needs the full setup and 100 s of job time
+// before the inverse is even reached).  Everything here is O(N^2/P) per
+// rank; at N=138240 on 288 ranks it is the production problem shape
+// exactly, in a driver that runs in minutes.
+//
+//   S2D_N   : global dimension          (default 720, laptop friendly)
+//   S2D_NB  : block size                (default N/P rows-per-rank if that
+//                                        is exact, else 48)
+//
+// The matrix is diagonally dominant (the recursion does not pivot); its
+// conditioning is BENIGN, so this rehearses scale and speed, not the real
+// operator's numerics -- the production VERIFY does that.
+//
+// Certificate: Cert = A0 . Ainv by the (independently validated) SUMMA,
+// then every rank checks ITS OWN local elements against the identity.
+// One GlobalMax at the end to report; the pipeline itself is pure P2P.
+//
+//   T1 : max|A.Ainv - I|  <  1e-8
+//   T2 : round-trip redistribution of the INVERSE bitwise consistent
+//        (CyclicToRows then RowsToCyclic reproduces the device data).
+//////////////////////////////////////////////////////////////////////////////
+
+#include <Grid/Grid.h>
+#include <Grid/algorithms/multigrid/BlockCyclicSchurInverse.h>
+#include <Grid/algorithms/multigrid/BlockCyclicRedistribute.h>
+
+using namespace Grid;
+
+static ComplexD Fill(int64_t i, int64_t j, int64_t N)
+{
+  double x = std::sin(0.7*i + 1.3*j);
+  double y = std::cos(1.9*i - 0.4*j);
+  if ( i==j ) return ComplexD(3.0*64 + x, 0.5);   // dominance independent of N
+  // band-limit the off-diagonal so row sums stay bounded as N grows:
+  // only |i-j| <= 64 entries are non-zero => sum |offdiag| <= 128*1.42 < 3*64
+  if ( std::abs((double)(i-j)) > 64.0 ) return ComplexD(0.0,0.0);
+  return ComplexD(x,y);
+}
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+                                                       GridDefaultSimd(Nd, vComplexD::Nsimd()),
+                                                       GridDefaultMpi());
+  const int P  = grid->ProcessorCount();
+  const int me = grid->ThisRank();
+
+  int64_t N  = getenv("S2D_N")  ? atol(getenv("S2D_N"))  : 720;
+  int64_t nb;
+  if      ( getenv("S2D_NB") ) nb = atol(getenv("S2D_NB"));
+  else if ( N % P == 0 )       nb = N/P;
+  else                         nb = 48;
+  GRID_ASSERT( N >= 1 ); GRID_ASSERT( nb >= 1 );
+
+  int Pr,Pc;
+  BlockCyclicLayout::ChooseProcessGrid(P,Pr,Pc);
+
+  // uniform-as-possible 1D ownership, as the production import produces
+  std::vector<int64_t> rowStart(P+1); rowStart[0]=0;
+  for(int r=0;r<P;r++) rowStart[r+1] = rowStart[r] + N/P + ( r < (int)(N%P) ? 1 : 0 );
+  int64_t myrows = rowStart[me+1]-rowStart[me];
+  int64_t row0   = rowStart[me];
+
+  std::cout << GridLogMessage << "Test_schur2d_scale: N=" << N << " nb=" << nb
+            << "  grid " << Pr << "x" << Pc << "  rows/rank ~" << myrows
+            << "  matrix " << (double)N*N*16.0/1.0e9 << " GB global, "
+            << (double)myrows*N*16.0/1.0e9 << " GB/rank rows" << std::endl;
+
+  ////////////////////////////////////////////////////////////////////////
+  // Fill MY rows only: O(N^2/P) host work, no global matrix anywhere.
+  ////////////////////////////////////////////////////////////////////////
+  double t0 = usecond();
+  std::vector<ComplexD> h((uint64_t)std::max<int64_t>(myrows,1)*N);
+  thread_for(jj, N, {
+    for(int64_t i=0;i<myrows;i++) h[i + jj*myrows] = Fill(row0+i, jj, N);
+  });
+  deviceVector<ComplexD> rows1d(h.size());
+  acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD));
+  double t1 = usecond();
+
+  ////////////////////////////////////////////////////////////////////////
+  // The DENSE_SCHUR2D pipeline, phase-timed.  A0 keeps the original for
+  // the certificate.
+  ////////////////////////////////////////////////////////////////////////
+  BlockCyclicMatrix A (grid,N,nb,Pr,Pc);
+  BlockCyclicMatrix A0(grid,N,nb,Pr,Pc);
+  BlockCyclicSchurInverse RSI2;
+
+  BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,A);
+  double t2 = usecond();
+  if ( A.data.size() )
+    acceleratorCopyDeviceToDevice((void *)&A.data[0],(void *)&A0.data[0],
+                                  A.data.size()*sizeof(ComplexD));
+  double t3 = usecond();
+  RSI2.Invert(A);
+  double t4 = usecond();
+  BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&rows1d[0],myrows);
+  double t5 = usecond();
+  RSI2.ReportTelemetry(grid);
+
+  ////////////////////////////////////////////////////////////////////////
+  // T1 : certificate by SUMMA, checked locally, reported by one GlobalMax.
+  ////////////////////////////////////////////////////////////////////////
+  int failures = 0;
+  {
+    BlockCyclicMatrix Cert(grid,N,nb,Pr,Pc);
+    BlockCyclicSumma  SUMMA;
+    SUMMA.Multiply(ComplexD(1.0,0.0),A0,A,ComplexD(0.0,0.0),Cert, 0,N,0,N,0,N);
+    double t6 = usecond();
+
+    BlockCyclicLayout &L = Cert.layout;
+    std::vector<ComplexD> hc((uint64_t)std::max<int64_t>(L.mloc*L.nloc,1));
+    if ( L.mloc*L.nloc )
+      acceleratorCopyFromDevice(&Cert.data[0], &hc[0],
+                                (uint64_t)L.mloc*L.nloc*sizeof(ComplexD));
+    double mx = 0.0;
+    for(int64_t lj=0;lj<L.nloc;lj++){
+      int64_t gj = BlockCyclicLayout::LocalToGlobal(lj, nb, L.pcol, Pc);
+      for(int64_t li=0;li<L.mloc;li++){
+        int64_t gi = BlockCyclicLayout::LocalToGlobal(li, nb, L.prow, Pr);
+        ComplexD id = (gi==gj) ? ComplexD(1.0,0.0) : ComplexD(0.0,0.0);
+        mx = std::max(mx, std::abs(hc[li+lj*L.mloc]-id));
+      }
+    }
+    RealD gmx = mx;  grid->GlobalMax(gmx);
+    if ( gmx > 1.0e-8 ) failures++;
+
+    std::cout << GridLogMessage << "Test_schur2d_scale phases (s):"
+              << "  fill "         << (t1-t0)/1.0e6
+              << "  redist->2D "   << (t2-t1)/1.0e6
+              << "  invert "       << (t4-t3)/1.0e6
+              << "  redist->rows " << (t5-t4)/1.0e6
+              << "  certify "      << (t6-t5)/1.0e6 << std::endl;
+    std::cout << GridLogMessage << "Test_schur2d_scale CERTIFICATE max|A.Ainv - I| = "
+              << gmx << (gmx > 1.0e-8 ? "   ** FAIL **" : "   PASS") << std::endl;
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T2 : the rows now hold the inverse; push them back out and compare
+  // against A on device -- redistribution must be bitwise invertible on
+  // real (non-synthetic-import) data too.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    BlockCyclicMatrix B(grid,N,nb,Pr,Pc);
+    BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,B);
+    std::vector<ComplexD> x((uint64_t)std::max<int64_t>(A.layout.mloc*A.layout.nloc,1));
+    std::vector<ComplexD> y(x.size());
+    if ( A.layout.mloc*A.layout.nloc ){
+      acceleratorCopyFromDevice(&A.data[0], &x[0], x.size()*sizeof(ComplexD));
+      acceleratorCopyFromDevice(&B.data[0], &y[0], y.size()*sizeof(ComplexD));
+    }
+    int bad = 0;
+    for(uint64_t i=0;i<x.size();i++) if ( !(x[i]==y[i]) ) bad++;
+    uint64_t gbad = bad;  grid->GlobalSum(gbad);
+    if ( gbad ) failures++;
+    std::cout << GridLogMessage << "Test_schur2d_scale ROUND TRIP mismatches = "
+              << gbad << (gbad ? "   ** FAIL **" : "   PASS") << std::endl;
+  }
+
+  {
+    uint64_t f = failures;  grid->GlobalSum(f);
+    failures = (int)f;
+  }
+  std::cout << GridLogMessage << (failures ? "Test_schur2d_scale: FAILURES"
+                                           : "Test_schur2d_scale: ALL PASS") << std::endl;
+  Grid_finalize();
+  return failures ? 1 : 0;
+}

tests/debug/Test_schur_dense_coarse.cc

@@ -262,8 +262,11 @@ int main (int argc, char ** argv)
   setenv("DENSE_PANEL_BYTES","65536",1);
   unsetenv("SLAB_FILE");
 
-  typedef DenseCoarseMatrix<vSpinColourVector,vTComplex,nbasis> DenseCC;
-  DenseCC dcm(LittleDiracOp, Coarse5d);
+  // Restructured interface: DenseCoarseMatrix<CComplex,nbasis>, constructed
+  // on the grid and fed by Import (which runs the certificate chain).
+  typedef DenseCoarseMatrix<vTComplex,nbasis> DenseCC;
+  DenseCC dcm(Coarse5d);
+  dcm.Import(LittleDiracOp);
 
   std::cout << GridLogMessage << "T6  audit relative slab difference (schur vs single) = "
             << dcm.schurAuditRel << std::endl;

tests/debug/Test_summa.cc

@@ -0,0 +1,270 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_summa.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//////////////////////////////////////////////////////////////////////////////
+// Regression gate for BlockCyclicSumma -- stage 2 of the 2D distributed
+// dense inverse.  CPU build under mpirun:
+//
+//   mpirun -n 1 ./Test_summa --grid 8.8.8.8 --mpi 1.1.1.1
+//   mpirun -n 2 ./Test_summa --grid 8.8.8.8 --mpi 1.1.1.2
+//   mpirun -n 3 ./Test_summa --grid 8.8.8.12 --mpi 1.1.1.3
+//   mpirun -n 4 ./Test_summa --grid 8.8.8.8 --mpi 1.1.1.4
+//
+// Every stage sweeps all process-grid factorisations of P (including the
+// degenerate 1xP and Px1 rings) and a battery of (N,nb) with ragged
+// trailing blocks, nb>N, and more processes than blocks.  Reference is a
+// host triple loop on the replicated global matrix; tolerance 1e-11 on
+// max element error (the distributed and reference summation orders
+// differ, so bitwise equality is not expected AGAINST THE REFERENCE --
+// but IS expected between repeated distributed runs, which is T5).
+//
+//   T1 : C = A.B, full range, all layouts x all grids.
+//   T2 : C = beta C + alpha A.B, preloaded C, complex alpha/beta.
+//   T3 : windowed products, block-aligned sub-ranges incl. ragged N end.
+//   T4 : in-place windows: same matrix as A, B and C on disjoint windows.
+//   T5 : determinism: repeated product bitwise identical.
+//////////////////////////////////////////////////////////////////////////////
+
+#include <Grid/Grid.h>
+#include <Grid/algorithms/multigrid/BlockCyclicSumma.h>
+
+using namespace Grid;
+
+static int failures = 0;
+
+static void Report(const std::string &name, bool pass, const std::string &detail="")
+{
+  std::cout << GridLogMessage << "  " << name << (pass ? "  PASS" : "  ** FAIL **");
+  if ( detail.size() ) std::cout << "   " << detail;
+  std::cout << std::endl;
+  if ( !pass ) failures++;
+}
+
+// Deterministic pseudo-random fill: identical on every rank, no RNG state.
+static ComplexD Fill(int64_t i, int64_t j, int salt)
+{
+  double x = std::sin(0.7*i + 1.3*j + 0.31*salt) ;
+  double y = std::cos(1.9*i - 0.4*j + 0.77*salt) ;
+  return ComplexD(x,y);
+}
+
+// Host reference:  C[i0:i1,j0:j1] = beta C + alpha A[i-,k-].B[k-,j-]
+static void RefGemm(ComplexD alpha, const std::vector<ComplexD> &A,
+                    const std::vector<ComplexD> &B,
+                    ComplexD beta, std::vector<ComplexD> &C, int64_t N,
+                    int64_t i0,int64_t i1,int64_t j0,int64_t j1,int64_t k0,int64_t k1)
+{
+  for(int64_t j=j0;j<j1;j++){
+    for(int64_t i=i0;i<i1;i++){
+      ComplexD acc(0.0,0.0);
+      for(int64_t k=k0;k<k1;k++) acc += A[i+k*N]*B[k+j*N];
+      C[i+j*N] = beta*C[i+j*N] + alpha*acc;
+    }
+  }
+}
+
+static double MaxDiff(const std::vector<ComplexD> &X, const std::vector<ComplexD> &Y)
+{
+  double m = 0.0;
+  for(uint64_t i=0;i<X.size();i++) m = std::max(m, std::abs(X[i]-Y[i]));
+  return m;
+}
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+                                                       GridDefaultSimd(Nd, vComplexD::Nsimd()),
+                                                       GridDefaultMpi());
+  const int P = grid->ProcessorCount();
+
+  // All factorisations Pr*Pc == P: squarest, plus both degenerate rings.
+  std::vector<std::pair<int,int>> grids;
+  for(int r=1;r<=P;r++) if ( P%r==0 ) grids.push_back({r,P/r});
+
+  struct Cfg { int64_t N; int64_t nb; };
+  std::vector<Cfg> cfgs = { {24,4}, {26,4}, {17,5}, {8,8}, {30,7}, {3,4}, {33,3} };
+
+  const double tol = 1.0e-11;
+  BlockCyclicSumma SUMMA;
+
+  std::cout << GridLogMessage << "BlockCyclicSumma regression: P=" << P
+            << ", " << grids.size() << " process grids, " << cfgs.size()
+            << " layouts" << std::endl;
+
+  ////////////////////////////////////////////////////////////////////////
+  // T1 : full product, alpha=1 beta=0.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true; double worst = 0.0;
+    for(auto &g : grids){
+      for(auto &c : cfgs){
+        int64_t N = c.N;
+        std::vector<ComplexD> Ag(N*N), Bg(N*N), Cg(N*N, ComplexD(0.0,0.0)), Cd;
+        for(int64_t j=0;j<N;j++) for(int64_t i=0;i<N;i++){
+          Ag[i+j*N]=Fill(i,j,1); Bg[i+j*N]=Fill(i,j,2);
+        }
+        BlockCyclicMatrix A(grid,N,c.nb,g.first,g.second);
+        BlockCyclicMatrix B(grid,N,c.nb,g.first,g.second);
+        BlockCyclicMatrix C(grid,N,c.nb,g.first,g.second);
+        A.ImportGlobal(Ag); B.ImportGlobal(Bg);
+        SUMMA.Multiply(ComplexD(1.0,0.0),A,B,ComplexD(0.0,0.0),C, 0,N, 0,N, 0,N);
+        C.ExportGlobal(Cd);
+        RefGemm(ComplexD(1.0,0.0),Ag,Bg,ComplexD(0.0,0.0),Cg,N, 0,N,0,N,0,N);
+        double d = MaxDiff(Cd,Cg); worst = std::max(worst,d);
+        if ( d > tol ) ok = false;
+      }
+    }
+    Report("T1  full C = A.B, all grids x layouts", ok,
+           "max err "+std::to_string(worst));
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T2 : beta C + alpha A.B with preloaded C, complex coefficients.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true; double worst = 0.0;
+    ComplexD alpha(0.5,-0.25), beta(-1.0,0.75);
+    for(auto &g : grids){
+      for(auto &c : cfgs){
+        int64_t N = c.N;
+        std::vector<ComplexD> Ag(N*N), Bg(N*N), Cg(N*N), Cd;
+        for(int64_t j=0;j<N;j++) for(int64_t i=0;i<N;i++){
+          Ag[i+j*N]=Fill(i,j,3); Bg[i+j*N]=Fill(i,j,4); Cg[i+j*N]=Fill(i,j,5);
+        }
+        BlockCyclicMatrix A(grid,N,c.nb,g.first,g.second);
+        BlockCyclicMatrix B(grid,N,c.nb,g.first,g.second);
+        BlockCyclicMatrix C(grid,N,c.nb,g.first,g.second);
+        A.ImportGlobal(Ag); B.ImportGlobal(Bg); C.ImportGlobal(Cg);
+        SUMMA.Multiply(alpha,A,B,beta,C, 0,N, 0,N, 0,N);
+        C.ExportGlobal(Cd);
+        RefGemm(alpha,Ag,Bg,beta,Cg,N, 0,N,0,N,0,N);
+        double d = MaxDiff(Cd,Cg); worst = std::max(worst,d);
+        if ( d > tol ) ok = false;
+      }
+    }
+    Report("T2  C = beta C + alpha A.B, preloaded C", ok,
+           "max err "+std::to_string(worst));
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T3 : windowed products.  Block-aligned sub-ranges, including the
+  // ragged top end g1==N, on a layout with a partial trailing block.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true; double worst = 0.0;
+    int64_t N = 26, nb = 4;                    // 6 full blocks + ragged 2
+    struct Rng { int64_t i0,i1,j0,j1,k0,k1; };
+    std::vector<Rng> rngs = {
+      { 0,8,   8,16,  16,24 },                 // interior windows
+      { 4,12,  0,4,   12,26 },                 // ragged k end
+      { 16,26, 20,26, 0,8   },                 // ragged i and j ends
+      { 0,4,   0,4,   4,8   },                 // minimal one-block windows
+      { 0,26,  0,26,  8,12  },                 // full ij, thin k
+    };
+    for(auto &g : grids){
+      std::vector<ComplexD> Ag(N*N), Bg(N*N), Cg(N*N), Cd;
+      for(int64_t j=0;j<N;j++) for(int64_t i=0;i<N;i++){
+        Ag[i+j*N]=Fill(i,j,6); Bg[i+j*N]=Fill(i,j,7);
+      }
+      for(auto &r : rngs){
+        for(int64_t j=0;j<N;j++) for(int64_t i=0;i<N;i++) Cg[i+j*N]=Fill(i,j,8);
+        BlockCyclicMatrix A(grid,N,nb,g.first,g.second);
+        BlockCyclicMatrix B(grid,N,nb,g.first,g.second);
+        BlockCyclicMatrix C(grid,N,nb,g.first,g.second);
+        A.ImportGlobal(Ag); B.ImportGlobal(Bg); C.ImportGlobal(Cg);
+        SUMMA.Multiply(ComplexD(1.0,0.0),A,B,ComplexD(1.0,0.0),C,
+                       r.i0,r.i1, r.j0,r.j1, r.k0,r.k1);
+        C.ExportGlobal(Cd);
+        RefGemm(ComplexD(1.0,0.0),Ag,Bg,ComplexD(1.0,0.0),Cg,N,
+                r.i0,r.i1, r.j0,r.j1, r.k0,r.k1);
+        double d = MaxDiff(Cd,Cg); worst = std::max(worst,d);
+        if ( d > tol ) ok = false;
+      }
+    }
+    Report("T3  windowed products, ragged ends", ok,
+           "max err "+std::to_string(worst));
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T4 : in-place windows of ONE matrix, exactly the stage-3 usage:
+  //   M[0:b, 2b:3b] = M[0:b, b:2b] . M[b:2b, 2b:3b]
+  // C window disjoint from both operand windows (asserted in Multiply).
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true; double worst = 0.0;
+    int64_t nb = 4, N = 4*nb;
+    for(auto &g : grids){
+      std::vector<ComplexD> Mg(N*N), Md, Mr;
+      for(int64_t j=0;j<N;j++) for(int64_t i=0;i<N;i++) Mg[i+j*N]=Fill(i,j,9);
+      BlockCyclicMatrix M(grid,N,nb,g.first,g.second);
+      M.ImportGlobal(Mg);
+      SUMMA.Multiply(ComplexD(1.0,0.0),M,M,ComplexD(0.0,0.0),M,
+                     0,nb,  2*nb,3*nb,  nb,2*nb);
+      M.ExportGlobal(Md);
+      Mr = Mg;
+      RefGemm(ComplexD(1.0,0.0),Mg,Mg,ComplexD(0.0,0.0),Mr,N,
+              0,nb, 2*nb,3*nb, nb,2*nb);
+      double d = MaxDiff(Md,Mr); worst = std::max(worst,d);
+      if ( d > tol ) ok = false;
+    }
+    Report("T4  in-place disjoint windows (stage-3 usage)", ok,
+           "max err "+std::to_string(worst));
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // T5 : determinism.  The summation order is fixed (ascending k-block),
+  // so repeated distributed products must agree BITWISE -- the property
+  // the P2P design buys and a collective reduce cannot promise.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    bool ok = true;
+    for(auto &g : grids){
+      int64_t N = 30, nb = 7;
+      std::vector<ComplexD> Ag(N*N), Bg(N*N), C1, C2;
+      for(int64_t j=0;j<N;j++) for(int64_t i=0;i<N;i++){
+        Ag[i+j*N]=Fill(i,j,10); Bg[i+j*N]=Fill(i,j,11);
+      }
+      BlockCyclicMatrix A(grid,N,nb,g.first,g.second);
+      BlockCyclicMatrix B(grid,N,nb,g.first,g.second);
+      BlockCyclicMatrix C(grid,N,nb,g.first,g.second);
+      A.ImportGlobal(Ag); B.ImportGlobal(Bg);
+      SUMMA.Multiply(ComplexD(1.0,0.0),A,B,ComplexD(0.0,0.0),C, 0,N,0,N,0,N);
+      C.ExportGlobal(C1);
+      SUMMA.Multiply(ComplexD(1.0,0.0),A,B,ComplexD(0.0,0.0),C, 0,N,0,N,0,N);
+      C.ExportGlobal(C2);
+      for(uint64_t i=0;i<C1.size();i++)
+        if ( !(C1[i]==C2[i]) ) ok = false;     // BITWISE, not toleranced
+    }
+    Report("T5  repeated product bitwise identical", ok);
+  }
+
+  {
+    uint64_t f = failures;  grid->GlobalSum(f);
+    if ( f && !failures )
+      std::cout << GridLogMessage << "  ** failures on OTHER ranks: " << f << " **" << std::endl;
+    failures = (int)f;
+  }
+  std::cout << GridLogMessage << (failures ? "Test_summa: FAILURES"
+                                           : "Test_summa: ALL PASS") << std::endl;
+  Grid_finalize();
+  return failures ? 1 : 0;
+}