Better reduction

Grid/algorithms/iterative/GCRCoefficients.h

@@ -28,29 +28,57 @@ NAMESPACE_BEGIN(Grid);
 //////////////////////////////////////////////////////////////////////////////
 struct GCRCoefficients {
   int mmax = 0;
-  std::vector<ComplexD>               a_sum;   // [k]
-  std::vector<int>                    a_n;
-  std::vector<std::vector<ComplexD> > b_sum;   // [k][j]
-  std::vector<std::vector<int> >      b_n;
+  // Every recorded call is kept: calls[c] = list of (a_k, [b_kj]) per step.
+  // A(k)/B(k,j) return the coefficients of the SELECTED call: by default the
+  // last complete one.  Selection "mean" averages coefficients over calls --
+  // kept for comparison only: the mean of the coefficients of a nonlinear
+  // recurrence is not the mean of the polynomials, and in practice (Frontier
+  // M3, 2026-08-26) it was worse than every individual call.
+  enum Select { Last=0, First=1, Index=2, Mean=3 };
+  Select select = Last;
+  int    index  = 0;
+  typedef std::vector<std::pair<ComplexD,std::vector<ComplexD> > > Call;
+  std::vector<Call> calls;
+  Call current;
   void RecordA(int k, ComplexD a){
-    if ( (int)a_sum.size() <= k ) { a_sum.resize(k+1,ComplexD(0.0)); a_n.resize(k+1,0); }
-    a_sum[k] += a; a_n[k]++;
+    if ( k==0 && current.size() ) { calls.push_back(current); current.clear(); }
+    if ( (int)current.size() <= k ) current.resize(k+1);
+    current[k].first = a;
   }
   void RecordB(int k, const std::vector<ComplexD> &b){
-    if ( (int)b_sum.size() <= k ) { b_sum.resize(k+1); b_n.resize(k+1); }
-    if ( b_sum[k].size() < b.size() ) { b_sum[k].resize(b.size(),ComplexD(0.0)); b_n[k].resize(b.size(),0); }
-    for(int j=0;j<(int)b.size();j++){ b_sum[k][j] += b[j]; b_n[k][j]++; }
+    if ( (int)current.size() <= k ) current.resize(k+1);
+    current[k].second = b;
   }
-  int  Steps(void) const { return a_sum.size(); }
-  int  Calls(void) const { return a_n.size() ? a_n[0] : 0; }
-  ComplexD A(int k) const { return a_sum[k]/(double)a_n[k]; }
-  int  NB(int k) const { return (k<(int)b_sum.size()) ? b_sum[k].size() : 0; }
-  ComplexD B(int k,int j) const { return b_sum[k][j]/(double)b_n[k][j]; }
+  void Flush(void){ if ( current.size() ) { calls.push_back(current); current.clear(); } }
+  int  Calls(void) const { return calls.size() + (current.size() ? 1 : 0); }
+  const Call & Chosen(void) const {
+    GRID_ASSERT( calls.size() || current.size() );
+    if ( calls.empty() ) return current;
+    if ( select==First ) return calls.front();
+    if ( select==Index ) { GRID_ASSERT(index>=0 && index<(int)calls.size()); return calls[index]; }
+    return calls.back();
+  }
+  int  Steps(void) const { return select==Mean ? MeanSteps() : Chosen().size(); }
+  int  NB(int k) const   { return select==Mean ? MeanNB(k)   : Chosen()[k].second.size(); }
+  ComplexD A(int k) const    { return select==Mean ? MeanA(k)   : Chosen()[k].first; }
+  ComplexD B(int k,int j) const { return select==Mean ? MeanB(k,j) : Chosen()[k].second[j]; }
+  // mean over calls (comparison only)
+  int MeanSteps(void) const { int m=0; for(auto &c:calls) m = std::max(m,(int)c.size()); return m; }
+  int MeanNB(int k) const { int m=0; for(auto &c:calls) if(k<(int)c.size()) m = std::max(m,(int)c[k].second.size()); return m; }
+  ComplexD MeanA(int k) const { ComplexD s(0.0); int n=0; for(auto &c:calls) if(k<(int)c.size()){ s+=c[k].first; n++; } return s/(double)n; }
+  ComplexD MeanB(int k,int j) const { ComplexD s(0.0); int n=0; for(auto &c:calls) if(k<(int)c.size() && j<(int)c[k].second.size()){ s+=c[k].second[j]; n++; } return s/(double)n; }
   void Report(const std::string &name) const {
-    std::cout << GridLogMessage << "GCRCoefficients " << name << ": " << Calls() << " calls, " << Steps() << " steps, mmax " << mmax << std::endl;
+    const char *sel[4]={"last","first","index","mean"};
+    std::cout << GridLogMessage << "GCRCoefficients " << name << ": " << Calls() << " calls, " << Steps() << " steps, mmax " << mmax
+              << ", selection " << sel[select] << std::endl;
     for(int k=0;k<Steps();k++){
       std::cout << GridLogMessage << "  step " << k << " a=(" << real(A(k)) << "," << imag(A(k)) << ")";
       for(int j=0;j<NB(k);j++) std::cout << " b[" << j << "]=(" << real(B(k,j)) << "," << imag(B(k,j)) << ")";
+      // spread of a_k across calls: how different the individual polynomials are
+      if ( calls.size()>1 ) {
+        RealD lo=1e300, hi=0; for(auto &c:calls) if(k<(int)c.size()){ RealD x=real(c[k].first), y=imag(c[k].first); RealD m=std::sqrt(x*x+y*y); lo=std::min(lo,m); hi=std::max(hi,m); }
+        std::cout << "   |a| over calls [" << lo << "," << hi << "]";
+      }
       std::cout << std::endl;
     }
   }

Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h

@@ -85,6 +85,9 @@ public:
   // the iteration.
   GCRCoefficients *Recorder = nullptr;
   void SetCoefficientRecorder(GCRCoefficients *r) { Recorder = r; if(r) r->mmax = mmax; };
+  // Free the persistent history (e.g. when this solver is replaced by a
+  // replayed polynomial): re-made on the next call if ever needed again.
+  void ReleaseHistory(void) { q.clear(); p.clear(); qq.clear(); hist_grid = nullptr; };
 
   PrecGeneralisedConjugateResidualNonHermitian(RealD tol,Integer maxit,LinearOperatorBase<Field> &_Linop,LinearFunction<Field> &Prec,int _mmax,int _nstep) : 
     Tolerance(tol), 

Grid/algorithms/multigrid/BlockCyclic.h

@@ -62,13 +62,19 @@ public:
   // 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.
   ///////////////////////////////////////////////////////////////////////////
+  // Nearest-to-square factorisation with Pr >= Pc.  The orientation matters
+  // for SUMMA: B panels travel along process columns in slots of
+  // S = ceil(Pc/Pr) panels, padded; with Pr >= Pc, S == 1 and there is no
+  // padding.  (288 as 16x18 had S=2 -- half the B-ring bytes were zeros;
+  // 18x16 moves ~30% fewer bytes for the same inverse.)
   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;
+    int r = 1;
+    for(int f=1; (int64_t)f*f <= (int64_t)P; f++)
+      if ( P % f == 0 ) r = f;
+    Pc = r;
+    Pr = P / r;
   }
 
   ///////////////////////////////////////////////////////////////////////////

Grid/algorithms/multigrid/BlockCyclicSchurInverse.h

@@ -256,11 +256,103 @@ public:
   // PUBLIC ENTRY.  In-place inverse of the whole matrix.  Scratch (4x the
   // matrix footprint) is allocated here and released on return.
   ///////////////////////////////////////////////////////////////////////////
