From 167753aa3f77168e73132c5d047a146148c3e087 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 25 Aug 2026 12:07:11 -0400 Subject: [PATCH] Adding ring all reduce for large vectors --- Grid/GridCore.h | 3 +- Grid/algorithms/multigrid/DenseCoarseMatrix.h | 16 +- Grid/threads/Threads.h | 5 +- tests/debug/Test_multi_reduction.cc | 117 +++++++++ tests/debug/Test_pgcr_history.cc | 109 ++++++++ tests/debug/Test_ring_allreduce.cc | 122 +++++++++ tests/debug/Test_schur2d_vs_slate.cc | 241 ++++++++++++++++++ 7 files changed, 607 insertions(+), 6 deletions(-) create mode 100644 tests/debug/Test_multi_reduction.cc create mode 100644 tests/debug/Test_pgcr_history.cc create mode 100644 tests/debug/Test_ring_allreduce.cc create mode 100644 tests/debug/Test_schur2d_vs_slate.cc diff --git a/Grid/GridCore.h b/Grid/GridCore.h index 90d01f2e4..51bcecd3a 100644 --- a/Grid/GridCore.h +++ b/Grid/GridCore.h @@ -54,7 +54,8 @@ Author: paboyle #include #include #include -#include +#include +#include #include #include #include diff --git a/Grid/algorithms/multigrid/DenseCoarseMatrix.h b/Grid/algorithms/multigrid/DenseCoarseMatrix.h index f581f71b4..11ee52388 100644 --- a/Grid/algorithms/multigrid/DenseCoarseMatrix.h +++ b/Grid/algorithms/multigrid/DenseCoarseMatrix.h @@ -251,11 +251,13 @@ public: for(int j=0;j=0 && devSum<=3); std::cout << GridLogMessage << "DenseCoarseMatrix: slab resident on device (" << sbytes/1024./1024. << " MB/rank), split-K NK=" << NK << " (Kc=" << Kc << "); " - << (devSum ? "DEVICE-buffer allreduce (GPU-aware MPI)" : "host allreduce") - << std::endl; + << sumName[devSum] << std::endl; } //////////////////////////////////////////////////////////////////// @@ -961,7 +963,13 @@ public: } t2 = usecond(); { GRID_TRACE("DenseAllreduce"); - grid->GlobalSumVector((ComplexF *)&dX[0], (int)nX); + // DENSE_DEVICE_SUM=1 : device-buffer MPI_Allreduce (Cray MPICH aborts + // above ~8 MB: 12 RHS at N=138240 is 13.3 MB) + // DENSE_DEVICE_SUM=2 : CartesianRingAllReduce, P2P only, no size cliff + // DENSE_DEVICE_SUM=3 : flat RingAllReduce, P2P only + if (devSum==2) CartesianRingAllReduce(grid,(ComplexF *)&dX[0],nX); + else if (devSum==3) RingAllReduce(grid,(ComplexF *)&dX[0],nX); + else grid->GlobalSumVector((ComplexF *)&dX[0], (int)nX); } t3 = usecond(); } else { diff --git a/Grid/threads/Threads.h b/Grid/threads/Threads.h index c517ca13c..9d767c48d 100644 --- a/Grid/threads/Threads.h +++ b/Grid/threads/Threads.h @@ -77,11 +77,14 @@ inline void thread_bcopy(const void *from, void *to,size_t bytes) { const uint64_t *ufrom = (const uint64_t *)from; uint64_t *uto = (uint64_t *)to; - GRID_ASSERT(bytes%8==0); uint64_t words=bytes/8; thread_for(w,words,{ uto[w] = ufrom[w]; }); + // Tail: byte counts that are not a multiple of 8 (e.g. odd-length float + // buffers) copy the remainder serially instead of asserting. + uint64_t tail = bytes%8; + if ( tail ) bcopy((const char *)from + words*8, (char *)to + words*8, tail); } #else inline void thread_bcopy(const void *from, void *to,size_t bytes) diff --git a/tests/debug/Test_multi_reduction.cc b/tests/debug/Test_multi_reduction.cc new file mode 100644 index 000000000..973ae9809 --- /dev/null +++ b/tests/debug/Test_multi_reduction.cc @@ -0,0 +1,117 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/debug/Test_multi_reduction.cc + + Copyright (C) 2026 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////// +// Gate for the batched Krylov linear algebra in Lattice_reduction.h: +// innerProductMulti vs innerProduct, one j at a time +// axpyMulti vs sequential axpy +// axpyMultiNorm vs sequential axpy + norm2 +// on TWO field types through the SAME code path: +// fine : LatticeFermionD on the 4D grid (vSpinColourVector) +// coarse: Lattice> on a 5D grid with the extra +// dimension standing in for nrhs, as the coarse mRHS solver does. +// Batch sizes cover every template width and the >16 chunking path. +// +// mpirun -n 1 ./Test_multi_reduction --grid 8.8.8.8 --mpi 1.1.1.1 +////////////////////////////////////////////////////////////////////////////// + +#include + +using namespace Grid; + +static int failures = 0; +static void Report(const std::string &name, bool pass, const std::string &detail="") +{ + std::cout << GridLogMessage << " " << name << (pass ? " PASS" : " ** FAIL **"); + if ( detail.size() ) std::cout << " " << detail; + std::cout << std::endl; + if ( !pass ) failures++; +} +static double Cabs(const ComplexD &z){ double re=z.real(), im=z.imag(); return std::sqrt(re*re+im*im); } + +template +void Check(const std::string &tag, GridBase *grid, GridParallelRNG &RNG, const std::vector &batches) +{ + const int Mmax = 20; + std::vector x; for(int j=0;j win(m); + std::vector b(m); + for(int j=0;j ip; + innerProductMulti(ip,win,r); + double worst=0.0; + for(int j=0;j seeds({1,2,3,4}); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds); + GridParallelRNG RNG5(CGrid); RNG5.SeedFixedIntegers(seeds); + + std::vector batches({1,2,3,4,5,8,9,16,17,20}); + + Check("fine ", UGrid, RNG4, batches); + + const int nbasis = 8; + typedef iVector>>,nbasis> CoarseSiteVector; + typedef Lattice CoarseField; + Check("coarse", CGrid, RNG5, batches); + + std::cout << GridLogMessage << (failures ? "Test_multi_reduction: FAILURES" + : "Test_multi_reduction: ALL PASS") << std::endl; + Grid_finalize(); + return failures ? 1 : 0; +} diff --git a/tests/debug/Test_pgcr_history.cc b/tests/debug/Test_pgcr_history.cc new file mode 100644 index 000000000..8adc47c43 --- /dev/null +++ b/tests/debug/Test_pgcr_history.cc @@ -0,0 +1,109 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/debug/Test_pgcr_history.cc + + Copyright (C) 2026 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////// +// Regression gate for the PERSISTENT history arrays in +// PrecGeneralisedConjugateResidualNonHermitian (q, p, qq allocated once per +// grid instead of per GCRnStep call). +// +// mpirun -n 1 ./Test_pgcr_history --grid 8.8.8.8 --mpi 1.1.1.1 +// +// Persistence is only safe if no stale history is ever read. The test: +// 1. solve src1 -> x1 (fresh allocation) +// 2. solve src2 -> y (DIRTIES the persistent history) +// 3. solve src1 -> x2 (must reuse arrays, must not see src2) +// T1 : x1 == x2 BITWISE -- any consumption of stale history breaks this +// T2 : true residual of x1 below tolerance -- the solver still solves +// T3 : a second solver object with a different mmax on the same grid +// (exercises the (re)allocation path) also converges. +// +// Wilson at mass 0.5 on a hot configuration: the Hermitian part is positive +// definite, so unpreconditioned GCR converges without breakdown. +////////////////////////////////////////////////////////////////////////////// + +#include + +using namespace Grid; + +static int failures = 0; +static void Report(const std::string &name, bool pass, const std::string &detail="") +{ + std::cout << GridLogMessage << " " << name << (pass ? " PASS" : " ** FAIL **"); + if ( detail.size() ) std::cout << " " << detail; + std::cout << std::endl; + if ( !pass ) failures++; +} + +int main(int argc, char **argv) +{ + Grid_init(&argc, &argv); + + GridCartesian *UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + + std::vector seeds({1,2,3,4}); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds); + + LatticeGaugeFieldD Umu(UGrid); + SU::HotConfiguration(RNG4, Umu); + + RealD mass = 0.5; + WilsonFermionD Dw(Umu, *UGrid, *UrbGrid, mass); + NonHermitianLinearOperator Op(Dw); + TrivialPrecon simple; + + PrecGeneralisedConjugateResidualNonHermitian PGCR(1.0e-10, 2000, Op, simple, 4, 8); + PGCR.SetZeroGuess(1); + + LatticeFermionD src1(UGrid), src2(UGrid), x1(UGrid), x2(UGrid), y(UGrid), r(UGrid); + gaussian(RNG4, src1); + gaussian(RNG4, src2); + + x1 = Zero(); PGCR(src1, x1); + y = Zero(); PGCR(src2, y); // dirty the history with an unrelated solve + x2 = Zero(); PGCR(src1, x2); + + { + r = x1 - x2; + RealD dd = norm2(r); + Report("T1 repeat solve after a different source is BITWISE identical", dd == 0.0, + "|x1-x2|^2 = "+std::to_string(dd)); + } + { + Op.Op(x1, r); r = r - src1; + RealD res = std::sqrt(norm2(r)/norm2(src1)); + Report("T2 true residual of the solve", res < 1.0e-9, "rel "+std::to_string(res)); + } + { + PrecGeneralisedConjugateResidualNonHermitian PGCR2(1.0e-10, 2000, Op, simple, 2, 8); + PGCR2.SetZeroGuess(1); + PGCR2.LogCoefficients(1); // exercise the coefficient log format + LatticeFermionD x3(UGrid); x3 = Zero(); PGCR2(src1, x3); + Op.Op(x3, r); r = r - src1; + RealD res = std::sqrt(norm2(r)/norm2(src1)); + Report("T3 second solver object, different mmax, converges", res < 1.0e-9, "rel "+std::to_string(res)); + } + + std::cout << GridLogMessage << (failures ? "Test_pgcr_history: FAILURES" + : "Test_pgcr_history: ALL PASS") << std::endl; + Grid_finalize(); + return failures ? 1 : 0; +} diff --git a/tests/debug/Test_ring_allreduce.cc b/tests/debug/Test_ring_allreduce.cc new file mode 100644 index 000000000..ff2e5ab1d --- /dev/null +++ b/tests/debug/Test_ring_allreduce.cc @@ -0,0 +1,122 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/debug/Test_ring_allreduce.cc + + Copyright (C) 2026 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////// +// Gate for RingAllReduce / CartesianRingAllReduce (P2P-only vector +// all-reduce) against GlobalSumVector (MPI_Allreduce): +// T1 flat ring == GlobalSumVector, RealD/ComplexD/ComplexF/RealF, +// n in {1, P-1, P, P+1, 1000003, 2^20} (n + +using namespace Grid; + +static int failures = 0; +static void Report(const std::string &name, bool pass, const std::string &detail="") +{ + std::cout << GridLogMessage << " " << name << (pass ? " PASS" : " ** FAIL **"); + if ( detail.size() ) std::cout << " " << detail; + std::cout << std::endl; + if ( !pass ) failures++; +} + +template double Mag(const T &x){ return std::fabs((double)x); } +template<> double Mag(const ComplexD &z){ return std::sqrt(z.real()*z.real()+z.imag()*z.imag()); } +template<> double Mag(const ComplexF &z){ return std::sqrt((double)z.real()*z.real()+(double)z.imag()*z.imag()); } + +template T Fill(uint64_t i, int rank){ return T( 0.5*std::sin(0.01*i + 0.7*rank) + 1.0e-3*rank ); } +template<> ComplexD Fill(uint64_t i,int rank){ return ComplexD(0.5*std::sin(0.01*i+0.7*rank), 0.3*std::cos(0.02*i-0.1*rank)); } +template<> ComplexF Fill(uint64_t i,int rank){ return ComplexF(0.5*std::sin(0.01*i+0.7*rank), 0.3*std::cos(0.02*i-0.1*rank)); } + +template +void Check(const std::string &tname, GridCartesian *grid, double tol) +{ + int P = grid->ProcessorCount(); + int me = grid->ThisRank(); + std::vector sizes({1,(uint64_t)std::max(P-1,1),(uint64_t)P,(uint64_t)P+1,1000003,1<<20}); + for(uint64_t n : sizes){ + std::vector h(n); for(uint64_t i=0;i(i,me); + std::vector ref(h); grid->GlobalSumVector(&ref[0],(int)n); + + deviceVector d(n); + for(int variant=0;variant<2;variant++){ + acceleratorCopyToDevice(&h[0],&d[0],n*sizeof(T)); + if(variant==0) RingAllReduce(grid,&d[0],n); + else CartesianRingAllReduce(grid,&d[0],n); + std::vector out(n); acceleratorCopyFromDevice(&d[0],&out[0],n*sizeof(T)); + double worst=0.0, scale=0.0; + for(uint64_t i=0;i(out[i]-ref[i])); scale=std::max(scale,Mag(ref[i])); } + RealD w=worst; grid->GlobalMax(w); + std::ostringstream os; os<<"n="< h(n); for(uint64_t i=0;i(i,me); + deviceVector d(n); std::vector a(n),b(n); + acceleratorCopyToDevice(&h[0],&d[0],n*sizeof(T)); RingAllReduce(grid,&d[0],n); acceleratorCopyFromDevice(&d[0],&a[0],n*sizeof(T)); + acceleratorCopyToDevice(&h[0],&d[0],n*sizeof(T)); RingAllReduce(grid,&d[0],n); acceleratorCopyFromDevice(&d[0],&b[0],n*sizeof(T)); + RealD diff = (memcmp(&a[0],&b[0],n*sizeof(T))!=0) ? 1.0 : 0.0; + grid->GlobalSum(diff); + Report("T3 bitwise repeat "+tname, diff==0.0); + } +} + +int main(int argc, char **argv) +{ + Grid_init(&argc, &argv); + GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi()); + std::cout << GridLogMessage << "Ring allreduce test: P=" << grid->ProcessorCount() + << " processor grid " << grid->_processors << std::endl; + + Check ("RealD ", grid, 1.0e-13); + Check("ComplexD", grid, 1.0e-13); + Check ("RealF ", grid, 1.0e-5); + Check("ComplexF", grid, 1.0e-5); + + // T4 timing at 16 MB of ComplexF (the dense-apply size at 12 RHS is 13.3 MB) + { + uint64_t n = 2*1024*1024; + deviceVector d(n); std::vector h(n,ComplexF(1.0,0.0)); + for(int rep=0;rep<2;rep++){ + acceleratorCopyToDevice(&h[0],&d[0],n*sizeof(ComplexF)); + double t0=usecond(); RingAllReduce(grid,&d[0],n); double t1=usecond(); + acceleratorCopyToDevice(&h[0],&d[0],n*sizeof(ComplexF)); + double t2=usecond(); CartesianRingAllReduce(grid,&d[0],n); double t3=usecond(); + double t4=usecond(); grid->GlobalSumVector(&h[0],(int)n); double t5=usecond(); + if(rep) std::cout << GridLogMessage << "T4 timing 16 MB ComplexF: flat ring " << (t1-t0)/1000. + << " ms, cartesian ring " << (t3-t2)/1000. << " ms, GlobalSumVector(host) " << (t5-t4)/1000. << " ms" << std::endl; + } + } + + std::cout << GridLogMessage << (failures ? "Test_ring_allreduce: FAILURES" : "Test_ring_allreduce: ALL PASS") << std::endl; + Grid_finalize(); + return failures ? 1 : 0; +} diff --git a/tests/debug/Test_schur2d_vs_slate.cc b/tests/debug/Test_schur2d_vs_slate.cc new file mode 100644 index 000000000..d76d2d366 --- /dev/null +++ b/tests/debug/Test_schur2d_vs_slate.cc @@ -0,0 +1,241 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/debug/Test_schur2d_vs_slate.cc + + Copyright (C) 2026 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////// +// LIKE-FOR-LIKE: Grid's 2D block-cyclic Schur inverse versus SLATE's +// getrf+getri, on the SAME synthetic matrix, same ranks, same job, with EVERY +// layout cost on the clock. +// +// Grid : rows -> RowsToCyclic -> BlockCyclicSchurInverse -> CyclicToRows +// SLATE : rows -> RowsToCyclic -> D2H -> fromScaLAPACK -> getrf -> getri +// -> H2D -> CyclicToRows +// +// Both legs pay the identical rows<->2D redistribution (the production import +// produces rank-major rows). SLATE's EXTRA layout cost is the host round +// trip: fromScaLAPACK wraps host memory in exactly the layout BlockCyclicLayout +// already uses (csrc=0 cyclic, column-major local, ld=mloc, mb=nb, row-major +// process grid), so no other relayout exists to charge. +// +// Both inverses are certified by the SAME instrument: Grid's SUMMA product +// against an untouched copy, max |A.Ainv - I| checked locally, one GlobalMax. +// +// Precision inventory (both): fp64 throughout. Algorithms differ: SLATE is +// LU with partial pivoting then getri; Grid is pivot-free recursive Schur. +// +// Build: needs SLATE (spack load slate) and +// CXXFLAGS += -DGRID_HAVE_SLATE -I$SLATE_ROOT/include +// LDFLAGS += -L$SLATE_ROOT/lib -lslate -lblaspp -llapackpp +// Without GRID_HAVE_SLATE the SLATE leg reports itself as not built and the +// Grid leg still runs, so the binary is always usable. +// +// MPI threading: the SLATE build initialises MPI at MPI_THREAD_MULTIPLE before +// Grid_init (see main); Grid's own SERIALIZED level gives Bus errors in +// slate::listBcast with >1 OpenMP thread. On an oversubscribed laptop set +// OMP_WAIT_POLICY=passive and OPENBLAS_NUM_THREADS=1, or SLATE's spinning +// tasks and OpenMPI's polling livelock each other (0.03 s -> minutes). +// +// Validated on the laptop (CPU, HostTask) 2026-08-25: 1x2 and 2x2 process +// grids, nb dividing N (48|720) and ragged (N=730, nb=50); both legs certify +// max|A.Ainv-I| ~ 1e-15. +// +// S2D_N, S2D_NB as in Test_schur2d_scale (default nb = N/P). +////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#ifdef GRID_HAVE_SLATE +#include +#endif + +using namespace Grid; + +static double Cabs(const ComplexD &z){ double re=z.real(), im=z.imag(); return std::sqrt(re*re+im*im); } + +static ComplexD Fill(int64_t i, int64_t j) +{ + double x = std::sin(0.7*i + 1.3*j); + double y = std::cos(1.9*i - 0.4*j); + if ( i==j ) return ComplexD(3.0*64 + x, 0.5); + if ( std::fabs((double)(i-j)) > 64.0 ) return ComplexD(0.0,0.0); + return ComplexD(x,y); +} + +// max |A0 . Ainv - I| over my local elements, reduced once. +static double Certify(GridBase *grid, BlockCyclicMatrix &A0, BlockCyclicMatrix &Ainv, int64_t nb, int Pr, int Pc) +{ + BlockCyclicMatrix Cert(grid, A0.layout.N, nb, Pr, Pc); + BlockCyclicSumma SUMMA; + int64_t N = A0.layout.N; + SUMMA.Multiply(ComplexD(1.0,0.0),A0,Ainv,ComplexD(0.0,0.0),Cert, 0,N,0,N,0,N); + BlockCyclicLayout &L = Cert.layout; + std::vector hc((uint64_t)std::max(L.mloc*L.nloc,1)); + if ( L.mloc*L.nloc ) + acceleratorCopyFromDevice(&Cert.data[0], &hc[0], (uint64_t)L.mloc*L.nloc*sizeof(ComplexD)); + double mx = 0.0; + for(int64_t lj=0;ljGlobalMax(gmx); + return gmx; +} + +int main(int argc, char **argv) +{ +#ifdef GRID_HAVE_SLATE + // SLATE issues MPI calls from concurrent OpenMP tasks and needs + // MPI_THREAD_MULTIPLE. Grid (GridStd.h #undef GRID_COMMS_THREADS) only asks + // for MPI_THREAD_SERIALIZED, which reproducibly gives Bus errors inside + // slate::BaseMatrix::listBcast with >1 OpenMP thread. Grid_init honours an + // already-initialised MPI, so initialise it here at the level SLATE needs. + { + int provided = 0; + MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &provided); + if ( provided != MPI_THREAD_MULTIPLE ) { + fprintf(stderr, "MPI_THREAD_MULTIPLE not provided (got %d); SLATE leg unsafe\n", provided); + GRID_ASSERT(provided == MPI_THREAD_MULTIPLE); + } + } +#endif + Grid_init(&argc, &argv); + GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd, vComplexD::Nsimd()), + GridDefaultMpi()); + const int P = grid->ProcessorCount(); + const int me = grid->ThisRank(); + + int64_t N = getenv("S2D_N") ? atol(getenv("S2D_N")) : 720; + int64_t nb = getenv("S2D_NB") ? atol(getenv("S2D_NB")) : ( (N%P==0) ? N/P : 48 ); + int Pr,Pc; BlockCyclicLayout::ChooseProcessGrid(P,Pr,Pc); + + std::vector rowStart(P+1); rowStart[0]=0; + for(int r=0;r h((uint64_t)std::max(myrows,1)*N); + thread_for(jj, N, { for(int64_t i=0;i rows1d(h.size()); + + //////////////////////////////////////////////////////////////////////// + // LEG 1: Grid + //////////////////////////////////////////////////////////////////////// + { + acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD)); + BlockCyclicMatrix A(grid,N,nb,Pr,Pc), A0(grid,N,nb,Pr,Pc); + BlockCyclicSchurInverse RSI2; + 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)); + double t2=usecond(); + RSI2.Invert(A); + double t3=usecond(); + BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&rows1d[0],myrows); + double t4=usecond(); + double cert = Certify(grid,A0,A,nb,Pr,Pc); + std::cout << GridLogMessage << "GRID : redist->2D " << (t1-t0)/1e6 + << " invert " << (t3-t2)/1e6 << " redist->rows " << (t4-t3)/1e6 + << " TOTAL " << ((t1-t0)+(t3-t2)+(t4-t3))/1e6 << " s" + << " certificate " << cert << std::endl; + } + + //////////////////////////////////////////////////////////////////////// + // LEG 2: SLATE, every layout step timed and charged. + //////////////////////////////////////////////////////////////////////// +#ifdef GRID_HAVE_SLATE + { + 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); // same as Grid leg + double t1=usecond(); + if ( A.data.size() ) + acceleratorCopyDeviceToDevice((void *)&A.data[0],(void *)&A0.data[0],A.data.size()*sizeof(ComplexD)); + + // SLATE's extra layout cost: host copy in the ScaLAPACK layout we already hold. + BlockCyclicLayout &L = A.layout; + uint64_t nloc = (uint64_t)L.mloc*L.nloc; + std::vector hA(nloc ? nloc : 1); + double t2=usecond(); + if ( nloc ) acceleratorCopyFromDevice(&A.data[0], (void *)&hA[0], nloc*sizeof(ComplexD)); + double t3=usecond(); + + // Wrap: same block size, csrc=0 cyclic ownership, column-major local + // storage with lld = mloc, ROW-major process grid (rank = p*Pc + q), + // exactly BlockCyclicLayout's conventions. No further relayout exists. + // + // The communicator MUST be Grid's cartesian one, not MPI_COMM_WORLD: + // BlockCyclicLayout numbers processes by grid->ThisRank(), and Grid's + // OptimalCommunicator permutes ranks relative to the world communicator, + // so under MPI_COMM_WORLD SLATE and Grid disagree on tile ownership + // (observed: Bus error inside listBcast on a 2x2 grid). + auto S = slate::Matrix::fromScaLAPACK(N, N, &hA[0], (int64_t)std::max(L.mloc,1), + nb, nb, slate::GridOrder::Row, Pr, Pc, grid->communicator); + // Target follows the Grid build: devices on GPU builds, host tasks on a + // CPU build (laptop validation of the SLATE calls at small N). +#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; + double t4=usecond(); + slate::getrf(S, pivots, opts); // LU, partial pivoting + double t5=usecond(); + 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) + if ( nloc ) acceleratorCopyToDevice((void *)&hA[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 : redist->2D " << (t1-t0)/1e6 + << " D2H " << (t3-t2)/1e6 << " wrap " << (t4-t3)/1e6 + << " getrf " << (t5-t4)/1e6 << " getri " << (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 -DGRID_HAVE_SLATE and link -lslate -lblaspp -llapackpp)" << std::endl; +#endif + + Grid_finalize(); + return 0; +}