SUMMA telemetry: ring time histogram by message size

Per log2-size bucket: messages, GB, seconds inside SendToRecvFrom, GB/s,
% of ring time.  Decomposes the 2 GB/s/rank average (probe: 11-20 GB/s
at 8 MB) into latency-bound small messages vs slow large ones vs wait.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RAhdQHrkzKmzfxpxW1whLn

Grid/algorithms/multigrid/BlockCyclicSchurInverse.h

@@ -397,6 +397,25 @@ public:
               << " gemm min/max " << gmin << "/" << gmax
               << " | ring bytes/rank " << gb << " GB -> " << (ring>0 ? gb/ring : 0.0) << " GB/s/rank (boss)"
               << std::endl;
+    // Ring time decomposed by message size (boss rank; every rank sends the
+    // same sequence of sizes).  Time is wall time inside SendToRecvFrom, so it
+    // includes waiting for the partner -- a bucket whose GB/s is far below the
+    // probe's for the same size is wait, not wire.
+    std::cout << GridLogMessage << "BlockCyclicSumma ring histogram (boss):  size-bucket  msgs  GB  secs  GB/s  %time" << std::endl;
+    std::streamsize oldprec = std::cout.precision();
+    for(int b=0;b<SUMMA.NHIST;b++){
+      if ( !SUMMA.histN[b] ) continue;
+      double sec = SUMMA.histUs[b]/1.0e6, g = SUMMA.histBytes[b]/1.0e9;
+      double lo = (double)(1ull<<b);
+      char sz[32]; if (lo>=1048576) snprintf(sz,32,"%6.1f MB",lo/1048576.); else if (lo>=1024) snprintf(sz,32,"%6.1f KB",lo/1024.); else snprintf(sz,32,"%6.0f B ",lo);
+      std::cout << GridLogMessage << "   >=" << sz
+                << std::setw(8) << SUMMA.histN[b]
+                << std::setw(10) << std::setprecision(3) << g
+                << std::setw(9) << std::setprecision(3) << sec
+                << std::setw(9) << std::setprecision(3) << (sec>0 ? g/sec : 0.0)
+                << std::setw(8) << std::setprecision(3) << (ring>0 ? 100.0*sec/ring : 0.0) << std::endl;
+    }
+    std::cout.precision(oldprec);
   }
 };
 

Grid/algorithms/multigrid/BlockCyclicSumma.h

@@ -143,7 +143,14 @@ public:
   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; }
+  // Per-message-size histogram (bucket = floor(log2 bytes)): count, bytes,
+  // microseconds -- decomposes the ring time by packet size so a low average
+  // GB/s can be attributed (many small latency-bound messages vs slow large
+  // ones vs partner-wait).  The 8 MB probe runs at 11-20 GB/s; SUMMA averaged 2.
+  static const int NHIST=48;
+  uint64_t histN[NHIST]={0}, histBytes[NHIST]={0}; double histUs[NHIST]={0};
+  void HistAdd(uint64_t bytes, double us){ int b=0; while((bytes>>b)>1) b++; histN[b]++; histBytes[b]+=bytes; histUs[b]+=us; }
+  void ResetTelemetry(void){ tAlloc=tPack=tRingA=tRingB=tGemm=0; bytesRing=nRingMsg=nMultiply=nGemm=0; for(int b=0;b<NHIST;b++){histN[b]=histBytes[b]=0; histUs[b]=0;} }
 
   static int Overlap(int64_t a0,int64_t a1,int64_t b0,int64_t b1)
   { return (a0 < b1) && (b0 < a1); }
@@ -273,9 +280,11 @@ public:
         for(int t=1;t<Pc;t++){
           int cs = (pcol - t + 1 + Pc*Pc) % Pc;
           int cr = (pcol - t     + Pc*Pc) % Pc;
+          double tm = usecond();
           grid->SendToRecvFrom((void *)(&Abuf[0]+slotA*cs), dest,
                                (void *)(&Abuf[0]+slotA*cr), src,
                                slotA*sizeof(ComplexD));
+          HistAdd(slotA*sizeof(ComplexD), usecond()-tm);
           bytesRing += slotA*sizeof(ComplexD); nRingMsg++;
         }
         tRingA += usecond();
