diff --git a/Grid/lattice/Lattice_reduction.h b/Grid/lattice/Lattice_reduction.h index d3b612a0e..4fd86c674 100644 --- a/Grid/lattice/Lattice_reduction.h +++ b/Grid/lattice/Lattice_reduction.h @@ -809,6 +809,21 @@ static void sliceInnerProductMatrix( Eigen::MatrixXcd &mat, const Lattice // Same code path for every Lattice: fine fermion fields and coarse // multi-RHS fields (nrhs folded into the grid) alike. ///////////////////////////////////////////////////////////////////////////// +// The batch of views and coefficients is passed to the kernel BY VALUE as +// lambda-captured kernel arguments (a ViewPack), so there is no per-call +// deviceVector allocation and no synchronous host->device memcpy of pointer +// tables: the only sync points are the reductions themselves. +// LatticeView has no default constructor, so the pack holds raw aligned +// storage and views are copied in bytewise (as basisRotateJ does through +// acceleratorPut); the struct is trivially copyable as a kernel argument. +template +struct ViewPack { + alignas(View) unsigned char raw[B*sizeof(View)]; + ComplexD b[B]; + accelerator_inline const View & v(int j) const { return reinterpret_cast(raw)[j]; } + void set(int j,const View &view){ memcpy(raw+j*sizeof(View),&view,sizeof(View)); } +}; + template void rankInnerProductMultiChunk(ComplexD *out,int m, const std::vector*> &left, @@ -822,14 +837,14 @@ void rankInnerProductMultiChunk(ComplexD *out,int m, GridBase *grid = right.Grid(); const uint64_t sites = grid->oSites(); - hostVector h_left_v(m); - deviceVector d_left_v(m); + std::vector h_v; h_v.reserve(m); + ViewPack pack; for(int j=0;jView(AcceleratorRead); + h_v.push_back(left[j]->View(AcceleratorRead)); + pack.set(j,h_v[j]); } - acceleratorCopyToDevice(&h_left_v[0],&d_left_v[0],m*sizeof(View)); - View *left_vp = &d_left_v[0]; + for(int j=m;j partial(sites); batch_t *partial_v = &partial[0]; @@ -839,13 +854,13 @@ void rankInnerProductMultiChunk(ComplexD *out,int m, auto r = right_v[ss]; batch_t acc; for(int j=0;j &out, if ( out.size() ) right.Grid()->GlobalSumVector(&out[0],(int)out.size()); } -// z = z + sum_j b[j] x[j]; if do_norm, returns global |z|^2 from the same pass. -template -RealD axpyMultiNormImpl(Lattice &z,const std::vector &b, - const std::vector*> &x,int do_norm) +// z = z + sum_{j +void axpyMultiChunk(Lattice &z,const ComplexD *b, + const std::vector*> &x,int m, + int do_norm, + decltype(innerProduct(vobj(),vobj())) *inner_tmp_v) { typedef decltype(z.View(AcceleratorRead)) View; - - int m = x.size(); - GRID_ASSERT((int)b.size()>=m); + GRID_ASSERT(m>=1 && m<=B); GridBase *grid = z.Grid(); const uint64_t nsimd = grid->Nsimd(); const uint64_t sites = grid->oSites(); - hostVector h_x_v(std::max(m,1)); - deviceVector d_x_v(std::max(m,1)); - hostVector h_b(std::max(m,1)); - deviceVector d_b(std::max(m,1)); + std::vector h_v; h_v.reserve(m); + ViewPack pack; for(int j=0;jView(AcceleratorRead); - h_b[j] = b[j]; + h_v.push_back(x[j]->View(AcceleratorRead)); + pack.set(j,h_v[j]); + pack.b[j] = b[j]; } - if ( m ) { - acceleratorCopyToDevice(&h_x_v[0],&d_x_v[0],m*sizeof(View)); - acceleratorCopyToDevice(&h_b[0],&d_b[0],m*sizeof(ComplexD)); - } - View *x_vp = &d_x_v[0]; - ComplexD *b_p = &d_b[0]; + for(int j=m;j inner_tmp(do_norm ? sites : 1); - inner_t *inner_tmp_v = &inner_tmp[0]; - accelerator_for(ss,sites,nsimd,{ auto acc = coalescedRead(z_v[ss]); - for(int j=0;j +RealD axpyMultiNormImpl(Lattice &z,const std::vector &b, + const std::vector*> &x,int do_norm) +{ + typedef decltype(innerProduct(vobj(),vobj())) inner_t; + int m = x.size(); + GRID_ASSERT((int)b.size()>=m); + GridBase *grid = z.Grid(); + const uint64_t sites = grid->oSites(); + + deviceVector inner_tmp(do_norm ? sites : 1); + inner_t *inner_tmp_v = &inner_tmp[0]; + + if ( m==0 ) { + if ( do_norm ) return norm2(z); + return 0.0; + } + for(int j0=0;j0=m); + std::vector*> sub(x.begin()+j0,x.begin()+j0+mm); + int dn = do_norm && last; + if ( mm<=2 ) axpyMultiChunk<2> (z,&b[j0],sub,mm,dn,inner_tmp_v); + else if ( mm<=4 ) axpyMultiChunk<4> (z,&b[j0],sub,mm,dn,inner_tmp_v); + else if ( mm<=8 ) axpyMultiChunk<8> (z,&b[j0],sub,mm,dn,inner_tmp_v); + else axpyMultiChunk<16>(z,&b[j0],sub,mm,dn,inner_tmp_v); + } RealD nrm = 0.0; if ( do_norm ) { diff --git a/systems/Frontier/slate_vs_schur2d.job b/systems/Frontier/slate_vs_schur2d.job index 4873f800f..d04a2d330 100644 --- a/systems/Frontier/slate_vs_schur2d.job +++ b/systems/Frontier/slate_vs_schur2d.job @@ -52,15 +52,24 @@ chmod +x ./select_gpu root=$HOME/ParallelIO/systems/Frontier source $root/sourceme-rocm7.2.sh -export OMP_NUM_THREADS=7 +# ONE host thread per rank, for both codes. Grid uses no host OpenMP on the +# GPU build; SLATE with 7 threads issuing device-buffer MPI concurrently +# deadlocks in getrf's tile broadcast under Cray MPICH (slate_debug.job V1 +# hung, V2 with OMP_NUM_THREADS=1 completed, 2026-08-25). Each rank has a +# whole GCD; the parallelism is on the device for both. +export OMP_NUM_THREADS=1 export MPICH_GPU_SUPPORT_ENABLED=1 export MPICH_SMP_SINGLE_COPY_MODE=CMA export MPICH_OFI_NIC_POLICY=GPU -export MPICH_MAX_THREAD_SAFETY=multiple # SLATE calls MPI from OpenMP tasks +export MPICH_MAX_THREAD_SAFETY=multiple # SLATE still requests MPI_THREAD_MULTIPLE module load libfabric BIN=$root/tests/debug/Test_schur2d_vs_slate +# Three legs per stage: Grid 2D Schur, SLATE getrf+getri, SLATE getrf+getrs(I). +# getri is a host loop in SLATE (minutes at N=138240; measured once already): +# S2D_SKIP_GETRI=1 drops it so S3 is ~3 min. Unset to measure it again. +export S2D_SKIP_GETRI=${S2D_SKIP_GETRI:-1} OPTS1="--accelerator-threads 8 --shm 4096 --shm-mpi 1 --device-mem 32000" echo "=========================================================" @@ -80,7 +89,7 @@ S2D_N=4096 srun -N1 -n8 ./select_gpu $BIN --mpi 1.1.2.4 --grid 16.16.16.16 $OPTS 2>&1 | tee s1.out S1RC=${PIPESTATUS[0]} S1CERT=$(grep -c "certificate" s1.out) -echo "S1 exit code $S1RC, certificates printed $S1CERT (expect 2)" +echo "S1 exit code $S1RC, certificates printed $S1CERT (expect 2, or 3 with getri)" ############################################################################## echo "=========================================================" @@ -94,7 +103,7 @@ echo "=========================================================" echo "S3: THE comparison, 36 nodes, 288 ranks, N=138240 (nb=480, grid 16x18)" echo "=========================================================" ############################################################################## -if [ "$S1RC" -eq 0 ] && [ "$S1CERT" -eq 2 ] +if [ "$S1RC" -eq 0 ] && [ "$S1CERT" -ge 2 ] then S2D_N=138240 srun -N36 -n288 ./select_gpu $BIN --mpi 3.6.4.4 --grid 48.48.48.96 $OPTS1 echo "S3 exit code $?" @@ -105,6 +114,6 @@ fi echo "=========================================================" echo "summary (both legs, all stages)" echo "=========================================================" -grep -h -E "Grid-vs-SLATE|GRID :|SLATE :" slurm-$SLURM_JOB_ID.out 2>/dev/null +grep -h -E "Grid-vs-SLATE|GRID :|SLATE :|SLATE-getrs" slurm-$SLURM_JOB_ID.out 2>/dev/null grep -h "git commit hash" slurm-$SLURM_JOB_ID.out 2>/dev/null | sort -u echo "=========================================================" diff --git a/tests/debug/Test_schur2d_vs_slate.cc b/tests/debug/Test_schur2d_vs_slate.cc index df720d7cd..da3175450 100644 --- a/tests/debug/Test_schur2d_vs_slate.cc +++ b/tests/debug/Test_schur2d_vs_slate.cc @@ -55,7 +55,19 @@ Author: Peter Boyle // grids, nb dividing N (48|720) and ragged (N=730, nb=50); both legs certify // max|A.Ainv-I| ~ 1e-15. // +// Frontier (HIP, Target::Devices, Cray MPICH, 2026-08-25): run with +// OMP_NUM_THREADS=1. With 7 threads per rank SLATE's getrf deadlocks in its +// tile broadcast (concurrent device-buffer MPI from OpenMP tasks); with one +// thread it completes (8 GCDs, N=4096: certificate 2.7e-15). One host +// thread per rank is the like-for-like anyway: Grid's GPU build uses none. +// Also: the site `slate` module (cpu env) must NOT be loaded -- its host-only +// libblaspp shadows the ROCm one via LD_LIBRARY_PATH and throws +// "device BLAS not available" from host_malloc_pinned. +// // S2D_N, S2D_NB as in Test_schur2d_scale (default nb = N/P). +// S2D_SKIP_GETRI=1 skips the getri leg (host loop; ~4 min at N=138240). +// S2D_NOWARM=1 skips the warm-up. +// A third leg, getrf+getrs(I), is SLATE's device-resident inverse route. ////////////////////////////////////////////////////////////////////////////// #include @@ -153,13 +165,22 @@ int main(int argc, char **argv) // a hang shows WHERE even through block-buffered stdout. auto Stage = [&](const char *s){ std::cout << GridLogMessage << "stage: " << s << std::endl << std::flush; }; if ( !getenv("S2D_NOWARM") ) { - int64_t Nw = 8*P; int64_t nbw = 8; - std::vector rs(P+1); for(int r=0;r<=P;r++) rs[r]=8*r; - std::vector hw(8*Nw); for(int64_t jj=0;jj a 122 s SLATE + // warm-up dominated by 288-way tile broadcasts.) Ranks beyond the first + // Nw/nbw own nothing in the rows layout, which RowsToCyclic handles. + int64_t Nw = 64; int64_t nbw = 8; + // same partition rule as the main leg: Nw/P rows each, remainder to the + // first ranks; at P > Nw most ranks contribute zero rows. + std::vector rs(P+1); rs[0]=0; + for(int r=0;r hw(std::max(rw,1)*Nw); + for(int64_t jj=0;jj dw(hw.size()); acceleratorCopyToDevice(&hw[0],&dw[0],hw.size()*sizeof(ComplexD)); BlockCyclicMatrix W(grid,Nw,nbw,Pr,Pc); Stage("warm-up Grid RowsToCyclic"); - BlockCyclicRedistribute::RowsToCyclic(grid,rs,&dw[0],8,W); + BlockCyclicRedistribute::RowsToCyclic(grid,rs,&dw[0],rw,W); Stage("warm-up Grid Invert"); BlockCyclicSchurInverse RSIw; RSIw.Invert(W); #ifdef HAVE_SLATE @@ -203,7 +224,7 @@ int main(int argc, char **argv) acceleratorCopyDeviceToDevice((void *)&A.data[0],(void *)&A0.data[0],A.data.size()*sizeof(ComplexD)); double t2=usecond(); Stage("Grid Invert"); - RSI2.Invert(A); + { GRID_TRACE("GridInvert"); RSI2.Invert(A); } double t3=usecond(); BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&rows1d[0],myrows); double t4=usecond(); @@ -218,7 +239,7 @@ int main(int argc, char **argv) // LEG 2: SLATE, every layout step timed and charged. //////////////////////////////////////////////////////////////////////// #ifdef HAVE_SLATE - { + if ( !getenv("S2D_SKIP_GETRI") ) { // S2D_SKIP_GETRI=1: getri is a host loop, minutes at N=138240 typedef std::complex scalar_t; acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD)); BlockCyclicMatrix A(grid,N,nb,Pr,Pc), A0(grid,N,nb,Pr,Pc); @@ -262,10 +283,10 @@ int main(int argc, char **argv) slate::Pivots pivots; Stage("SLATE getrf"); double t4=usecond(); - slate::getrf(S, pivots, opts); // LU, partial pivoting + { GRID_TRACE("SLATE_getrf"); slate::getrf(S, pivots, opts); } // LU, partial pivoting double t5=usecond(); Stage("SLATE getri"); - slate::getri(S, pivots, opts); // in-place inverse from the factor + { GRID_TRACE("SLATE_getri"); slate::getri(S, pivots, opts); } // in-place inverse from the factor double t6=usecond(); // back onto the device, in our layout (getri applies the pivots itself) @@ -281,6 +302,72 @@ int main(int argc, char **argv) << " TOTAL " << ((t1-t0)+(t8-t2))/1e6 << " s" << " certificate " << cert << std::endl; } + + //////////////////////////////////////////////////////////////////////// + // LEG 3: SLATE getrf + getrs(I) -- SLATE's device-resident route to an + // explicit inverse. getri's L-side loop is hard-coded Target::HostTask + // (src/getri.cc: copy/gemmA/trsm/permuteRows all ), so under + // Target::Devices it runs as a 288-step host loop on one thread per rank. + // getrs is two trsm's that honour the target (plus host row permutes). + // The identity RHS is built on the host in the same ScaLAPACK layout and + // its construction is on SLATE's clock. + //////////////////////////////////////////////////////////////////////// + { + typedef std::complex scalar_t; + acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD)); + BlockCyclicMatrix A(grid,N,nb,Pr,Pc), A0(grid,N,nb,Pr,Pc); + double t0=usecond(); + BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,A); + double t1=usecond(); + if ( A.data.size() ) + acceleratorCopyDeviceToDevice((void *)&A.data[0],(void *)&A0.data[0],A.data.size()*sizeof(ComplexD)); + BlockCyclicLayout &L = A.layout; + uint64_t nloc = (uint64_t)L.mloc*L.nloc; + std::vector hA(nloc ? nloc : 1), hB(nloc ? nloc : 1); + double t2=usecond(); + if ( nloc ) acceleratorCopyFromDevice(&A.data[0], (void *)&hA[0], nloc*sizeof(ComplexD)); + for(int64_t lj=0;lj::fromScaLAPACK(N, N, &hA[0], (int64_t)std::max(L.mloc,1), + nb, nb, slate::GridOrder::Row, Pr, Pc, grid->communicator); + auto B = slate::Matrix::fromScaLAPACK(N, N, &hB[0], (int64_t)std::max(L.mloc,1), + nb, nb, slate::GridOrder::Row, Pr, Pc, grid->communicator); +#if defined(GRID_HIP) || defined(GRID_CUDA) || defined(GRID_SYCL) + slate::Target target = slate::Target::Devices; +#else + slate::Target target = slate::Target::HostTask; +#endif + slate::Options opts = { + { slate::Option::Target, target }, + { slate::Option::Lookahead, 1 }, + { slate::Option::InnerBlocking, 16 }, + }; + slate::Pivots pivots; + Stage("SLATE getrf (getrs leg)"); + double t4=usecond(); + { GRID_TRACE("SLATE_getrf"); slate::getrf(S, pivots, opts); } + double t5=usecond(); + Stage("SLATE getrs(I)"); + { GRID_TRACE("SLATE_getrs"); slate::getrs(S, pivots, B, opts); } // B <- A^{-1} I + double t6=usecond(); + if ( nloc ) acceleratorCopyToDevice((void *)&hB[0], &A.data[0], nloc*sizeof(ComplexD)); + double t7=usecond(); + BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&rows1d[0],myrows); + double t8=usecond(); + double cert = Certify(grid,A0,A,nb,Pr,Pc); + std::cout << GridLogMessage << "SLATE-getrs : redist->2D " << (t1-t0)/1e6 + << " D2H+I " << (t3-t2)/1e6 << " wrap " << (t4-t3)/1e6 + << " getrf " << (t5-t4)/1e6 << " getrs " << (t6-t5)/1e6 + << " H2D " << (t7-t6)/1e6 << " redist->rows " << (t8-t7)/1e6 + << " TOTAL " << ((t1-t0)+(t8-t2))/1e6 << " s" + << " certificate " << cert << std::endl; + } #else std::cout << GridLogMessage << "SLATE : leg not built (compile with -DHAVE_SLATE and link -lslate -lblaspp -llapackpp)" << std::endl; #endif