+  ///////////////////////////////////////////////////////////////////////////
+  // POWER - CLOCK - GROUND.  Before the inverse runs, print the preconditions
+  // that have differed between fast (62 s) and slow (133-141 s) runs of the
+  // SAME inverse, and measure ONE SendToRecvFrom to the ACTUAL ring partners
+  // at three sizes, device and host buffers.  Same code in the harness and in
+  // production, so the two processes are compared on the identical primitive
+  // before any explanation of the SUMMA rings is attempted.
+  // SCHUR2D_PROBE=0 disables (costs ~0.1-0.5 s).
+  ///////////////////////////////////////////////////////////////////////////
+  void Probe(BlockCyclicMatrix &A)
+  {
+    BlockCyclicLayout &L = A.layout;
+    GridBase *grid = A.grid;
+    int me = grid->ThisRank();
+    // --- banner ---
+    int thr = -1;
+#ifdef GRID_COMMS_MPI3
+    MPI_Query_thread(&thr);
+#endif
+    const char *omp = getenv("OMP_NUM_THREADS");
+    MemoryStatus ms = MemoryManager::GetFootprint();
+    std::cout << GridLogMessage << "Schur2D PROBE banner: MPI thread level " << thr
+              << " (0 single,1 funneled,2 serialized,3 multiple)  OMP_NUM_THREADS=" << (omp?omp:"unset")
+              << "  MemoryManager device bytes " << ms.DeviceBytes/1.0e9 << " GB (LRU " << ms.DeviceLRUBytes/1.0e9
+              << " GB, cap " << ms.DeviceMaxBytes/1.0e9 << " GB)"
+              << "  grid " << L.Pr << "x" << L.Pc << " nb " << L.nb << std::endl;
+#ifdef GRID_HIP
+    if ( me==0 ) acceleratorMem();
+#endif
+    // --- ring partners exactly as SUMMA uses them ---
+    int prow=L.prow, pcol=L.pcol, Pr=L.Pr, Pc=L.Pc;
+    struct Ring { const char *name; int dest, src; };
+    Ring rings[2] = { {"ringA(row, q+-1)", prow*Pc + (pcol+1)%Pc,       prow*Pc + (pcol-1+Pc)%Pc},
+                      {"ringB(col, p+-1)", ((prow+1)%Pr)*Pc + pcol,     ((prow-1+Pr)%Pr)*Pc + pcol} };
+    uint64_t sizes[3] = { 64ull*1024, 1024ull*1024, 8ull*1024*1024 };
+    uint64_t maxb = sizes[2];
+    deviceVector<char> dsend(maxb), drecv(maxb);
+    std::vector<char>  hsend(maxb), hrecv(maxb);
+    for(int r=0;r<2;r++){
+      if ( (r==0 && Pc==1) || (r==1 && Pr==1) ) continue;
+      int off = grid->IsOffNode(rings[r].dest);
+      for(int si=0;si<3;si++){
+        uint64_t bytes = sizes[si];
+        // warm one, time five, both memory spaces
+        grid->SendToRecvFrom(&dsend[0], rings[r].dest, &drecv[0], rings[r].src, bytes);
+        double t0=usecond();
+        for(int i=0;i<5;i++) grid->SendToRecvFrom(&dsend[0], rings[r].dest, &drecv[0], rings[r].src, bytes);
+        double td=(usecond()-t0)/5.0;
+        grid->SendToRecvFrom(&hsend[0], rings[r].dest, &hrecv[0], rings[r].src, bytes);
+        t0=usecond();
+        for(int i=0;i<5;i++) grid->SendToRecvFrom(&hsend[0], rings[r].dest, &hrecv[0], rings[r].src, bytes);
+        double th=(usecond()-t0)/5.0;
+        // spread over ranks
+        RealD dmax=td, dmin=-td; grid->GlobalMax(dmax); grid->GlobalMax(dmin); dmin=-dmin;
+        std::cout << GridLogMessage << "Schur2D PROBE " << rings[r].name << (off?" OFF-node":" on-node")
+                  << " " << bytes/1024 << " KB: device " << td << " us (" << bytes/td/1.0e3 << " GB/s) [min/max over ranks " << dmin << "/" << dmax << " us]"
+                  << "  host " << th << " us (" << bytes/th/1.0e3 << " GB/s)" << std::endl;
+      }
+    }
+    ///////////////////////////////////////////////////////////////////////
+    // The SUMMA's conditions, one at a time, at 8 MB on ring B:
+    //  (a) LARGE persistent buffers (the rings use ~0.5 GB Abuf/Bbuf), sending
+    //      from offset 0 and from deep inside the region;
+    //  (b) a pack kernel + accelerator_barrier immediately before each
+    //      message, as the SUMMA does.
+    // Isolated 8 MB messages ran at 11-20 GB/s while the SUMMA averaged 2.1;
+    // whichever variant drops to ~2 GB/s names the condition.
+    ///////////////////////////////////////////////////////////////////////
+    {
+      int r = (Pr>1) ? 1 : 0;
+      uint64_t bytes = sizes[2];
+      uint64_t big = 512ull*1024*1024;
+      deviceVector<char> bsend(big), brecv(big);
+      for(int variant=0; variant<3; variant++){
+        uint64_t so = (variant==1) ? big-bytes : 0;               // deep offset in the large region
+        char *sp=&bsend[so], *rp=&brecv[so];
+        grid->SendToRecvFrom(sp, rings[r].dest, rp, rings[r].src, bytes);
+        double t0=usecond();
+        for(int i=0;i<5;i++){
+          if ( variant==2 ) { accelerator_for(k, bytes/8, 1, { ((uint64_t *)sp)[k] = (uint64_t)k; }); accelerator_barrier(); }
+          grid->SendToRecvFrom(sp, rings[r].dest, rp, rings[r].src, bytes);
+        }
+        double t=(usecond()-t0)/5.0;
+        RealD tmax=t, tmin=-t; grid->GlobalMax(tmax); grid->GlobalMax(tmin); tmin=-tmin;
+        const char *vn[3]={"512MB buffer, offset 0","512MB buffer, offset 504MB","pack kernel + barrier before each send"};
+        std::cout << GridLogMessage << "Schur2D PROBE " << rings[r].name << " 8192 KB device, " << vn[variant] << ": "
+                  << t << " us (" << bytes/t/1.0e3 << " GB/s) [min/max over ranks " << tmin << "/" << tmax << " us]" << std::endl;
+      }
+    }
+  }
+
   void Invert(BlockCyclicMatrix &A)
   {
     BlockCyclicLayout &L = A.layout;
     GRID_ASSERT( L.N >= 1 );
     int64_t nblocks = (L.N + L.nb - 1)/L.nb;
+    if ( !(getenv("SCHUR2D_PROBE") && atoi(getenv("SCHUR2D_PROBE"))==0) ) Probe(A);
 
     BlockCyclicMatrix Bt(A.grid, L.N, L.nb, L.Pr, L.Pc);
     BlockCyclicMatrix Ct(A.grid, L.N, L.nb, L.Pr, L.Pc);

Grid/algorithms/multigrid/BlockCyclicSumma.h

@@ -132,6 +132,15 @@ public:
   // row, ring B along the process column, then the local GEMMs.  The rings
   // are synchronous SendToRecvFrom, so tRing is time the GPU is idle unless
   // a future version overlaps them with the GEMMs.
+  // PERSISTENT ring buffers: allocated once (grow-only) and reused by every
+  // Multiply, so the device addresses handed to MPI never change.  Fresh
+  // per-call buffers rotated through the caching allocator's blocks, and a
+  // GPU-direct RDMA registration cache keyed on address then re-registers
+  // per message: measured 1.26 GB/s/rank on 13 MB ring messages (production,
+  // GRID_ALLOC_NCACHE_LARGE=64) against 62 s total for the same inverse when
+  // hipMalloc returned a stable address.
+  deviceVector<ComplexD> Abuf;
+  deviceVector<ComplexD> Bbuf;
   double   tAlloc=0, tPack=0, tRingA=0, tRingB=0, tGemm=0;
   uint64_t bytesRing=0, nRingMsg=0, nMultiply=0, nGemm=0;
   void ResetTelemetry(void){ tAlloc=tPack=tRingA=tRingB=tGemm=0; bytesRing=nRingMsg=nMultiply=nGemm=0; }
@@ -199,8 +208,8 @@ public:
     const uint64_t slotB  = (uint64_t)S*slotB1;
     nMultiply++;
     tAlloc -= usecond();
-    deviceVector<ComplexD> Abuf( slotA*Pc ? slotA*Pc : 1 );
-    deviceVector<ComplexD> Bbuf( slotB*Pr ? slotB*Pr : 1 );
+    if ( Abuf.size() < std::max<uint64_t>(slotA*Pc,1) ) Abuf.resize( std::max<uint64_t>(slotA*Pc,1) );
+    if ( Bbuf.size() < std::max<uint64_t>(slotB*Pr,1) ) Bbuf.resize( std::max<uint64_t>(slotB*Pr,1) );
     tAlloc += usecond();
 
     deviceVector<ComplexD *> ap(1), bp(1), cp(1);

Grid/algorithms/multigrid/DenseCoarseMatrix.h

@@ -119,6 +119,9 @@ public:
   deviceVector<ComplexF>  dSlab;
   deviceVector<ComplexF>  dX;       // N x MRHS_MAX
   deviceVector<ComplexF>  dY;       // nrows x MRHS_MAX
+  deviceVector<ComplexF>  dG;       // N x MRHS_MAX lex-major staging for the allgather (devSum==4)
+  deviceVector<int>       dLex2Rank;// lex index of a process coordinate -> its rank (allgather block order -> row-block order)
+  int                     myLex;
   deviceVector<ComplexF>  dPartial; // NK x (nrows x MRHS_MAX)
   deviceVector<ComplexF*> aptrs;    // slab K-chunk pointers   (lda = N)
   deviceVector<ComplexF*> xptrs;    // X    K-chunk pointers   (ldb = N)
@@ -252,9 +255,22 @@ public:
       acceleratorCopyToDevice(&h[0],&cptrs[0],NK*sizeof(ComplexF*));
 
       devSum = getenv("DENSE_DEVICE_SUM") ? atoi(getenv("DENSE_DEVICE_SUM")) : 0;
-      const char *sumName[4] = {"host allreduce","DEVICE-buffer allreduce (GPU-aware MPI)",
-                                "DEVICE cartesian ring allreduce (P2P)","DEVICE flat ring allreduce (P2P)"};
-      GRID_ASSERT(devSum>=0 && devSum<=3);
+      const char *sumName[5] = {"host allreduce","DEVICE-buffer allreduce (GPU-aware MPI)",
+                                "DEVICE cartesian ring allreduce (P2P)","DEVICE flat ring allreduce (P2P)",
+                                "DEVICE cartesian ring ALLGATHER (P2P, ~8x fewer bytes than the padded allreduce)"};
+      GRID_ASSERT(devSum>=0 && devSum<=4);
+      if ( devSum==4 ) {
+        dG.resize((uint64_t)N*MRHS_MAX);
+        // allgather delivers blocks in lexicographic-coordinate order; the row
+        // blocks of x are in RANK order.  Same table as BuildRankMajorMap.
+        int P = grid->ProcessorCount();
+        std::vector<int> l2r(P);
+        for(int lp=0; lp<P; lp++){ Coordinate pc(nd); Lexicographic::CoorFromIndex(pc, lp, grid->_processors); l2r[lp] = grid->RankFromProcessorCoor(pc); }
+        dLex2Rank.resize(P);
+        acceleratorCopyToDevice(&l2r[0], &dLex2Rank[0], P*sizeof(int));
+        myLex = CartesianLexIndex(grid);
+        GRID_ASSERT( l2r[myLex] == grid->ThisRank() );
+      }
       std::cout << GridLogMessage << "DenseCoarseMatrix: slab resident on device ("
                 << sbytes/1024./1024. << " MB/rank), split-K NK=" << NK << " (Kc=" << Kc << "); "
                 << sumName[devSum] << std::endl;
@@ -957,7 +973,34 @@ public:
     int64_t  Kc = N / NK;
     double t1 = usecond();
     double t2, t3;
-    if (devSum) {
+    if (devSum==4) {
+      // ALLGATHER: x is not a reduction -- every rank owns rows
+      // [me*nrows,(me+1)*nrows) of x and needs all of it.  Only MY rows go
+      // host->device (nrows x nr, ~15 KB at nr=1), rank-major staged, gathered
+      // along the process grid, then scattered into the column-major dX
+      // (ld = N) the split-K GEMM reads.
+      const int me = grid->ThisRank();
+      const uint64_t chunk = (uint64_t)nrows*nr;              // my block: [r][i]
+      { GRID_TRACE("DenseH2D");
+        std::vector<ComplexF> hG(chunk);
+        for(int r=0;r<nr;r++)
+          memcpy(&hG[(uint64_t)r*nrows], &hX[(uint64_t)r*N + (uint64_t)me*nrows], nrows*sizeof(ComplexF));
+        acceleratorCopyToDevice(&hG[0], &dG[(uint64_t)myLex*chunk], chunk*sizeof(ComplexF));   // my slot is my LEX index
+      }
+      t2 = usecond();
+      { GRID_TRACE("DenseAllgather");
+        CartesianRingAllGather(grid, (ComplexF *)&dG[0], chunk);
+        // scatter lex block L=[r][i] -> dX[r*N + rank(L)*nrows + i]
+        ComplexF *g = &dG[0]; ComplexF *x = &dX[0]; int *l2r = &dLex2Rank[0];
+        const int64_t nrw = nrows; const int64_t NN = N; const int nrr = nr;
+        accelerator_for(idx, (uint64_t)N*nr, 1, {
+          int64_t r = idx / NN;  int64_t gi = idx - r*NN;
+          int64_t L = gi / nrw;  int64_t i  = gi - L*nrw;
+          x[r*NN + (int64_t)l2r[L]*nrw + i] = g[L*(nrw*nrr) + r*nrw + i];
+        });
+      }
+      t3 = usecond();
+    } else if (devSum) {
       { GRID_TRACE("DenseH2D");
         acceleratorCopyToDevice(&hX[0],&dX[0],nX*sizeof(ComplexF));
       }
@@ -967,6 +1010,7 @@ public:
         //                      above ~8 MB: 12 RHS at N=138240 is 13.3 MB)
         // DENSE_DEVICE_SUM=2 : CartesianRingAllReduce, P2P only, no size cliff
         // DENSE_DEVICE_SUM=3 : flat RingAllReduce, P2P only
+        // DENSE_DEVICE_SUM=4 : CartesianRingAllGather (branch above)
         if      (devSum==2) CartesianRingAllReduce(grid,(ComplexF *)&dX[0],nX);
         else if (devSum==3) RingAllReduce(grid,(ComplexF *)&dX[0],nX);
         else                grid->GlobalSumVector((ComplexF *)&dX[0], (int)nX);

Grid/algorithms/multigrid/Smoothers.h

@@ -91,6 +91,8 @@ public:
   LinearOperatorBase<Field> &Linop;
   RealD lo, hi; int order;
   std::vector<RealD> Coeffs;
+  int Verbose = 0;
+  std::string name = "cheb";
   ChebyshevNonHermitianSmoother(RealD _lo,RealD _hi,int _order,LinearOperatorBase<Field> &Op)
     : Linop(Op), lo(_lo), hi(_hi), order(_order)
   {
@@ -126,6 +128,7 @@ public:
       if ( Coeffs[n] != 0.0 ) axpy(out,Coeffs[n],*Tnp,out);
       Field *swizzle=Tnm; Tnm=Tn; Tn=Tnp; Tnp=swizzle;
     }
+    if ( Verbose ) { Linop.Op(out,y); y = y - in; std::cout << GridLogMessage << " " << name << " cheb |r|/|r0| = " << std::sqrt(norm2(y)/norm2(in)) << std::endl; }
   }
 };
 
@@ -205,16 +208,32 @@ public:
   std::vector<std::vector<ComplexD> > b;
   GridBase *hist_grid = nullptr;
   std::vector<Field> p;
+  int Verbose = 0;                       // 1: print |r_m|/|r_0| per call (one extra reduction)
+  std::string name = "replay";
   GCRReplaySmoother(LinearOperatorBase<Field> &Op, const GCRCoefficients &c) : Linop(Op)
   {
     mmax  = c.mmax;  GRID_ASSERT(mmax>=1);
     nstep = c.Steps(); GRID_ASSERT(nstep>=1);
     a.resize(nstep); b.resize(nstep);
+    // No std::abs / std::isfinite on ComplexD: it is thrust::complex on HIP
+    // builds and neither overload exists.  Work on the real and imaginary
+    // parts explicitly.
+    auto finite = [](ComplexD z){ RealD x=real(z), y=imag(z); return (x==x) && (y==y) && (x-x==0.0) && (y-y==0.0); };
+    auto cabs   = [](ComplexD z){ RealD x=real(z), y=imag(z); return std::sqrt(x*x+y*y); };
+    RealD amax=0.0, bmax=0.0;
     for(int k=0;k<nstep;k++){
       a[k] = c.A(k);
+      GRID_ASSERT( finite(a[k]) );
+      amax = std::max(amax, cabs(a[k]));
       b[k].resize(c.NB(k));
-      for(int j=0;j<c.NB(k);j++) b[k][j] = c.B(k,j);
+      GRID_ASSERT( c.NB(k) <= mmax-1 );
+      for(int j=0;j<c.NB(k);j++){
+        b[k][j] = c.B(k,j);
+        GRID_ASSERT( finite(b[k][j]) );
+        bmax = std::max(bmax, cabs(b[k][j]));
+      }
     }
+    std::cout << GridLogMessage << " GCRReplaySmoother: max|a| " << amax << " max|b| " << bmax << std::endl;
     std::cout << GridLogMessage << " GCRReplaySmoother: " << nstep << " steps, mmax " << mmax
               << ", from " << c.Calls() << " recorded calls" << std::endl;
   }
@@ -231,11 +250,15 @@ public:
     r = src;
     psi = Zero();
     p[0] = r;
+    RealD r0 = Verbose ? norm2(src) : 0.0;
     for(int k=0;k<nstep;k++){
       int kp=k+1, peri_k=k%mmax, peri_kp=kp%mmax;
       Linop.Op(p[peri_k],q);                 // q_k = A p_k
       axpy(psi,  a[k], p[peri_k], psi);
-      if ( k==nstep-1 ) break;
+      if ( k==nstep-1 ) {
+        if ( Verbose ) { axpy(r,-a[k],q,r); std::cout << GridLogMessage << " " << name << " replay |r|/|r0| = " << std::sqrt(norm2(r)/r0) << std::endl; }
+        break;
+      }
       axpy(r,   -a[k], q,         r);
       p[peri_kp] = r;
       for(int j=0;j<(int)b[k].size();j++){

Grid/communicator/RingAllReduce.h

@@ -116,4 +116,68 @@ void CartesianRingAllReduce(CartesianCommunicator *comm, T *buf, uint64_t n)
   }
 }
 
+/////////////////////////////////////////////////////////////////////////////
+// Cartesian ring ALLGATHER, point-to-point only.
+//
+//   CartesianRingAllGather(comm, buf, chunk)
+//     buf holds P*chunk elements of T.  Block index = the Grid LEXICOGRAPHIC
+//     index of the owning process coordinate (dimension 0 fastest,
+//     Lexicographic::CoorFromIndex convention), NOT the MPI rank: on entry my
+//     chunk is at buf[CartesianLexIndex(comm)*chunk]; on exit block L is the
+//     chunk of the process at coordinate CoorFromIndex(L).  Map to ranks with
+//     comm->RankFromProcessorCoor.  (Ranks and coordinates are NOT related
+//     lexicographically on Frontier -- the OptimalCommunicator relabels ranks
+//     for shared-memory locality; assuming rank order gave a wrong inverse,
+//     VERIFY 0.9965, 2026-08-26.)
+//
+// Dimension by dimension from dimension 0 (fastest) upward: each stage is a
+// ring over the P_d ranks of that line, after which the held block is the
+// concatenation over that coordinate in the lexicographic nesting.  Bytes
+// sent per rank ~N = P*chunk in total (dominated by the last stage): 8x less
+// than a zero-padded CartesianRingAllReduce.  Steps: sum_d (P_d-1).  Exact.
+//
+// Written for the dense coarse-coarse apply (every rank owns rows of A^{-1}
+// and needs the whole x): 1.86 ms with the allreduce ring at 288 ranks was
+// wire speed but moved 35 MB per rank to deliver 4.4.
+/////////////////////////////////////////////////////////////////////////////
+inline int CartesianLexIndex(CartesianCommunicator *comm)
+{
+  int idx=0, stride=1;
+  for(int d=0; d<(int)comm->_ndimension; d++){ idx += comm->_processor_coor[d]*stride; stride *= comm->_processors[d]; }
+  return idx;
+}
+
+template<class T>
+void CartesianRingAllGather(CartesianCommunicator *comm, T *buf, uint64_t chunk)
+{
+  int P  = comm->ProcessorCount();
+  if ( P==1 || chunk==0 ) return;
+  int Nd = comm->_ndimension;
+  int mylex = CartesianLexIndex(comm);
+  deviceVector<T> work((uint64_t)P*chunk);
+  // ping-pong between buf and work; the held block lives at offset `off` in `cur`
+  T *cur = buf;       uint64_t off = (uint64_t)mylex*chunk;
+  T *oth = &work[0];
+  uint64_t blk = chunk;                        // elements in the held block
+  for(int d=0; d<Nd; d++){                     // dimension 0 first: it is the fastest lex index
+    int Pd = comm->_processors[d];
+    if ( Pd==1 ) continue;
+    int med = comm->_processor_coor[d];
+    int next, prev;
+    comm->ShiftedRanks(d, 1, prev, next);      // (dim, shift, source, dest)
+    GRID_ASSERT( (blk*sizeof(T))%4 == 0 );
+    acceleratorCopyDeviceToDevice((void *)(cur+off), (void *)(oth+(uint64_t)med*blk), blk*sizeof(T));
+    for(int t=1;t<Pd;t++){
+      int sendslot = (med - t + 1 + Pd) % Pd;
+      int recvslot = (med - t     + Pd) % Pd;
+      comm->SendToRecvFrom((void *)(oth+(uint64_t)sendslot*blk), next,
+                           (void *)(oth+(uint64_t)recvslot*blk), prev, blk*sizeof(T));
+    }
+    T *tmp = cur; cur = oth; oth = tmp; off = 0;
+    blk *= Pd;
+  }
+  GRID_ASSERT( blk == (uint64_t)P*chunk );
+  if ( cur != buf ) acceleratorCopyDeviceToDevice((void *)cur, (void *)buf, blk*sizeof(T));
+}
+
 NAMESPACE_END(Grid);

Grid/perfmon/Timer.h

@@ -78,7 +78,12 @@ inline std::ostream& operator<< (std::ostream & stream, const GridUsecs & now)
   auto     seconds    = now/second ; 
   auto     subseconds = now%second ;
   auto     fill       = stream.fill();
-  stream << seconds<<"."<<std::setw(6)<<std::setfill('0')<<subseconds.count()<<" s";
+  auto     flags      = stream.flags();
+  // The logger leaves the stream in std::left (channel-name field); under
+  // std::left setw(6)/setfill('0') pads on the RIGHT, printing 403.059986 as
+  // 403.599860 -- timestamps with <6 sub-second digits jumped by up to 10x.
+  stream << seconds<<"."<<std::right<<std::setw(6)<<std::setfill('0')<<subseconds.count()<<" s";
+  stream.flags(flags);
   stream.fill(fill);
   return stream;
 }

examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc

@@ -103,6 +103,10 @@ int   PowerIterations      = 0;     // >0: power-iterate the smoother operators
 std::string FineSmootherMode   = "gcr";
 std::string CoarseSmootherMode = "gcr";
 int   PolyRecordIters      = 4;
+int   PolyRecordStart      = 0;     // outer step at which recording begins (0: from the first step)
+std::string PolyRecordSelect = "last"; // which recorded call to replay: last|first|mean (mean is the bad one)
+int   PolyRefresh          = 0;     // >0: every PolyRefresh outer steps, one adaptive step re-records the polynomial (HDCG: every 10)
+int   PolyVerbose          = 0;     // 1: fixed-polynomial smoothers print |r_m|/|r_0| per call
 RealD FineChebLo   = 3.0,  FineChebHi   = 137.0;   // from the harvested polynomial and PowerIterations edge
 RealD CoarseChebLo = 8.0,  CoarseChebHi = 45.0;
 int   CoarseSmootherNstep  = 2;
@@ -152,6 +156,10 @@ void ParseEnvironment(void)
   if(getenv("FineSmootherMode"))   FineSmootherMode   = getenv("FineSmootherMode");
   if(getenv("CoarseSmootherMode")) CoarseSmootherMode = getenv("CoarseSmootherMode");
   if(getenv("PolyRecordIters"))    PolyRecordIters    = atoi(getenv("PolyRecordIters"));
