diff --git a/Grid/algorithms/blas/BatchedBlas.h b/Grid/algorithms/blas/BatchedBlas.h index c0be5dd91..a4de5d696 100644 --- a/Grid/algorithms/blas/BatchedBlas.h +++ b/Grid/algorithms/blas/BatchedBlas.h @@ -995,7 +995,7 @@ public: /////////////////////////////////////////////////////////////////////////////////// // Explicit-leading-dimension complex double GEMM. Mirror of the ComplexF // overload above; motivating use is the fp64 distributed recursive Schur - // inversion (RecursiveSchurInverse), whose operands are column windows of + // inversion (the distributed Schur recursion), whose operands are column windows of // larger row-slab allocations. /////////////////////////////////////////////////////////////////////////////////// void gemmBatched(GridBLASOperation_t OpA, diff --git a/Grid/algorithms/multigrid/BlockCyclicRedistribute.h b/Grid/algorithms/multigrid/BlockCyclicRedistribute.h index 9fb3e55aa..69eb24d45 100644 --- a/Grid/algorithms/multigrid/BlockCyclicRedistribute.h +++ b/Grid/algorithms/multigrid/BlockCyclicRedistribute.h @@ -23,11 +23,47 @@ Author: Peter Boyle NAMESPACE_BEGIN(Grid); +/////////////////////////////////////////////////////////////////////////////// +// My rows of a distributed dense matrix in the 1D rank-major row layout: +// rank r owns contiguous global rows [rowStart[r], rowStart[r+1]) of an +// N x N matrix, stored rows x cols column major, ld = rows; element (i,j) +// at data[i + j*ld]. This is the layout the stencil->dense import +// produces and the apply slab consumes; the 2D inverse slots between them +// through the redistribution below. +/////////////////////////////////////////////////////////////////////////////// +class BlockRows +{ +public: + deviceVector data; + int64_t rows; + int64_t cols; + int64_t ld; + + BlockRows() + { + rows = 0; + cols = 0; + ld = 0; + } + void Resize(int64_t r, int64_t c) + { + rows = r; + cols = c; + ld = r; + data.resize((uint64_t)r*c); + } + ComplexD *ColumnWindow(int64_t col0) + { + GRID_ASSERT( col0 >= 0 ); + GRID_ASSERT( col0 <= cols ); + return &data[(uint64_t)col0*ld]; + } +}; + /////////////////////////////////////////////////////////////////////////////// // Stage 4 of the 2D distributed dense inverse: redistribution between the -// 1D rank-major row layout (BlockRows: rank r owns contiguous global rows -// [rowStart[r], rowStart[r+1]) of an N x N matrix, stored rows x N column -// major with ld = rows) and the 2D block-cyclic layout. +// 1D rank-major row layout (BlockRows above) and the 2D block-cyclic +// layout. // // This is what lets the EXISTING stencil->dense import, its certificate, // the fp32 slab conversion and the apply path all remain byte-for-byte diff --git a/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h b/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h index 5c1a96382..14344d708 100644 --- a/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h +++ b/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h @@ -28,12 +28,11 @@ NAMESPACE_BEGIN(Grid); // Stage 3 of the 2D distributed dense inverse: the recursive Schur // complement on a block-cyclic matrix, in place. // -// The nine-step algebra is IDENTICAL to RecursiveSchurInverse (1D); what -// changes is the decomposition. The recursion splits the GLOBAL INDEX -// RANGE at the block boundary nearest the midpoint -- not the rank range -- -// so every rank owns part of every sub-block at every depth, and the -// ownership gating (inI/inJ, dummy operands, zero-width rank ranges) of the -// 1D scheme has no analogue here: it is simply gone. +// The recursion splits the GLOBAL INDEX RANGE at the block boundary +// nearest the midpoint -- not the rank range as the retired 1D +// RecursiveSchurInverse did -- so every rank owns part of every sub-block +// at every depth, and the ownership gating (inI/inJ, dummy operands, +// zero-width rank ranges) a 1D scheme needs has no analogue here. // // I = [c0,m) J = [m,c1) (block-aligned, m the mid block boundary) // 1. recurse I : A11 -> A11inv (in place) @@ -64,8 +63,7 @@ NAMESPACE_BEGIN(Grid); // LEAF. A leaf is a single diagonal block, and block (b,b) of a // block-cyclic layout lives ENTIRELY on rank (b%Pr, b%Pc). The leaf // inversion is therefore purely local -- pack the strided block dense, -// GridBLASInverse, unpack -- with NO communication and no assembly. The -// leaf-assembly transport question of the 1D scheme does not arise. +// GridBLASInverse, unpack -- with NO communication and no assembly. // Successive leaves cycle over ranks, so leaf work is naturally spread. // // COMMUNICATION. Every transfer in the whole inversion is a @@ -75,12 +73,9 @@ NAMESPACE_BEGIN(Grid); // one optional exception: it performs reductions, and is only ever called // explicitly by a caller who wants the numbers. // -// NUMERICS. No pivoting, exactly as the 1D scheme: every A11 and every -// Schur complement met on the way down must be non-singular. The growth -// telemetry stands in for pivoting; note the recursion splits differently -// from the 1D rank-range tree, so DIFFERENT sub-blocks are inverted and -// telemetry values are NOT comparable with the 1D implementation's -- -// re-baseline, do not compare. +// NUMERICS. No pivoting: every A11 and every Schur complement met on the +// way down must be non-singular. The growth telemetry stands in for +// pivoting. /////////////////////////////////////////////////////////////////////////////// class BlockCyclicSchurInverse @@ -92,14 +87,18 @@ public: // Telemetry: accumulated LOCALLY, no comms unless ReportTelemetry(). double telLeafMaxInv; uint64_t nLeaf; - // BIG LEAVES (SCHUR2D_LEAF_SPAN=s, default 1 = off). Below span s blocks a - // sub-block lives on <= s of the Pr process rows / s of the Pc columns; the - // SUMMA rings then run on a few ranks while the rest block in their next - // SendToRecvFrom (histogram 2026-08-27: 93% of ring time in the 3.7 MB - // single-block panels of exactly these levels). Instead: gather the - // (s*nb)^2 sub-block to one rank, invert locally, scatter back. - int leafSpan = -1; - int leafLU = -1; // SCHUR2D_LEAF_LU=1: big-leaf inverse via GridBLASInverse::inverseLU (blocked getrf_64 + identity getrs_64) instead of getri_batched + // BIG LEAVES. Below span s blocks a sub-block lives on <= s of the Pr + // process rows / s of the Pc columns; when s is small RELATIVE TO THE + // GRID the SUMMA rings run on a few ranks while the rest block in their + // next SendToRecvFrom (histogram 2026-08-27: 93% of ring time in the + // 3.7 MB single-block panels of exactly these levels). Instead: gather + // the (s*nb)^2 sub-block to one rank, invert locally (inverseLU), scatter + // back. Fires only while span < min(Pr,Pc) -- when the whole grid + // participates the rings are not degenerate and gathering would only + // concentrate memory (and at the top of the tree, gather the whole + // matrix). Default 9: banked on Frontier 2026-08-27 (invert 132 -> 27.6 s + // at N=138240 on a 16x18 grid; s=18 gave 30.6, s=4 37.0). + int leafSpan = 9; uint64_t nBigLeaf = 0; int64_t maxBigW = 0; double tBigGather = 0, tBigInv = 0, tBigScatter = 0; uint64_t nNode; @@ -284,18 +283,9 @@ public: } tBigGather += usecond(); - // ---- invert on root ---- + // ---- invert on root: blocked getrf_64 + identity getrs_64 ---- tBigInv -= usecond(); - if ( leafLU < 0 ) leafLU = getenv("SCHUR2D_LEAF_LU") ? atoi(getenv("SCHUR2D_LEAF_LU")) : 0; - if ( me == root ) { - if ( leafLU ) { - INV.inverseLU(W, &dense[0]); - } else { - deviceVector bp(1); std::vector ptr(1); ptr[0] = &dense[0]; - acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexD*)); - INV.inverseBatched(W, bp); - } - } + if ( me == root ) INV.inverseLU(W, &dense[0]); tBigInv += usecond(); // ---- scatter ---- @@ -330,8 +320,7 @@ public: int64_t span = b1-b0; GRID_ASSERT( span >= 1 ); if ( span == 1 ) { Leaf(A, b0); return; } - if ( leafSpan < 0 ) leafSpan = getenv("SCHUR2D_LEAF_SPAN") ? atoi(getenv("SCHUR2D_LEAF_SPAN")) : 1; - if ( span <= leafSpan ) { BigLeaf(A, b0, b1); return; } + if ( span <= leafSpan && span < std::min(L.Pr,L.Pc) ) { BigLeaf(A, b0, b1); return; } GRID_TRACE("SchurNode"); nNode++; @@ -376,150 +365,11 @@ public: // PUBLIC ENTRY. In-place inverse of the whole matrix. Scratch (4x the // matrix footprint) is allocated here and released on return. /////////////////////////////////////////////////////////////////////////// - /////////////////////////////////////////////////////////////////////////// - // POWER - CLOCK - GROUND. Before the inverse runs, print the preconditions - // that have differed between fast (62 s) and slow (133-141 s) runs of the - // SAME inverse, and measure ONE SendToRecvFrom to the ACTUAL ring partners - // at three sizes, device and host buffers. Same code in the harness and in - // production, so the two processes are compared on the identical primitive - // before any explanation of the SUMMA rings is attempted. - // SCHUR2D_PROBE=0 disables (costs ~0.1-0.5 s). - /////////////////////////////////////////////////////////////////////////// - void Probe(BlockCyclicMatrix &A) - { - BlockCyclicLayout &L = A.layout; - GridBase *grid = A.grid; - int me = grid->ThisRank(); - // --- banner --- - int thr = -1; -#ifdef GRID_COMMS_MPI3 - MPI_Query_thread(&thr); -#endif - const char *omp = getenv("OMP_NUM_THREADS"); - MemoryStatus ms = MemoryManager::GetFootprint(); - std::cout << GridLogMessage << "Schur2D PROBE banner: MPI thread level " << thr - << " (0 single,1 funneled,2 serialized,3 multiple) OMP_NUM_THREADS=" << (omp?omp:"unset") - << " MemoryManager device bytes " << ms.DeviceBytes/1.0e9 << " GB (LRU " << ms.DeviceLRUBytes/1.0e9 - << " GB, cap " << ms.DeviceMaxBytes/1.0e9 << " GB)" - << " grid " << L.Pr << "x" << L.Pc << " nb " << L.nb << std::endl; -#ifdef GRID_HIP - if ( me==0 ) acceleratorMem(); -#endif - // --- ring partners exactly as SUMMA uses them --- - int prow=L.prow, pcol=L.pcol, Pr=L.Pr, Pc=L.Pc; - struct Ring { const char *name; int dest, src; }; - Ring rings[2] = { {"ringA(row, q+-1)", prow*Pc + (pcol+1)%Pc, prow*Pc + (pcol-1+Pc)%Pc}, - {"ringB(col, p+-1)", ((prow+1)%Pr)*Pc + pcol, ((prow-1+Pr)%Pr)*Pc + pcol} }; - // 2/3/4 MB added 2026-08-27: the SUMMA histogram put 93% of ring time in - // [2,4) MB messages at 0.3 GB/s while >=4 MB ran at 11-13 GB/s. - const int NSZ = 6; - uint64_t sizes[NSZ] = { 64ull*1024, 1024ull*1024, 2048ull*1024, 3072ull*1024, 4096ull*1024, 8ull*1024*1024 }; - uint64_t maxb = sizes[NSZ-1]; - deviceVector dsend(maxb), drecv(maxb); - std::vector hsend(maxb), hrecv(maxb); - for(int r=0;r<2;r++){ - if ( (r==0 && Pc==1) || (r==1 && Pr==1) ) continue; - int off = grid->IsOffNode(rings[r].dest); - for(int si=0;siSendToRecvFrom(&dsend[0], rings[r].dest, &drecv[0], rings[r].src, bytes); - double t0=usecond(); - for(int i=0;i<5;i++) grid->SendToRecvFrom(&dsend[0], rings[r].dest, &drecv[0], rings[r].src, bytes); - double td=(usecond()-t0)/5.0; - grid->SendToRecvFrom(&hsend[0], rings[r].dest, &hrecv[0], rings[r].src, bytes); - t0=usecond(); - for(int i=0;i<5;i++) grid->SendToRecvFrom(&hsend[0], rings[r].dest, &hrecv[0], rings[r].src, bytes); - double th=(usecond()-t0)/5.0; - // spread over ranks - RealD dmax=td, dmin=-td; grid->GlobalMax(dmax); grid->GlobalMax(dmin); dmin=-dmin; - std::cout << GridLogMessage << "Schur2D PROBE " << rings[r].name << (off?" OFF-node":" on-node") - << " " << bytes/1024 << " KB: device " << td << " us (" << bytes/td/1.0e3 << " GB/s) [min/max over ranks " << dmin << "/" << dmax << " us]" - << " host " << th << " us (" << bytes/th/1.0e3 << " GB/s)" << std::endl; - } - } - /////////////////////////////////////////////////////////////////////// - // ONE-GCD LEAF MICROBENCHMARK (boss only). The big-leaf inverse at - // W=4320 measured 0.53 s with BOTH getri_batched and getrf_64+getrs_64 - // (2026-08-27) -- ~0.4 TF/s on a GCD that runs zgemm at ~15. Time the - // three primitives in isolation on a well-conditioned dense matrix so the - // leaf's rate can be compared with the machine's, and getrf split from - // getrs. Sizes: the W of span 4 / 9 / 18 leaves on accelerator builds; - // tiny on CPU builds (Eigen would take minutes at 4320). - /////////////////////////////////////////////////////////////////////// - if ( me == 0 ) { -#if defined(GRID_HIP) || defined(GRID_CUDA) || defined(GRID_SYCL) - std::vector Ws({1920, 4320, 8640}); -#else - std::vector Ws({240, 480}); -#endif - for(int64_t W : Ws){ - deviceVector M((uint64_t)W*W), C((uint64_t)W*W); - { ComplexD *m = &M[0]; const int64_t WW = W; // diagonally dominant: (i==j ? W : 0) + cos/sin noise - accelerator_for(idx,(uint64_t)W*W,1,{ int64_t j=idx/WW, i=idx-j*WW; double x=0.37*i+0.61*j; - m[idx] = ComplexD((i==j)?(double)WW:0.0,0.0) + ComplexD(std::cos(x),std::sin(1.3*x)); }); - accelerator_barrier(); } - double flopLU = 8.0/3.0*(double)W*W*W; // complex LU ~ (4 real flops per complex mult-add) * (2/3 n^3) - double flopGEMM = 8.0*(double)W*W*W; // complex GEMM - // 1. getri_batched (batch 1) - double tb; - { deviceVector bp(1); std::vector ptr(1); ptr[0]=&M[0]; - acceleratorCopyToDevice(&ptr[0],&bp[0],sizeof(ComplexD*)); - double t0=usecond(); INV.inverseBatched(W,bp); tb=usecond()-t0; } - // 2. inverseLU (getrf_64 + identity getrs_64), timed inside - double tl; { double t0=usecond(); INV.inverseLU(W,&M[0]); tl=usecond()-t0; } - // 3. one zgemm W x W x W for the machine rate - double tg; - { deviceVector ap(1),bp(1),cp(1); std::vector ptr(1); - ptr[0]=&M[0]; acceleratorCopyToDevice(&ptr[0],&ap[0],sizeof(ComplexD*)); acceleratorCopyToDevice(&ptr[0],&bp[0],sizeof(ComplexD*)); - ptr[0]=&C[0]; acceleratorCopyToDevice(&ptr[0],&cp[0],sizeof(ComplexD*)); - double t0=usecond(); - SUMMA.BLAS.gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N,(int)W,(int)W,(int)W,ComplexD(1.0,0.0),ap,(int)W,bp,(int)W,ComplexD(0.0,0.0),cp,(int)W); - SUMMA.BLAS.synchronise(); tg=usecond()-t0; } - std::cout << GridLogMessage << "Schur2D PROBE leaf W=" << W - << ": getri_batched " << tb/1.0e6 << " s (" << flopLU/tb/1.0e6 << " TF/s)" - << " getrf_64+getrs_64 " << tl/1.0e6 << " s (getrf " << INV.lastGetrfUs/1.0e6 << " getrs " << INV.lastGetrsUs/1.0e6 << ")" - << " zgemm " << tg/1.0e6 << " s (" << flopGEMM/tg/1.0e6 << " TF/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[NSZ-1]; - 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) { BlockCyclicLayout &L = A.layout; GRID_ASSERT( L.N >= 1 ); int64_t nblocks = (L.N + L.nb - 1)/L.nb; - if ( !(getenv("SCHUR2D_PROBE") && atoi(getenv("SCHUR2D_PROBE"))==0) ) Probe(A); BlockCyclicMatrix Bt(A.grid, L.N, L.nb, L.Pr, L.Pc); BlockCyclicMatrix Ct(A.grid, L.N, L.nb, L.Pr, L.Pc); @@ -550,8 +400,8 @@ public: << std::endl; if ( nBigLeaf ) { RealD ti = tBigInv/1.0e6; grid->GlobalMax(ti); // inverse runs on the root of each leaf: report the max over ranks - std::cout << GridLogMessage << "BlockCyclicSchurInverse: BIG LEAVES (SCHUR2D_LEAF_SPAN=" << leafSpan - << (leafLU>0 ? ", SCHUR2D_LEAF_LU=1: getrf_64+getrs_64" : ", getri_batched") << "): " << nBigLeaf + std::cout << GridLogMessage << "BlockCyclicSchurInverse: BIG LEAVES (span " << leafSpan + << ", inverseLU): " << nBigLeaf << " leaves, max W " << maxBigW << " boss secs: gather " << tBigGather/1.0e6 << " scatter " << tBigScatter/1.0e6 << " inverse (max over ranks) " << ti << std::endl; @@ -576,8 +426,7 @@ public: // same sequence of sizes). Time is wall time inside SendToRecvFrom, so it // includes waiting for the partner -- a bucket whose GB/s is far below the // probe's for the same size is wait, not wire. - std::cout << GridLogMessage << "BlockCyclicSumma ring histogram (boss): size-bucket msgs GB xfer-secs GB/s %time" - << (SUMMA.handshake>0 ? " handshake-secs (partner wait, excluded from xfer)" : "") << std::endl; + std::cout << GridLogMessage << "BlockCyclicSumma ring histogram (boss): size-bucket msgs GB xfer-secs GB/s %time" << std::endl; std::streamsize oldprec = std::cout.precision(); for(int b=0;b0 ? g/sec : 0.0) << std::setw(8) << std::setprecision(3) << (ring>0 ? 100.0*sec/ring : 0.0); - if ( SUMMA.handshake>0 ) std::cout << std::setw(12) << std::setprecision(3) << SUMMA.histHsUs[b]/1.0e6; std::cout << std::endl; } std::cout.precision(oldprec); diff --git a/Grid/algorithms/multigrid/BlockCyclicSumma.h b/Grid/algorithms/multigrid/BlockCyclicSumma.h index da3dc7bb5..a8f92dac9 100644 --- a/Grid/algorithms/multigrid/BlockCyclicSumma.h +++ b/Grid/algorithms/multigrid/BlockCyclicSumma.h @@ -146,16 +146,12 @@ public: // Per-message-size histogram (bucket = floor(log2 bytes)): count, bytes, // microseconds -- decomposes the ring time by packet size so a low average // GB/s can be attributed (many small latency-bound messages vs slow large - // ones vs partner-wait). The 8 MB probe runs at 11-20 GB/s; SUMMA averaged 2. + // ones vs partner-wait). Time is wall time inside SendToRecvFrom, so a + // bucket far below the wire rate for its size is wait, not wire. static const int NHIST=48; - uint64_t histN[NHIST]={0}, histBytes[NHIST]={0}; double histUs[NHIST]={0}, histHsUs[NHIST]={0}; - // SUMMA_HANDSHAKE=1: a 4-byte SendToRecvFrom with the same partner - // immediately before each ring message, timed separately (histHsUs). - // Handshake time = partner-arrival skew; the remainder = transfer. Splits - // the [2,4) MB bucket's 12 ms/msg (2026-08-27) into wait vs wire. - int handshake = -1; int hsTx=0, hsRx=0; - void HistAdd(uint64_t bytes, double us, double hs=0.0){ int b=0; while((bytes>>b)>1) b++; histN[b]++; histBytes[b]+=bytes; histUs[b]+=us; histHsUs[b]+=hs; } - void ResetTelemetry(void){ tAlloc=tPack=tRingA=tRingB=tGemm=0; bytesRing=nRingMsg=nMultiply=nGemm=0; for(int b=0;b>b)>1) b++; histN[b]++; histBytes[b]+=bytes; histUs[b]+=us; } + void ResetTelemetry(void){ tAlloc=tPack=tRingA=tRingB=tGemm=0; bytesRing=nRingMsg=nMultiply=nGemm=0; for(int b=0;b(slotA*Pc,1) ) Abuf.resize( std::max(slotA*Pc,1) ); if ( Bbuf.size() < std::max(slotB*Pr,1) ) Bbuf.resize( std::max(slotB*Pr,1) ); @@ -286,12 +281,11 @@ public: for(int t=1;tSendToRecvFrom((void *)&hsTx, dest, (void *)&hsRx, src, sizeof(int)); ths = usecond()-tm; tm = usecond(); } + double tm = usecond(); grid->SendToRecvFrom((void *)(&Abuf[0]+slotA*cs), dest, (void *)(&Abuf[0]+slotA*cr), src, slotA*sizeof(ComplexD)); - HistAdd(slotA*sizeof(ComplexD), usecond()-tm, ths); + HistAdd(slotA*sizeof(ComplexD), usecond()-tm); bytesRing += slotA*sizeof(ComplexD); nRingMsg++; } tRingA += usecond(); @@ -307,12 +301,11 @@ public: for(int t=1;tSendToRecvFrom((void *)&hsTx, dest, (void *)&hsRx, src, sizeof(int)); ths = usecond()-tm; tm = usecond(); } + double tm = usecond(); grid->SendToRecvFrom((void *)(&Bbuf[0]+slotB*rs), dest, (void *)(&Bbuf[0]+slotB*rr), src, slotB*sizeof(ComplexD)); - HistAdd(slotB*sizeof(ComplexD), usecond()-tm, ths); + HistAdd(slotB*sizeof(ComplexD), usecond()-tm); bytesRing += slotB*sizeof(ComplexD); nRingMsg++; } tRingB += usecond(); diff --git a/Grid/algorithms/multigrid/DenseCoarseMatrix.h b/Grid/algorithms/multigrid/DenseCoarseMatrix.h index 5b0873ec0..51bbbecc3 100644 --- a/Grid/algorithms/multigrid/DenseCoarseMatrix.h +++ b/Grid/algorithms/multigrid/DenseCoarseMatrix.h @@ -29,7 +29,6 @@ Author: Peter Boyle #include #include -#include #include #include @@ -39,46 +38,33 @@ NAMESPACE_BEGIN(Grid); ////////////////////////////////////////////////////////////////////////////////////// // DenseCoarseMatrix: a coarsened operator treated as a DENSE matrix -- explicit, -// row-distributed A^{-1} of a GeneralCoarsenedMatrix. Library-grade successor of -// the example-local DistributedDenseInverse (Example_pvdagm_mrhs_3level_dense.cc, -// FROZEN as the regression baseline). +// row-distributed A^{-1} of a GeneralCoarsenedMatrix. // -// What is new versus the example class: // - Stencil -> dense DIRECT IMPORT. The coarse operator IS the dense matrix // unrolled: Dense[(s,a),(s+shift_p,b)] += A[p][s]_{a,b}. Rows of my sites are -// assembled from purely LOCAL _A[p] data: no operator applies, no comms -- the -// O(N) probe assembly (93 s at N=69120) is retired. ACCUMULATE (+=) because on -// short axes distinct shifts wrap to the same neighbour. An IMPORT CERTIFICATE -// compares the dense apply against Op.M on a NON-CONSTANT vector (a constant one -// cannot see a shift-sign error); DENSE_IMPORT_SIGN=-1 flips the convention -// without recompiling. +// assembled from purely LOCAL _A[p] data: no operator applies, no comms. +// ACCUMULATE (+=) because on short axes distinct shifts wrap to the same +// neighbour. An IMPORT CERTIFICATE compares the dense apply against Op.M on a +// NON-CONSTANT vector (a constant one cannot see a shift-sign error). +// +// - Inversion is END-TO-END fp64 through the 2D block-cyclic recursive Schur +// complement (BlockCyclicSchurInverse): fp64 rank-major import -> +// RowsToCyclic -> in-place recursion (pure point-to-point SUMMA rings and +// local leaves; bitwise reproducible) -> CyclicToRows -> ONE terminal +// rounding into the fp32 apply slab. Distributed at every N and P; banked +// 3.87x against its retired 1D predecessor and >10x against SLATE at +// N=138240 on 288 GCDs. +// // - Split-K apply through GridBLAS.gemmBatched with EXPLICIT leading dimensions // (arXiv:2409.03904 fig 11): the tiny-output/huge-K GEMM Y = slab^T X becomes -// DENSE_SPLITK chunk-GEMMs by pointer offset into the resident slab (lda = N), -// partials reduced in one accelerator_for. Platform-agnostic: deviceVector + -// GridBLAS run the SAME code on HIP/CUDA/SYCL and CPU(Eigen). -// - deviceVector everywhere in the apply path; the ONE surviving naked-HIP block -// is the boss inversion buffer (quarantined below, documented). +// SPLITK chunk-GEMMs by pointer offset into the resident slab (lda = N), +// partials reduced in one accelerator_for. The source vector is assembled by +// a cartesian ring ALLGATHER of device buffers (pure P2P; ~8x fewer bytes +// than a padded allreduce, and no collective size cliffs). Platform-agnostic: +// deviceVector + GridBLAS run the SAME code on HIP/CUDA/SYCL and CPU(Eigen). // -// Setup: SLAB_FILE= loads per-rank . (header-guarded N/nrows/ -// nbasis -- the interchange format shared with the frozen example; the STEM must -// encode cfg/mass/blocking/nbasis, only the header is guarded). Absent: direct -// import -> import certificate -> chunked zero-fill+GlobalSum gather streamed to -// the boss GCD -> cgetrf_64 (ILP64) -> rows of A^{-1} via blocked identity -// cgetrs_64 + broadcast, each rank keeping the rows of its own sites -> save. -// VERIFY ||A Ainv x - x||/||x|| runs in BOTH paths (and now certifies the DEVICE -// slab + split-K path, since the single-RHS apply routes through the same core). -// -// Env: SLAB_FILE DENSE_SPLITK (default 32, snapped to a divisor of N) -// DENSE_DEVICE_SUM DENSE_IMPORT_SIGN DENSE_APPLY_PROFILE DENSE_CC_CHECK -// DENSE_SCHUR (0/absent: single-GCD gather-invert; 1: distributed -// recursive Schur; 2: AUDIT -- run BOTH on the same imported A, report -// the slab difference, keep the Schur result) DENSE_PANEL_BYTES -// -// The DENSE_SCHUR=1 path is the RecursiveSchurInverse distributed -// factorisation: it lifts the fp32 N ~ 90k boss-HBM ceiling (the CC-grid -// 256-rank SIMD cap remains -- separate issue). Internal only: slab layout, -// apply path, SLAB_FILE format and VERIFY are identical in every mode. +// VERIFY ||A Ainv x - x||/||x|| certifies the DEVICE slab + split-K path at the +// end of Import, since the single-RHS apply routes through the same core. // // Tensor-depth agnostic: site scalar objects treated as contiguous ComplexD // (iScalar wrappers add no data), so any MG level's coarse operator imports. @@ -111,15 +97,15 @@ public: std::vector myGsite; // global lex site index of my site ss std::vector slab; // nrows x N row-major: A during setup, rows of A^{-1} after - static const int64_t CHUNKROWS = 1024; // getrs harvest block (trsm efficiency + fewer broadcasts) static const int MRHS_MAX = 32; + static const int SPLITK = 32; // requested split-K chunk count, snapped DOWN to a divisor of N // Apply machinery: resident slab + persistent buffers + AOT split-K pointers. GridBLAS BLAS; deviceVector dSlab; deviceVector dX; // N x MRHS_MAX deviceVector dY; // nrows x MRHS_MAX - deviceVector dG; // N x MRHS_MAX lex-major staging for the allgather (devSum==4) + deviceVector dG; // N x MRHS_MAX lex-major staging for the allgather deviceVector dLex2Rank;// lex index of a process coordinate -> its rank (allgather block order -> row-block order) deviceVector dRm2G; // rank-major index (rank*nrows + ss*nbasis + b) -> global column (gsite*nbasis + b) of x / the slab int myLex; @@ -130,8 +116,6 @@ public: std::vector hX; std::vector hY; int NK; // split-K chunk count (divides N) - int devSum; - double schurAuditRel; // DENSE_SCHUR=2: rel slab diff single-vs-schur (-1 = not run) DenseCoarseMatrix(GridBase *g) : grid(g) @@ -142,7 +126,6 @@ public: N = grid->gSites() * nbasis; lsites = grid->lSites(); nrows = (int64_t)lsites * nbasis; - schurAuditRel = -1.0; std::cout << GridLogMessage << "DenseCoarseMatrix: N = " << N << " (" << grid->gSites() << " sites x " << nbasis << ")" @@ -180,47 +163,11 @@ public: void Import(CoarseOp &Op) { double t0 = usecond(); - //////////////////////////////////////////////////////////////////// - // 0. Slab cache: SLAB_FILE= -> per-rank raw file .. - // SAME format as the frozen example (interchange compatible). - //////////////////////////////////////////////////////////////////// - bool loaded = false; - char *sfile = getenv("SLAB_FILE"); - std::string slabfile; - if (sfile) { - slabfile = std::string(sfile) + "." + std::to_string(grid->ThisRank()); - FILE *f = fopen(slabfile.c_str(),"rb"); - if (f) { - int64_t hdr[4] = {0,0,0,0}; - GRID_ASSERT( fread(hdr,sizeof(int64_t),4,f) == 4 ); - GRID_ASSERT( hdr[0] == (int64_t)0x44454E5345 ); // magic "DENSE" - GRID_ASSERT( hdr[1] == N && hdr[2] == (int64_t)nrows && hdr[3] == (int64_t)nbasis ); - uint64_t nelem = (uint64_t)nrows * N; - GRID_ASSERT( fread(&slab[0], sizeof(ComplexF), nelem, f) == nelem ); - fclose(f); - loaded = true; - std::cout << GridLogMessage << "DenseCoarseMatrix: slab loaded from " - << slabfile << " -- skipping import/factor/solve" << std::endl; - } else { - std::cout << GridLogMessage << "DenseCoarseMatrix: slab cache " << slabfile - << " absent -- full setup, will write it" << std::endl; - } - } - if (!loaded) { + { ImportDense(Op); // slab <- my rows of A (LOCAL, no comms) ImportCertificate(Op); // dense apply == Op.M, before inversion InvertDense(Op); // slab <- my rows of A^{-1} double t1 = usecond(); - if (sfile) { - FILE *f = fopen(slabfile.c_str(),"wb"); - GRID_ASSERT(f != nullptr); - int64_t hdr[4] = { (int64_t)0x44454E5345, N, (int64_t)nrows, (int64_t)nbasis }; - GRID_ASSERT( fwrite(hdr,sizeof(int64_t),4,f) == 4 ); - uint64_t nelem = (uint64_t)nrows * N; - GRID_ASSERT( fwrite(&slab[0], sizeof(ComplexF), nelem, f) == nelem ); - fclose(f); - std::cout << GridLogMessage << "DenseCoarseMatrix: slab written to " << slabfile << std::endl; - } std::cout << GridLogMessage << "DenseCoarseMatrix: import+invert took " << (t1-t0)/1.0e6 << " s" << std::endl; } @@ -233,11 +180,9 @@ public: dSlab.resize((uint64_t)nrows*N); acceleratorCopyToDevice(&slab[0],&dSlab[0],sbytes); - // DENSE_SPLITK: requested chunk count, snapped DOWN to a divisor of N. - int req = getenv("DENSE_SPLITK") ? atoi(getenv("DENSE_SPLITK")) : 32; - if (req < 1) req = 1; + // Split-K chunk count, snapped DOWN to a divisor of N. NK = 1; - for(int j=1;j<=req;j++) if ( (N % j) == 0 ) NK = j; + for(int j=1;j<=SPLITK;j++) if ( (N % j) == 0 ) NK = j; int64_t Kc = N / NK; dX.resize((uint64_t)N*MRHS_MAX); @@ -255,12 +200,7 @@ public: for(int j=0;j=0 && devSum<=4); - 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. @@ -282,8 +222,8 @@ public: acceleratorCopyToDevice(&rm2g[0], &dRm2G[0], N*sizeof(int64_t)); } std::cout << GridLogMessage << "DenseCoarseMatrix: slab resident on device (" - << sbytes/1024./1024. << " MB/rank), split-K NK=" << NK << " (Kc=" << Kc << "); " - << sumName[devSum] << std::endl; + << sbytes/1024./1024. << " MB/rank), split-K NK=" << NK << " (Kc=" << Kc + << "); DEVICE cartesian ring ALLGATHER (P2P)" << std::endl; } //////////////////////////////////////////////////////////////////// @@ -355,8 +295,6 @@ public: { double t = -usecond(); Coordinate gdims = grid->GlobalDimensions(); - int sign = getenv("DENSE_IMPORT_SIGN") ? atoi(getenv("DENSE_IMPORT_SIGN")) : 1; - GRID_ASSERT( sign==1 || sign==-1 ); uint64_t nelem = (uint64_t)nrows * N; thread_for(i, nelem, { slab[i] = ComplexF(0.0,0.0); }); @@ -367,28 +305,11 @@ public: // extract the unpadded field before peeking with unpadded coordinates // (exactly as MultiGeneralCoarsenedMatrix::CopyMatrix does). CoarseMatrix Aun(grid); Op.ExtractMatrix(p,Aun); - if ( getenv("DENSE_IMPORT_DEBUG") ) { - // Peek-path vs field-norm audit: sum |peekLocalSite|^2 must match norm2 - double pk = 0.0; - autoView(Adbg, Aun, CpuRead); - for(int ss=0; ssGlobalSumVector(&gpk, 1); - std::cout << GridLogMessage << "DenseCoarseMatrix: DEBUG p=" << p - << " norm2(_A[p]) " << norm2(Aun) - << " norm2(Extract) " << norm2(Aun) - << " sum|peek|^2 " << gpk << std::endl; - } autoView(Av, Aun, CpuRead); thread_for(ss, lsites, { Coordinate ncoor(nd); for(int d=0; d_lstart[d] + myLcoor[ss][d] + sign*shift[d]; + int64_t g = grid->_lstart[d] + myLcoor[ss][d] + shift[d]; ncoor[d] = (int)((g % gdims[d] + gdims[d]) % gdims[d]); } int64_t nsite; @@ -402,8 +323,8 @@ public: // column (nbr,a). BUG LEDGER 2026-08-14: the original mapping // wrote (s,a),(nbr,b) -- caught by the IMPORT CERTIFICATE on its // FIRST fresh-import exercise (Test_schur_dense_coarse); every - // production slab predates this path (probe-import SLAB_FILEs), - // so no production output is suspect. + // production slab predated this path (probe-import era), so no + // production output is suspect. for(int b=0; b= 1.0e-3 ) { std::cout << GridLogMessage << "DenseCoarseMatrix: IMPORT CERTIFICATE FAILED. If O(1), the " - << "stencil shift-sign convention is opposite: rerun with DENSE_IMPORT_SIGN=-1" + << "stencil shift-sign convention of the coarse operator has changed: " + << "the import in ImportDense/ImportDenseFP64 must change with it" << std::endl; } GRID_ASSERT(rel < 1.0e-3); } //////////////////////////////////////////////////////////////////// - // 3. Invert dispatcher. slab holds my rows of A on entry, my rows - // of A^{-1} on exit, in every mode. - // DENSE_SCHUR absent/0 : single-GCD gather-invert (the oracle) - // DENSE_SCHUR=1 : distributed recursive Schur - // DENSE_SCHUR=2 : AUDIT -- both on the same A; report the - // slab difference; keep the Schur result - // (so VERIFY certifies the new path). - //////////////////////////////////////////////////////////////////// - template - void InvertDense(CoarseOp &Op) - { - char *sc = getenv("DENSE_SCHUR"); - int mode = sc ? atoi(sc) : 0; - - if ( mode == 0 ) - { - InvertDenseSingle(); - return; - } - if ( mode == 1 ) - { - InvertDenseSchur(Op); - return; - } - GRID_ASSERT( mode == 2 ); - std::vector Aimp(slab); // imported A - InvertDenseSingle(); - std::vector ref(slab); // Ainv, single path - slab = Aimp; - InvertDenseSchur(Op); // slab = Ainv, Schur path - - // NaN-PROOF comparison: max() masks NaN, so count non-finite - // entries in each result explicitly. - double mx = 0.0; - double mr = 0.0; - int64_t badschur = 0; - int64_t badsingle = 0; - for(uint64_t i=0; i<(uint64_t)nrows*N; i++) - { - double as = abs(ComplexD(slab[i])); - double ar = abs(ComplexD(ref[i])); - if ( !std::isfinite(as) ) badschur++; - if ( !std::isfinite(ar) ) badsingle++; - if ( std::isfinite(as) && std::isfinite(ar) ) - { - mx = std::max(mx, (double)abs(ComplexD(slab[i]) - ComplexD(ref[i]))); - mr = std::max(mr, ar); - } - } - RealD gmx = mx; - RealD gmr = mr; - RealD gbs = (RealD)badschur; - RealD gbr = (RealD)badsingle; - grid->GlobalMax(gmx); - grid->GlobalMax(gmr); - grid->GlobalSumVector(&gbs, 1); - grid->GlobalSumVector(&gbr, 1); - schurAuditRel = gmx/gmr; - std::cout << GridLogMessage << "DenseCoarseMatrix: DENSE_SCHUR=2 AUDIT " - << "max|Ainv_schur - Ainv_single| = " << gmx - << " relative " << schurAuditRel - << " non-finite: schur " << (int64_t)gbs << " single " << (int64_t)gbr - << " (two fp32 roundings of the same inverse; expect ~ growth * eps32)" - << std::endl; - GRID_ASSERT( gbs == 0 ); - GRID_ASSERT( gbr == 0 ); - } - - //////////////////////////////////////////////////////////////////// - // 3a. Single-GCD invert: chunked zero-fill+GlobalSum gather of A - // streamed to the boss GCD, cgetrf_64 (ILP64), rows of A^{-1} by - // blocked identity cgetrs_64 + broadcast; each rank keeps its own - // rows (in `slab`, overwriting A). Proven path; the SCHUR oracle. - //////////////////////////////////////////////////////////////////// - void InvertDenseSingle(void) - { - double t1 = usecond(); - int boss = grid->IsBoss(); - std::vector Afull; -#ifdef GRID_HIP - // QUARANTINED naked HIP: the boss-only N^2 inversion buffer (34GB at - // N=65536) must come from raw HBM; EvictAll flushes the device-copy - // layer to make the window. (FreePool of the allocator free-list - // awaits the type-dispatched fix.) Confined to setup; the apply path - // is pure Grid primitives. - rocblas_float_complex *dA = nullptr; - rocblas_float_complex *dB = nullptr; - int64_t *dIpiv = nullptr; - uint64_t Abytes = (uint64_t)N * N * sizeof(ComplexF); - MemoryManager::EvictAll(); - if (boss) { - auto aerr = hipMalloc((void **)&dA, Abytes); - if (aerr != hipSuccess) { - std::cout << GridLogMessage << "DenseCoarseMatrix: hipMalloc of " - << Abytes/1024./1024./1024. << " GB FAILED -- reduce --device-mem" << std::endl; - GRID_ASSERT(aerr == hipSuccess); - } - std::cout << GridLogMessage << "DenseCoarseMatrix: device inversion buffer allocated (" - << Abytes/1024./1024./1024. << " GB)" << std::endl; - } -#else - if (boss) Afull.resize((uint64_t)N * N); -#endif - { - std::unordered_map rowmap; // global row -> my slab row - for(int ss=0; ss chunk((uint64_t)CHUNKROWS * N); - for(int64_t row0=0; row0second) * N; - uint64_t dst = (uint64_t)(r-row0) * N; - for(int64_t j=0;jGlobalSumVector(&chunk[0], (int)nelem); - if (boss) { -#ifdef GRID_HIP - GRID_ASSERT( hipMemcpy((char *)dA + (uint64_t)row0*N*sizeof(ComplexF), - &chunk[0], nelem*sizeof(ComplexF), - hipMemcpyHostToDevice) == hipSuccess ); -#else - uint64_t dst = (uint64_t)row0 * N; - for(uint64_t i=0;i LU of A^T. - //////////////////////////////////////////////////////////////////// - if (boss) { -#ifdef GRID_HIP - std::cout << GridLogMessage << "DenseCoarseMatrix: rocSOLVER cgetrf_64 (ILP64 LU) N=" << N - << " in place on resident device buffer" << std::endl; - rocblas_handle handle = GridBLASInverse::Handle(); - int64_t *dInfo; - 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); - GRID_ASSERT( hipDeviceSynchronize() == hipSuccess ); - int64_t info_h = -1; - 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); - 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 - // A^{-1}, in exactly the linear layout the harvest expects. -#else - // Eigen fallback: small local CPU tests only. - std::cout << GridLogMessage << "DenseCoarseMatrix: Eigen fallback inversion N=" << N - << (N > 10000 ? " (WARNING: SLOW; use the HIP/rocSOLVER path)" : "") - << std::endl; - typedef Eigen::Matrix,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> MatF; - Eigen::Map A(reinterpret_cast*>(&Afull[0]), N, N); - MatF Ainv = A.inverse(); - A = Ainv; -#endif - } - double t3 = usecond(); - std::cout << GridLogMessage << "DenseCoarseMatrix: factorisation took " - << (t3-t2)/1.0e6 << " s" << std::endl; - - //////////////////////////////////////////////////////////////////// - // Blocked solve + broadcast: rows of A^{-1} chunk by chunk; each - // rank keeps the rows of its own sites (ownership-aligned). - //////////////////////////////////////////////////////////////////// - { - std::unordered_map rowmap; - for(int ss=0; ss chunk((uint64_t)CHUNKROWS * N); - for(int64_t row0=0; row0Broadcast(0, &chunk[0], nelem*sizeof(ComplexF)); - for(int64_t r=row0; rsecond) * N; - uint64_t src = (uint64_t)(r-row0) * N; - for(int64_t j=0;j rank-major column map, computed LOCALLY. + // 3a. Global column -> rank-major column map, computed LOCALLY. // Rank-major ordering: rank q's rows/columns are the contiguous // block [q*nrows, (q+1)*nrows), ordered by q's local site index // (uniform local volumes make ownership arithmetic exact). @@ -789,7 +475,7 @@ public: } //////////////////////////////////////////////////////////////////// - // 3c. Direct stencil -> fp64 rank-major import of MY ROWS of A (the + // 3b. Direct stencil -> fp64 rank-major import of MY ROWS of A (the // end-to-end fp64 path: the stencil source IS ComplexD; nothing is // rounded through fp32 on the way into the inversion). Same // loop/sign/accumulate/transposed-contraction discipline as @@ -805,8 +491,6 @@ public: void ImportDenseFP64(CoarseOp &Op, BlockRows &S, std::vector &g2rm) { Coordinate gdims = grid->GlobalDimensions(); - int sign = getenv("DENSE_IMPORT_SIGN") ? atoi(getenv("DENSE_IMPORT_SIGN")) : 1; - GRID_ASSERT( sign==1 || sign==-1 ); std::vector h((uint64_t)nrows*N, ComplexD(0.0,0.0)); for(int p=0; p_lstart[d] + myLcoor[ss][d] + sign*shift[d]; + int64_t g = grid->_lstart[d] + myLcoor[ss][d] + shift[d]; ncoor[d] = (int)((g % gdims[d] + gdims[d]) % gdims[d]); } int64_t nsite; @@ -868,14 +552,14 @@ public: } //////////////////////////////////////////////////////////////////// - // 3d. Distributed recursive Schur invert, END-TO-END fp64 (decision + // 3c. The inverse: distributed recursive Schur, END-TO-END fp64 (decision // 2026-08-14): stencil (ComplexD) -> fp64 rank-major import -> // fp64 recursion -> ONE terminal rounding into the fp32 apply // slab. Everything downstream (device residency, split-K apply, - // VERIFY, SLAB_FILE) is untouched. + // VERIFY) is untouched. //////////////////////////////////////////////////////////////////// template - void InvertDenseSchur(CoarseOp &Op) + void InvertDense(CoarseOp &Op) { double t1 = usecond(); int P = grid->ProcessorCount(); @@ -905,26 +589,17 @@ public: ImportDenseFP64(Op, S, g2rm); //////////////////////////////////////////////////////////////// - // DENSE_SCHUR2D=1 : invert via the 2D block-cyclic recursion - // (BlockCyclicSchurInverse) instead of the 1D rank-range one. - // The SAME imported rank-major rows S go in and come back, so the - // import certificate above and the slab rounding / VERIFY below are - // identical for both paths: a clean A/B on one imported operator. - // - // Everything in the 2D path -- redistribution, SUMMA rings, leaf -- - // is point-to-point SendToRecvFrom; no collectives at all. - // DENSE_NB overrides the block size (default: rows-per-rank, which - // makes the redistribution edges maximally regular). + // The 2D block-cyclic recursion (BlockCyclicSchurInverse). + // Everything in it -- redistribution, SUMMA rings, leaves -- is + // point-to-point SendToRecvFrom; no collectives at all. Block size + // nb = rows-per-rank makes the redistribution edges maximally + // regular. //////////////////////////////////////////////////////////////// - int use2d = getenv("DENSE_SCHUR2D") ? atoi(getenv("DENSE_SCHUR2D")) : 0; - int64_t panelBytes = getenv("DENSE_PANEL_BYTES") ? atol(getenv("DENSE_PANEL_BYTES")) - : (int64_t)1024*1024*1024; // 1D path only double t2, t3; - if ( use2d ) { int Pr,Pc; BlockCyclicLayout::ChooseProcessGrid(P, Pr, Pc); - int64_t nb = getenv("DENSE_NB") ? atol(getenv("DENSE_NB")) : nrows; + int64_t nb = nrows; GRID_ASSERT( nb >= 1 ); std::cout << GridLogMessage << "DenseCoarseMatrix: 2D SCHUR invert, process grid " << Pr << " x " << Pc << " nb " << nb @@ -938,14 +613,6 @@ public: t3 = usecond(); RSI2.ReportTelemetry(grid); } - else - { - RecursiveSchurInverse RSI(grid, N, rowStart, panelBytes); - t2 = usecond(); - RSI.Invert(S); - t3 = usecond(); - RSI.ReportTelemetry(); - } // The single terminal rounding: fp64 inverse -> fp32 apply slab // (row-major, global columns) @@ -962,28 +629,26 @@ public: } double t4 = usecond(); std::cout << GridLogMessage << "DenseCoarseMatrix: SCHUR fp64 distributed invert took " - << (t4-t1)/1.0e6 << " s (recursion " << (t3-t2)/1.0e6 << " s), panelBytes " - << panelBytes << std::endl; + << (t4-t1)/1.0e6 << " s (recursion " << (t3-t2)/1.0e6 << " s)" << std::endl; } //////////////////////////////////////////////////////////////////// // CORE apply on packed data: hX[N x nr] (zero-filled, my sites only) - // -> allreduce -> split-K GEMM against the resident slab -> reduce - // partials -> hY[nrows x nr] (column major). Platform-agnostic: + // -> ring allgather -> split-K GEMM against the resident slab -> + // reduce partials -> hY[nrows x nr] (column major). Platform-agnostic: // deviceVector + GridBLAS (Eigen fallback on CPU builds). - // fp32 allreduce is EXACT: zero-fill assembly gives every element - // exactly one contributing rank. + // tprof (optional): per-phase microseconds {allgather, H2D, gemm+reduce, + // D2H}, printed by the caller on GridLogPerformance. //////////////////////////////////////////////////////////////////// void SlabApplyPacked(int nr, double *tprof) { GRID_TRACE("DenseSlabApply"); GRID_ASSERT(nr <= MRHS_MAX); - uint64_t nX = (uint64_t)N * nr; uint64_t nY = (uint64_t)nrows * nr; int64_t Kc = N / NK; double t1 = usecond(); double t2, t3; - if (devSum==4) { + { // ALLGATHER: x is not a reduction -- every rank owns the rows of x at // global columns myGsite[ss]*nbasis+b (scattered by site coordinate, NOT // a contiguous block) and needs all of it. Only MY rows go host->device @@ -1011,31 +676,6 @@ public: }); } t3 = usecond(); - } else if (devSum) { - { GRID_TRACE("DenseH2D"); - acceleratorCopyToDevice(&hX[0],&dX[0],nX*sizeof(ComplexF)); - } - t2 = usecond(); - { GRID_TRACE("DenseAllreduce"); - // 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 - // DENSE_DEVICE_SUM=4 : CartesianRingAllGather (branch above) - 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 { - { GRID_TRACE("DenseAllreduce"); - grid->GlobalSumVector(&hX[0], (int)nX); - } - t2 = usecond(); - { GRID_TRACE("DenseH2D"); - acceleratorCopyToDevice(&hX[0],&dX[0],nX*sizeof(ComplexF)); - } - t3 = usecond(); } // Y = op(slab,T) . X : row-major slab (nrows x N) == col-major A^T // (N x nrows, lda=N) => transpose gives the nrows x N operator. @@ -1065,8 +705,8 @@ public: } double t5 = usecond(); if (tprof) { - tprof[0] = devSum ? (t3-t2) : (t2-t1); // allreduce - tprof[1] = devSum ? (t2-t1) : (t3-t2); // H2D + tprof[0] = t3-t2; // allgather + tprof[1] = t2-t1; // H2D tprof[2] = t4-t3; // gemm+reduce tprof[3] = t5-t4; // D2H } @@ -1208,16 +848,14 @@ public: double t6 = usecond(); std::cout << GridLogMessage << "DenseCoarseMatrix: apply6D " << nr << " rhs took " << (t6-t0)/1000.0 << " ms" << std::endl; - if ( getenv("DENSE_APPLY_PROFILE") ) { - std::cout << GridLogMessage << "DenseCoarseMatrix: apply6D profile:" - << " pack " << (t1-t0)/1000.0 - << (devSum==4 ? " allgather " : " allreduce ") << tprof[0]/1000.0 - << " H2D " << tprof[1]/1000.0 - << " gemm+reduce "<< tprof[2]/1000.0 - << " D2H " << tprof[3]/1000.0 - << " unpack " << (t6-t5)/1000.0 - << " ms" << std::endl; - } + std::cout << GridLogPerformance << "DenseCoarseMatrix: apply6D profile:" + << " pack " << (t1-t0)/1000.0 + << " allgather " << tprof[0]/1000.0 + << " H2D " << tprof[1]/1000.0 + << " gemm+reduce "<< tprof[2]/1000.0 + << " D2H " << tprof[3]/1000.0 + << " unpack " << (t6-t5)/1000.0 + << " ms" << std::endl; } }; diff --git a/Grid/algorithms/multigrid/RecursiveSchurInverse.h b/Grid/algorithms/multigrid/RecursiveSchurInverse.h deleted file mode 100644 index dd62d4f0c..000000000 --- a/Grid/algorithms/multigrid/RecursiveSchurInverse.h +++ /dev/null @@ -1,874 +0,0 @@ -/************************************************************************************* - - Grid physics library, www.github.com/paboyle/Grid - - Source file: RecursiveSchurInverse.h - - 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. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License along - with this program; if not, write to the Free Software Foundation, Inc., - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - - See the full license in the file "LICENSE" in the top level distribution directory -*************************************************************************************/ -/* END LEGAL */ -#pragma once - -#include -#include - -#include - -NAMESPACE_BEGIN(Grid); - -/////////////////////////////////////////////////////////////////////////////// -// RecursiveSchurInverse: distributed dense inversion by recursive Schur -// complement over a binary rank-range tree. -// -// Contract: rank r owns global rows [rowStart[r], rowStart[r+1]) of an -// N x N matrix in rank-major ordering, and receives its rows of the -// inverse in the same layout. Consumes GridBase collectives, GridBLAS -// and GridBLASInverse only; Eigen reference backends permit CPU unit -// testing under mpirun (Test_schur_inverse). -// -// Arithmetic is fp64 throughout; the caller rounds once into fp32 storage. -// -// Execution: SPMD full-tree walk. Every rank makes the identical call -// sequence; data participation is ownership-gated; all collectives are -// world-wide, so no deadlock surface exists. -// -// Storage: BlockRows is column-major, ld = rows; element (i,j) at -// data[i + j*ld]; a column window is the contiguous slice at data[col0*ld]. -/////////////////////////////////////////////////////////////////////////////// - -/////////////////////////////////////////////////////////////////////////////// -// My rows of a distributed dense matrix: rows x cols, column major, ld = rows. -/////////////////////////////////////////////////////////////////////////////// -class BlockRows -{ -public: - deviceVector data; - int64_t rows; - int64_t cols; - int64_t ld; - - BlockRows() - { - rows = 0; - cols = 0; - ld = 0; - } - void Resize(int64_t r, int64_t c) - { - rows = r; - cols = c; - ld = r; - data.resize((uint64_t)r*c); - } - ComplexD *ColumnWindow(int64_t col0) - { - GRID_ASSERT( col0 >= 0 ); - GRID_ASSERT( col0 <= cols ); - return &data[(uint64_t)col0*ld]; - } -}; - -class RecursiveSchurInverse -{ -public: - GridBase *grid; // world collectives only - int64_t N; // global matrix dimension - int P; // ranks - int me; // this rank - std::vector rowStart; // P+1 entries: rank-major row ownership - int64_t myRow0; - int64_t myNrows; - int64_t panelBytes; // gather panel budget (DENSE_PANEL_BYTES) - - GridBLAS BLAS; - GridBLASInverse INV; - - // Growth telemetry: one entry per merge node, walk order - std::vector telNormB; // ||B||_F = ||A11inv A12||_F - std::vector telSratio; // ||S||_F / ||A22||_F - double telLeafMaxInv; // max |(leaf inverse)_ij| over leaves - - // Phase timers/counters, reported by ReportTelemetry - double tMemset; // device panel zero-fill - double tDeposit; // owner rows -> panel (device kernel) - double tAllreduce; // panel collective (see tBarrier) - double tBarrier; // arrival skew, when DENSE_BARRIER_PROBE - double tRepack; // rank-major -> panel (gather path) - double tGemm; // strided gemm + synchronise - double tLeaf; // leaf inversions - double tARmin; // fastest single panel collective - double tARmax; // slowest single panel collective - std::vector tARall; // every panel collective, for percentiles - uint64_t bytesAllreduce; - uint64_t nAllreduce; // panel collectives - uint64_t nGatherGemm; // GatherGemm calls - uint64_t nGather; // gather-path collectives (debug gate) - uint64_t nGatherV; // ... of which MPI_Allgatherv - uint64_t nBcast; // ... of which Bcast-assembled - - // Persistent grow-only device panel; assembly and collectives are - // device-resident. Device builds require GPU-aware MPI. - deviceVector dPanelBuf; - // Rank-major receive staging for the AllGatherV path, plus the per-panel-row - // owner maps that drive the repack. Grow-only, same discipline as dPanelBuf: - // a fresh device allocation per collective is catastrophic (measured). - deviceVector dRecvBuf; - deviceVector dOwnerOff; // panel row -> owner's first panel row - deviceVector dOwnerRows; // panel row -> owner's row count - - // DENSE_GATHER : transport for the panel assembly. - // 0 = zero-fill + GlobalSumVector (default; moves the payload twice) - // 1 = MPI_Allgatherv (KNOWN BROKEN here, see guard below) - // 2 = C sequential MPI_Bcast, one per contributing rank. Bcast has no - // count vector, so the zero-count shape that breaks Allgatherv - // cannot arise. Costs C collectives instead of 1 and the roots do - // not transmit concurrently, so the byte cost is ~2*panel per rank - // -- the same as the allreduce. It wins only if Bcast is faster - // PER BYTE than Allreduce, which is why it must be measured. - // DENSE_GATHER_MIN_BYTES : panels below this stay on the allreduce path. - // Default 0 (always gather). The Schur tree puts - // 94% of its bytes in panels with ~3.7 MB per-rank - // messages, i.e. squarely bandwidth bound; only the - // two deepest levels (1.2% of bytes, 16-65 KB per - // 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 : 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). - // Static and communicator-free so synthetic tables unit-test directly. - /////////////////////////////////////////////////////////////////////////// - static void CheckRowStart(const std::vector &table, int64_t N) - { - int P = (int)table.size() - 1; - GRID_ASSERT( P >= 1 ); - GRID_ASSERT( table[0] == 0 ); - GRID_ASSERT( table[P] == N ); - for(int r=0; r= table[r] ); // zero-row ranks permitted - } - } - - /////////////////////////////////////////////////////////////////////////// - // Build the ownership table from each rank's local row count: zero-fill - // allgather (the standing comms idiom) then prefix sum. Every rank - // returns the identical table. - /////////////////////////////////////////////////////////////////////////// - static std::vector MakeRowStart(GridBase *g, int64_t myNrows) - { - int P = g->ProcessorCount(); - int me = g->ThisRank(); - - std::vector counts(P, 0); - counts[me] = (uint64_t)myNrows; - g->GlobalSumVector(&counts[0], P); - - std::vector table(P+1); - table[0] = 0; - for(int r=0; r &rowStart_, - int64_t panelBytes_) - { - grid = g; - N = N_; - P = g->ProcessorCount(); - me = g->ThisRank(); - rowStart = rowStart_; - panelBytes = panelBytes_; - - GRID_ASSERT( (int)rowStart.size() == P+1 ); - CheckRowStart(rowStart, N); - - myRow0 = rowStart[me]; - myNrows = rowStart[me+1] - rowStart[me]; - - telLeafMaxInv = 0.0; - - 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") ? 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==1) && !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"); - } - } - - /////////////////////////////////////////////////////////////////////////// - // The communication primitive. - // - // C(:, colC : colC+widthB) <- beta * C(:, colC : colC+widthB) - // + alpha * A(:, colA : colA+widthA) * Bsub - // - // Bsub is the widthA x widthB sub-block of a row-distributed operand - // owned by ranks [rB0, rB1): owner r contributes its rows of - // B(:, colB : colB+widthB) at sub-block row offset - // rowStart[r] - rowStart[rB0], gathered in panelBytes row-chunks by - // device zero-fill + deposit kernel + GlobalSumVector. - // - // Every rank calls; non-owners of B add zeros; ranks with A.rows == 0 - // skip local compute but make every collective call. Column offsets - // are local buffer offsets -- non-participants pass 0. - /////////////////////////////////////////////////////////////////////////// - void GatherGemm(ComplexD alpha, - BlockRows &A, int64_t colA, int64_t widthA, - int rB0, int rB1, - BlockRows &B, int64_t colB, int64_t widthB, - ComplexD beta, - BlockRows &C, int64_t colC) - { - GRID_ASSERT( rB0 >= 0 ); - GRID_ASSERT( rB1 > rB0 ); - GRID_ASSERT( rB1 <= P ); - - int64_t k = rowStart[rB1] - rowStart[rB0]; - int64_t m = A.rows; - int64_t n = widthB; - GRID_ASSERT( widthA == k ); - GRID_ASSERT( n >= 1 ); - - int owner = ( me >= rB0 ) && ( me < rB1 ) && ( B.rows > 0 ); - int64_t myOff = 0; - if ( owner ) - { - myOff = rowStart[me] - rowStart[rB0]; - } - - if ( m > 0 ) - { - GRID_ASSERT( colA + widthA <= A.cols ); - GRID_ASSERT( colC + widthB <= C.cols ); - GRID_ASSERT( C.rows == m ); - } - - nGatherGemm++; - GRID_TRACE("GatherGemm"); - - if ( owner ) - { - GRID_ASSERT( colB + widthB <= B.cols ); - } - - // COLUMN chunking, not row chunking. Each rank's contribution to a - // column chunk is rows_r x nchunk at ld = B.rows -- i.e. exactly the - // contiguous window B.ColumnWindow(colB+j0) -- so the AllGatherV send - // needs no pack. Row chunking would slice each column and force one. - // It also makes every chunk a full-k product, so beta applies directly - // and the cross-chunk accumulation disappears. - int64_t bufs = useGather ? 2 : 1; // panel, plus receive staging - int64_t nc = panelBytes / ( bufs * (int64_t)sizeof(ComplexD) * k ); - if ( nc < 1 ) nc = 1; - if ( nc > n ) nc = n; - GRID_ASSERT( k*nc < 2147483647L ); // collective counts are int - - deviceVector &dPanel = dPanelBuf; - if ( dPanel.size() < (uint64_t)k*nc ) dPanel.resize((uint64_t)k*nc); - deviceVector ap(1); - deviceVector bp(1); - deviceVector cp(1); - std::vector ptr(1); - - // Panel-row -> owner maps for the repack. Depend only on [rB0,rB1), - // so they are built once per call rather than once per chunk. - if ( useGather ) - { - if ( dRecvBuf.size() < (uint64_t)k*nc ) dRecvBuf.resize((uint64_t)k*nc); - if ( dOwnerOff.size() < (uint64_t)k ) dOwnerOff.resize((uint64_t)k); - if ( dOwnerRows.size() < (uint64_t)k ) dOwnerRows.resize((uint64_t)k); - std::vector hoff(k), hrows(k); - for(int r=rB0; r= gatherMinBytes ) ? useGather : 0; - - // Arrival skew is charged to the barrier, so that tAllreduce measures - // transfer alone. Diagnostic only; off by default. - if ( barrierProbe ) { double tb = -usecond(); grid->Barrier(); tBarrier += tb+usecond(); } - - double tar = -usecond(); - if ( gatherThis ) - { - ////////////////////////////////////////////////////////////////// - // AllGatherV: move the payload once, no arithmetic, no zero fill. - // Receive is rank major -- a concatenation of rows_r x nchunk - // column-major blocks -- which is the correct ROW order but the - // wrong LAYOUT, so one device repack follows. - ////////////////////////////////////////////////////////////////// - std::vector counts(P,0), displs(P,0); - for(int r=rB0; rAllGatherV(send, counts[me], - (void *)&dRecvBuf[0], counts, displs, sizeof(ComplexD)); - } - else - { - //////////////////////////////////////////////////////////////// - // C sequential broadcasts into the SAME rank-major staging - // buffer, so the repack below is shared with the AllGatherV arm. - // The root must already hold its own block: B.ColumnWindow is - // contiguous rows_r x nchunk, and so is its slot in dRecvBuf. - // Every rank issues all C broadcasts in the same order, so - // collective matching stays positional and safe. - //////////////////////////////////////////////////////////////// - nBcast++; - if ( owner ) - { - int64_t off_me = rowStart[me] - rowStart[rB0]; - acceleratorCopyDeviceToDevice((void *)B.ColumnWindow(colB+j0), - (void *)&dRecvBuf[off_me*nchunk], - (uint64_t)B.rows*nchunk*sizeof(ComplexD)); - } - for(int r=rB0;rBroadcast(r,(void *)&dRecvBuf[off_r*nchunk], - (uint64_t)rows_r*nchunk*sizeof(ComplexD)); - } - } - tar += usecond(); - - tRepack -= usecond(); - { - ComplexD *pan = &dPanel[0]; - ComplexD *rcv = &dRecvBuf[0]; - int64_t *ooff = &dOwnerOff[0]; - int64_t *orows = &dOwnerRows[0]; - int64_t kk = k, ncw = nchunk; - accelerator_for(idx, panelWords, 1, { - int64_t j = idx / kk; - int64_t p = idx - j*kk; - int64_t o = ooff[p]; - int64_t rr= orows[p]; - pan[idx] = rcv[(uint64_t)(o*ncw + (p-o) + j*rr)]; - }); - } - tRepack += usecond(); - } - else - { - ////////////////////////////////////////////////////////////////// - // Zero fill + deposit + GlobalSumVector. Exact because the zero - // fill leaves exactly one contributing rank per element, but it - // moves the payload twice and reduces over zeros. - ////////////////////////////////////////////////////////////////// - tar += usecond(); - tMemset -= usecond(); - acceleratorMemSet(&dPanel[0], 0, panelBytesThis); - tMemset += usecond(); - if ( owner ) - { - int64_t brows = B.rows; - int64_t kk = k; - int64_t dof = myOff; - ComplexD *src = B.ColumnWindow(colB+j0); - ComplexD *dst = &dPanel[0]; - tDeposit -= usecond(); - accelerator_for(idx, (uint64_t)(brows*nchunk), 1, { - int64_t j = idx / brows; - int64_t i = idx - j*brows; - dst[(uint64_t)(dof + i + j*kk)] = src[(uint64_t)(i + j*brows)]; - }); - tDeposit += usecond(); - } - tar = -usecond(); - grid->GlobalSumVector(&dPanel[0], (int)panelWords); - tar += usecond(); - } - tAllreduce += tar; - tARmin = std::min(tARmin, tar); - tARmax = std::max(tARmax, tar); - tARall.push_back(tar); - bytesAllreduce += panelBytesThis; - nAllreduce++; - - if ( m > 0 ) - { - // Full-k product into this column chunk of C: beta applies directly. - ptr[0] = A.ColumnWindow(colA); - acceleratorCopyToDevice(&ptr[0], &ap[0], sizeof(ComplexD*)); - ptr[0] = &dPanel[0]; - acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexD*)); - ptr[0] = C.ColumnWindow(colC+j0); - acceleratorCopyToDevice(&ptr[0], &cp[0], sizeof(ComplexD*)); - - tGemm -= usecond(); - BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N, - (int)m, (int)nchunk, (int)k, - alpha, ap, (int)A.ld, - bp, (int)k, - beta, cp, (int)C.ld); - BLAS.synchronise(); - tGemm += usecond(); - } - } - } - - /////////////////////////////////////////////////////////////////////////// - // Local Frobenius norm-squared of a full-height column window. - // NO comms; callers GlobalSum the result. Host staging, setup-scale. - /////////////////////////////////////////////////////////////////////////// - double FrobNorm2Local(BlockRows &X, int64_t col0, int64_t w) - { - if ( X.rows == 0 ) return 0.0; - GRID_ASSERT( col0 + w <= X.cols ); - uint64_t len = (uint64_t)X.rows*w; - std::vector h(len); - acceleratorCopyFromDevice(X.ColumnWindow(col0), &h[0], len*sizeof(ComplexD)); - // Member real()/imag(): portable across std::complex (CPU) and - // thrust::complex (HIP), where std::norm does not resolve. - double s = 0.0; - for(uint64_t i=0; i bp(1); - std::vector ptr(1); - ptr[0] = Arows.ColumnWindow(col0); - acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexD*)); - INV.inverseBatched(w, bp); - - // Telemetry: max |element| of the leaf inverse - { - std::vector h(len); - acceleratorCopyFromDevice(Arows.ColumnWindow(col0), &h[0], len*sizeof(ComplexD)); - double mx = 0.0; - for(uint64_t i=0; i A11inv - // 2. B = A11inv.A12 (I rows) - // 3. C = A21.A11inv (J rows) - // 4. S = A22 - A21.B in place (J rows) [alpha=-1, beta=1] - // 5. recurse J: S -> Sinv - // 6. T = Sinv.C (J rows) - // 7. U = B.Sinv (I rows) - // 8. X11 = A11inv + U.C in place (I rows) [beta=1] - // 9. X12 = -U, X21 = -T local negates; X22 = Sinv already in place - /////////////////////////////////////////////////////////////////////////// - void SchurNode(int r0, int r1, int64_t col0, int64_t width, BlockRows &Arows) - { - int span = r1 - r0; - GRID_ASSERT( span >= 1 ); - GRID_ASSERT( width == rowStart[r1] - rowStart[r0] ); - - if ( span == 1 ) - { - if ( ( me == r0 ) && ( myNrows > 0 ) ) - { - LeafInvert(col0, width, Arows); - } - return; - } - - int mid = ( r0 + r1 ) / 2; - int64_t WI = rowStart[mid] - rowStart[r0]; - int64_t WJ = rowStart[r1] - rowStart[mid]; - - // Zero-width child ranges (all ranks of a half owning no rows) are a - // KNOWN LIMITATION: fail loudly rather than divide mysteriously. - GRID_ASSERT( WI > 0 ); - GRID_ASSERT( WJ > 0 ); - - int inI = ( me >= r0 ) && ( me < mid ); - int inJ = ( me >= mid ) && ( me < r1 ); - - ComplexD one ( 1.0,0.0); - ComplexD mone (-1.0,0.0); - ComplexD zero ( 0.0,0.0); - - BlockRows dummy; - - // 1. A11 -> A11inv - SchurNode(r0, mid, col0, WI, Arows); - - // 2. B = A11inv . A12 (I rows; gather A12 from I owners) - BlockRows Bbuf; - if ( inI ) Bbuf.Resize(myNrows, WJ); - { - BlockRows &Aop = inI ? Arows : dummy; - BlockRows &Cop = inI ? Bbuf : dummy; - int64_t cA = inI ? col0 : 0; - GatherGemm(one, Aop, cA, WI, - r0, mid, - Arows, col0+WI, WJ, - zero, Cop, 0); - } - double nB = FrobNorm2Local(Bbuf, 0, inI ? WJ : 0); - grid->GlobalSumVector(&nB, 1); - telNormB.push_back(std::sqrt(nB)); - - // 3. C = A21 . A11inv (J rows; gather A11inv from I owners) - BlockRows Cbuf; - if ( inJ ) Cbuf.Resize(myNrows, WI); - { - BlockRows &Aop = inJ ? Arows : dummy; - BlockRows &Cop = inJ ? Cbuf : dummy; - int64_t cA = inJ ? col0 : 0; - GatherGemm(one, Aop, cA, WI, - r0, mid, - Arows, col0, WI, - zero, Cop, 0); - } - - // 4. S = A22 - A21 . B in place on my A22 window (J rows) - double nA22 = FrobNorm2Local( inJ ? Arows : dummy, inJ ? col0+WI : 0, inJ ? WJ : 0 ); - grid->GlobalSumVector(&nA22, 1); - { - BlockRows &Aop = inJ ? Arows : dummy; - BlockRows &Cop = inJ ? Arows : dummy; - int64_t cA = inJ ? col0 : 0; - int64_t cC = inJ ? col0+WI : 0; - GatherGemm(mone, Aop, cA, WI, - r0, mid, - Bbuf, 0, WJ, - one, Cop, cC); - } - double nS = FrobNorm2Local( inJ ? Arows : dummy, inJ ? col0+WI : 0, inJ ? WJ : 0 ); - grid->GlobalSumVector(&nS, 1); - telSratio.push_back( std::sqrt(nS) / ( std::sqrt(nA22) + 1.0e-300 ) ); - - // 5. S -> Sinv - SchurNode(mid, r1, col0+WI, WJ, Arows); - - // 6. T = Sinv . C (J rows; gather C from J owners) - BlockRows Tbuf; - if ( inJ ) Tbuf.Resize(myNrows, WI); - { - BlockRows &Aop = inJ ? Arows : dummy; - BlockRows &Cop = inJ ? Tbuf : dummy; - int64_t cA = inJ ? col0+WI : 0; - GatherGemm(one, Aop, cA, WJ, - mid, r1, - Cbuf, 0, WI, - zero, Cop, 0); - } - - // 7. U = B . Sinv (I rows; gather Sinv from J owners) - BlockRows Ubuf; - if ( inI ) Ubuf.Resize(myNrows, WJ); - { - BlockRows &Aop = inI ? Bbuf : dummy; - BlockRows &Cop = inI ? Ubuf : dummy; - GatherGemm(one, Aop, 0, WJ, - mid, r1, - Arows, col0+WI, WJ, - zero, Cop, 0); - } - - // 8. X11 = A11inv + U . C in place (I rows; gather C from J owners) - { - BlockRows &Aop = inI ? Ubuf : dummy; - BlockRows &Cop = inI ? Arows : dummy; - int64_t cC = inI ? col0 : 0; - GatherGemm(one, Aop, 0, WJ, - mid, r1, - Cbuf, 0, WI, - one, Cop, cC); - } - - // 9. Off-diagonal signs, local - if ( inI ) NegateCopy(Ubuf, Arows, col0+WI, WJ); - if ( inJ ) NegateCopy(Tbuf, Arows, col0, WI); - } - - /////////////////////////////////////////////////////////////////////////// - // PUBLIC ENTRY. Arows: my rows of the rank-major N x N matrix (fp64). - // On exit Arows holds my rows of the inverse, still fp64; the caller - // owns the single terminal rounding into its fp32 apply storage. - /////////////////////////////////////////////////////////////////////////// - void Invert(BlockRows &Arows) - { - GRID_ASSERT( Arows.rows == myNrows ); - GRID_ASSERT( Arows.cols == N ); - - telNormB.resize(0); - telSratio.resize(0); - telLeafMaxInv = 0.0; - - tMemset = 0.0; - tDeposit = 0.0; - tAllreduce = 0.0; - tBarrier = 0.0; - tRepack = 0.0; - tGemm = 0.0; - tLeaf = 0.0; - tARmin = 1.0e30; - tARmax = 0.0; - tARall.clear(); - bytesAllreduce = 0; - nAllreduce = 0; - nGatherGemm = 0; - nGather = 0; - nGatherV = 0; - nBcast = 0; - - SchurNode(0, P, 0, N, Arows); - - RealD mx = telLeafMaxInv; - grid->GlobalMax(mx); - telLeafMaxInv = mx; - } - - // All telemetry values are globally reduced or boss-local; safe to - // stream on every rank (Grid quiesces stdout to the boss unless - // --debug-stdout). NOTE: tAllreduce INCLUDES wait/imbalance -- a rank - // arriving early books its wait here; the min/max spread across ranks - // separates true wire time (~min) from skew (max-min). - void ReportTelemetry(void) - { - for(uint64_t i=0; iGlobalMax(armax); - grid->GlobalMax(armin); - armin = -armin; - - std::cout << GridLogMessage << "Schur phases (boss rank, seconds):" - << " memset " << tMemset/1.0e6 - << " deposit " << tDeposit/1.0e6 - << " allreduce " << tAllreduce/1.0e6 - << " barrier " << tBarrier/1.0e6 - << " repack " << tRepack/1.0e6 - << " gemm " << tGemm/1.0e6 - << " leaf " << tLeaf/1.0e6 - << std::endl; - std::cout << GridLogMessage << "Schur comms:" - << " [transports run: allreduce " << (nAllreduce - nGatherV - nBcast) - << " allgatherv " << nGatherV << " bcast " << nBcast << "]" - << ( barrierProbe ? " [barrier probe on]" : "" ) - << " GatherGemm calls " << nGatherGemm - << " panel allreduces " << nAllreduce - << " allreduce GB " << bytesAllreduce/1024./1024./1024. - << " allreduce s min/max over ranks " << armin/1.0e6 - << " / " << armax/1.0e6 - << std::endl; - std::cout << GridLogMessage << "Schur comms per-call (boss):" - << " min " << (nAllreduce ? tARmin/1.0e3 : 0.0) << " ms" - << " avg " << (nAllreduce ? tAllreduce/nAllreduce/1.0e3 : 0.0) << " ms" - << " max " << tARmax/1.0e3 << " ms" - << " 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 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 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; - } - } -}; - -NAMESPACE_END(Grid); diff --git a/Grid/parallelIO/BinaryIO.h b/Grid/parallelIO/BinaryIO.h index d586472a2..3d9d8e240 100644 --- a/Grid/parallelIO/BinaryIO.h +++ b/Grid/parallelIO/BinaryIO.h @@ -266,9 +266,10 @@ class BinaryIO { // GRID_BINARYIO_NOAGGREGATE falls back to plain lexicographic I/O. static int DefaultControl(void) { - static int ctrl = getenv("GRID_BINARYIO_NOAGGREGATE") - ? BINARYIO_LEXICOGRAPHIC - : BINARYIO_LEXICOGRAPHIC|BINARYIO_AGGREGATE; + // static int ctrl = getenv("GRID_BINARYIO_NOAGGREGATE") + // ? BINARYIO_LEXICOGRAPHIC + // : BINARYIO_LEXICOGRAPHIC|BINARYIO_AGGREGATE; + static int ctrl = BINARYIO_LEXICOGRAPHIC|BINARYIO_AGGREGATE; return ctrl; } diff --git a/tests/debug/Test_allgather.cc b/tests/debug/Test_allgather.cc index 2709450fa..5b79a24c4 100644 --- a/tests/debug/Test_allgather.cc +++ b/tests/debug/Test_allgather.cc @@ -244,7 +244,7 @@ int main(int argc, char **argv) // T7 : the SAME shape as T6, but assembled by C sequential MPI_Bcast -- // one broadcast per contributing rank -- instead of one MPI_Allgatherv. // - // This is the transport of DENSE_GATHER=2. Bcast takes no count vector, + // This is the transport of the retired chunked-Bcast gather. Bcast takes no count vector, // so the zero-count asymmetry that makes T6 run at ~0.18 MB/s and trip // mpir_request.h:508 cannot arise. It costs C collectives rather than 1 // and the roots do not transmit concurrently, so the byte cost is about diff --git a/tests/debug/Test_schur2d.cc b/tests/debug/Test_schur2d.cc index 36553b609..e82bfc06b 100644 --- a/tests/debug/Test_schur2d.cc +++ b/tests/debug/Test_schur2d.cc @@ -49,8 +49,8 @@ using namespace Grid; static int failures = 0; // Portable |z|: ComplexD is std::complex on CPU builds and thrust::complex -// under HIP, where std::abs does not resolve (same trap RecursiveSchurInverse -// documents at FrobNorm2Local). Member real()/imag() work on both. +// under HIP, where std::abs does not resolve. Member real()/imag() work +// on both. static double Cabs(const ComplexD &z) { double re = z.real(), im = z.imag(); diff --git a/tests/debug/Test_schur2d_redist.cc b/tests/debug/Test_schur2d_redist.cc index 784fe2c56..2a9c43330 100644 --- a/tests/debug/Test_schur2d_redist.cc +++ b/tests/debug/Test_schur2d_redist.cc @@ -23,7 +23,7 @@ Author: Peter Boyle // // 1D rank-major rows -> block cyclic -> Invert -> back to 1D rows // -// which is exactly what DENSE_SCHUR2D runs inside DenseCoarseMatrix. +// which is exactly what the dense inverse runs inside DenseCoarseMatrix. // CPU build under mpirun at n = 1,2,3,4. // // T1 : RowsToCyclic against a direct ImportGlobal of the same matrix -- @@ -31,14 +31,9 @@ Author: Peter Boyle // T2 : round trip rows -> 2D -> rows -- BITWISE, uniform AND non-uniform // rowStart, layouts with ragged trailing blocks. // T3 : full pipeline inverse against a host Gauss-Jordan reference. -// T4 : CROSS-IMPLEMENTATION: the same matrix inverted by the 1D -// RecursiveSchurInverse and by the 2D pipeline; results compared -// element-wise. Two independent implementations, two independent -// decompositions, one answer. ////////////////////////////////////////////////////////////////////////////// #include -#include #include #include @@ -47,8 +42,8 @@ using namespace Grid; static int failures = 0; // Portable |z|: ComplexD is std::complex on CPU builds and thrust::complex -// under HIP, where std::abs does not resolve (same trap RecursiveSchurInverse -// documents at FrobNorm2Local). Member real()/imag() work on both. +// under HIP, where std::abs does not resolve. Member real()/imag() work +// on both. static double Cabs(const ComplexD &z) { double re = z.real(), im = z.imag(); @@ -190,12 +185,11 @@ int main(int argc, char **argv) } //////////////////////////////////////////////////////////////////////// - // T3 + T4 : the DENSE_SCHUR2D pipeline against the host reference and - // against the INDEPENDENT 1D RecursiveSchurInverse. + // T3 : the 2D pipeline against the host reference. //////////////////////////////////////////////////////////////////////// { - bool ok3 = true, ok4 = true; - double worst3 = 0.0, worst4 = 0.0; + bool ok3 = true; + double worst3 = 0.0; BlockCyclicSchurInverse RSI2; for(auto &g : grids){ for(auto &c : cfgs){ @@ -232,27 +226,9 @@ int main(int argc, char **argv) worst3 = std::max(worst3,d); if ( d > 1.0e-9 ) ok3 = false; } - - // ---- 1D RecursiveSchurInverse on the same matrix ---- - { - BlockRows Ar; Ar.Resize(myrows, N); - acceleratorCopyToDevice(&h[0], &Ar.data[0], h.size()*sizeof(ComplexD)); - std::vector rs = rowStart; - RecursiveSchurInverse RSI1(grid, N, rs, 1<<20); - RSI1.Invert(Ar); - std::vector h1d(h.size()); - acceleratorCopyFromDevice(&Ar.data[0], &h1d[0], h1d.size()*sizeof(ComplexD)); - for(int64_t j=0;j 1.0e-9 ) ok4 = false; - } - } } } Report("T3 2D pipeline vs host reference", ok3, "worst "+std::to_string(worst3)); - Report("T4 2D pipeline vs 1D RecursiveSchurInverse", ok4, "worst "+std::to_string(worst4)); } { diff --git a/tests/debug/Test_schur2d_scale.cc b/tests/debug/Test_schur2d_scale.cc index 1e05696d7..dad050a41 100644 --- a/tests/debug/Test_schur2d_scale.cc +++ b/tests/debug/Test_schur2d_scale.cc @@ -20,7 +20,7 @@ Author: Peter Boyle ////////////////////////////////////////////////////////////////////////////// // SCALE rehearsal for the 2D distributed dense inverse: the full -// DENSE_SCHUR2D pipeline -- 1D rows -> redistribute -> invert -> +// 2D block-cyclic pipeline -- 1D rows -> redistribute -> invert -> // redistribute back -> certificate -- on a SYNTHETIC matrix of any size, // with no multigrid machinery, no configuration and no subspace file. // @@ -31,8 +31,8 @@ Author: Peter Boyle // rank; at N=138240 on 288 ranks it is the production problem shape // exactly, in a driver that runs in minutes. // -// S2D_N : global dimension (default 720, laptop friendly) -// S2D_NB : block size (default N/P rows-per-rank if that +// --schur2d-global-dimension : N (default 720, laptop friendly) +// --schur2d-block-size : nb (default N/P rows-per-rank if that // is exact, else 48) // // The matrix is diagonally dominant (the recursion does not pivot); its @@ -55,8 +55,8 @@ Author: Peter Boyle using namespace Grid; // Portable |z|: ComplexD is std::complex on CPU builds and thrust::complex -// under HIP, where std::abs does not resolve (same trap RecursiveSchurInverse -// documents at FrobNorm2Local). Member real()/imag() work on both. +// under HIP, where std::abs does not resolve. Member real()/imag() work +// on both. static double Cabs(const ComplexD &z) { double re = z.real(), im = z.imag(); @@ -76,13 +76,6 @@ static ComplexD Fill(int64_t i, int64_t j, int64_t N) int main(int argc, char **argv) { - // Environment walk (2026-08-27, systems/Frontier/schur2d_env.job): knobs to - // reproduce the example's environment here, one at a time. Result: thread - // level, OMP_NUM_THREADS, device residency and sustained load all NIL; only - // "first job step on fresh nodes" (+3 s) is real. - // S2D_BALLAST_GB=x x GB of Lattice fields made device-resident before the invert - // S2D_PREHEAT_S=x x seconds of back-to-back zgemm before the invert - // OMP_NUM_THREADS set in the job, read by nothing here but the runtime Grid_init(&argc, &argv); GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), @@ -91,9 +84,12 @@ int main(int argc, char **argv) const int P = grid->ProcessorCount(); const int me = grid->ThisRank(); - int64_t N = getenv("S2D_N") ? atol(getenv("S2D_N")) : 720; + int64_t N = 720; + if ( GridCmdOptionExists(argv,argv+argc,"--schur2d-global-dimension") ) + N = atol(GridCmdOptionPayload(argv,argv+argc,"--schur2d-global-dimension").c_str()); int64_t nb; - if ( getenv("S2D_NB") ) nb = atol(getenv("S2D_NB")); + if ( GridCmdOptionExists(argv,argv+argc,"--schur2d-block-size") ) + nb = atol(GridCmdOptionPayload(argv,argv+argc,"--schur2d-block-size").c_str()); else if ( N % P == 0 ) nb = N/P; else nb = 48; GRID_ASSERT( N >= 1 ); GRID_ASSERT( nb >= 1 ); @@ -125,56 +121,13 @@ int main(int argc, char **argv) double t1 = usecond(); //////////////////////////////////////////////////////////////////////// - // The DENSE_SCHUR2D pipeline, phase-timed. A0 keeps the original for + // The 2D pipeline, phase-timed. A0 keeps the original for // the certificate. //////////////////////////////////////////////////////////////////////// BlockCyclicMatrix A (grid,N,nb,Pr,Pc); BlockCyclicMatrix A0(grid,N,nb,Pr,Pc); BlockCyclicSchurInverse RSI2; - // Pre-heat: drive the GCD with back-to-back zgemm for S2D_PREHEAT_S seconds - // before the invert. The example calls the inverse after ~100 s of full - // load on all 288 GCDs and its LOCAL kernels run 20-40% slower than the - // idle-start harness (GEMM 3.1 vs 2.5 s, leaf 0.76 vs 0.24 s) with the - // wires unchanged; thread level / OMP / residency (E1-E5) did not reproduce - // that. If sustained load does, it is clock/power management, not code. - if ( getenv("S2D_PREHEAT_S") ) { - double secs = atof(getenv("S2D_PREHEAT_S")); - const int64_t W = 4320; - deviceVector M((uint64_t)W*W), C((uint64_t)W*W); - { ComplexD *m = &M[0]; accelerator_for(idx,(uint64_t)W*W,1,{ m[idx] = ComplexD(1.0e-3*(idx%97),1.0e-3*(idx%89)); }); accelerator_barrier(); } - deviceVector ap(1),bp(1),cp(1); std::vector ptr(1); - ptr[0]=&M[0]; acceleratorCopyToDevice(&ptr[0],&ap[0],sizeof(ComplexD*)); acceleratorCopyToDevice(&ptr[0],&bp[0],sizeof(ComplexD*)); - ptr[0]=&C[0]; acceleratorCopyToDevice(&ptr[0],&cp[0],sizeof(ComplexD*)); - double t0=usecond(); int n=0; double tlast=0; - while ( (usecond()-t0)/1.0e6 < secs ) { - double t1=usecond(); - RSI2.SUMMA.BLAS.gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N,(int)W,(int)W,(int)W,ComplexD(1.0,0.0),ap,(int)W,bp,(int)W,ComplexD(0.0,0.0),cp,(int)W); - RSI2.SUMMA.BLAS.synchronise(); tlast=usecond()-t1; n++; - } - double tfirst = 0; (void)tfirst; - std::cout << GridLogMessage << "Test_schur2d_scale: pre-heat " << (usecond()-t0)/1.0e6 << " s, " << n << " zgemm W=" << W - << ", last zgemm " << tlast/1.0e6 << " s (" << 8.0*W*W*W/tlast/1.0e6 << " TF/s; idle-start rate 23.7)" << std::endl; - } - - // Device ballast: Lattice fields written on the accelerator so they sit in - // the MemoryManager's device LRU exactly as the example's fine-grid state does. - typedef Lattice,Ns> > BallastField; - std::vector ballast; - if ( getenv("S2D_BALLAST_GB") ) { - double gb = atof(getenv("S2D_BALLAST_GB")); - uint64_t fbytes = (uint64_t)grid->oSites()*sizeof(BallastField::vector_object); - int nf = (int)(gb*1.0e9/(double)fbytes + 0.5); - ballast.reserve(nf); - for(int i=0;ioSites(), 1, { v[ss] = Zero(); }); - } - std::cout << GridLogMessage << "Test_schur2d_scale: device ballast " << nf << " fields x " << fbytes/1.0e6 - << " MB = " << nf*fbytes/1.0e9 << " GB resident (S2D_BALLAST_GB=" << gb << ")" << std::endl; - } - BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,A); double t2 = usecond(); if ( A.data.size() ) diff --git a/tests/debug/Test_schur2d_vs_slate.cc b/tests/debug/Test_schur2d_vs_slate.cc index 46765104d..9990aafe9 100644 --- a/tests/debug/Test_schur2d_vs_slate.cc +++ b/tests/debug/Test_schur2d_vs_slate.cc @@ -64,9 +64,10 @@ Author: Peter Boyle // 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. +// --schur2d-global-dimension, --schur2d-block-size as in +// Test_schur2d_scale (default nb = N/P). +// --schur2d-skip-getri skips the getri leg (host loop; ~4 min at N=138240). +// --schur2d-nowarm skips the warm-up. // A third leg, getrf+getrs(I), is SLATE's device-resident inverse route. ////////////////////////////////////////////////////////////////////////////// @@ -138,8 +139,12 @@ int main(int argc, char **argv) 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 ); + int64_t N = 720; + if ( GridCmdOptionExists(argv,argv+argc,"--schur2d-global-dimension") ) + N = atol(GridCmdOptionPayload(argv,argv+argc,"--schur2d-global-dimension").c_str()); + int64_t nb = (N%P==0) ? N/P : 48; + if ( GridCmdOptionExists(argv,argv+argc,"--schur2d-block-size") ) + nb = atol(GridCmdOptionPayload(argv,argv+argc,"--schur2d-block-size").c_str()); int Pr,Pc; BlockCyclicLayout::ChooseProcessGrid(P,Pr,Pc); std::vector rowStart(P+1); rowStart[0]=0; @@ -161,10 +166,10 @@ int main(int argc, char **argv) // first. Run a small throwaway inverse through BOTH paths so the timed // legs below measure hot code. Not reported. //////////////////////////////////////////////////////////////////////// - // S2D_NOWARM=1 skips it (hang localisation). Stage markers are flushed so + // --schur2d-nowarm skips it (hang localisation). Stage markers are flushed so // 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") ) { + if ( !GridCmdOptionExists(argv,argv+argc,"--schur2d-nowarm") ) { // Fixed tiny size independent of P: the purpose is handle creation and // kernel loading, not work. (8*P at P=288 was N=2304 -> a 122 s SLATE // warm-up dominated by 288-way tile broadcasts.) Ranks beyond the first @@ -207,7 +212,7 @@ int main(int argc, char **argv) #endif std::cout << GridLogMessage << "warm-up done (both paths, N=" << Nw << ")" << std::endl << std::flush; } else { - std::cout << GridLogMessage << "warm-up SKIPPED (S2D_NOWARM)" << std::endl << std::flush; + std::cout << GridLogMessage << "warm-up SKIPPED (--schur2d-nowarm)" << std::endl << std::flush; } //////////////////////////////////////////////////////////////////////// @@ -240,7 +245,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 + if ( !GridCmdOptionExists(argv,argv+argc,"--schur2d-skip-getri") ) { // 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); diff --git a/tests/debug/Test_schur_dense_coarse.cc b/tests/debug/Test_schur_dense_coarse.cc index 1e750bffe..506b30eb4 100644 --- a/tests/debug/Test_schur_dense_coarse.cc +++ b/tests/debug/Test_schur_dense_coarse.cc @@ -27,25 +27,19 @@ Author: Peter Boyle /* END LEGAL */ // -// T6 of the RecursiveSchurInverse regression chain -// (schur_recursive_inverse_plan.txt 4B.5): the DenseCoarseMatrix GLUE, -// on a real (tiny) lattice coarse operator, CPU laptop build. +// The DenseCoarseMatrix GLUE test, on a real (tiny) lattice coarse +// operator, CPU laptop build. // // Builds a genuine GeneralCoarsenedMatrix (DWF MdagM + 0.5 shift for a -// guaranteed-invertible Galerkin coarse op, random aggregation basis, -// nbasis=8, 4^4 x Ls/1 blocking) and constructs DenseCoarseMatrix in -// DENSE_SCHUR=2 AUDIT mode with small DENSE_PANEL_BYTES (multi-panel -// gathers exercised through the glue). The constructor then runs, in -// order, all the certificates this stage exists to check: -// - fresh ImportDense (no SLAB_FILE) + IMPORT CERTIFICATE vs Op.M -// - InvertDenseSingle (the oracle) -// - InvertDenseSchur: self-certifying rank-major map, fp64 diagonal -// import certificate vs the fp32 slab, distributed recursion, -// growth telemetry -// - AUDIT: max|Ainv_schur - Ainv_single| over the full slab -// - VERIFY ||A Ainv x - x||/||x|| through the SCHUR result -// This program adds asserts on the audit number and a random-vector -// round trip. +// guaranteed-invertible Galerkin coarse op, random aggregation basis) +// and runs the whole Import certificate chain through the glue: +// - fresh ImportDense + IMPORT CERTIFICATE vs Op.M +// - fp64 rank-major import certificate vs the fp32 slab +// - the 2D block-cyclic recursion + growth telemetry +// - VERIFY ||A Ainv x - x||/||x|| through the device split-K apply +// This program adds an INDEPENDENT Eigen fp64 host-inverse oracle at +// small N (built from applies of M to unit vectors, so it shares no +// code with the import) and a random-vector round trip. // // Uniform local volume 12.12.12.12 (fine), per-dim blocks {4,4,3,3}, // coarse 3.3.4.4/rank, nbasis 4 (N = 576n): @@ -193,13 +187,18 @@ int main (int argc, char ** argv) /////////////////////////////////////////////////////////////////////// // Full-matrix conditioning probe at small N: dense columns by - // applying M to unit vectors, fp64 Eigen SVD. + // applying M to unit vectors, fp64 Eigen SVD. eA is kept: it is the + // INDEPENDENT oracle for the inverse below (built from applies of M, + // sharing no code with the stencil->dense import). /////////////////////////////////////////////////////////////////////// + int64_t Nprobe = Coarse5d->gSites() * nbasis; + Eigen::MatrixXcd eA; + bool haveOracle = false; { - int64_t Nprobe = Coarse5d->gSites() * nbasis; if ( Nprobe <= 700 ) { - Eigen::MatrixXcd eA(Nprobe, Nprobe); + eA.resize(Nprobe, Nprobe); + haveOracle = true; CoarseVector e(Coarse5d); CoarseVector Me(Coarse5d); for(int64_t j=0; j, constructed - // on the grid and fed by Import (which runs the certificate chain). typedef DenseCoarseMatrix DenseCC; DenseCC dcm(Coarse5d); dcm.Import(LittleDiracOp); - std::cout << GridLogMessage << "T6 audit relative slab difference (schur vs single) = " - << dcm.schurAuditRel << std::endl; - GRID_ASSERT( dcm.schurAuditRel >= 0.0 ); // audit actually ran - GRID_ASSERT( dcm.schurAuditRel < 1.0e-3 ); - /////////////////////////////////////////////////////////////////////// - // Random-vector round trip through the SCHUR inverse + // Random-vector round trip through the inverse /////////////////////////////////////////////////////////////////////// CoarseVector x(Coarse5d); CoarseVector y(Coarse5d); @@ -288,6 +276,38 @@ int main (int argc, char ** argv) << rel << std::endl; GRID_ASSERT( rel < 1.0e-2 ); + /////////////////////////////////////////////////////////////////////// + // Independent oracle at small N: y from dcm must match the Eigen fp64 + // solve of eA (dense columns of M itself) on the same x. Catches an + // inverse that is self-consistent with a WRONG import, which the round + // trip above cannot (dense and M would share the error). + /////////////////////////////////////////////////////////////////////// + if ( haveOracle ) + { + Eigen::VectorXcd xs(Nprobe), ys(Nprobe); + typedef typename CoarseVector::vector_object::scalar_object csobj; + for(int64_t j=0; j_ndimension); + Lexicographic::CoorFromIndex(gcoor, gsite, Coarse5d->GlobalDimensions()); + csobj s; + peekSite(s, x, gcoor); + ComplexD zz = ((ComplexD *)&s)[b]; + xs(j) = std::complex(zz.real(), zz.imag()); + peekSite(s, y, gcoor); + zz = ((ComplexD *)&s)[b]; + ys(j) = std::complex(zz.real(), zz.imag()); + } + Eigen::VectorXcd yref = eA.fullPivLu().solve(xs); + double dev = (ys - yref).cwiseAbs().maxCoeff(); + double ymax = yref.cwiseAbs().maxCoeff(); + std::cout << GridLogMessage << "T6 oracle max|Ainv x - eigen solve|/max|y| = " + << dev/ymax << " (fp32 slab vs fp64 host solve)" << std::endl; + GRID_ASSERT( dev/ymax < 1.0e-4 ); + } + std::cout << GridLogMessage << "Test_schur_dense_coarse: T6 ALL PASS" << std::endl; Grid_finalize(); diff --git a/tests/debug/Test_schur_inverse.cc b/tests/debug/Test_schur_inverse.cc deleted file mode 100644 index e0a2651e8..000000000 --- a/tests/debug/Test_schur_inverse.cc +++ /dev/null @@ -1,769 +0,0 @@ -/************************************************************************************* - - Grid physics library, www.github.com/paboyle/Grid - - Source file: Test_schur_inverse.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. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License along - with this program; if not, write to the Free Software Foundation, Inc., - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - - See the full license in the file "LICENSE" in the top level distribution directory -*************************************************************************************/ -/* END LEGAL */ - -// -// Staged regression gate for RecursiveSchurInverse (distributed dense -// inversion by recursive Schur complement) -- the laptop-side certificate -// chain of schur_recursive_inverse_plan.txt section 4B.5. Runs on a -// CPU-only build (Eigen BLAS backends) under mpirun: -// -// mpirun -n 1 ./Test_schur_inverse --grid 8.8.8.8 --mpi 1.1.1.1 -// mpirun -n 2 ./Test_schur_inverse --grid 8.8.8.8 --mpi 1.1.1.2 -// mpirun -n 3 ./Test_schur_inverse --grid 8.8.8.12 --mpi 1.1.1.3 -// mpirun -n 4 ./Test_schur_inverse --grid 8.8.8.8 --mpi 1.1.1.4 -// -// (n=3 exercises uneven row splits throughout.) The lattice exists only -// to furnish the communicator; no field is ever constructed. -// -// PRECISION: the inversion runs ENTIRELY in fp64 (decision 2026-08-14, -// superseding the fp32-merge design); certificates are eps64-scaled. -// The single terminal fp32 rounding belongs to the caller (tested at the -// glue level, Test_schur_dense_coarse). -// -// Stages present (cumulative -- earlier tests are never removed): -// T1a : ownership tables -- CheckRowStart on synthetic uneven partitions, -// MakeRowStart allgather vs closed form on the live communicator. -// T1b : STORAGE-CONVENTION PIN -- column-major + ld + window-offset -// semantics fixed once via identity multiplies through the -// explicit-ld gemmBatched, on INTEGER-VALUED data so all three -// cases below are EXACT (values well within the mantissa): -// (1) alpha=1,beta=0 read from an input column window -// (2) alpha=-1,beta=1 accumulate (the S-formation case) -// (3) write INTO an output column window, neighbours untouched -// No later failure can be a transposition/convention ambiguity. -// T2 : GatherGemm vs naive fp64 oracle (owner sub-ranges, alpha-beta -// cases, tiny+huge panels, half-participation call shape). -// T3 : LeafInvert in-place residual certificate. -// T4 : full recursive Invert vs Eigen fp64 oracle, growth-scaled -// certification, adversarial near-singular-A11 family with -// telemetry-spike assertion. -// -// Hard asserts throughout; thresholds pre-registered in the plan. -// -#include -#include -#include - -using namespace std; -using namespace Grid; - -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); - - GridCartesian Comm(GridDefaultLatt(), - GridDefaultSimd(Nd,vComplex::Nsimd()), - GridDefaultMpi()); - GridBase *grid = &Comm; - - //////////////////////////////////////////////////////////////// - // T1a : ownership tables - //////////////////////////////////////////////////////////////// - { - // Synthetic partitions of N=97 (prime: every P>1 is uneven) - const int64_t N = 97; - for(int P=1; P<=4; P++) - { - std::vector table(P+1); - table[0] = 0; - for(int r=0; rProcessorCount(); - int me = grid->ThisRank(); - - int64_t myNrows = 3 + me; - std::vector table = RecursiveSchurInverse::MakeRowStart(grid, myNrows); - - std::vector expect(P+1); - expect[0] = 0; - for(int r=0; r ComplexD - { - return ComplexD( (RealD)(1 + i + 10*j), (RealD)(i - j) ); - }; - - BlockRows A; - A.Resize(rows, cols); - { - std::vector Ahost((uint64_t)rows*cols); - for(int64_t j=0; j Idev((uint64_t)w*w); - { - std::vector Ihost((uint64_t)w*w, ComplexD(0.0,0.0)); - for(int64_t d=0; d Ap(1); - deviceVector Bp(1); - deviceVector Cp(1); - std::vector ptr_h(1); - - auto setptr = [&](deviceVector &d, ComplexD *p) - { - ptr_h[0] = p; - acceleratorCopyToDevice(&ptr_h[0], &d[0], sizeof(ComplexD*)); - }; - - //////////////////////////////////////////////////////////// - // Case 1: C = A(:, col0:col0+w) . I_w (alpha=1, beta=0) - //////////////////////////////////////////////////////////// - { - deviceVector Cdev((uint64_t)rows*w); - setptr(Ap, A.ColumnWindow(col0)); - setptr(Bp, &Idev[0]); - setptr(Cp, &Cdev[0]); - - BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N, - (int)rows, (int)w, (int)w, - one, Ap, (int)A.ld, - Bp, (int)w, - zero, Cp, (int)rows); - BLAS.synchronise(); - - std::vector Chost((uint64_t)rows*w); - acceleratorCopyFromDevice(&Cdev[0], &Chost[0], (uint64_t)rows*w*sizeof(ComplexD)); - for(int64_t j=0; j ComplexD - { - return ComplexD( (RealD)(100 + i + j), (RealD)7 ); - }; - deviceVector Cdev((uint64_t)rows*w); - { - std::vector Chost((uint64_t)rows*w); - for(int64_t j=0; j Chost((uint64_t)rows*w); - acceleratorCopyFromDevice(&Cdev[0], &Chost[0], (uint64_t)rows*w*sizeof(ComplexD)); - for(int64_t j=0; j Chost((uint64_t)rows*ccols, ComplexD(-999.0, 999.0)); - acceleratorCopyToDevice(&Chost[0], &C.data[0], (uint64_t)rows*ccols*sizeof(ComplexD)); - } - setptr(Ap, A.ColumnWindow(col0)); - setptr(Bp, &Idev[0]); - setptr(Cp, C.ColumnWindow(cw0)); - - BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N, - (int)rows, (int)w, (int)w, - one, Ap, (int)A.ld, - Bp, (int)w, - zero, Cp, (int)C.ld); - BLAS.synchronise(); - - std::vector Chost((uint64_t)rows*ccols); - acceleratorCopyFromDevice(&C.data[0], &Chost[0], (uint64_t)rows*ccols*sizeof(ComplexD)); - for(int64_t j=0; j= cw0) && (j < cw0+w) ) - { - GRID_ASSERT( got == f(i, col0 + (j-cw0)) ); - } - else - { - GRID_ASSERT( got == ComplexD(-999.0, 999.0) ); - } - } - } - } - - std::cout << GridLogMessage - << "T1b storage-convention pin (window read / S-accumulate / window write, exact) PASS" - << std::endl; - } - - //////////////////////////////////////////////////////////////// - // T2 : GatherGemm vs naive double-precision oracle. - // - // Every rank generates the SAME full N x N random fp64 operands - // from a fixed seed (no comms needed for the oracle), keeps only - // its own rows in BlockRows form, and after each GatherGemm call - // checks its output window element-by-element against a plain - // triple-loop ComplexD accumulation over the same entries. - // - // Sweep: N in {8, 96, 97}; owner ranges full/upper-half/single; - // (alpha,beta) in {(1,0), (-1,1)}; panelBytes tiny (ragged - // many-chunk gathers) and huge (single panel). Sentinel columns - // outside the output window must be untouched. Finally, a - // HALF-PARTICIPATION case rehearses the recursion call pattern: - // lower ranks own B but pass EMPTY A/C (collectives only). - //////////////////////////////////////////////////////////////// - { - int P = grid->ProcessorCount(); - int me = grid->ThisRank(); - - std::mt19937 rng(777); - std::uniform_real_distribution dist(-1.0,1.0); - - const int64_t nout = 5; // output width - const int64_t colB = 3; // B window offset - const int64_t colC = 2; // C window offset - - for(int64_t N : {8L, 96L, 97L}) - { - // Ownership: uneven for any P not dividing N - std::vector table(P+1); - table[0] = 0; - for(int r=0; r Aglob((uint64_t)N*N); - std::vector Bglob((uint64_t)N*N); - for(uint64_t i=0; i<(uint64_t)N*N; i++) Aglob[i] = ComplexD(dist(rng),dist(rng)); - for(uint64_t i=0; i<(uint64_t)N*N; i++) Bglob[i] = ComplexD(dist(rng),dist(rng)); - - // My rows of a full-matrix operand as a BlockRows - auto fillRows = [&](BlockRows &X, std::vector &glob, - int64_t row0, int64_t nr) - { - X.Resize(nr, N); - if ( nr == 0 ) return; - std::vector h((uint64_t)nr*N); - for(int64_t j=0; j > ranges; - ranges.push_back(std::make_pair(0, P)); - if ( P > 1 ) ranges.push_back(std::make_pair(P/2, P)); - if ( P > 1 ) ranges.push_back(std::make_pair(1, 2)); - - for(auto range : ranges) - { - int rB0 = range.first; - int rB1 = range.second; - int64_t ka0 = table[rB0]; // A-column window start = B row span - int64_t k = table[rB1] - table[rB0]; - - for(int acase=0; acase<2; acase++) - { - ComplexD alpha = ( acase==0 ) ? ComplexD( 1.0,0.0) : ComplexD(-1.0,0.0); - ComplexD beta = ( acase==0 ) ? ComplexD( 0.0,0.0) : ComplexD( 1.0,0.0); - - for(int64_t panelBytes : {64L, 1L<<30}) - { - RecursiveSchurInverse RSI(grid, N, table, panelBytes); - - BlockRows A; - BlockRows B; - BlockRows C; - fillRows(A, Aglob, r0, myNr); - fillRows(B, Bglob, r0, myNr); - - // Output: sentinel-filled, window at colC - const ComplexD sentinel(-999.0, 999.0); - const int64_t ccols = colC + nout + 2; - C.Resize(myNr, ccols); - std::vector C0((uint64_t)myNr*ccols, sentinel); - if ( acase == 1 ) - { - // beta=1 needs defined window content: g(i,j), integer-valued - for(int64_t j=0; j 0 ) - { - acceleratorCopyToDevice(&C0[0], &C.data[0], (uint64_t)myNr*ccols*sizeof(ComplexD)); - } - - RSI.GatherGemm(alpha, A, ka0, k, - rB0, rB1, - B, colB, nout, - beta, C, colC); - - std::vector Chost((uint64_t)myNr*ccols); - if ( myNr > 0 ) - { - acceleratorCopyFromDevice(&C.data[0], &Chost[0], (uint64_t)myNr*ccols*sizeof(ComplexD)); - } - - double tol = 1.0e-14 * (double)k; - for(int64_t j=0; j= colC) && (j < colC+nout) ) - { - int64_t jj = j - colC; - ComplexD acc(0.0,0.0); - if ( acase == 1 ) - { - acc = C0[(uint64_t)(i + j*myNr)]; - } - for(int64_t t=0; t 1 ) - { - int ph = ( P+1 ) / 2; - int64_t ka0 = table[0]; - int64_t k = table[ph] - table[0]; - int participant = ( me >= ph ); - - RecursiveSchurInverse RSI(grid, N, table, 64); - - BlockRows A; - BlockRows B; - BlockRows C; - fillRows(B, Bglob, r0, myNr); - if ( participant ) - { - fillRows(A, Aglob, r0, myNr); - C.Resize(myNr, nout); - } - - ComplexD one (1.0,0.0); - ComplexD zero(0.0,0.0); - int64_t cA = participant ? ka0 : 0; - RSI.GatherGemm(one, A, cA, k, - 0, ph, - B, colB, nout, // owners deposit from their B window - zero, C, 0); - - if ( participant ) - { - std::vector Chost((uint64_t)myNr*nout); - acceleratorCopyFromDevice(&C.data[0], &Chost[0], (uint64_t)myNr*nout*sizeof(ComplexD)); - double tol = 1.0e-14 * (double)k; - for(int64_t j=0; jProcessorCount(); - int me = grid->ThisRank(); - - int64_t w = 17 + 3*me; - uint64_t len = (uint64_t)w*w; - - std::vector table = RecursiveSchurInverse::MakeRowStart(grid, w); - RecursiveSchurInverse RSI(grid, table[P], table, 1<<20); - - // A = w I + R : well conditioned - std::mt19937 rng(31 + me); - std::uniform_real_distribution dist(-1.0,1.0); - std::vector Ahost(len); - for(uint64_t i=0; i X(len); - acceleratorCopyFromDevice(&Ar.data[0], &X[0], len*sizeof(ComplexD)); - - double maxdev = 0.0; - for(int64_t j=0; j 10 when P > 1); at fp64 the certificate - // barely notices it: that insensitivity IS the point - // of the fp64 conversion. - // - // N=64 runs with panelBytes=128 (ragged many-chunk gathers inside - // the recursion); larger N with 1 MB panels. - //////////////////////////////////////////////////////////////// - { - int P = grid->ProcessorCount(); - int me = grid->ThisRank(); - - std::mt19937 rng(2026); - std::uniform_real_distribution dist(-1.0,1.0); - - for(int64_t N : {64L, 200L, 513L}) - { - std::vector table(P+1); - table[0] = 0; - for(int r=0; r Aglob((uint64_t)N*N); - for(uint64_t i=0; i<(uint64_t)N*N; i++) Aglob[i] = ComplexD(dist(rng),dist(rng)); - for(int64_t d=0; d 1 ) ? table[1] : N/4; - for(int64_t j=0; j std::complex - { - return std::complex(z.real(), z.imag()); - }; - Eigen::MatrixXcd eA(N,N); - for(int64_t j=0; j h((uint64_t)myNr*N); - for(int64_t j=0; j Xfull((uint64_t)N*N, ComplexD(0.0,0.0)); - { - std::vector h((uint64_t)myNr*N); - acceleratorCopyFromDevice(&Arows.data[0], &h[0], (uint64_t)myNr*N*sizeof(ComplexD)); - for(int64_t j=0; jGlobalSumVector(&Xfull[0], (int)(N*N)); - - // cert1 = ||A X - I||_max - double cert1 = 0.0; - for(int64_t j=0; j 1) ) - { - GRID_ASSERT( maxNormB > 10.0 ); // the spike must REGISTER - } - } - } - - std::cout << GridLogMessage - << "T4 recursive Invert vs Eigen oracle (N=64/200/513, 3 families) PASS" - << std::endl; - } - - std::cout << GridLogMessage - << "Test_schur_inverse: ALL STAGES PASS" << std::endl; - - Grid_finalize(); -} diff --git a/tests/debug/Test_summa.cc b/tests/debug/Test_summa.cc index 3068bd9ed..c72bcff22 100644 --- a/tests/debug/Test_summa.cc +++ b/tests/debug/Test_summa.cc @@ -50,8 +50,8 @@ using namespace Grid; static int failures = 0; // Portable |z|: ComplexD is std::complex on CPU builds and thrust::complex -// under HIP, where std::abs does not resolve (same trap RecursiveSchurInverse -// documents at FrobNorm2Local). Member real()/imag() work on both. +// under HIP, where std::abs does not resolve. Member real()/imag() work +// on both. static double Cabs(const ComplexD &z) { double re = z.real(), im = z.imag(); diff --git a/tests/solver/Test_multigrid_common.h b/tests/solver/Test_multigrid_common.h index c1570a381..7f4c14be0 100644 --- a/tests/solver/Test_multigrid_common.h +++ b/tests/solver/Test_multigrid_common.h @@ -37,9 +37,9 @@ Zero zero; // serialization strategy of Grid? // clang-format off -struct MultiGridParams : Serializable { +struct WilsonMGParams : Serializable { public: - GRID_SERIALIZABLE_CLASS_MEMBERS(MultiGridParams, + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonMGParams, int, nLevels, std::vector>, blockSizes, // size == nLevels - 1 std::vector, smootherTol, // size == nLevels - 1 @@ -54,7 +54,7 @@ public: int, coarseSolverMaxInnerIter); // constructor with default values - MultiGridParams(int _nLevels = 2, + WilsonMGParams(int _nLevels = 2, std::vector> _blockSizes = {{4, 4, 4, 4}}, std::vector _smootherTol = {1e-14}, std::vector _smootherMaxOuterIter = {4}, @@ -82,7 +82,7 @@ public: }; // clang-format on -void checkParameterValidity(MultiGridParams const ¶ms) { +void checkParameterValidity(WilsonMGParams const ¶ms) { auto correctSize = params.nLevels - 1; @@ -101,7 +101,7 @@ public: std::vector Grids; std::vector PRNGs; - LevelInfo(GridCartesian *FineGrid, MultiGridParams const &mgParams) { + LevelInfo(GridCartesian *FineGrid, WilsonMGParams const &mgParams) { auto nCoarseLevels = mgParams.blockSizes.size(); @@ -176,7 +176,7 @@ public: int _CurrentLevel; int _NextCoarserLevel; - MultiGridParams &_MultiGridParams; + WilsonMGParams &_WilsonMGParams; LevelInfo & _LevelInfo; FineDiracMatrix & _FineMatrix; @@ -201,10 +201,10 @@ public: // Member Functions ///////////////////////////////////////////// - MultiGridPreconditioner(MultiGridParams &mgParams, LevelInfo &LvlInfo, FineDiracMatrix &FineMat, FineDiracMatrix &SmootherMat) + MultiGridPreconditioner(WilsonMGParams &mgParams, LevelInfo &LvlInfo, FineDiracMatrix &FineMat, FineDiracMatrix &SmootherMat) : _CurrentLevel(mgParams.nLevels - (nCoarserLevels + 1)) // _Level = 0 corresponds to finest , _NextCoarserLevel(_CurrentLevel + 1) // incremented for instances on coarser levels - , _MultiGridParams(mgParams) + , _WilsonMGParams(mgParams) , _LevelInfo(LvlInfo) , _FineMatrix(FineMat) , _SmootherMatrix(SmootherMat) @@ -212,7 +212,7 @@ public: , _CoarseMatrix(*_LevelInfo.Grids[_NextCoarserLevel]) { _NextPreconditionerLevel - = std::unique_ptr(new NextPreconditionerLevel(_MultiGridParams, _LevelInfo, _CoarseMatrix, _CoarseMatrix)); + = std::unique_ptr(new NextPreconditionerLevel(_WilsonMGParams, _LevelInfo, _CoarseMatrix, _CoarseMatrix)); resetTimers(); } @@ -261,7 +261,7 @@ public: conformable(in, out); // TODO: implement a W-cycle - if(_MultiGridParams.kCycle) + if(_WilsonMGParams.kCycle) kCycle(in, out); else vCycle(in, out); @@ -279,13 +279,13 @@ public: FineVector fineTmp(in.Grid()); - auto maxSmootherIter = _MultiGridParams.smootherMaxOuterIter[_CurrentLevel] * _MultiGridParams.smootherMaxInnerIter[_CurrentLevel]; + auto maxSmootherIter = _WilsonMGParams.smootherMaxOuterIter[_CurrentLevel] * _WilsonMGParams.smootherMaxInnerIter[_CurrentLevel]; TrivialPrecon fineTrivialPreconditioner; - FlexibleGeneralisedMinimalResidual fineFGMRES(_MultiGridParams.smootherTol[_CurrentLevel], + FlexibleGeneralisedMinimalResidual fineFGMRES(_WilsonMGParams.smootherTol[_CurrentLevel], maxSmootherIter, fineTrivialPreconditioner, - _MultiGridParams.smootherMaxInnerIter[_CurrentLevel], + _WilsonMGParams.smootherMaxInnerIter[_CurrentLevel], false); MdagMLinearOperator fineMdagMOp(_FineMatrix); @@ -336,19 +336,19 @@ public: FineVector fineTmp(in.Grid()); - auto smootherMaxIter = _MultiGridParams.smootherMaxOuterIter[_CurrentLevel] * _MultiGridParams.smootherMaxInnerIter[_CurrentLevel]; - auto kCycleMaxIter = _MultiGridParams.kCycleMaxOuterIter[_CurrentLevel] * _MultiGridParams.kCycleMaxInnerIter[_CurrentLevel]; + auto smootherMaxIter = _WilsonMGParams.smootherMaxOuterIter[_CurrentLevel] * _WilsonMGParams.smootherMaxInnerIter[_CurrentLevel]; + auto kCycleMaxIter = _WilsonMGParams.kCycleMaxOuterIter[_CurrentLevel] * _WilsonMGParams.kCycleMaxInnerIter[_CurrentLevel]; TrivialPrecon fineTrivialPreconditioner; - FlexibleGeneralisedMinimalResidual fineFGMRES(_MultiGridParams.smootherTol[_CurrentLevel], + FlexibleGeneralisedMinimalResidual fineFGMRES(_WilsonMGParams.smootherTol[_CurrentLevel], smootherMaxIter, fineTrivialPreconditioner, - _MultiGridParams.smootherMaxInnerIter[_CurrentLevel], + _WilsonMGParams.smootherMaxInnerIter[_CurrentLevel], false); - FlexibleGeneralisedMinimalResidual coarseFGMRES(_MultiGridParams.kCycleTol[_CurrentLevel], + FlexibleGeneralisedMinimalResidual coarseFGMRES(_WilsonMGParams.kCycleTol[_CurrentLevel], kCycleMaxIter, *_NextPreconditionerLevel, - _MultiGridParams.kCycleMaxInnerIter[_CurrentLevel], + _WilsonMGParams.kCycleMaxInnerIter[_CurrentLevel], false); MdagMLinearOperator fineMdagMOp(_FineMatrix); @@ -581,7 +581,7 @@ public: int _CurrentLevel; - MultiGridParams &_MultiGridParams; + WilsonMGParams &_WilsonMGParams; LevelInfo & _LevelInfo; FineDiracMatrix &_FineMatrix; @@ -594,9 +594,9 @@ public: // Member Functions ///////////////////////////////////////////// - MultiGridPreconditioner(MultiGridParams &mgParams, LevelInfo &LvlInfo, FineDiracMatrix &FineMat, FineDiracMatrix &SmootherMat) + MultiGridPreconditioner(WilsonMGParams &mgParams, LevelInfo &LvlInfo, FineDiracMatrix &FineMat, FineDiracMatrix &SmootherMat) : _CurrentLevel(mgParams.nLevels - (0 + 1)) - , _MultiGridParams(mgParams) + , _WilsonMGParams(mgParams) , _LevelInfo(LvlInfo) , _FineMatrix(FineMat) , _SmootherMatrix(SmootherMat) { @@ -613,12 +613,12 @@ public: conformable(_LevelInfo.Grids[_CurrentLevel], in.Grid()); conformable(in, out); - auto coarseSolverMaxIter = _MultiGridParams.coarseSolverMaxOuterIter * _MultiGridParams.coarseSolverMaxInnerIter; + auto coarseSolverMaxIter = _WilsonMGParams.coarseSolverMaxOuterIter * _WilsonMGParams.coarseSolverMaxInnerIter; // On the coarsest level we only have what I above call the fine level, no coarse one TrivialPrecon fineTrivialPreconditioner; FlexibleGeneralisedMinimalResidual fineFGMRES( - _MultiGridParams.coarseSolverTol, coarseSolverMaxIter, fineTrivialPreconditioner, _MultiGridParams.coarseSolverMaxInnerIter, false); + _WilsonMGParams.coarseSolverTol, coarseSolverMaxIter, fineTrivialPreconditioner, _WilsonMGParams.coarseSolverMaxInnerIter, false); MdagMLinearOperator fineMdagMOp(_FineMatrix); @@ -651,7 +651,7 @@ using NLevelMGPreconditioner = MultiGridPreconditioner std::unique_ptr>> -createMGInstance(MultiGridParams &mgParams, LevelInfo &levelInfo, Matrix &FineMat, Matrix &SmootherMat) { +createMGInstance(WilsonMGParams &mgParams, LevelInfo &levelInfo, Matrix &FineMat, Matrix &SmootherMat) { #define CASE_FOR_N_LEVELS(nLevels) \ case nLevels: \ diff --git a/tests/solver/Test_wilson_mg.cc b/tests/solver/Test_wilson_mg.cc index 6ab2295e9..4e9c73670 100644 --- a/tests/solver/Test_wilson_mg.cc +++ b/tests/solver/Test_wilson_mg.cc @@ -51,7 +51,7 @@ int main(int argc, char **argv) { RealD mass = -0.25; - MultiGridParams mgParams; + WilsonMGParams mgParams; std::string inputXml{"./mg_params.xml"}; if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) { diff --git a/tests/solver/Test_wilson_mg_mp.cc b/tests/solver/Test_wilson_mg_mp.cc index 9161593cb..e4646cea0 100644 --- a/tests/solver/Test_wilson_mg_mp.cc +++ b/tests/solver/Test_wilson_mg_mp.cc @@ -58,7 +58,7 @@ int main(int argc, char **argv) { RealD mass = -0.25; - MultiGridParams mgParams; + WilsonMGParams mgParams; std::string inputXml{"./mg_params.xml"}; if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) { diff --git a/tests/solver/Test_wilsonclover_mg.cc b/tests/solver/Test_wilsonclover_mg.cc index 116a379c5..e528e60fd 100644 --- a/tests/solver/Test_wilsonclover_mg.cc +++ b/tests/solver/Test_wilsonclover_mg.cc @@ -54,7 +54,7 @@ int main(int argc, char **argv) { RealD csw_r = 1.0; RealD csw_t = 1.0; - MultiGridParams mgParams; + WilsonMGParams mgParams; std::string inputXml{"./mg_params.xml"}; if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) { diff --git a/tests/solver/Test_wilsonclover_mg_lime.cc b/tests/solver/Test_wilsonclover_mg_lime.cc index 12d84f4de..bbe6b7f58 100644 --- a/tests/solver/Test_wilsonclover_mg_lime.cc +++ b/tests/solver/Test_wilsonclover_mg_lime.cc @@ -84,7 +84,7 @@ int main(int argc, char **argv) { RealD csw_r = 1.0; RealD csw_t = 1.0; - MultiGridParams mgParams; + WilsonMGParams mgParams; std::string inputXml{"./mg_params.xml"}; if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) { diff --git a/tests/solver/Test_wilsonclover_mg_mp.cc b/tests/solver/Test_wilsonclover_mg_mp.cc index 02322a825..6e994417c 100644 --- a/tests/solver/Test_wilsonclover_mg_mp.cc +++ b/tests/solver/Test_wilsonclover_mg_mp.cc @@ -60,7 +60,7 @@ int main(int argc, char **argv) { RealD csw_r = 1.0; RealD csw_t = 1.0; - MultiGridParams mgParams; + WilsonMGParams mgParams; std::string inputXml{"./mg_params.xml"}; if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) {