Changes to dense inverse AND improved timing on comms

This commit is contained in:
Peter Boyle
2026-08-22 11:43:13 -04:00
parent 152e9ecf3b
commit 1e287b6b30
8 changed files with 274 additions and 52 deletions
+4 -3
View File
@@ -148,14 +148,15 @@ public:
#endif
#ifdef GRID_HIP
std::cout << "hipblasCreate"<<std::endl;
hipblasCreate(&gridblasHandle);
GRID_ASSERT( hipblasCreate(&gridblasHandle) == HIPBLAS_STATUS_SUCCESS );
// Explicit HOST mode: the hipBLAS default is UNDOCUMENTED in the
// headers (enum 0 == HOST by cuBLAS-mirroring convention only);
// set it and print it so every log carries the ground truth.
hipblasSetPointerMode(gridblasHandle, HIPBLAS_POINTER_MODE_HOST);
GRID_ASSERT( hipblasSetPointerMode(gridblasHandle, HIPBLAS_POINTER_MODE_HOST)
== HIPBLAS_STATUS_SUCCESS );
{
hipblasPointerMode_t pm;
hipblasGetPointerMode(gridblasHandle,&pm);
GRID_ASSERT( hipblasGetPointerMode(gridblasHandle,&pm) == HIPBLAS_STATUS_SUCCESS );
std::cout << "GridBLAS: hipBLAS pointer mode "
<< ((pm==HIPBLAS_POINTER_MODE_DEVICE)?"DEVICE":"HOST") <<std::endl;
}
@@ -155,6 +155,7 @@ public:
}
void ImportFineGridVectors(std::vector <Field > &vecs, deviceVector<scalar> &blas)
{
GRID_TRACE("ImportFineGridVectors");
int nvec = vecs.size();
typedef typename Field::vector_object vobj;
// std::cout << GridLogMessage <<" BlockProjector importing "<<nvec<< " fine grid vectors" <<std::endl;
@@ -547,6 +548,7 @@ public:
void ExportFineGridVectors(std::vector <Field> &vecs, deviceVector<scalar> &blas)
{
GRID_TRACE("ExportFineGridVectors");
typedef typename Field::vector_object vobj;
int nvec = vecs.size();
@@ -622,6 +624,7 @@ public:
template<class vobj>
void ImportCoarseGridVectors(std::vector <Lattice<vobj> > &vecs, deviceVector<scalar> &blas)
{
GRID_TRACE("ImportCoarseGridVectors");
int nvec = vecs.size();
typedef typename vobj::scalar_object coarse_scalar_object;
@@ -688,6 +691,7 @@ public:
template<class vobj>
void ExportCoarseGridVectors(std::vector <Lattice<vobj> > &vecs, deviceVector<scalar> &blas)
{
GRID_TRACE("ExportCoarseGridVectors");
int nvec = vecs.size();
typedef typename vobj::scalar_object coarse_scalar_object;
// std::cout << GridLogMessage<<" BlockProjector exporting "<<nvec<< " coarse grid vectors" <<std::endl;
@@ -775,17 +779,20 @@ public:
deviceVector<scalar *> Cd(coarse_vol);
// std::cout << "BlockProject pointers"<<std::endl;
for(int c=0;c<coarse_vol;c++){
// BLAS_V[coarse_vol][nbasis][block_vol][words]
// BLAS_F[coarse_vol][nrhs][block_vol][words]
// BLAS_C[coarse_vol][nrhs][nbasis]
scalar * Vh = & BLAS_V[c*nbasis*block_vol*words];
scalar * Fh = & BLAS_F[c*nrhs*block_vol*words];
scalar * Ch = & BLAS_C[c*nrhs*nbasis];
acceleratorPut(Vd[c],Vh);
acceleratorPut(Fd[c],Fh);
acceleratorPut(Cd[c],Ch);
// ONE bulk transfer per table. acceleratorPut is a *synchronous* 8-byte
// hipMemcpy, so the elementwise form emitted 3*coarse_vol of them per call:
// a traced run showed 272k hipMemcpy calls costing 4.5 s of API time to move
// 0.7 s worth of bytes. Same fix as BatchedBlas.h's staging rewrite.
if ( coarse_vol ) {
std::vector<scalar *> hVd(coarse_vol), hFd(coarse_vol), hCd(coarse_vol);
for(int c=0;c<coarse_vol;c++){
hVd[c] = & BLAS_V[c*nbasis*block_vol*words];
hFd[c] = & BLAS_F[c*nrhs*block_vol*words];
hCd[c] = & BLAS_C[c*nrhs*nbasis];
}
acceleratorCopyToDevice(&hVd[0],&Vd[0],coarse_vol*sizeof(scalar *));
acceleratorCopyToDevice(&hFd[0],&Fd[0],coarse_vol*sizeof(scalar *));
acceleratorCopyToDevice(&hCd[0],&Cd[0],coarse_vol*sizeof(scalar *));
}
GridBLAS BLAS;
@@ -827,16 +834,20 @@ public:
deviceVector<scalar *> Fd(coarse_vol);
deviceVector<scalar *> Cd(coarse_vol);
for(int c=0;c<coarse_vol;c++){
// BLAS_V[coarse_vol][nbasis][block_vol][words]
// BLAS_F[coarse_vol][nrhs][block_vol][words]
// BLAS_C[coarse_vol][nrhs][nbasis]
scalar * Vh = & BLAS_V[c*nbasis*block_vol*words];
scalar * Fh = & BLAS_F[c*nrhs*block_vol*words];
scalar * Ch = & BLAS_C[c*nrhs*nbasis];
acceleratorPut(Vd[c],Vh);
acceleratorPut(Fd[c],Fh);
acceleratorPut(Cd[c],Ch);
// ONE bulk transfer per table. acceleratorPut is a *synchronous* 8-byte
// hipMemcpy, so the elementwise form emitted 3*coarse_vol of them per call:
// a traced run showed 272k hipMemcpy calls costing 4.5 s of API time to move
// 0.7 s worth of bytes. Same fix as BatchedBlas.h's staging rewrite.
if ( coarse_vol ) {
std::vector<scalar *> hVd(coarse_vol), hFd(coarse_vol), hCd(coarse_vol);
for(int c=0;c<coarse_vol;c++){
hVd[c] = & BLAS_V[c*nbasis*block_vol*words];
hFd[c] = & BLAS_F[c*nrhs*block_vol*words];
hCd[c] = & BLAS_C[c*nrhs*nbasis];
}
acceleratorCopyToDevice(&hVd[0],&Vd[0],coarse_vol*sizeof(scalar *));
acceleratorCopyToDevice(&hFd[0],&Fd[0],coarse_vol*sizeof(scalar *));
acceleratorCopyToDevice(&hCd[0],&Cd[0],coarse_vol*sizeof(scalar *));
}
/////////////////////////////////////////
@@ -971,13 +982,20 @@ public:
deviceVector<scalar *> &Fd,
deviceVector<scalar *> &Cd)
{
for(int c=0;c<coarse_vol;c++){
scalar * Vh = & BLAS_V[c*nbasis*block_vol*words];
scalar * Fh = & BLAS_F[c*nrhs*block_vol*words];
scalar * Ch = & BLAS_C[c*nrhs*nbasis];
acceleratorPut(Vd[c],Vh);
acceleratorPut(Fd[c],Fh);
acceleratorPut(Cd[c],Ch);
// ONE bulk transfer per table. acceleratorPut is a *synchronous* 8-byte
// hipMemcpy, so the elementwise form emitted 3*coarse_vol of them per call:
// a traced run showed 272k hipMemcpy calls costing 4.5 s of API time to move
// 0.7 s worth of bytes. Same fix as BatchedBlas.h's staging rewrite.
if ( coarse_vol ) {
std::vector<scalar *> hVd(coarse_vol), hFd(coarse_vol), hCd(coarse_vol);
for(int c=0;c<coarse_vol;c++){
hVd[c] = & BLAS_V[c*nbasis*block_vol*words];
hFd[c] = & BLAS_F[c*nrhs*block_vol*words];
hCd[c] = & BLAS_C[c*nrhs*nbasis];
}
acceleratorCopyToDevice(&hVd[0],&Vd[0],coarse_vol*sizeof(scalar *));
acceleratorCopyToDevice(&hFd[0],&Fd[0],coarse_vol*sizeof(scalar *));
acceleratorCopyToDevice(&hCd[0],&Cd[0],coarse_vol*sizeof(scalar *));
}
}
@@ -628,14 +628,15 @@ public:
GRID_ASSERT( hipMalloc((void **)&dIpiv, N*sizeof(int64_t)) == hipSuccess );
GRID_ASSERT( hipMalloc((void **)&dInfo, sizeof(int64_t)) == hipSuccess );
auto st1 = rocsolver_cgetrf_64(handle, (int64_t)N, (int64_t)N, dA, (int64_t)N, dIpiv, dInfo);
hipDeviceSynchronize();
GRID_ASSERT( hipDeviceSynchronize() == hipSuccess );
int64_t info_h = -1;
hipMemcpy(&info_h, dInfo, sizeof(int64_t), hipMemcpyDeviceToHost);
GRID_ASSERT( hipMemcpy(&info_h, dInfo, sizeof(int64_t),
hipMemcpyDeviceToHost) == hipSuccess );
std::cout << GridLogMessage << "DenseCoarseMatrix: cgetrf_64 status " << (int)st1
<< " info = " << (int)info_h << std::endl;
GRID_ASSERT(st1 == rocblas_status_success);
GRID_ASSERT(info_h == 0);
hipFree(dInfo);
GRID_ASSERT( hipFree(dInfo) == hipSuccess );
GRID_ASSERT( hipMalloc((void **)&dB, (uint64_t)CHUNKROWS*N*sizeof(ComplexF)) == hipSuccess );
// dA holds the LU of A^T; rows of A^{-1} are produced blockwise below via
// cgetrs_64 on identity-column blocks: A^T X = E => X columns = rows of
@@ -680,7 +681,7 @@ public:
(int64_t)N, (int64_t)nrow,
dA, (int64_t)N, dIpiv, dB, (int64_t)N);
GRID_ASSERT(strs == rocblas_status_success);
hipDeviceSynchronize();
GRID_ASSERT( hipDeviceSynchronize() == hipSuccess );
GRID_ASSERT( hipMemcpy(&chunk[0], dB, nelem*sizeof(ComplexF), hipMemcpyDeviceToHost) == hipSuccess );
#else
uint64_t src = (uint64_t)row0 * N;
@@ -700,9 +701,9 @@ public:
}
#ifdef GRID_HIP
if (boss) {
if (dA) hipFree(dA);
if (dB) hipFree(dB);
if (dIpiv) hipFree(dIpiv);
if (dA) GRID_ASSERT( hipFree(dA) == hipSuccess );
if (dB) GRID_ASSERT( hipFree(dB) == hipSuccess );
if (dIpiv) GRID_ASSERT( hipFree(dIpiv) == hipSuccess );
}
#endif
double t4 = usecond();
@@ -30,6 +30,8 @@ Author: Peter Boyle <pboyle@bnl.gov>
#include <Grid/algorithms/blas/BatchedBlas.h>
#include <Grid/algorithms/blas/BatchedInverse.h>
#include <algorithm>
NAMESPACE_BEGIN(Grid);
///////////////////////////////////////////////////////////////////////////////
@@ -114,6 +116,7 @@ public:
double tLeaf; // leaf inversions
double tARmin; // fastest single panel collective
double tARmax; // slowest single panel collective
std::vector<double> tARall; // every panel collective, for percentiles
uint64_t bytesAllreduce;
uint64_t nAllreduce; // panel collectives
uint64_t nGatherGemm; // GatherGemm calls
@@ -137,12 +140,24 @@ public:
// rank) are latency bound. This exists so that
// tail can be excluded with an env var rather than
// a rebuild, if it ever proves to matter.
// DENSE_BARRIER_PROBE : 1 => Barrier() before each collective, so that
// arrival skew lands in tBarrier and tAllreduce
// measures transfer alone.
// DENSE_BARRIER_PROBE : Barrier() before each panel collective.
// DEFAULTS ON, and it is an optimisation rather than instrumentation.
// Measured at 288 ranks, N=138240: comms 380.3 s -> 314.7 s and the whole
// invert 394.3 s -> 327.7 s, with the worst single collective falling
// 2474.9 ms -> 256.2 ms and effective rate 2.40 -> 4.27 GB/s. Convoy
// accounting alone predicts NO saving (barrier and collective both
// complete at max-over-ranks), so the mechanism is that the collective
// itself degrades under skewed arrival, not that skew is rebooked.
// Set DENSE_BARRIER_PROBE=0 to recover the old behaviour. It also
// separates arrival skew (tBarrier) from transfer (tAllreduce).
int useGather;
int64_t gatherMinBytes;
int barrierProbe;
// DENSE_GATHER_DEBUG=N : dump the descriptor of the first N AllGatherV
// calls from rank 0. The tiling invariant below is checked ALWAYS -- it
// costs one O(P) host loop against a collective, and a counts array that
// does not tile the panel is exactly the failure we are hunting.
int gatherDebug;
///////////////////////////////////////////////////////////////////////////
// Ownership-table validation: a proper partition of [0,N).
@@ -207,7 +222,42 @@ public:
useGather = getenv("DENSE_GATHER") ? atoi(getenv("DENSE_GATHER")) : 0;
gatherMinBytes = getenv("DENSE_GATHER_MIN_BYTES")
? atol(getenv("DENSE_GATHER_MIN_BYTES")) : 0;
barrierProbe = getenv("DENSE_BARRIER_PROBE") ? 1 : 0;
barrierProbe = getenv("DENSE_BARRIER_PROBE") ? atoi(getenv("DENSE_BARRIER_PROBE")) : 1;
gatherDebug = getenv("DENSE_GATHER_DEBUG") ? atoi(getenv("DENSE_GATHER_DEBUG")) : 0;
///////////////////////////////////////////////////////////////////////
// DENSE_GATHER=1 is KNOWN BROKEN on Cray MPICH and refuses to start.
//
// The gather asks MPI_Allgatherv to assemble a panel to which only the
// rank sub-range [rB0,rB1) contributes; every other rank passes count 0.
// At Schur depth d that is P/2^(d+1) contributors out of P, i.e. 9 of 288
// at depth 4 and ~1 of 288 at depth 7. Measured consequences at 288
// ranks (tests/debug/Test_allgather):
// T5 288/288 contributing, 18.4 MB on device : sub-second.
// T6 144/288 contributing, 9.2 MB on device : ~53 s.
// and in production, 9/288 contributing at 299 MB hangs and then aborts
// inside PMPI_Allgatherv ("req != NULL", mpir_request.h:508).
//
// The primitive itself is sound -- Test_allgather T1-T6 pass, and T3/T4
// verify it BITWISE against the zero-fill+GlobalSum idiom. What fails is
// MPI's handling of a collective in which almost every rank contributes
// nothing. Fixing it needs either P2P (a bisection allgather over
// [r0,r1)) or the 2D block-cyclic layout, which removes the shape
// altogether by giving every rank part of every sub-block.
//
// Fail here, in the first second, rather than 100 s and 36 nodes into a
// job. DENSE_GATHER_FORCE=1 proceeds anyway for debugging.
///////////////////////////////////////////////////////////////////////
if ( useGather && !getenv("DENSE_GATHER_FORCE") )
{
std::cout << GridLogError
<< "DENSE_GATHER=1 is disabled: MPI_Allgatherv on this MPI is"
<< " pathological when most ranks contribute count 0 (see the"
<< " note at RecursiveSchurInverse.h, and Test_allgather T5 vs"
<< " T6). Unset DENSE_GATHER, or set DENSE_GATHER_FORCE=1 to"
<< " proceed anyway." << std::endl;
GRID_ASSERT(0 && "DENSE_GATHER=1 known broken: see comment above");
}
}
///////////////////////////////////////////////////////////////////////////
@@ -330,6 +380,27 @@ public:
counts[r] = (int)((rowStart[r+1]-rowStart[r])*nchunk);
displs[r] = (int)((rowStart[r] -rowStart[rB0])*nchunk);
}
//////////////////////////////////////////////////////////////////
// The counts MUST tile the panel exactly: every one of the
// k*nchunk words has exactly one owner, or MPI is handed an
// inconsistent descriptor. Always checked; O(P) against a
// collective is free.
//////////////////////////////////////////////////////////////////
int64_t csum = 0; int nz = 0;
for(int r=0;r<P;r++){ csum += counts[r]; if ( counts[r] ) nz++; }
GRID_ASSERT( csum == (int64_t)panelWords );
GRID_ASSERT( displs[rB1-1] + counts[rB1-1] == (int)panelWords );
if ( gatherDebug && ((int)nAllreduce < gatherDebug) && (me==0) ) {
std::cout << GridLogMessage << "GATHER["<<nAllreduce<<"]"
<< " ranks ["<<rB0<<","<<rB1<<")"
<< " k "<<k<<" n "<<n<<" nchunk "<<nchunk
<< " panelWords "<<panelWords
<< " ("<<panelBytesThis/1024/1024<<" MB)"
<< " contributors "<<nz<<"/"<<P
<< " maxcount "<<*std::max_element(counts.begin(),counts.end())
<< " rank0: count "<<counts[me]<<(owner?" owner":" non-owner")
<< std::endl;
}
void *send = owner ? (void *)B.ColumnWindow(colB+j0) : (void *)&dRecvBuf[0];
grid->AllGatherV(send, counts[me],
(void *)&dRecvBuf[0], counts, displs, sizeof(ComplexD));
@@ -385,6 +456,7 @@ public:
tAllreduce += tar;
tARmin = std::min(tARmin, tar);
tARmax = std::max(tARmax, tar);
tARall.push_back(tar);
bytesAllreduce += panelBytesThis;
nAllreduce++;
@@ -657,6 +729,7 @@ public:
tLeaf = 0.0;
tARmin = 1.0e30;
tARmax = 0.0;
tARall.clear();
bytesAllreduce = 0;
nAllreduce = 0;
nGatherGemm = 0;
@@ -718,6 +791,29 @@ public:
<< " effective " << (tAllreduce>0 ? bytesAllreduce/tAllreduce*1.0e6/1.0e9 : 0.0)
<< " GB/s"
<< std::endl;
//////////////////////////////////////////////////////////////////////
// Distribution, not just min/avg/max. A barrier that merely REBOOKS
// skew shifts the median; a barrier that suppresses pathological
// collectives shortens the tail. p50 vs p99 separates the two from a
// single run, which min/avg/max cannot.
//////////////////////////////////////////////////////////////////////
if ( tARall.size() ) {
std::vector<double> v = tARall;
std::sort(v.begin(),v.end());
size_t n = v.size();
auto pc = [&](double f){ size_t i=(size_t)(f*(n-1)); return v[i]/1.0e3; };
double med = pc(0.50);
// how much time is spent in calls that are gross outliers
double tail=0.0; size_t ntail=0;
for(size_t i=0;i<n;i++) if ( v[i] > 4.0*med*1.0e3 ) { tail+=v[i]; ntail++; }
std::cout << GridLogMessage << "Schur comms distribution (boss):"
<< " p50 " << med
<< " p90 " << pc(0.90)
<< " p99 " << pc(0.99)
<< " ms calls > 4x median: " << ntail << "/" << n
<< " carrying " << tail/1.0e6 << " s"
<< std::endl;
}
}
};
+27
View File
@@ -315,6 +315,27 @@ public:
}
int traceID;
// Delivered-comms instrumentation, per CommunicateBegin/Complete pair.
//
// OffNodeBytes : bytes handed to MPI, i.e. EXCLUDING intranode traffic --
// StencilSendToRecvFrom* return off-node bytes only, which is the
// differentiation we want. It is BIDIRECTIONAL (send + receive) under
// ACCELERATOR_AWARE_MPI (Communicator_mpi3.cc:463,481). On the
// host-staged path the send is deferred to PollDtoH and its bytes are
// NOT currently counted, so figures from the two paths are not
// comparable. Prepare() returns 0.0 on the accelerator-aware path.
//
// CommTimer : microseconds between traceStart and traceStop, i.e. exactly
// the "Stencil::CommunicateBegin" roctx range -- transfer only, with the
// StencilBarrier and the compress kernels already excluded. It spans the
// window in which the interior kernel runs, so the derived rate is
// bandwidth delivered CONCURRENT WITH COMPUTE, which is the quantity a
// comms-only benchmark cannot see.
//
// InterNodeBandwidthMBps : the per-call rate. For a reportable figure
// accumulate bytes and time separately across calls and divide once --
// averaging per-call rates over-weights the fast calls. See
// benchmarks/Benchmark_dwf.cc.
double OffNodeBytes;
double CommTimer;
double InterNodeBandwidthMBps;
@@ -955,6 +976,12 @@ public:
bool preserve_shm=false)
{
SloppyComms = 0;
// Never leave the delivered-comms counters uninitialised: they are read
// from outside (benchmarks, drivers) and a stencil that has not yet
// exchanged would otherwise return denormal garbage.
OffNodeBytes = 0;
CommTimer = 0;
InterNodeBandwidthMBps = 0;
face_table_computed=0;
_grid = grid;
this->parameters=p;
+45 -3
View File
@@ -276,10 +276,15 @@ void Benchmark(int Ls, Coordinate Dirichlet,bool sloppy)
Dw.Dhop(src,result,0);
std::cout<<GridLogMessage<<"Called warmup"<<std::endl;
double t0=usecond();
double aveBW = 0;
// Accumulate BYTES and TIME separately and divide once. Averaging the
// per-call rate instead would over-weight the fast calls: one anomalously
// short exchange contributes an enormous rate to the mean.
double locBytes = 0; // off-node bytes moved by THIS rank, over all calls
double locTime = 0; // microseconds THIS rank spent in the exchange
for(int i=0;i<ncall;i++){
Dw.Dhop(src,result,0);
aveBW +=Dw.Stencil.InterNodeBandwidthMBps/ncall;
locBytes += Dw.Stencil.OffNodeBytes;
locTime += Dw.Stencil.CommTimer;
}
double t1=usecond();
FGrid->Barrier();
@@ -300,7 +305,34 @@ void Benchmark(int Ls, Coordinate Dirichlet,bool sloppy)
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t1-t0)<<std::endl;
std::cout<<GridLogMessage << "mflop/s per rank = "<< flops/(t1-t0)/NP<<std::endl;
std::cout<<GridLogMessage << "mflop/s per node = "<< flops/(t1-t0)/NN<<std::endl;
std::cout<<GridLogMessage << "average bandwidth = "<< aveBW <<std::endl;
//////////////////////////////////////////////////////////////////////
// Delivered inter-node bandwidth, measured CONCURRENT WITH COMPUTE.
//
// Bytes are summed over the whole machine; the exchange time is averaged
// over ranks; the result is reported per node. Reduced ONCE here rather
// than inside CommunicateComplete -- a collective in the exchange would
// perturb the quantity being measured.
//
// bytes / microsecond == MB/s exactly (decimal MB).
//
// NB OffNodeBytes is bidirectional (send + receive) only under
// ACCELERATOR_AWARE_MPI; on the host-staged path the send side is
// deferred to PollDtoH and is not currently counted, so numbers from the
// two paths are not comparable.
//////////////////////////////////////////////////////////////////////
{
double sumBytes = locBytes; FGrid->GlobalSum(sumBytes); // whole machine
double sumTime = locTime; FGrid->GlobalSum(sumTime);
double avgTime = sumTime/NP; // mean over ranks
double BWnode = sumBytes/avgTime/NN; // MB/s per node
std::cout<<GridLogMessage << "Comms: off-node bytes/call/node = "
<< sumBytes/ncall/NN/1.0e6 <<" MB"<<std::endl;
std::cout<<GridLogMessage << "Comms: in-exchange time = "
<< avgTime/ncall <<" us/call (mean over "<<NP<<" ranks)"<<std::endl;
std::cout<<GridLogMessage << "Comms: delivered inter-node bandwidth = "
<< BWnode <<" MB/s per node ("<< BWnode/1000.0
<<" GB/s per node, "<<NN<<" nodes)"<<std::endl;
}
err = ref-result;
n2e = norm2(err);
std::cout<<GridLogMessage << "norm diff "<< n2e<< " Line "<<__LINE__ <<std::endl;
@@ -408,8 +440,11 @@ void Benchmark(int Ls, Coordinate Dirichlet,bool sloppy)
FGrid->Barrier();
Dw.DhopEO(src_o,r_e,DaggerNo);
double t0=usecond();
double locBytes = 0, locTime = 0;
for(int i=0;i<ncall;i++){
Dw.DhopEO(src_o,r_e,DaggerNo);
locBytes += Dw.StencilOdd.OffNodeBytes; // DhopEO drives StencilOdd,
locTime += Dw.StencilOdd.CommTimer; // NOT the full-grid Stencil
}
double t1=usecond();
FGrid->Barrier();
@@ -420,6 +455,13 @@ void Benchmark(int Ls, Coordinate Dirichlet,bool sloppy)
std::cout<<GridLogMessage << "Deo mflop/s = "<< flops/(t1-t0)<<std::endl;
std::cout<<GridLogMessage << "Deo mflop/s per rank "<< flops/(t1-t0)/NP<<std::endl;
std::cout<<GridLogMessage << "Deo mflop/s per node "<< flops/(t1-t0)/NN<<std::endl;
{
double sumBytes = locBytes; FGrid->GlobalSum(sumBytes);
double sumTime = locTime; FGrid->GlobalSum(sumTime);
double avgTime = sumTime/NP;
std::cout<<GridLogMessage << "Deo comms: delivered inter-node bandwidth = "
<< sumBytes/avgTime/NN <<" MB/s per node"<<std::endl;
}
}
Dw.DhopEO(src_o,r_e,DaggerNo);
Dw.DhopOE(src_e,r_o,DaggerNo);
@@ -90,6 +90,14 @@ RealD OuterTol = 1.0e-8;
int OuterMmax = 8;
int OuterNstep = 8;
// Halo exchange in reduced precision on the fine stencils. Stencil::SloppyComms
// is a plain runtime setter (Stencil.h:303) applied to the full/even/odd stencils
// of an operator that is already built, so this costs no extra storage and no
// second operator. Benchmark_dwf at this decomposition measures ~1.6x on the
// fine Dhop, consistent with the uncompressed halo being ~1.6x the interior
// compute time. Set FineSloppyComms=0 to recover the exact-comms behaviour.
int FineSloppyComms = 1;
void ParseEnvironment(void)
{
if(getenv("MASS")) mass = atof(getenv("MASS"));
@@ -104,6 +112,7 @@ void ParseEnvironment(void)
if(getenv("CoarseSolverOrder")) CoarseSolverOrder = atoi(getenv("CoarseSolverOrder"));
if(getenv("OuterTol")) OuterTol = atof(getenv("OuterTol"));
if(getenv("OuterMmax")) OuterMmax = atoi(getenv("OuterMmax"));
if(getenv("FineSloppyComms")) FineSloppyComms = atoi(getenv("FineSloppyComms"));
if(getenv("OuterNstep")) OuterNstep = atoi(getenv("OuterNstep"));
if(getenv("LATT")){
Coordinate l;
@@ -475,6 +484,16 @@ int main (int argc, char ** argv)
MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c);
// Reduced-precision halo on both fine operators. PVdagM and ShiftedPVdagM are
// thin wrappers over these same objects, so this covers coarsening, the fine
// smoother, the V-cycle and the outer Krylov alike. The final residual check
// below turns it off again: a verification computed with a sloppy operator
// would certify the wrong matrix.
Ddwf.SloppyComms(FineSloppyComms);
Dpv .SloppyComms(FineSloppyComms);
std::cout << GridLogMessage << "Fine stencils: SloppyComms = " << FineSloppyComms
<< (FineSloppyComms ? " (reduced-precision halo)" : " (exact halo)") << std::endl;
typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
PVdagM_t PVdagM(Ddwf,Dpv);
@@ -815,13 +834,20 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage << "V2 3-level solve Nrhs "<<nr<<" total " << w.Elapsed()
<< " (per RHS: " << w.useconds()/1.0e6/nr << " s)" << std::endl;
// Verify against the EXACT operator: the solve may have used a reduced
// precision halo, but the residual we report must not.
Ddwf.SloppyComms(0);
Dpv .SloppyComms(0);
{ LatticeFermionD Ax(FGrid); RealD worst=0.0;
for(int r=0;r<nr;r++){ PVdagM.Op(sol[r],Ax); Ax=Ax-src[r];
RealD rn=std::sqrt(norm2(Ax)/norm2(src[r]));
std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": rhs["<<r<<"] true residual = " << rn << std::endl;
worst=std::max(worst,rn); }
std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": worst-case residual = " << worst << std::endl;
std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": worst-case residual = " << worst
<< " (exact-halo verification)" << std::endl;
}
Ddwf.SloppyComms(FineSloppyComms);
Dpv .SloppyComms(FineSloppyComms);
// The operators borrow these grids and build a PaddedCell on them, so
// they must let go before the grids are destroyed.
+17 -6
View File
@@ -201,8 +201,16 @@ int main(int argc, char **argv)
// 288 ranks, while the host-buffer stages above passed. Reproducing that
// here costs one node and a second instead of a 36-node job.
////////////////////////////////////////////////////////////////////////
// AG_WORDS : ComplexD words each CONTRIBUTING rank sends (default 4096)
// AG_NPART : how many ranks contribute at all in T6 (default P/2).
// The GatherGemm shape at Schur depth d has only P/2^(d+1) contributors,
// so depth 7 at 288 ranks is AG_NPART=2 AG_WORDS=518400 (8.3 MB each,
// 18.7 MB gathered). Depth 0 is AG_NPART=144 AG_WORDS=232800.
const int64_t AGW = getenv("AG_WORDS") ? atol(getenv("AG_WORDS")) : 4096;
const int AGN = getenv("AG_NPART") ? atoi(getenv("AG_NPART")) : (P+1)/2;
{
const int64_t W = 4096; // words per contributing rank
const int64_t W = AGW; // words per contributing rank
std::vector<int> counts(P,0), displs(P,0);
int64_t total=0;
for(int r=0;r<P;r++){ counts[r]=(int)W; displs[r]=(int)total; total+=W; }
@@ -220,13 +228,15 @@ int main(int argc, char **argv)
for(int r=0;r<P;r++)
for(int64_t i=0;i<W;i++)
if ( hall[displs[r]+i] != Stamp(r,i) ) ok=false;
Report("T5 AllGatherV on DEVICE buffers", ok);
Report("T5 AllGatherV on DEVICE, all "+std::to_string(P)+" contributing, "+
std::to_string(W*16/1024)+" KB each, "+
std::to_string(total*16/1048576)+" MB gathered", ok);
}
{
// Only the first half of the ranks contribute; the rest send count 0.
const int64_t W = 4096;
int half = (P+1)/2;
const int64_t W = AGW;
int half = AGN < 1 ? 1 : (AGN > P ? P : AGN);
std::vector<int> counts(P,0), displs(P,0);
int64_t total=0;
for(int r=0;r<P;r++){
@@ -247,8 +257,9 @@ int main(int argc, char **argv)
for(int r=0;r<half;r++)
for(int64_t i=0;i<W;i++)
if ( hall[displs[r]+i] != Stamp(r,i) ) ok=false;
Report("T6 AllGatherV on DEVICE, "+std::to_string(P-half)+"/"+std::to_string(P)+
" ranks sending count 0", ok);
Report("T6 AllGatherV on DEVICE, "+std::to_string(half)+"/"+std::to_string(P)+
" contributing, "+std::to_string(W*16/1024)+" KB each, "+
std::to_string(total*16/1048576)+" MB gathered", ok);
}
std::cout << GridLogMessage << (failures ? "AllGather regression: FAILURES"