+  if(getenv("PolyRecordStart"))    PolyRecordStart    = atoi(getenv("PolyRecordStart"));
+  if(getenv("PolyRecordSelect"))   PolyRecordSelect   = getenv("PolyRecordSelect");
+  if(getenv("PolyRefresh"))        PolyRefresh        = atoi(getenv("PolyRefresh"));
+  if(getenv("PolyVerbose"))        PolyVerbose        = atoi(getenv("PolyVerbose"));
   if(getenv("FineChebLo"))         FineChebLo         = atof(getenv("FineChebLo"));
   if(getenv("FineChebHi"))         FineChebHi         = atof(getenv("FineChebHi"));
   if(getenv("CoarseChebLo"))       CoarseChebLo       = atof(getenv("CoarseChebLo"));
@@ -548,6 +556,17 @@ public:
 
 int main (int argc, char ** argv)
 {
+  // GRID_MPI_THREAD_MULTIPLE=1: initialise MPI at MPI_THREAD_MULTIPLE before
+  // Grid_init (Grid asks for SERIALIZED).  The SLATE harness does this and its
+  // copy of the 2D inverse ran at ~2x the production ring rate (62 s vs
+  // 133-141 s).  Second hypothesis behind OMP_NUM_THREADS; test one at a time.
+  // Pair with MPICH_MAX_THREAD_SAFETY=multiple.
+  if ( getenv("GRID_MPI_THREAD_MULTIPLE") ) {
+    int provided = 0;
+    MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &provided);
+    std::cout << "GRID_MPI_THREAD_MULTIPLE: requested MPI_THREAD_MULTIPLE, provided " << provided
+              << (provided==MPI_THREAD_MULTIPLE ? " (MULTIPLE)" : " (NOT multiple)") << std::endl;
+  }
   Grid_init(&argc,&argv);
   ParseEnvironment();
 
