From e5f1553232365383d020685aef4d279ff0f8342e Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 26 Aug 2026 23:36:31 -0400 Subject: [PATCH] Better reduction --- .../multigrid/BlockCyclicSchurInverse.h | 30 +++++++++++++ Grid/algorithms/multigrid/DenseCoarseMatrix.h | 27 ++++++++--- Grid/communicator/RingAllReduce.h | 45 +++++++++++-------- tests/debug/Test_ring_allreduce.cc | 32 ++++++++----- 4 files changed, 98 insertions(+), 36 deletions(-) diff --git a/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h b/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h index a530e54fe..c459abfd7 100644 --- a/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h +++ b/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h @@ -315,6 +315,36 @@ public: << " 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 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) diff --git a/Grid/algorithms/multigrid/DenseCoarseMatrix.h b/Grid/algorithms/multigrid/DenseCoarseMatrix.h index d91873fd5..72c9991ff 100644 --- a/Grid/algorithms/multigrid/DenseCoarseMatrix.h +++ b/Grid/algorithms/multigrid/DenseCoarseMatrix.h @@ -119,7 +119,9 @@ public: deviceVector dSlab; deviceVector dX; // N x MRHS_MAX deviceVector dY; // nrows x MRHS_MAX - deviceVector dG; // N x MRHS_MAX rank-major staging for the allgather (devSum==4) + deviceVector dG; // N x MRHS_MAX lex-major staging for the allgather (devSum==4) + deviceVector dLex2Rank;// lex index of a process coordinate -> its rank (allgather block order -> row-block order) + int myLex; deviceVector dPartial; // NK x (nrows x MRHS_MAX) deviceVector aptrs; // slab K-chunk pointers (lda = N) deviceVector xptrs; // X K-chunk pointers (ldb = N) @@ -257,7 +259,18 @@ public: "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); + 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 l2r(P); + for(int lp=0; lp_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; @@ -972,18 +985,18 @@ public: std::vector hG(chunk); for(int r=0;r dX[r*N + q*nrows + i] - ComplexF *g = &dG[0]; ComplexF *x = &dX[0]; + // 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 q = gi / nrw; int64_t i = gi - q*nrw; - x[idx] = g[q*(nrw*nrr) + r*nrw + i]; + 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(); diff --git a/Grid/communicator/RingAllReduce.h b/Grid/communicator/RingAllReduce.h index 7761347c7..310f28567 100644 --- a/Grid/communicator/RingAllReduce.h +++ b/Grid/communicator/RingAllReduce.h @@ -120,44 +120,52 @@ 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. On entry rank r's chunk is at -// buf[r*chunk]; on exit every rank holds all P chunks in RANK order. +// 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 the fastest-varying process coordinate -// (dim Nd-1) to the slowest: each stage is a ring over the P_d ranks of that -// line, after which the held block is the concatenation over that -// coordinate; because MPI Cartesian ranks are lexicographic with the last -// coordinate fastest, the final concatenation IS rank order -- no -// permutation. Bytes sent per rank ~ chunk*(P-1) ... dominated by the last -// stage, i.e. ~N = P*chunk total: 8x less than a zero-padded -// CartesianRingAllReduce of the same vector (which reduce-scatters AND -// gathers along every dimension). Steps: sum_d (P_d-1). Exact (no -// arithmetic): the result is bitwise the same as the padded allreduce. +// 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), measured 1.86 ms for 4.4 MB at 288 ranks with the -// allreduce ring -- at wire speed, but moving 35 MB per rank to deliver 4.4. +// 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 void CartesianRingAllGather(CartesianCommunicator *comm, T *buf, uint64_t chunk) { int P = comm->ProcessorCount(); - int me = comm->ThisRank(); if ( P==1 || chunk==0 ) return; int Nd = comm->_ndimension; + int mylex = CartesianLexIndex(comm); deviceVector 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)me*chunk; + 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=Nd-1; d>=0; d--){ + for(int d=0; d_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 ); - // place my block in slot med of the staging area (oth[0 .. Pd*blk)) acceleratorCopyDeviceToDevice((void *)(cur+off), (void *)(oth+(uint64_t)med*blk), blk*sizeof(T)); for(int t=1;tSendToRecvFrom((void *)(oth+(uint64_t)sendslot*blk), next, (void *)(oth+(uint64_t)recvslot*blk), prev, blk*sizeof(T)); } - // the staging area is the new held block T *tmp = cur; cur = oth; oth = tmp; off = 0; blk *= Pd; } diff --git a/tests/debug/Test_ring_allreduce.cc b/tests/debug/Test_ring_allreduce.cc index dc23156fb..9158d3bfa 100644 --- a/tests/debug/Test_ring_allreduce.cc +++ b/tests/debug/Test_ring_allreduce.cc @@ -102,30 +102,42 @@ int main(int argc, char **argv) Check ("RealF ", grid, 1.0e-5); Check("ComplexF", grid, 1.0e-5); - // T5: CartesianRingAllGather == zero-padded GlobalSumVector, BITWISE - // (no arithmetic in either path for disjoint chunks), all types, chunk - // sizes including 1 element and non-multiples of anything. + // 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(); + int P=grid->ProcessorCount(), me=grid->ThisRank(), mylex=CartesianLexIndex(grid); + std::vector l2r(P); + for(int lp=0; lp_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({1,3,64,1000,65537})){ uint64_t n=chunk*P; std::vector h(n,ComplexD(0.0,0.0)), ref; - for(uint64_t i=0;i(me*chunk+i,me); + for(uint64_t i=0;i(me*chunk+i,me); // rank-order reference ref=h; grid->GlobalSumVector(&ref[0],(int)n); - deviceVector d(n); acceleratorCopyToDevice(&h[0],&d[0],n*sizeof(ComplexD)); + std::vector hin(n,ComplexD(0.0,0.0)); + for(uint64_t i=0;i d(n); acceleratorCopyToDevice(&hin[0],&d[0],n*sizeof(ComplexD)); CartesianRingAllGather(grid,&d[0],chunk); std::vector out(n); acceleratorCopyFromDevice(&d[0],&out[0],n*sizeof(ComplexD)); - RealD diff=(memcmp(&out[0],&ref[0],n*sizeof(ComplexD))!=0)?1.0:0.0; grid->GlobalSum(diff); - Report("T5 CartesianRingAllGather bitwise == padded GlobalSumVector, ComplexD chunk="+std::to_string(chunk), diff==0.0); + int bad=0; for(int lp=0;lpGlobalSum(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 h(n,ComplexF(0.0,0.0)), ref; for(uint64_t i=0;i(me*chunk+i,me); ref=h; grid->GlobalSumVector(&ref[0],(int)n); - deviceVector d(n); acceleratorCopyToDevice(&h[0],&d[0],n*sizeof(ComplexF)); + std::vector hin(n,ComplexF(0.0,0.0)); for(uint64_t i=0;i d(n); acceleratorCopyToDevice(&hin[0],&d[0],n*sizeof(ComplexF)); CartesianRingAllGather(grid,&d[0],chunk); std::vector out(n); acceleratorCopyFromDevice(&d[0],&out[0],n*sizeof(ComplexF)); - RealD diff=(memcmp(&out[0],&ref[0],n*sizeof(ComplexF))!=0)?1.0:0.0; grid->GlobalSum(diff); + int bad=0; for(int lp=0;lpGlobalSum(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