From 357ede36642223f1e605bbf583e87d6e8f20bc83 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 3 Sep 2026 17:33:32 -0400 Subject: [PATCH] Improved FFT -- 1.5-2x when there are 2-6 ranks in a given axis of the cartesian communicator. Barrel shift -> all to all (x2) and distributed FFT work fully load balanced without redundant work. There is little more I can do now on FFT. Comms dominated and running distributed work dividing bandwidth optimal RingAllToAll /ccs/home/paboyle/ParallelIO/systems/Frontier/tests/core/Test_fft_prop --mpi 3.6.4.4 --grid 48.48.48.96 --accelerator-threads 8 --shm 4096 --shm-mpi 1 --device-mem 32000 --log Error,Warning,Message,Performance ************************************************* Benchmarking FFT of LatticeFermionD on plane wave ************************************************* Grid : Performance : 0.501524 s : FFT took 0.001311 s (transpose P=3) Grid : Performance : 0.501531 s : FFT pack 5.9e-05 s Grid : Performance : 0.501533 s : FFT alltoall 0.000828 s Grid : Performance : 0.501534 s : FFT reorder 0.000204 s Grid : Performance : 0.501535 s : FFT kernels 1e-05 s Grid : Performance : 0.501536 s : FFT unpack 5e-05 s Grid : Performance : 0.509992 s : FFT took 0.001829 s (transpose P=6) Grid : Performance : 0.510000 s : FFT pack 6e-05 s Grid : Performance : 0.510002 s : FFT alltoall 0.001436 s Grid : Performance : 0.510003 s : FFT reorder 0.000202 s Grid : Performance : 0.510005 s : FFT kernels 9e-06 s Grid : Performance : 0.510006 s : FFT unpack 5.3e-05 s Grid : Performance : 0.517690 s : FFT took 0.001599 s (transpose P=4) Grid : Performance : 0.517698 s : FFT pack 6e-05 s Grid : Performance : 0.517700 s : FFT alltoall 0.001258 s Grid : Performance : 0.517701 s : FFT reorder 0.0002 s Grid : Performance : 0.517702 s : FFT kernels 9e-06 s Grid : Performance : 0.517703 s : FFT unpack 4.9e-05 s Grid : Performance : 0.524858 s : FFT took 0.001561 s (transpose P=4) Grid : Performance : 0.524865 s : FFT pack 5.8e-05 s Grid : Performance : 0.524867 s : FFT alltoall 0.001213 s Grid : Performance : 0.524868 s : FFT reorder 0.000209 s Grid : Performance : 0.524869 s : FFT kernels 8e-06 s Grid : Performance : 0.524870 s : FFT unpack 4.9e-05 s ************************************************* FFT of [48 48 48 96] LatticeFermionD took 0.030916 s ************************************************* --- Grid/algorithms/FFT.h | 214 ++++++++++++++++++++++++++++- Grid/allocator/MemoryManager.cc | 2 +- Grid/communicator/RingAllReduce.h | 50 +++++++ systems/Frontier/nrhs_fault_36.job | 11 +- tests/debug/Test_ring_allreduce.cc | 24 ++++ 5 files changed, 295 insertions(+), 6 deletions(-) diff --git a/Grid/algorithms/FFT.h b/Grid/algorithms/FFT.h index c9d9d6a89..5b5ba223a 100644 --- a/Grid/algorithms/FFT.h +++ b/Grid/algorithms/FFT.h @@ -233,6 +233,15 @@ static void FFT_dim_execute( typedef typename vobj::vector_type vector_type; typedef typename FFTW::FFTW_scalar FFTW_scalar; +#if 0 + // ======================= ORIGINAL barrel-shift path ======================= + // Preserved for reference. Superseded by the transpose / all-to-all path + // below (the single active path for ALL P): the barrel is a P-fold redundant + // all-gather -- every rank assembles and transforms all Nperp lines of length + // G, keeping only its own L points. The transpose partitions the Nperp + // perpendicular lines across the P ranks along dim, so each rank transforms + // only ceil(Nperp/P) lines and moves (P-1)/P of its data instead of P-1 + // redundant copies. See CartesianRingAllToAll in communicator/RingAllReduce.h. const int Ndim = grid->Nd(); int L = grid->_ldimensions[dim]; int G = grid->_fdimensions[dim]; @@ -363,6 +372,203 @@ static void FFT_dim_execute( std::cout << GridLogPerformance << " of which shift" << t_shift/1.0e6 << " s" << std::endl; std::cout << GridLogPerformance << " FFT kernels " << t_fft/1.0e6 << " s" << std::endl; std::cout << GridLogPerformance << " FFT insert " << t_insert/1.0e6 << " s" << std::endl; +#endif + + // ==================== transpose / all-to-all pencil FFT ==================== + // + // P ranks lie along dim; each holds a local L-slab (L=ldim) of every global + // line (length G=fdim=L*P) and Nperp perpendicular lines. We partition the + // Nperp lines across the P ranks: the rank at coord `own` owns the lines with + // olin/Oloc == own, gathers their full G points via a cartesian all-to-all, + // transforms only Oloc = ceil(Nperp/P) of them, then scatters the result back + // with the inverse all-to-all. + // + // The owned-line count is CEIL-padded to a multiple of P (Oloc*P >= Nperp) so + // the all-to-all stays SYMMETRIC (one uniform chunk) for ANY (P,Nperp) -- in + // particular a P carrying a factor absent from Nperp, e.g. P=3 on the T axis + // of 128^3x288 whose perpendicular volume is a pure power of two. This keeps + // the transpose as total over decompositions as the barrel it replaces (no + // new geometry constraint), at a cost of <= (P-1) padded lines out of Nperp. + // Padding slots olin in [Nperp, Oloc*P) are never packed and never unpacked. + // + // The load-bearing identity: the all-to-all block index == cartesian coord + // along dim == which L-slab [c*L, c*L+L) of the global line -- so the block a + // rank receives carries the global-x tag needed to order the FFT input. + // + // Degenerate P: P=1 -> both all-to-alls are self-copies and the reorders are + // the identity (G=L), i.e. a pure local FFT. P=2 -> each all-to-all moves + // half a field, two of them one field, equal to the barrel's single Cshift. + { + const int Ndim = grid->Nd(); + int L = grid->_ldimensions[dim]; + int G = grid->_fdimensions[dim]; + int Ncomp = sizeof(sobj) / sizeof(scalar); + int P = grid->_processors[dim]; + int64_t Nperp = 1; + for (int d = 0; d < Ndim; d++) + if (d != dim) Nperp *= grid->_ldimensions[d]; + + int64_t Oloc = (Nperp + P - 1) / P; // ceil: owned (padded) lines/rank + int64_t chunk = (int64_t)L * Oloc * Ncomp; // one all-to-all block (uniform) + int64_t nbuf = (int64_t)P * chunk; // == Ncomp*Oloc*G, one field's worth + int64_t howmany_local = (int64_t)Ncomp * Oloc; + + scalar div; + if (sign == FFTW_BACKWARD) div = 1.0 / G; + else if (sign == FFTW_FORWARD) div = 1.0; + else GRID_ASSERT(0); + + double t_total = -usecond(); + double t_pack = 0, t_a2a = 0, t_reorder = 0, t_fft = 0, t_unpack = 0; + + deviceVector sbuf(nbuf); + deviceVector rbuf(nbuf); + deviceVector pgbuf(nbuf); // FFTW pencil buffer, Ncomp*Oloc lines of G + scalar *sbuf_v = &sbuf[0]; + scalar *rbuf_v = &rbuf[0]; + scalar *pgbuf_v = &pgbuf[0]; + + // deterministic ceil-pad slots (never read back, but keeps padded FFT lines finite) + acceleratorMemSet(sbuf_v, 0, nbuf*sizeof(scalar)); + + const Coordinate ldims = grid->_ldimensions; + const Coordinate rdims = grid->_rdimensions; + const Coordinate sdims = grid->_simd_layout; + const int Nsimd = vobj::Nsimd(); + + // ---- 1. pack: source -> sbuf. block = owner coord; payload xloc + L*(slot + Oloc*w) + t_pack -= usecond(); + { + autoView(s_v, source, AcceleratorRead); + accelerator_for(idx, grid->oSites(), Nsimd, { +#ifdef GRID_SIMT + { + int lane = acceleratorSIMTlane(Nsimd); +#else + for (int lane = 0; lane < Nsimd; lane++) { +#endif + Coordinate icoor(Ndim), ocoor(Ndim); + Lexicographic::CoorFromIndex(icoor, lane, sdims); + Lexicographic::CoorFromIndex(ocoor, idx, rdims); + int64_t xloc = ocoor[dim] + icoor[dim]*rdims[dim]; + int64_t olin = 0, str = 1; + for (int d = 0; d < Ndim; d++) { + if (d == dim) continue; + int64_t c = ocoor[d] + icoor[d]*rdims[d]; + olin += str * c; + str *= ldims[d]; + } + int64_t own = olin / Oloc; + int64_t slot = olin - own*Oloc; + vector_type *from = (vector_type *)&s_v[idx]; + for (int w = 0; w < Ncomp; w++) { + scalar_type stmp = getlane(from[w], lane); + sbuf_v[ own*chunk + xloc + L*(slot + Oloc*w) ] = stmp; + } +#ifdef GRID_SIMT + } +#else + } +#endif + }); + } + t_pack += usecond(); + + // ---- 2. forward all-to-all: gather my owned lines' L-slabs from every rank + t_a2a -= usecond(); + CartesianRingAllToAll(grid, sbuf_v, rbuf_v, (uint64_t)chunk, dim); + t_a2a += usecond(); + + // ---- 3. reorder rbuf -> pgbuf: contiguous G-lines (w,slot), xpos = src*L + xloc + t_reorder -= usecond(); + accelerator_for(q, nbuf, 1, { + int64_t xpos = q % G; + int64_t t = q / G; // = w*Oloc + slot + int64_t slot = t % Oloc; + int64_t w = t / Oloc; + int64_t src = xpos / L; + int64_t xloc = xpos % L; + pgbuf_v[q] = rbuf_v[ src*chunk + xloc + L*(slot + Oloc*w) ]; + }); + t_reorder += usecond(); + + // ---- 4. FFT: Ncomp*Oloc contiguous lines of length G (istride 1, idist G) + { + FFTW_scalar *in = (FFTW_scalar *)pgbuf_v; + FFTW_scalar *out = (FFTW_scalar *)pgbuf_v; + t_fft -= usecond(); + FFTW::fftw_execute_dft(p, in, out, sign); + t_fft += usecond(); + } + flops_call = 5.0 * (double)howmany_local * G * log2(G); + usec = (uint64_t)t_fft; + flops = flops_call; + + // ---- 5a. reorder pgbuf -> sbuf: block = destination coord; xpos = dst*L + xloc + t_reorder -= usecond(); + accelerator_for(j, nbuf, 1, { + int64_t dst = j / chunk; + int64_t r = j % chunk; + int64_t xloc = r % L; + int64_t u = r / L; // = slot + Oloc*w + int64_t slot = u % Oloc; + int64_t w = u / Oloc; + int64_t xpos = dst*L + xloc; + sbuf_v[j] = pgbuf_v[ w*Oloc*G + slot*G + xpos ]; + }); + t_reorder += usecond(); + + // ---- 5b. inverse all-to-all: scatter transformed L-slabs back + t_a2a -= usecond(); + CartesianRingAllToAll(grid, sbuf_v, rbuf_v, (uint64_t)chunk, dim); + t_a2a += usecond(); + + // ---- 5c. unpack rbuf -> result (x div); block = owner coord of each line + t_unpack -= usecond(); + { + autoView(r_v, result, AcceleratorWrite); + accelerator_for(idx, grid->oSites(), Nsimd, { +#ifdef GRID_SIMT + { + int lane = acceleratorSIMTlane(Nsimd); +#else + for (int lane = 0; lane < Nsimd; lane++) { +#endif + Coordinate icoor(Ndim), ocoor(Ndim); + Lexicographic::CoorFromIndex(icoor, lane, sdims); + Lexicographic::CoorFromIndex(ocoor, idx, rdims); + int64_t xloc = ocoor[dim] + icoor[dim]*rdims[dim]; + int64_t olin = 0, str = 1; + for (int d = 0; d < Ndim; d++) { + if (d == dim) continue; + int64_t c = ocoor[d] + icoor[d]*rdims[d]; + olin += str * c; + str *= ldims[d]; + } + int64_t own = olin / Oloc; + int64_t slot = olin - own*Oloc; + vector_type *to = (vector_type *)&r_v[idx]; + for (int w = 0; w < Ncomp; w++) { + scalar_type stmp = div * rbuf_v[ own*chunk + xloc + L*(slot + Oloc*w) ]; + putlane(to[w], stmp, lane); + } +#ifdef GRID_SIMT + } +#else + } +#endif + }); + } + t_unpack += usecond(); + t_total += usecond(); + + std::cout << GridLogPerformance << " FFT took " << t_total/1.0e6 << " s (transpose P=" << P << ")" << std::endl; + std::cout << GridLogPerformance << " FFT pack " << t_pack/1.0e6 << " s" << std::endl; + std::cout << GridLogPerformance << " FFT alltoall " << t_a2a/1.0e6 << " s" << std::endl; + std::cout << GridLogPerformance << " FFT reorder " << t_reorder/1.0e6 << " s" << std::endl; + std::cout << GridLogPerformance << " FFT kernels " << t_fft/1.0e6 << " s" << std::endl; + std::cout << GridLogPerformance << " FFT unpack " << t_unpack/1.0e6 << " s" << std::endl; + } } class FFT : public FFTbase { @@ -405,8 +611,10 @@ public: int64_t Nperp = 1; for (int d = 0; d < Ndim; d++) if (d != dim) Nperp *= _grid->_ldimensions[d]; + int P = _grid->_processors[dim]; + int64_t Oloc = (Nperp + P - 1) / P; // ceil-padded owned lines/rank (see FFT_dim_execute) int n[] = {G}; - int howmany = Ncomp * Nperp; + int howmany = Ncomp * (int)Oloc; deviceVector dummy(2); FFTW_scalar *buf = (FFTW_scalar *)&dummy[0]; @@ -442,7 +650,9 @@ private: int64_t Nperp = 1; for (int dd = 0; dd < Ndim; dd++) if (dd != d) Nperp *= _grid->_ldimensions[dd]; - int howmany = Ncomp * (int)Nperp; + int P = _grid->_processors[d]; + int64_t Oloc = (Nperp + P - 1) / P; // ceil-padded owned lines/rank (see FFT_dim_execute) + int howmany = Ncomp * (int)Oloc; int n[] = {G}; deviceVector dummy(2); diff --git a/Grid/allocator/MemoryManager.cc b/Grid/allocator/MemoryManager.cc index c8f4b27d6..db84cde41 100644 --- a/Grid/allocator/MemoryManager.cc +++ b/Grid/allocator/MemoryManager.cc @@ -366,7 +366,7 @@ void *MemoryManager::Insert(void *ptr,size_t bytes,AllocationCacheEntry *entries GRID_ASSERT(omp_in_parallel()==0); #endif - if (ncache == 0) return ptr; + if (ncache == 0) { if(freed) *freed = bytes; return ptr; } // uncached: the incoming block is the one freed void * ret = NULL; int v = -1; diff --git a/Grid/communicator/RingAllReduce.h b/Grid/communicator/RingAllReduce.h index b78a28039..54d503d17 100644 --- a/Grid/communicator/RingAllReduce.h +++ b/Grid/communicator/RingAllReduce.h @@ -201,4 +201,54 @@ void CartesianRingAllGather(CartesianCommunicator *comm, T *buf, uint64_t chunk, if ( cur != buf ) acceleratorCopyDeviceToDevice((void *)cur, (void *)buf, blk*sizeof(T)); } +///////////////////////////////////////////////////////////////////////////// +// Cartesian ALL-TO-ALL along one processor dimension, point-to-point only. +// +// CartesianRingAllToAll(comm, sbuf, rbuf, chunk, dim) +// sbuf, rbuf each hold P_dim*chunk elements of T. Block index is the +// process COORDINATE along dim (not the MPI rank -- ShiftedRanks resolves +// the OptimalCommunicator relabelling): sbuf[c*chunk] is my data destined +// for the rank at coordinate c along dim; on exit rbuf[c*chunk] is the +// block that rank sent me. All other coordinates equal mine. +// +// Direct-pairwise, NOT ring-forwarding: P_dim-1 symmetric exchanges, step k +// with the +/-k neighbours (ShiftedRanks == MPI_Cart_shift, arbitrary shift), +// each moving one chunk with no store-and-forward. Per-rank traffic +// (P_dim-1)*chunk -- a factor P below the AllGather's replicated P*chunk, which +// is the point for the pencil-FFT transpose (partition the orthogonal coords +// rather than replicate the whole transform line on every rank). +// +// On a switched fabric a distance-k send is one logical hop, so direct-pairwise +// moves the minimum; ring-forwarding would only win on a true physical ring. +// Every step is one SendToRecvFrom of one chunk: no collective, no size cliff. +// Send/recv-rank binding follows Cshift_mpi.h: ShiftedRanks(dim,k,xmit,recv) +// sends to coord-k, receives from coord+k. Steps: P_dim-1. Exact. +///////////////////////////////////////////////////////////////////////////// +template +void CartesianRingAllToAll(CartesianCommunicator *comm, + T *sbuf, T *rbuf, uint64_t chunk, int dim) +{ + int Nd = comm->_ndimension; + GRID_ASSERT( dim >= 0 && dim < Nd ); + int P = comm->_processors[dim]; + int me = comm->_processor_coor[dim]; + if ( chunk==0 ) return; + GRID_ASSERT( (chunk*sizeof(T))%4 == 0 ); // SendToRecvFrom counts int32 words + uint64_t bytes = chunk*sizeof(T); + + // my own block never goes on the wire + acceleratorCopyDeviceToDevice((void *)&sbuf[(uint64_t)me*chunk], + (void *)&rbuf[(uint64_t)me*chunk], bytes); + if ( P==1 ) return; + + for(int k=1;kShiftedRanks(dim, k, xmit_to_rank, recv_from_rank); // xmit=coord-k, recv=coord+k + uint64_t sidx = (uint64_t)((me - k + P) % P); // block destined for the coord-k rank + uint64_t ridx = (uint64_t)((me + k) % P); // block arriving from the coord+k rank + comm->SendToRecvFrom((void *)&sbuf[sidx*chunk], xmit_to_rank, + (void *)&rbuf[ridx*chunk], recv_from_rank, bytes); + } +} + NAMESPACE_END(Grid); diff --git a/systems/Frontier/nrhs_fault_36.job b/systems/Frontier/nrhs_fault_36.job index b80e1d473..23328c2f2 100644 --- a/systems/Frontier/nrhs_fault_36.job +++ b/systems/Frontier/nrhs_fault_36.job @@ -69,6 +69,7 @@ module load libfabric BIN=$root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix OPTS1="--accelerator-threads 8 --shm 4096 --shm-mpi 1 --device-mem 32000" +OVERLAP="--comms-overlap" # set empty for a no-overlap job vol=48.48.48.96 MPI_GEOM=3.6.4.4 @@ -99,7 +100,7 @@ run_cell () { name=$1; shift echo "----- $name : $* -----" export GRID_STDOUT_ROOT=$RUNDIR/fault36_$name - env "$@" srun -N36 -n288 --kill-on-bad-exit=1 ./select_gpu $BIN --mpi ${MPI_GEOM} --grid $vol $OPTS1 --comms-overlap \ + env "$@" srun -N36 -n288 --kill-on-bad-exit=1 ./select_gpu $BIN --mpi ${MPI_GEOM} --grid $vol $OPTS1 $OVERLAP \ --debug-stdout --log Error,Warning,Message,Memory > log.fault36.$name 2>&1 echo " exit $?"; sleep 30 f=$(grep -l "Memory access fault" $GRID_STDOUT_ROOT/*/Grid.stderr.* 2>/dev/null | head -1) @@ -117,6 +118,8 @@ run_cell () { #run_cell A_serialised AMD_SERIALIZE_KERNEL=3 AMD_SERIALIZE_COPY=3 AMD_LOG_LEVEL=3 # 5371703-era: OOM found; done #run_cell B_plain AMD_LOG_LEVEL=1 +# NOTE (2026-08-29): FineSloppyComms is NOT set in this job; the example defaults it to 1, so every +# cell below ran with SLOPPY comms ON (confirmed by the PARAM line). Only 5371826/5372414 set it to 0. # Ladder for the NO_TRANSLATION / hang fault (both at NRHS=12): # 5371703 sloppy ON, no kdreg2 -> NO_TRANSLATION (Dhop halo) at outer 36 # 5371826 sloppy OFF, no kdreg2 -> 12-RHS converged, HANG in single-RHS solve @@ -146,5 +149,7 @@ run_cell L_cache4 FI_MR_CACHE_MAX_COUNT=4 # PB: room for Pa # at the first exchange (as the reproducer predicts for the cache-off path); G-M all solved # Nrhs 6 + Nrhs 1 (NOTE: NRHS=6 above, the three failures in the ledger were at NRHS=12). # Missing control: the new PaddedCell, NO env knobs, at NRHS=12, twice (intermittent failure). -run_cell N_control NRHS=12 -run_cell N_control2 NRHS=12 +run_cell N_control NRHS=12 # 5374670: PASSED +run_cell N_control2 NRHS=12 # 5374670: NO_TRANSLATION, Dhop Waitall(16), 3 ranks; last logged event on the faulting rank = eviction 49 of the 62-eviction restart burst after outer step 48 (170 MB hipFree via the ring cache), the failing exchange followed it +# Overlap on/off is a WHOLE-JOB A/B: set OVERLAP empty above and submit again. Measured +# 2026-08-29/30: it does not change the failure rate. diff --git a/tests/debug/Test_ring_allreduce.cc b/tests/debug/Test_ring_allreduce.cc index c6f110fb0..d8410e9ad 100644 --- a/tests/debug/Test_ring_allreduce.cc +++ b/tests/debug/Test_ring_allreduce.cc @@ -204,6 +204,30 @@ int main(int argc, char **argv) } } + // T8: single-dimension ALL-TO-ALL. Fingerprint block j (destined for coord j) + // on the rank at coord m with a tag encoding the (source,dest) coordinate pair, + // src*Pd+dst; after the exchange block i (from coord i) must carry the tag the + // coord-i rank prepared for me, i.e. Fill(e, i*Pd + m). Analytic reference, so + // a coordinate<->rank relabelling bug (wrong ShiftedRanks partner) shows up as a + // mismatched source coordinate rather than silently passing. + { + int Nd=grid->_ndimension; + for(int d=0; d_processors[d]; if ( Pd==1 ) continue; + int m=grid->_processor_coor[d]; + uint64_t chunk=1009, n=chunk*(uint64_t)Pd; + std::vector sh(n), rh(n), ref(n); + for(int j=0;j(e, m*Pd+j); + for(int i=0;i(e, i*Pd+m); + deviceVector sd(n), rd(n); + acceleratorCopyToDevice(&sh[0],&sd[0],n*sizeof(ComplexD)); + CartesianRingAllToAll(grid,&sd[0],&rd[0],chunk,d); + acceleratorCopyFromDevice(&rd[0],&rh[0],n*sizeof(ComplexD)); + RealD diff=(memcmp(&rh[0],&ref[0],n*sizeof(ComplexD))!=0)?1.0:0.0; grid->GlobalSum(diff); + Report("T8 CartesianRingAllToAll(dim="+std::to_string(d)+") bitwise == analytic reference, P_d="+std::to_string(Pd), diff==0.0); + } + } + // T4 timing at 16 MB of ComplexF (the dense-apply size at 12 RHS is 13.3 MB) { uint64_t n = 2*1024*1024;