@@ -989,33 +1008,74 @@ int main (int argc, char ** argv)
     std::unique_ptr<GCRReplaySmoother<LatticeFermionD> > FineReplay;
     std::unique_ptr<GCRReplaySmoother<CoarseVector> >    CoarseReplay;
     GCRCoefficients recF, recC;
+    {
+      GCRCoefficients::Select sel = GCRCoefficients::Last;
+      if ( PolyRecordSelect=="first" ) sel = GCRCoefficients::First;
+      if ( PolyRecordSelect=="mean" )  sel = GCRCoefficients::Mean;
+      recF.select = sel; recC.select = sel;
+    }
     if ( FineSmootherMode == "cheb" ) {
       FineCheb.reset(new ChebyshevNonHermitianSmoother<LatticeFermionD>(FineChebLo,FineChebHi,FineSmootherOrder,ShiftedPVdagM));
+      FineCheb->Verbose = PolyVerbose; FineCheb->name = "Fsmoother";
       FineSmootherSlot.Set(*FineCheb,"Fsmoother Chebyshev");
     }
     if ( CoarseSmootherMode == "cheb" ) {
       CoarseCheb.reset(new ChebyshevNonHermitianSmoother<CoarseVector>(CoarseChebLo,CoarseChebHi,CoarseSmootherNstep,ShiftedC));
+      CoarseCheb->Verbose = PolyVerbose; CoarseCheb->name = "Csmoother";
       CoarseSmootherSlot.Set(*CoarseCheb,"Csmoother Chebyshev");
     }