@@ -291,9 +300,11 @@ public:
         for(int t=1;t<Pr;t++){
           int rs = (prow - t + 1 + Pr*Pr) % Pr;
           int rr = (prow - t     + Pr*Pr) % Pr;
+          double tm = usecond();
           grid->SendToRecvFrom((void *)(&Bbuf[0]+slotB*rs), dest,
                                (void *)(&Bbuf[0]+slotB*rr), src,
                                slotB*sizeof(ComplexD));
+          HistAdd(slotB*sizeof(ComplexD), usecond()-tm);
           bytesRing += slotB*sizeof(ComplexD); nRingMsg++;
         }
         tRingB += usecond();

Grid/algorithms/multigrid/DenseCoarseMatrix.h

@@ -121,6 +121,7 @@ public:
   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)
+  deviceVector<int64_t>   dRm2G;    // rank-major index (rank*nrows + ss*nbasis + b) -> global column (gsite*nbasis + b) of x / the slab
   int                     myLex;
   deviceVector<ComplexF>  dPartial; // NK x (nrows x MRHS_MAX)
   deviceVector<ComplexF*> aptrs;    // slab K-chunk pointers   (lda = N)
@@ -270,6 +271,15 @@ public:
         acceleratorCopyToDevice(&l2r[0], &dLex2Rank[0], P*sizeof(int));
         myLex = CartesianLexIndex(grid);
         GRID_ASSERT( l2r[myLex] == grid->ThisRank() );
+        // x and the slab columns are in GLOBAL-SITE order (hX[myGsite*nbasis+b]);
+        // the gathered blocks are in RANK-MAJOR order (rank*nrows + ss*nbasis + b).
+        // On one rank the two coincide, which is how the laptop passed while
+        // Frontier VERIFYed 0.9965 (2026-08-27).  Scatter through the inverse map.
+        std::vector<int64_t> g2rm; BuildRankMajorMap(g2rm);
+        std::vector<int64_t> rm2g(N); for(int64_t g=0; g<N; g++) rm2g[g2rm[g]] = g;
+        for(int ss=0; ss<lsites; ss++) GRID_ASSERT( rm2g[(int64_t)grid->ThisRank()*nrows + (int64_t)ss*nbasis] == myGsite[ss]*nbasis );
+        dRm2G.resize(N);
+        acceleratorCopyToDevice(&rm2g[0], &dRm2G[0], N*sizeof(int64_t));
       }
       std::cout << GridLogMessage << "DenseCoarseMatrix: slab resident on device ("
                 << sbytes/1024./1024. << " MB/rank), split-K NK=" << NK << " (Kc=" << Kc << "); "
@@ -974,29 +984,30 @@ public:
     double t1 = usecond();
     double t2, t3;
     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]
+      // ALLGATHER: x is not a reduction -- every rank owns the rows of x at
+      // global columns myGsite[ss]*nbasis+b (scattered by site coordinate, NOT
+      // a contiguous block) and needs all of it.  Only MY rows go host->device
+      // (nrows x nr, ~15 KB at nr=1), packed rank-major [r][ss*nbasis+b],
+      // gathered along the process grid, then scattered through rm2g into the
+      // column-major dX (ld = N, global-site columns) the split-K GEMM reads.
+      const uint64_t chunk = (uint64_t)nrows*nr;              // my block: [r][ss*nbasis+b]
       { 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));
+          for(int ss=0; ss<lsites; ss++)
+            memcpy(&hG[(uint64_t)r*nrows + (uint64_t)ss*nbasis], &hX[(uint64_t)r*N + (uint64_t)myGsite[ss]*nbasis], nbasis*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];
+        // scatter lex block L=[r][i] -> dX[r*N + rm2g[rank(L)*nrows + i]]
+        ComplexF *g = &dG[0]; ComplexF *x = &dX[0]; int *l2r = &dLex2Rank[0]; int64_t *rm2g = &dRm2G[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];
+          x[r*NN + rm2g[(int64_t)l2r[L]*nrw + i]] = g[L*(nrw*nrr) + r*nrw + i];
         });
       }
       t3 = usecond();