From effa31e761a48a2d271e0cfc905a5824c0c9cce5 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Mon, 24 Aug 2026 15:45:06 -0400
Subject: [PATCH] Drop redundant set zero's and optimise persistent history
 vectors for smoothers

---
 ...GeneralisedConjugateResidualNonHermitian.h | 22 ++++++++++++++++---
 .../multigrid/RecursiveSchurInverse.h         | 16 +++++++++++---
 2 files changed, 32 insertions(+), 6 deletions(-)

diff --git a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h
index f4ff8b410..5e24a5eb1 100644
--- a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h
+++ b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h
@@ -57,6 +57,11 @@ public:
   GridStopWatch LinalgTimer;
   std::string name;
   int ZeroGuess  = 0;  // caller contract: guess is always zero => first-cycle r0 = src, skip the apply
+  // persistent GCR history (see GCRnStep)
+  GridBase           *hist_grid = nullptr;
+  std::vector<Field>  q;
+  std::vector<Field>  p;
+  std::vector<RealD>  qq;
   int FirstCycle = 0;
 
   LinearFunction<Field>     &Preconditioner;
@@ -145,9 +150,20 @@ public:
     ////////////////////////////////
     // history for flexible orthog
     ////////////////////////////////
-    std::vector<Field> q(mmax,grid);
-    std::vector<Field> p(mmax,grid);
-    std::vector<RealD> qq(mmax);
+    // History arrays are PERSISTENT across calls (allocated once per grid,
+    // re-made only if the grid or mmax changes).  The per-call form
+    // std::vector<Field>(mmax,grid) built a temporary and copy-constructed it
+    // mmax times on every restart cycle -- measured ~5 ms per fine-smoother
+    // call.  Safe: every entry is written before it is read within a cycle
+    // (q[kp],p[kp] assigned before the northog loop can reach them), so no
+    // stale content is ever consumed.
+    if ( hist_grid != grid || (int)q.size() != mmax ) {
+      q.clear(); p.clear();
+      q.reserve(mmax); p.reserve(mmax);
+      for(int i=0;i<mmax;i++){ q.emplace_back(grid); p.emplace_back(grid); }
+      qq.assign(mmax, 0.0);
+      hist_grid = grid;
+    }
       
     GCRLogLevel<< "PGCR nStep("<<nstep<<")"<<std::endl;
 
diff --git a/Grid/algorithms/multigrid/RecursiveSchurInverse.h b/Grid/algorithms/multigrid/RecursiveSchurInverse.h
index 985b8dfee..dd62d4f0c 100644
--- a/Grid/algorithms/multigrid/RecursiveSchurInverse.h
+++ b/Grid/algorithms/multigrid/RecursiveSchurInverse.h
@@ -120,7 +120,9 @@ public:
   uint64_t              bytesAllreduce;
   uint64_t              nAllreduce;    // panel collectives
   uint64_t              nGatherGemm;   // GatherGemm calls
-  uint64_t              nGather;       // AllGatherV collectives (debug gate)
+  uint64_t              nGather;       // gather-path collectives (debug gate)
+  uint64_t              nGatherV;      // ... of which MPI_Allgatherv
+  uint64_t              nBcast;        // ... of which Bcast-assembled
 
   // Persistent grow-only device panel; assembly and collectives are
   // device-resident.  Device builds require GPU-aware MPI.
@@ -368,7 +370,10 @@ public:
       uint64_t panelWords = (uint64_t)k*nchunk;
       uint64_t panelBytesThis = panelWords*sizeof(ComplexD);
 
-      int gatherThis = useGather && ( (int64_t)panelBytesThis >= gatherMinBytes );
+      // The MODE (0/1/2), not a boolean: `useGather && cond` collapses 2 to 1
+      // and silently dispatched DENSE_GATHER=2 onto the known-broken
+      // AllGatherV path (Frontier hang, 2026-08-24).  Identical on all ranks.
+      int gatherThis = ( (int64_t)panelBytesThis >= gatherMinBytes ) ? useGather : 0;
 
       // Arrival skew is charged to the barrier, so that tAllreduce measures
       // transfer alone.  Diagnostic only; off by default.
@@ -417,6 +422,7 @@ public:
         }
         if ( gatherThis == 1 )
         {
+          nGatherV++;
           void *send = owner ? (void *)B.ColumnWindow(colB+j0) : (void *)&dRecvBuf[0];
           grid->AllGatherV(send, counts[me],
                            (void *)&dRecvBuf[0], counts, displs, sizeof(ComplexD));
@@ -431,6 +437,7 @@ public:
           // Every rank issues all C broadcasts in the same order, so
           // collective matching stays positional and safe.
           ////////////////////////////////////////////////////////////////
+          nBcast++;
           if ( owner )
           {
             int64_t off_me = rowStart[me] - rowStart[rB0];
@@ -777,6 +784,8 @@ public:
     nAllreduce     = 0;
     nGatherGemm    = 0;
     nGather        = 0;
+    nGatherV       = 0;
+    nBcast         = 0;
 
     SchurNode(0, P, 0, N, Arows);
 
@@ -820,7 +829,8 @@ public:
               << "  leaf "      << tLeaf/1.0e6
               << std::endl;
     std::cout << GridLogMessage << "Schur comms:"
-              << ( useGather ? "  [AllGatherV]" : "  [zero-fill+GlobalSum]" )
+              << "  [transports run: allreduce " << (nAllreduce - nGatherV - nBcast)
+              << " allgatherv " << nGatherV << " bcast " << nBcast << "]"
               << ( barrierProbe ? " [barrier probe on]" : "" )
               << "  GatherGemm calls " << nGatherGemm
               << "  panel allreduces " << nAllreduce