-    if ( FineSmootherMode   == "replay" ) SmootherGCR.SetCoefficientRecorder(&recF);
-    if ( CoarseSmootherMode == "replay" ) CoarseSmootherGCR.SetCoefficientRecorder(&recC);
-    L1PGCR.OnStep = [&](int step){
-      if ( step != PolyRecordIters ) return;
+    // Recording window [PolyRecordStart, PolyRecordStart+PolyRecordIters).
+    // M3 (2026-08-26): the GCR polynomial changes fast over the first outer
+    // steps (per-call |r|/|r0| 0.0034 -> 0.017 over steps 1-4) and the MEAN
+    // of those is a poor smoother (replay 0.033-0.049 per call); record a
+    // settled window instead.
+    if ( PolyRecordStart == 0 ) {
+      if ( FineSmootherMode   == "replay" ) SmootherGCR.SetCoefficientRecorder(&recF);
+      if ( CoarseSmootherMode == "replay" ) CoarseSmootherGCR.SetCoefficientRecorder(&recC);
+    }
+    // Record -> replay, with optional periodic re-recording ("re-record, not
+    // fade away": HDCG refreshed its polynomial every 10 steps, tracking the
+    // evolving spectral content of the residual).  Schedule on outer steps:
+    //   [PolyRecordStart, +PolyRecordIters)          adaptive GCR, recording
+    //   then replay of the selected recorded call;
+    //   if PolyRefresh>0: every PolyRefresh steps, ONE adaptive recording
+    //   step, then replay of that call.
+    auto BuildReplays = [&](void){
       if ( FineSmootherMode == "replay" ) {
         SmootherGCR.SetCoefficientRecorder(nullptr);
-        recF.Report("Fsmoother");
+        recF.Flush(); recF.Report("Fsmoother");
         FineReplay.reset(new GCRReplaySmoother<LatticeFermionD>(ShiftedPVdagM,recF));
+        FineReplay->Verbose = PolyVerbose; FineReplay->name = "Fsmoother";
         FineSmootherSlot.Set(*FineReplay,"Fsmoother replay");
+        SmootherGCR.ReleaseHistory();
       }
       if ( CoarseSmootherMode == "replay" ) {
         CoarseSmootherGCR.SetCoefficientRecorder(nullptr);
-        recC.Report("Csmoother");
+        recC.Flush(); recC.Report("Csmoother");
         CoarseReplay.reset(new GCRReplaySmoother<CoarseVector>(ShiftedC,recC));
+        CoarseReplay->Verbose = PolyVerbose; CoarseReplay->name = "Csmoother";
         CoarseSmootherSlot.Set(*CoarseReplay,"Csmoother replay");
+        CoarseSmootherGCR.ReleaseHistory();
+      }
+    };
+    auto StartRecording = [&](int step){
+      if ( FineSmootherMode   == "replay" ) { { auto sel=recF.select; recF = GCRCoefficients(); recF.select=sel; } SmootherGCR.SetCoefficientRecorder(&recF);       FineSmootherSlot.Set(SmootherGCR,"Fsmoother GCR (recording)"); }
+      if ( CoarseSmootherMode == "replay" ) { { auto sel=recC.select; recC = GCRCoefficients(); recC.select=sel; } CoarseSmootherGCR.SetCoefficientRecorder(&recC); CoarseSmootherSlot.Set(CoarseSmootherGCR,"Csmoother GCR (recording)"); }
+      std::cout << GridLogMessage << "Smoother coefficient recording starts at outer step " << step << std::endl;
+    };
+    int switchStep = PolyRecordStart + PolyRecordIters;
+    L1PGCR.OnStep = [&](int step){
+      if ( FineSmootherMode != "replay" && CoarseSmootherMode != "replay" ) return;
+      if ( step == PolyRecordStart && PolyRecordStart > 0 ) StartRecording(step);
+      if ( step == switchStep ) { BuildReplays(); return; }
+      if ( PolyRefresh > 0 && step > switchStep ) {
+        int since = step - switchStep;
+        if ( since % PolyRefresh == 0 )     { StartRecording(step); }       // one adaptive, recorded step
+        if ( since % PolyRefresh == 1 )     { BuildReplays(); }             // then replay it
       }
     };
     std::cout << GridLogMessage << "Smoother modes: fine " << FineSmootherMode << "  coarse " << CoarseSmootherMode
-              << (FineSmootherMode=="replay"||CoarseSmootherMode=="replay" ? "  (record for "+std::to_string(PolyRecordIters)+" outer steps)" : "")
+              << (FineSmootherMode=="replay"||CoarseSmootherMode=="replay" ? "  (record outer steps "+std::to_string(PolyRecordStart)+".."+std::to_string(PolyRecordStart+PolyRecordIters)+")" : "")
               << std::endl;
 
     std::vector<LatticeFermionD> src(nr,FGrid), sol(nr,FGrid);

systems/Frontier/schur2d_ladder.job

@@ -80,7 +80,7 @@ S2D_N=13824 srun -N1 -n8 ./select_gpu $root/tests/debug/Test_schur2d_scale \
 
 ##############################################################################
 echo "========================================================="
-echo "F3: scale rehearsal, 36 nodes, N=138240 (nb=480, grid 16x18)"
+echo "F3: scale rehearsal, 36 nodes, N=138240 (nb=480, grid 18x16)"
 echo "     THE number: invert phase vs the 1D baseline 328 s"
 echo "========================================================="
 ##############################################################################

systems/Frontier/slate_vs_schur2d.job

@@ -20,7 +20,7 @@
 #
 #   S1 : 1 node,  8 ranks, N=4096    -- API/device/thread-level shakeout
 #   S2 : 1 node,  8 ranks, N=13824   -- one-node production-shape rehearsal
-#   S3 : 36 nodes, 288 ranks, N=138240 (nb=480, grid 16x18) -- THE comparison
+#   S3 : 36 nodes, 288 ranks, N=138240 (nb=480, grid 18x16) -- THE comparison
 #
 # S3 is gated on S1 passing (rc=0 and both certificates printed): a broken
 # SLATE leg should not burn the 36-node allocation.
@@ -100,7 +100,7 @@ S2D_N=13824 srun -N1 -n8 ./select_gpu $BIN --mpi 1.1.2.4 --grid 16.16.16.16 $OPT
 
 ##############################################################################
 echo "========================================================="
-echo "S3: THE comparison, 36 nodes, 288 ranks, N=138240 (nb=480, grid 16x18)"
+echo "S3: THE comparison, 36 nodes, 288 ranks, N=138240 (nb=480, grid 18x16)"
 echo "========================================================="
 ##############################################################################
 if [ "$S1RC" -eq 0 ] && [ "$S1CERT" -ge 2 ]

systems/Frontier/smoother_modes.job

@@ -4,7 +4,7 @@
 #SBATCH --ntasks-per-node=8
 #SBATCH --cpus-per-task=7
 #SBATCH --gpus-per-node=8
-#SBATCH --time=1:00:00
+#SBATCH --time=1:45:00
 #SBATCH --account=phy157_dwf
 #SBATCH --gpu-bind=none
 #SBATCH --exclusive
@@ -23,6 +23,7 @@
 #                        heavy one at Nrhs=1)
 #   M4 cheb   / gcr      fine Chebyshev [FineChebLo,FineChebHi] order Fso
 #   M5 cheb   / cheb
+#   M6 gcr    / replay   coarse frozen only
 #
 # Laptop 8^4 findings (hot config, Ls=4, NBASIS=8): replay/replay converges
 # (28 vs 23 outer); cheb on the FINE level diverges there while cheb on the
@@ -88,38 +89,64 @@ export DENSE_CC=1
 export DENSE_APPLY_PROFILE=1
 unset  DENSE_CC_CHECK
 export DENSE_SPLITK=128
-export DENSE_DEVICE_SUM=2          # cartesian P2P ring: no 8 MB device-allreduce cliff (NRHS=12 abort)
+export DENSE_DEVICE_SUM=4          # cartesian P2P ring ALLGATHER: ~8x fewer bytes than the padded allreduce (=2); no collectives, no size cliff
 export GRID_ALLOC_NCACHE_LARGE=64
 export NRHS=4
 export PowerIterations=0
 export SmootherCoeffLog=0
 
 # frozen-polynomial controls
-export PolyRecordIters=4           # outer steps recorded before the switch
+export PolyRecordIters=8           # outer steps recorded
+export PolyRecordStart=8           # ...starting here: the early-step polynomials are unrepresentative (M3)
+export PolyRecordSelect=last       # replay ONE recorded call's polynomial (PB: every individual call beats the coefficient mean)
+export PolyRefresh=5               # re-record every 5 outer steps: BFM BfmHDCG.C:2243, k%5==1 -> LdopM1MirsPolyRecord, single call, replayed 4 steps
+# Inverse ring-rate hypotheses, ONE AT A TIME: (1) OMP_NUM_THREADS=1 (set above);
+# (2) if (1) fails, uncomment the two lines below (harness ran 62 s with these).
+#export MPICH_MAX_THREAD_SAFETY=multiple
+#export GRID_MPI_THREAD_MULTIPLE=1
+export PolyVerbose=1               # frozen smoothers print |r_m|/|r_0| per call: separates 'bad polynomial' from 'linear V-cycle stagnates the outer'
 export FineChebLo=3.0              # harvested |R|<0.1 edge / PowerIteration edge x1.05
 export FineChebHi=137.0
 export CoarseChebLo=8.0
 export CoarseChebHi=47.0           # shift 2.0: edge 43.3 x1.08
 
-run_mode () {
-  name=$1; export FineSmootherMode=$2; export CoarseSmootherMode=$3
-  echo "----- $name : FineSmootherMode=$FineSmootherMode CoarseSmootherMode=$CoarseSmootherMode -----"
-  fname=log.modes.$name
-  srun -N36 -n288 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \
+# Reference: the banked ADAPTIVE optimum, Fso6 / sm4 / Css2.0 / Nstep2 / svm8 ->
+# 28.57 s (Nrhs=1), ~14.9 s/RHS (Nrhs=4).  The stationary smoother converged at
+# the deliberate overshoot (order 12, fine shift 1.0, coarse Nstep 6); the
+# ladder below walks back towards the banked point.  A cell wins if it stays
+# convergent AND beats 28.57 s.  Each cell ~5 min.
+run_cell () {
+  name=$1; export FineSmootherOrder=$2; export FineSmootherShift=$3; export CoarseSmootherNstep=$4
+  export FineSmootherMode=$5; export CoarseSmootherMode=$6
+  echo "----- $name : Fso=$FineSmootherOrder Fss=$FineSmootherShift Csn=$CoarseSmootherNstep fine=$FineSmootherMode coarse=$CoarseSmootherMode -----"
+  fname=log.ladder.$name
+  srun -N36 -n288 --kill-on-bad-exit=1 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \
        --mpi ${MPI_GEOM} --grid $vol $OPTS1 --comms-overlap > $fname 2>&1
-  echo "  $(grep -h 'V2 3-level solve Nrhs' $fname | tr '\n' ' ')"
-  echo "  $(grep -h 'Fouter MrhsPGCR: Converged' $fname | sed 's/.*Converged/Converged/' | tr '\n' ' ')"
-  echo "  $(grep -h 'FINAL Nrhs .: worst' $fname | tr '\n' ' ')"
-  grep -h "SwitchableSmoother\|GCRCoefficients .*calls" $fname | head -6
+  echo "  exit $?"; sleep 60
+  echo "  $(grep -h 'V2 3-level solve Nrhs' $fname | sed 's/.*V2/V2/' | tr '\n' ' ')"
+  echo "  $(grep -h 'Fouter MrhsPGCR: Converged' $fname | sed 's/.*Converged/Converged/' | cut -c1-60 | tr '\n' ' ')"
+  echo "  replay per-call |r|/|r0| (Nrhs=1 solve): $(awk '/THREE-level solve, Nrhs = 1/{s=1} s && /Fsmoother replay \|r\|/{v=$NF; n++; t+=v; if(v>mx)mx=v} END{if(n) printf "mean %.4f max %.4f over %d calls", t/n, mx, n}' $fname)"
+  grep -h "SCHUR fp64 distributed invert took\|GB/s/rank" $fname | sed 's/^Grid : Message : [0-9.]* s : //' | cut -c1-120 | head -2
 }
 
-run_mode M1_gcr_gcr        gcr    gcr
-run_mode M2_replay_replay  replay replay
-run_mode M3_replay_gcr     replay gcr
-run_mode M4_cheb_gcr       cheb   gcr
-run_mode M5_cheb_cheb      cheb   cheb
+# Ladder 2 (after L0-L3 of the first ladder): the stationary polynomial
+# wants a LARGER shift than the adaptive one -- it must suppress the
+# high-mode deviation on average, not per call.  Measured (Csn=2):
+#   (12,1.0) 48 outer 31.4 s   (8,1.0) 58 outer 30.9 s   (8,0.5) 60 outer 34.3 s
+# Reference (adaptive, Fso6/Fss0.1): 56 outer, 28.57 s.  So walk the SHIFT UP
+# at orders 8 and 6, trading shift against order.  Nrhs=1 is the target.
+#          name          Fso  Fss   Csn  fine    coarse
+run_cell   S1_08_20       8    2.0   2    replay  gcr     # the promising cell: ~50 outer at 0.53 s/step would be ~26.5 s
+run_cell   S2_08_30       8    3.0   2    replay  gcr
+run_cell   S3_06_20       6    2.0   2    replay  gcr     # banked order, larger shift
+run_cell   S4_06_10       6    1.0   2    replay  gcr     # the L3 of the second table (never ran)
+run_cell   S5_10_20       10   2.0   2    replay  gcr     # between 8 and 12
+# coarse replay, recorded at a decent polynomial and backed off, at the best fine point
+run_cell   S6_coarse6     8    2.0   6    replay  replay
+run_cell   S7_coarse4     8    2.0   4    replay  replay
 
 echo "========================================================="
 echo "summary"
-grep -h "V2 3-level solve Nrhs 1" log.modes.* | sed 's/.*V2/V2/'
+for f in log.ladder.*; do echo "$f: $(grep -h "V2 3-level solve Nrhs 1" $f | sed "s/.*V2/V2/")   outer $(grep -h "Fouter MrhsPGCR: Converged" $f | tail -1 | grep -oE "iteration [0-9]+")"; done
+echo "reference (adaptive, banked): 28.57 s Nrhs=1, 14.9 s/RHS Nrhs=4"
 echo "========================================================="

tests/debug/Test_blockcyclic.cc

@@ -38,7 +38,7 @@ Author: Peter Boyle <pboyle@bnl.gov>
 //        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.
+//   T6 : ChooseProcessGrid: exact factorisation, Pr>=Pc, most-square.
 //////////////////////////////////////////////////////////////////////////////
 
 #include <Grid/Grid.h>
@@ -227,15 +227,15 @@ int main(int argc, char **argv)
       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 ( Pr < Pc )    ok = false;                 // Pr >= Pc (SUMMA B panels unpadded)
+      // most-square: no divisor r with Pc < r <= sqrt(P)
+      for(int r=Pc+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);
+    if ( !(Pr==18 && Pc==16) ) ok = false;   // Pr >= Pc: SUMMA B panels unpadded (S = ceil(Pc/Pr) = 1)
+    Report("T6  ChooseProcessGrid exact, Pr>=Pc, most-square (288 -> 18x16)", ok);
   }
 
   ////////////////////////////////////////////////////////////////////////

tests/debug/Test_poly_smoother.cc

@@ -0,0 +1,201 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_poly_smoother.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 */
+
+//////////////////////////////////////////////////////////////////////////////
+// Fixed-polynomial smoothers (Smoothers.h) against the adaptive GCR they
+// replace, on a shifted Wilson operator (non-Hermitian, spectrum in the right
+// half plane):
+//
+//   record : GCR(mmax=2, 8 steps) with a coefficient recorder on 16 sources
+//   T1     : GCRReplaySmoother on a FRESH source reduces the residual to
+//            within a factor 2 of the live GCR on the same source
+//   T2     : replay and GCR solutions agree to the coefficient spread (<10%)
+//   T3     : ChebyshevNonHermitianSmoother with Op=HermOp reproduces the legacy
+//            HermOp ChebyshevSmoother (real spectrum); on the complex-spectrum
+//            Wilson op its (poor) reduction is printed for information only
+//   T4     : the historical HermOp ChebyshevSmoother compiles with both
+//            constructor signatures and reduces the MdagM residual
+//   T5     : replay is bitwise repeatable (no reductions => no reordering)
+//
+//   mpirun -n 1 ./Test_poly_smoother --grid 8.8.8.8 --mpi 1.1.1.1
+//////////////////////////////////////////////////////////////////////////////
+
+#include <Grid/Grid.h>
+#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.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++;
+}
+
+template<class Field>
+class ShiftedOp : public LinearOperatorBase<Field> {
+  LinearOperatorBase<Field> &_Op; RealD shift;
+public:
+  ShiftedOp(RealD s, LinearOperatorBase<Field> &Op) : _Op(Op), shift(s) {}
+  void OpDiag  (const Field &in, Field &out) { GRID_ASSERT(0); }
+  void OpDir   (const Field &in, Field &out,int dir,int disp) { GRID_ASSERT(0); }
+  void OpDirAll(const Field &in, std::vector<Field> &out) { GRID_ASSERT(0); }
+  void Op      (const Field &in, Field &out) { _Op.Op(in,out);    out = out + shift*in; }
+  void AdjOp   (const Field &in, Field &out) { _Op.AdjOp(in,out); out = out + shift*in; }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); }
+  void HermOp  (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+template<class Field>
+RealD Residual(LinearOperatorBase<Field> &Op, const Field &src, const Field &x)
+{
+  Field r(src.Grid()); Op.Op(x,r); r = r - src;
+  return std::sqrt(norm2(r)/norm2(src));
+}
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  GridCartesian         *UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+                                       GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi());
+  GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+  std::vector<int> seeds({1,2,3,4});
+  GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds);
+
+  LatticeGaugeFieldD Umu(UGrid);
+  SU<Nc>::HotConfiguration(RNG4, Umu);
+
+  RealD mass = 0.5, shift = 0.1;
+  WilsonFermionD Dw(Umu, *UGrid, *UrbGrid, mass);
+  NonHermitianLinearOperator<WilsonFermionD, LatticeFermionD> Op(Dw);
+  ShiftedOp<LatticeFermionD> SOp(shift, Op);
+  TrivialPrecon<LatticeFermionD> simple;
+
+  const int mmax = 2, nstep = 8, ncal = 16;
+
+  //////////////////////////////////////////////////////////////////////
+  // record
+  //////////////////////////////////////////////////////////////////////
+  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> GCR(0.0, 1, SOp, simple, mmax, nstep);
+  GCR.SetZeroGuess(1); GCR.Name("smoother");
+  GCRCoefficients rec;
+  GCR.SetCoefficientRecorder(&rec);
+  LatticeFermionD src(UGrid), x(UGrid);
+  for(int c=0;c<ncal;c++){ gaussian(RNG4,src); x = Zero(); GCR(src,x); }
+  GCR.SetCoefficientRecorder(nullptr);
+  rec.Flush();
+  rec.Report("smoother");
+  Report("record: steps and calls", rec.Steps()==nstep && rec.Calls()==ncal,
+         std::to_string(rec.Steps())+" steps, "+std::to_string(rec.Calls())+" calls");
+
+  //////////////////////////////////////////////////////////////////////
+  // fresh source: live GCR vs replay
+  //////////////////////////////////////////////////////////////////////
+  gaussian(RNG4,src);
+  LatticeFermionD xg(UGrid), xr(UGrid), xr2(UGrid), xc(UGrid), d(UGrid);
+  xg = Zero(); GCR(src,xg);
+  RealD rg = Residual(SOp,src,xg);
+
+  GCRReplaySmoother<LatticeFermionD> Replay(SOp, rec);
+  Replay(src,xr);
+  RealD rr = Residual(SOp,src,xr);
+  Report("T1 replay residual within 2x of live GCR", rr < 2.0*rg,
+         "GCR |r|/|r0| = "+std::to_string(rg)+"  replay "+std::to_string(rr));
+
+  d = xr - xg;
+  RealD rel = std::sqrt(norm2(d)/norm2(xg));
+  Report("T2 replay solution vs GCR solution", rel < 0.1, "rel "+std::to_string(rel));
+
+  //////////////////////////////////////////////////////////////////////
+  // Chebyshev 1/x on [lo,hi], hi from a power iteration on SOp
+  //////////////////////////////////////////////////////////////////////
+  RealD hi;
+  {
+    LatticeFermionD v(UGrid), Av(UGrid); gaussian(RNG4,v);
+    RealD n = std::sqrt(norm2(v)); v = v*(1.0/n);
+    for(int i=0;i<60;i++){ SOp.Op(v,Av); hi = std::sqrt(norm2(Av)); v = Av*(1.0/hi); }
+  }
+  RealD lo = 0.5;
+  std::cout << GridLogMessage << "power iteration |lambda_max| ~ " << hi << "  Chebyshev range [" << lo << "," << 1.05*hi << "]" << std::endl;
+  ChebyshevNonHermitianSmoother<LatticeFermionD> Cheb(lo, 1.05*hi, nstep, SOp);
+  Cheb(src,xc);
+  RealD rc = Residual(SOp,src,xc);
+  // Not gated: Wilson's spectrum has O(1) imaginary parts and a real-interval
+  // Chebyshev fit to 1/x degrades exponentially off the axis (Bernstein
+  // ellipse).  Printed as the reminder of what a non-real spectrum does.
+  std::cout << GridLogMessage << "  (info) ChebyshevNonHermitian on the COMPLEX-spectrum Wilson op: |r|/|r0| = " << rc
+            << "   (GCR " << rg << ")  -- expected poor; PVdagM smoother ops have real coefficients" << std::endl;
+
+  //////////////////////////////////////////////////////////////////////
+  // legacy HermOp Chebyshev smoother, both constructor forms, on MdagM
+  //////////////////////////////////////////////////////////////////////
+  {
+    MdagMLinearOperator<WilsonFermionD, LatticeFermionD> HermOp(Dw);
+    LatticeFermionD hsrc(UGrid), hx(UGrid), hr(UGrid); gaussian(RNG4,hsrc);
+    RealD hhi = 0.0;
+    { LatticeFermionD v(UGrid), Av(UGrid); gaussian(RNG4,v); RealD n=std::sqrt(norm2(v)); v=v*(1.0/n);
+      for(int i=0;i<40;i++){ HermOp.HermOp(v,Av); hhi=std::sqrt(norm2(Av)); v=Av*(1.0/hhi); } }
+    ChebyshevSmoother<LatticeFermionD>                 S4(0.5, 1.05*hhi, 12, HermOp);
+    ChebyshevSmoother<LatticeFermionD, WilsonFermionD> S5(0.5, 1.05*hhi, 12, HermOp, Dw);   // historical 5-arg form
+    S4(hsrc,hx); HermOp.HermOp(hx,hr); hr = hr - hsrc;
+    RealD r4 = std::sqrt(norm2(hr)/norm2(hsrc));
+    S5(hsrc,hx); HermOp.HermOp(hx,hr); hr = hr - hsrc;
+    RealD r5 = std::sqrt(norm2(hr)/norm2(hsrc));
+    Report("T4 legacy ChebyshevSmoother (4- and 5-arg) reduces MdagM residual", r4 < 0.5 && r5 == r4,
+           "|r|/|r0| = "+std::to_string(r4)+" / "+std::to_string(r5));
+
+    // T3: the Op()-based Clenshaw on a REAL-spectrum operator must reproduce
+    // the legacy HermOp smoother: same coefficients, same recurrence.
+    struct HermAsOp : public LinearOperatorBase<LatticeFermionD> {
+      LinearOperatorBase<LatticeFermionD> &H;
+      HermAsOp(LinearOperatorBase<LatticeFermionD> &h) : H(h) {}
+      void OpDiag  (const LatticeFermionD &in, LatticeFermionD &out) { GRID_ASSERT(0); }
+      void OpDir   (const LatticeFermionD &in, LatticeFermionD &out,int dir,int disp) { GRID_ASSERT(0); }
+      void OpDirAll(const LatticeFermionD &in, std::vector<LatticeFermionD> &out) { GRID_ASSERT(0); }
+      void Op      (const LatticeFermionD &in, LatticeFermionD &out) { H.HermOp(in,out); }
+      void AdjOp   (const LatticeFermionD &in, LatticeFermionD &out) { H.HermOp(in,out); }
+      void HermOpAndNorm(const LatticeFermionD &in, LatticeFermionD &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); }
+      void HermOp  (const LatticeFermionD &in, LatticeFermionD &out) { H.HermOp(in,out); }
+    } HOp(HermOp);
+    ChebyshevNonHermitianSmoother<LatticeFermionD> C3(0.5, 1.05*hhi, 12, HOp);
+    LatticeFermionD hx3(UGrid), dd(UGrid);
+    C3(hsrc,hx3); HermOp.HermOp(hx3,hr); hr = hr - hsrc;
+    RealD r3 = std::sqrt(norm2(hr)/norm2(hsrc));
+    S4(hsrc,hx); dd = hx - hx3;
+    RealD reldiff = std::sqrt(norm2(dd)/norm2(hx));
+    Report("T3 ChebyshevNonHermitian(Op=HermOp) == legacy ChebyshevSmoother", r3 < 0.5 && reldiff < 1.0e-12,
+           "|r|/|r0| = "+std::to_string(r3)+"  rel diff of solutions "+std::to_string(reldiff));
+  }
+
+  //////////////////////////////////////////////////////////////////////
+  // determinism
+  //////////////////////////////////////////////////////////////////////
+  Replay(src,xr2);
+  d = xr - xr2;
+  Report("T5 replay bitwise repeatable", norm2(d)==0.0);
+
+  std::cout << GridLogMessage << (failures ? "Test_poly_smoother: FAILURES" : "Test_poly_smoother: ALL PASS") << std::endl;
+  Grid_finalize();
+  return failures ? 1 : 0;
+}

tests/debug/Test_ring_allreduce.cc

@@ -26,6 +26,7 @@ Author: Peter Boyle <pboyle@bnl.gov>
 //   T2 cartesian ring == GlobalSumVector, same sweep
 //   T3 bitwise repeatable (deterministic order)
 //   T4 timing at 16 MB, both rings vs GlobalSumVector
+//   T5 CartesianRingAllGather bitwise == padded GlobalSumVector; timing at the dense-apply shape
 //
 //   mpirun -n 4 ./Test_ring_allreduce --grid 16.16.16.32 --mpi 1.1.2.2
 // (2D process grid so the cartesian variant exercises more than one ring)
@@ -101,6 +102,54 @@ int main(int argc, char **argv)
   Check<RealF>   ("RealF   ", grid, 1.0e-5);
   Check<ComplexF>("ComplexF", grid, 1.0e-5);
 
+  // T5: CartesianRingAllGather delivers blocks in LEX-coordinate order; the
+  // reference is the zero-padded GlobalSumVector in RANK order, compared
+  // through the lex->rank table (bitwise: no arithmetic on either path).
+  // On a machine where ranks are relabelled (Frontier OptimalCommunicator)
+  // this is the test that catches a rank/coordinate confusion.
+  {
+    int P=grid->ProcessorCount(), me=grid->ThisRank(), mylex=CartesianLexIndex(grid);
+    std::vector<int> l2r(P);
+    for(int lp=0; lp<P; lp++){ Coordinate pc(grid->_ndimension); Lexicographic::CoorFromIndex(pc, lp, grid->_processors); l2r[lp]=grid->RankFromProcessorCoor(pc); }
+    Report("T5 lex->rank table consistent for my rank", l2r[mylex]==me);
+    int perm=0; for(int lp=0;lp<P;lp++) if(l2r[lp]!=lp) perm=1;
+    std::cout << GridLogMessage << "  (ranks " << (perm ? "ARE" : "are not") << " permuted relative to lex coordinates on this run)" << std::endl;
+    for(uint64_t chunk : std::vector<uint64_t>({1,3,64,1000,65537})){
+      uint64_t n=chunk*P;
+      std::vector<ComplexD> h(n,ComplexD(0.0,0.0)), ref;
+      for(uint64_t i=0;i<chunk;i++) h[me*chunk+i]=Fill<ComplexD>(me*chunk+i,me);   // rank-order reference
+      ref=h; grid->GlobalSumVector(&ref[0],(int)n);
+      std::vector<ComplexD> hin(n,ComplexD(0.0,0.0));
+      for(uint64_t i=0;i<chunk;i++) hin[mylex*chunk+i]=h[me*chunk+i];              // my slot is my lex index
+      deviceVector<ComplexD> d(n); acceleratorCopyToDevice(&hin[0],&d[0],n*sizeof(ComplexD));
+      CartesianRingAllGather(grid,&d[0],chunk);
+      std::vector<ComplexD> out(n); acceleratorCopyFromDevice(&d[0],&out[0],n*sizeof(ComplexD));
+      int bad=0; for(int lp=0;lp<P;lp++) if(memcmp(&out[lp*chunk],&ref[(uint64_t)l2r[lp]*chunk],chunk*sizeof(ComplexD))!=0) bad=1;
+      RealD diff=bad; grid->GlobalSum(diff);
+      Report("T5 CartesianRingAllGather (lex blocks) bitwise == rank-order padded GlobalSumVector, ComplexD chunk="+std::to_string(chunk), diff==0.0);
+    }
+    { uint64_t chunk=1001, n=chunk*P;
+      std::vector<ComplexF> h(n,ComplexF(0.0,0.0)), ref;
+      for(uint64_t i=0;i<chunk;i++) h[me*chunk+i]=Fill<ComplexF>(me*chunk+i,me);
+      ref=h; grid->GlobalSumVector(&ref[0],(int)n);
+      std::vector<ComplexF> hin(n,ComplexF(0.0,0.0)); for(uint64_t i=0;i<chunk;i++) hin[mylex*chunk+i]=h[me*chunk+i];
+      deviceVector<ComplexF> d(n); acceleratorCopyToDevice(&hin[0],&d[0],n*sizeof(ComplexF));
+      CartesianRingAllGather(grid,&d[0],chunk);
+      std::vector<ComplexF> out(n); acceleratorCopyFromDevice(&d[0],&out[0],n*sizeof(ComplexF));
+      int bad=0; for(int lp=0;lp<P;lp++) if(memcmp(&out[lp*chunk],&ref[(uint64_t)l2r[lp]*chunk],chunk*sizeof(ComplexF))!=0) bad=1;
+      RealD diff=bad; grid->GlobalSum(diff);
+      Report("T5 CartesianRingAllGather bitwise, ComplexF chunk=1001", diff==0.0);
+    }
+    // timing: the dense-apply shape, N=138240 x 4 rhs of ComplexF, chunk = N*4/P
+    { uint64_t chunk=(uint64_t)138240*4/P, n=chunk*P;
+      deviceVector<ComplexF> d(n); std::vector<ComplexF> h(n,ComplexF(1.0,0.0)); acceleratorCopyToDevice(&h[0],&d[0],n*sizeof(ComplexF));
+      double t0=usecond(); CartesianRingAllGather(grid,&d[0],chunk); double t1=usecond();
+      acceleratorCopyToDevice(&h[0],&d[0],n*sizeof(ComplexF));
+      double t2=usecond(); CartesianRingAllReduce(grid,&d[0],n); double t3=usecond();
+      std::cout << GridLogMessage << "T5 timing N=138240 x 4 ComplexF (" << n*8/1.0e6 << " MB): allgather " << (t1-t0)/1000. << " ms, cartesian allreduce " << (t3-t2)/1000. << " ms" << std::endl;
+    }
+  }
+
   // T4 timing at 16 MB of ComplexF (the dense-apply size at 12 RHS is 13.3 MB)
   {
     uint64_t n = 2*1024*1024;