diff --git a/tests/debug/Test_blockcyclic.cc b/tests/debug/Test_blockcyclic.cc new file mode 100644 index 000000000..f9e6388ad --- /dev/null +++ b/tests/debug/Test_blockcyclic.cc @@ -0,0 +1,280 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/debug/Test_blockcyclic.cc + + Copyright (C) 2026 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////// +// Regression gate for BlockCyclicLayout -- stage 1 of the 2D distributed +// dense inverse (documentation/DistributedDenseInverse2D.tex). +// +// The layout is pure index arithmetic, so this test is EXHAUSTIVE rather +// than statistical: every stage sweeps a battery of (N, nb, Pr, Pc) +// configurations chosen for their edge cases -- nb=1, nb>N, N%nb!=0, +// more processes than blocks, prime N, the production shape -- and checks +// every global index (T2/T3) or every element of an outer-product grid (T4) +// against a brute-force reference built by walking blocks. +// +// No communication: correct at mpirun -n 1, and identical at any rank count. +// +// T1 : NumLocal partitions N (sums over coords; against brute force). +// T2 : GlobalToLocal / LocalToGlobal round trip, every index. +// T3 : local indices are dense [0,mloc): a bijection, not just a cover. +// T4 : 2D ownership: every element has exactly one owner; per-rank counts +// equal mloc*nloc; LocalOffset is a bijection onto [0,mloc*nloc). +// T5 : block contiguity: within any owned global block, consecutive +// global rows are consecutive local rows (what SUMMA panels rely on). +// T6 : ChooseProcessGrid: exact factorisation, Pr<=Pc, most-square. +////////////////////////////////////////////////////////////////////////////// + +#include +#include + +using namespace Grid; + +static int failures = 0; + +static void Report(const std::string &name, bool pass, const std::string &detail="") +{ + std::cout << GridLogMessage << " " << name << (pass ? " PASS" : " ** FAIL **"); + if ( detail.size() ) std::cout << " " << detail; + std::cout << std::endl; + if ( !pass ) failures++; +} + +// One-dimensional configurations: {N, nb, Pg} +struct Cfg1d { int64_t N; int64_t nb; int Pg; }; +static const std::vector configs1d = { + { 0, 1, 1}, // empty matrix + { 1, 1, 1}, // minimal + { 16, 4, 4}, // exact tiling + { 17, 4, 4}, // trailing partial block + { 16, 1, 4}, // nb=1: pure cyclic + { 16, 32, 4}, // nb>N: one short block, most coords empty + { 16, 4, 7}, // more processes than blocks + { 97, 13, 5}, // prime N, awkward everything + {138240, 480, 16}, // production N, row dimension of 16x18 + {138240, 480, 18}, // production N, column dimension of 16x18 + {138240, 512, 16}, // production N, non-dividing block +}; + +int main(int argc, char **argv) +{ + Grid_init(&argc, &argv); + + std::cout << GridLogMessage << "BlockCyclicLayout regression (communicator-free index arithmetic)" << std::endl; + + //////////////////////////////////////////////////////////////////////// + // T1 : NumLocal partitions N, and agrees with brute force block-walking. + //////////////////////////////////////////////////////////////////////// + { + bool ok = true; + for(auto &c : configs1d){ + // brute force: walk global blocks, count elements per coordinate + std::vector counts(c.Pg, 0); + for(int64_t b=0; b*c.nb < c.N || (c.N==0 && b<0); b++){ + if ( b*c.nb >= c.N ) break; + int64_t lo = b*c.nb; + int64_t hi = std::min(c.N, lo+c.nb); + counts[b % c.Pg] += hi-lo; + } + int64_t sum = 0; + for(int p=0;p(c.N, 200000); // full N except production + for(int64_t g=0; g(0,c.N-1000); g LocalToGlobal round trip", ok); + } + + //////////////////////////////////////////////////////////////////////// + // T3 : per coordinate, the local indices hit [0, NumLocal) exactly once. + //////////////////////////////////////////////////////////////////////// + { + bool ok = true; + for(auto &c : configs1d){ + if ( c.N > 4096 ) continue; // dense bitmap check: small cfgs only + for(int p=0;p hit(n, 0); + for(int64_t g=0; g= n ) { ok = false; continue; } + hit[l]++; + } + for(int64_t l=0;l cfgs = { + { 16, 4, 2, 2}, + { 17, 4, 2, 3}, // partial block, non-square grid + { 23, 5, 3, 2}, + { 12, 2, 3, 4}, // 12 ranks + { 30, 7, 4, 2}, + }; + for(auto &c : cfgs){ + int P = c.Pr*c.Pc; + std::vector L; + for(int r=0;r owned(P,0); + std::vector> slot(P); + for(int r=0;r= P ) { ok=false; continue; } + // every layout instance must agree on the owner + if ( L[r].OwnerRank(i,j) != r ) ok = false; + if ( !L[r].Owns(i,j) ) ok = false; + owned[r]++; + int64_t off = L[r].LocalOffset(i,j); + if ( off < 0 || off >= (int64_t)slot[r].size() ) { ok=false; continue; } + slot[r][off]++; + } + } + int64_t tot=0; + for(int r=0;r(c.N-1, 100000); + for(int64_t g=0; g Pc ) ok = false; + // most-square: no divisor r with Pr < r <= sqrt(P) + for(int r=Pr+1; (int64_t)r*r <= (int64_t)P; r++) + if ( P % r == 0 ) ok = false; + } + int Pr,Pc; + BlockCyclicLayout::ChooseProcessGrid(288,Pr,Pc); + if ( !(Pr==16 && Pc==18) ) ok = false; + Report("T6 ChooseProcessGrid exact, Pr<=Pc, most-square (288 -> 16x18)", ok); + } + + //////////////////////////////////////////////////////////////////////// + // T7 : RangeToLocal. For every block-aligned range of every small + // configuration: the owned global indices of [g0,g1) map exactly onto + // local [l0,l1), contiguously and in order (brute force). + //////////////////////////////////////////////////////////////////////// + { + bool ok = true; + for(auto &c : configs1d){ + if ( c.N == 0 || c.N > 4096 ) continue; + int64_t nblocks = (c.N + c.nb - 1)/c.nb; + for(int p=0;p c.N ) continue; + int64_t l0,l1; + BlockCyclicLayout::RangeToLocal(g0,g1,c.N,c.nb,p,c.Pg,l0,l1); + // brute force: owned globals in [g0,g1) in ascending order + int64_t expect = l0; + for(int64_t g=g0; g + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////// +// Regression gate for BlockCyclicSchurInverse -- stage 3 of the 2D +// distributed dense inverse. CPU build under mpirun: +// +// mpirun -n 1 ./Test_schur2d --grid 8.8.8.8 --mpi 1.1.1.1 +// mpirun -n 2 ./Test_schur2d --grid 8.8.8.8 --mpi 1.1.1.2 +// mpirun -n 3 ./Test_schur2d --grid 8.8.8.12 --mpi 1.1.1.3 +// mpirun -n 4 ./Test_schur2d --grid 8.8.8.8 --mpi 1.1.1.4 +// +// Sweeps all process-grid factorisations of P and a battery of (N,nb) +// including ragged trailing blocks, a single-leaf matrix (nblocks==1), +// and nb=3 with many blocks. The matrices are diagonally dominant -- +// the recursion does not pivot, exactly like the 1D implementation, and +// the test respects that contract. +// +// T1 : certificate max|A . Ainv - I| with the product computed by the +// (independently validated) distributed SUMMA on an untouched copy. +// T2 : element-wise against a host Gauss-Jordan reference inverse +// (partial pivoting, fp64). +// T3 : repeated inversion bitwise identical -- the P2P determinism +// property, which a collective-reduce implementation cannot offer. +////////////////////////////////////////////////////////////////////////////// + +#include +#include + +using namespace Grid; + +static int failures = 0; + +static void Report(const std::string &name, bool pass, const std::string &detail="") +{ + std::cout << GridLogMessage << " " << name << (pass ? " PASS" : " ** FAIL **"); + if ( detail.size() ) std::cout << " " << detail; + std::cout << std::endl; + if ( !pass ) failures++; +} + +static ComplexD Fill(int64_t i, int64_t j, int salt) +{ + double x = std::sin(0.7*i + 1.3*j + 0.31*salt); + double y = std::cos(1.9*i - 0.4*j + 0.77*salt); + return ComplexD(x,y); +} + +// Diagonally dominant test matrix: no pivoting required at any depth. +static void MakeMatrix(std::vector &A, int64_t N, int salt) +{ + A.resize((uint64_t)N*N); + for(int64_t j=0;j A, std::vector &X, int64_t N) +{ + X.assign((uint64_t)N*N, ComplexD(0.0,0.0)); + for(int64_t i=0;i mx ){ mx=std::abs(A[r+c*N]); piv=r; } + GRID_ASSERT( mx > 0.0 ); + if ( piv != c ) + for(int64_t j=0;jProcessorCount(); + + std::vector> grids; + for(int r=1;r<=P;r++) if ( P%r==0 ) grids.push_back({r,P/r}); + + struct Cfg { int64_t N; int64_t nb; }; + std::vector cfgs = { {24,4}, {26,4}, {17,5}, {30,7}, {8,8}, {33,3}, {40,5} }; + + BlockCyclicSchurInverse RSI; + BlockCyclicSumma SUMMA; + + std::cout << GridLogMessage << "BlockCyclicSchurInverse regression: P=" << P + << ", " << grids.size() << " process grids, " << cfgs.size() + << " layouts" << std::endl; + + //////////////////////////////////////////////////////////////////////// + // T1 + T2, one pass: invert, certify with distributed SUMMA on an + // untouched copy, and compare against the host reference inverse. + //////////////////////////////////////////////////////////////////////// + { + bool okC = true, okR = true; + double worstC = 0.0, worstR = 0.0; + for(auto &g : grids){ + for(auto &c : cfgs){ + int64_t N = c.N; + std::vector Ag, Ref, Ainv, Cert; + MakeMatrix(Ag, N, 12); + HostInverse(Ag, Ref, N); + + BlockCyclicMatrix A (grid,N,c.nb,g.first,g.second); + BlockCyclicMatrix A0(grid,N,c.nb,g.first,g.second); + BlockCyclicMatrix Ce(grid,N,c.nb,g.first,g.second); + A.ImportGlobal(Ag); + A0.ImportGlobal(Ag); + + RSI.Invert(A); // in place: A now holds Ainv + + // certificate: Ce = A0 . Ainv, distributed + SUMMA.Multiply(ComplexD(1.0,0.0),A0,A,ComplexD(0.0,0.0),Ce, 0,N,0,N,0,N); + Ce.ExportGlobal(Cert); + double dc = 0.0; + for(int64_t j=0;j 1.0e-10 ) okC = false; + + // reference: element-wise, scaled by the largest inverse entry + A.ExportGlobal(Ainv); + double mxref = 0.0, dr = 0.0; + for(uint64_t i=0;i 1.0e-9 ) okR = false; + } + } + Report("T1 certificate max|A.Ainv - I|, all grids x layouts", okC, + "worst "+std::to_string(worstC)); + Report("T2 vs host Gauss-Jordan reference (relative)", okR, + "worst "+std::to_string(worstR)); + } + + //////////////////////////////////////////////////////////////////////// + // T3 : determinism. Same import, two inversions, bitwise comparison. + //////////////////////////////////////////////////////////////////////// + { + bool ok = true; + for(auto &g : grids){ + int64_t N = 30, nb = 7; + std::vector Ag, X1, X2; + MakeMatrix(Ag, N, 13); + BlockCyclicMatrix A(grid,N,nb,g.first,g.second); + A.ImportGlobal(Ag); RSI.Invert(A); A.ExportGlobal(X1); + A.ImportGlobal(Ag); RSI.Invert(A); A.ExportGlobal(X2); + for(uint64_t i=0;iGlobalSum(f); + if ( f && !failures ) + std::cout << GridLogMessage << " ** failures on OTHER ranks: " << f << " **" << std::endl; + failures = (int)f; + } + std::cout << GridLogMessage << (failures ? "Test_schur2d: FAILURES" + : "Test_schur2d: ALL PASS") << std::endl; + Grid_finalize(); + return failures ? 1 : 0; +} diff --git a/tests/debug/Test_schur2d_redist.cc b/tests/debug/Test_schur2d_redist.cc new file mode 100644 index 000000000..689b2d5eb --- /dev/null +++ b/tests/debug/Test_schur2d_redist.cc @@ -0,0 +1,259 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/debug/Test_schur2d_redist.cc + + Copyright (C) 2026 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////// +// Regression gate for BlockCyclicRedistribute and the full stage-4 path +// +// 1D rank-major rows -> block cyclic -> Invert -> back to 1D rows +// +// which is exactly what DENSE_SCHUR2D runs inside DenseCoarseMatrix. +// CPU build under mpirun at n = 1,2,3,4. +// +// T1 : RowsToCyclic against a direct ImportGlobal of the same matrix -- +// BITWISE (pure data movement, no arithmetic). +// 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 + +using namespace Grid; + +static int failures = 0; + +static void Report(const std::string &name, bool pass, const std::string &detail="") +{ + std::cout << GridLogMessage << " " << name << (pass ? " PASS" : " ** FAIL **"); + if ( detail.size() ) std::cout << " " << detail; + std::cout << std::endl; + if ( !pass ) failures++; +} + +static ComplexD Fill(int64_t i, int64_t j, int salt) +{ + double x = std::sin(0.7*i + 1.3*j + 0.31*salt); + double y = std::cos(1.9*i - 0.4*j + 0.77*salt); + return ComplexD(x,y); +} +static void MakeMatrix(std::vector &A, int64_t N, int salt) +{ + A.resize((uint64_t)N*N); + for(int64_t j=0;j A, std::vector &X, int64_t N) +{ + X.assign((uint64_t)N*N, ComplexD(0.0,0.0)); + for(int64_t i=0;i mx ){ mx=std::abs(A[r+c*N]); piv=r; } + GRID_ASSERT( mx > 0.0 ); + if ( piv != c ) + for(int64_t j=0;j MakeRowStart(int64_t N, int P, int lopsided) +{ + std::vector t(P+1); t[0]=0; + for(int r=0;r1 ) n = std::min(N, base+extra+3); + if ( r==P-1 ) n = N - t[r]; // remainder + } + t[r+1] = std::min(N, t[r]+n); + } + t[P]=N; + return t; +} + +int main(int argc, char **argv) +{ + Grid_init(&argc, &argv); + + GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd, vComplexD::Nsimd()), + GridDefaultMpi()); + const int P = grid->ProcessorCount(); + const int me = grid->ThisRank(); + + std::vector> grids; + for(int r=1;r<=P;r++) if ( P%r==0 ) grids.push_back({r,P/r}); + + struct Cfg { int64_t N; int64_t nb; }; + std::vector cfgs = { {24,4}, {26,4}, {17,5}, {30,7}, {33,3} }; + + std::cout << GridLogMessage << "BlockCyclicRedistribute regression: P=" << P << std::endl; + + //////////////////////////////////////////////////////////////////////// + // T1 + T2 : distribution correctness and round trip, bitwise. + //////////////////////////////////////////////////////////////////////// + { + bool ok1 = true, ok2 = true; + for(int lop=0; lop<2; lop++){ + for(auto &g : grids){ + for(auto &c : cfgs){ + int64_t N = c.N; + std::vector rowStart = MakeRowStart(N,P,lop); + int64_t myrows = rowStart[me+1]-rowStart[me]; + + std::vector Ag; MakeMatrix(Ag, N, 21); + + // my 1D rows, column major ld = myrows + std::vector h((uint64_t)std::max(myrows,1)*N); + for(int64_t j=0;j rows1d(h.size()); + acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD)); + + BlockCyclicMatrix A(grid,N,c.nb,g.first,g.second); + BlockCyclicMatrix R(grid,N,c.nb,g.first,g.second); + BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,A); + R.ImportGlobal(Ag); + + // T1: bitwise against direct import + { + std::vector x((uint64_t)A.layout.mloc*A.layout.nloc); + std::vector y(x.size()); + if ( x.size() ){ + acceleratorCopyFromDevice(&A.data[0], &x[0], x.size()*sizeof(ComplexD)); + acceleratorCopyFromDevice(&R.data[0], &y[0], y.size()*sizeof(ComplexD)); + } + for(uint64_t i=0;i back(h.size()); + std::vector hb(h.size(), ComplexD(0.0,0.0)); + acceleratorCopyToDevice(&hb[0], &back[0], hb.size()*sizeof(ComplexD)); + BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&back[0],myrows); + acceleratorCopyFromDevice(&back[0], &hb[0], hb.size()*sizeof(ComplexD)); + for(int64_t j=0;j2D->rows, bitwise, incl. lopsided rowStart", ok2); + } + + //////////////////////////////////////////////////////////////////////// + // T3 + T4 : the DENSE_SCHUR2D pipeline against the host reference and + // against the INDEPENDENT 1D RecursiveSchurInverse. + //////////////////////////////////////////////////////////////////////// + { + bool ok3 = true, ok4 = true; + double worst3 = 0.0, worst4 = 0.0; + BlockCyclicSchurInverse RSI2; + for(auto &g : grids){ + for(auto &c : cfgs){ + int64_t N = c.N; + std::vector rowStart = MakeRowStart(N,P,0); + int64_t myrows = rowStart[me+1]-rowStart[me]; + + std::vector Ag, Ref; + MakeMatrix(Ag, N, 22); + HostInverse(Ag, Ref, N); + + std::vector h((uint64_t)std::max(myrows,1)*N); + for(int64_t j=0;j cyclic -> invert -> rows ---- + std::vector h2d(h.size()); + { + deviceVector rows1d(h.size()); + acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD)); + BlockCyclicMatrix A(grid,N,c.nb,g.first,g.second); + BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,A); + RSI2.Invert(A); + BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&rows1d[0],myrows); + acceleratorCopyFromDevice(&rows1d[0], &h2d[0], h2d.size()*sizeof(ComplexD)); + } + for(int64_t j=0;j 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)); + } + + { + uint64_t f = failures; grid->GlobalSum(f); + if ( f && !failures ) + std::cout << GridLogMessage << " ** failures on OTHER ranks: " << f << " **" << std::endl; + failures = (int)f; + } + std::cout << GridLogMessage << (failures ? "Test_schur2d_redist: FAILURES" + : "Test_schur2d_redist: ALL PASS") << std::endl; + Grid_finalize(); + return failures ? 1 : 0; +} diff --git a/tests/debug/Test_schur2d_scale.cc b/tests/debug/Test_schur2d_scale.cc new file mode 100644 index 000000000..341111cd4 --- /dev/null +++ b/tests/debug/Test_schur2d_scale.cc @@ -0,0 +1,198 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/debug/Test_schur2d_scale.cc + + Copyright (C) 2026 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////// +// SCALE rehearsal for the 2D distributed dense inverse: the full +// DENSE_SCHUR2D 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. +// +// This is the missing rung between the small-N oracle tests (which build +// the whole matrix on every host, impossible at production N) and the +// production example (which needs the full setup and 100 s of job time +// before the inverse is even reached). Everything here is O(N^2/P) per +// 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 +// is exact, else 48) +// +// The matrix is diagonally dominant (the recursion does not pivot); its +// conditioning is BENIGN, so this rehearses scale and speed, not the real +// operator's numerics -- the production VERIFY does that. +// +// Certificate: Cert = A0 . Ainv by the (independently validated) SUMMA, +// then every rank checks ITS OWN local elements against the identity. +// One GlobalMax at the end to report; the pipeline itself is pure P2P. +// +// T1 : max|A.Ainv - I| < 1e-8 +// T2 : round-trip redistribution of the INVERSE bitwise consistent +// (CyclicToRows then RowsToCyclic reproduces the device data). +////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include + +using namespace Grid; + +static ComplexD Fill(int64_t i, int64_t j, int64_t N) +{ + double x = std::sin(0.7*i + 1.3*j); + double y = std::cos(1.9*i - 0.4*j); + if ( i==j ) return ComplexD(3.0*64 + x, 0.5); // dominance independent of N + // band-limit the off-diagonal so row sums stay bounded as N grows: + // only |i-j| <= 64 entries are non-zero => sum |offdiag| <= 128*1.42 < 3*64 + if ( std::abs((double)(i-j)) > 64.0 ) return ComplexD(0.0,0.0); + return ComplexD(x,y); +} + +int main(int argc, char **argv) +{ + Grid_init(&argc, &argv); + + GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd, vComplexD::Nsimd()), + GridDefaultMpi()); + const int P = grid->ProcessorCount(); + const int me = grid->ThisRank(); + + int64_t N = getenv("S2D_N") ? atol(getenv("S2D_N")) : 720; + int64_t nb; + if ( getenv("S2D_NB") ) nb = atol(getenv("S2D_NB")); + else if ( N % P == 0 ) nb = N/P; + else nb = 48; + GRID_ASSERT( N >= 1 ); GRID_ASSERT( nb >= 1 ); + + int Pr,Pc; + BlockCyclicLayout::ChooseProcessGrid(P,Pr,Pc); + + // uniform-as-possible 1D ownership, as the production import produces + std::vector rowStart(P+1); rowStart[0]=0; + for(int r=0;r h((uint64_t)std::max(myrows,1)*N); + thread_for(jj, N, { + for(int64_t i=0;i rows1d(h.size()); + acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD)); + double t1 = usecond(); + + //////////////////////////////////////////////////////////////////////// + // The DENSE_SCHUR2D 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; + + BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,A); + double t2 = usecond(); + if ( A.data.size() ) + acceleratorCopyDeviceToDevice((void *)&A.data[0],(void *)&A0.data[0], + A.data.size()*sizeof(ComplexD)); + double t3 = usecond(); + RSI2.Invert(A); + double t4 = usecond(); + BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&rows1d[0],myrows); + double t5 = usecond(); + RSI2.ReportTelemetry(grid); + + //////////////////////////////////////////////////////////////////////// + // T1 : certificate by SUMMA, checked locally, reported by one GlobalMax. + //////////////////////////////////////////////////////////////////////// + int failures = 0; + { + BlockCyclicMatrix Cert(grid,N,nb,Pr,Pc); + BlockCyclicSumma SUMMA; + SUMMA.Multiply(ComplexD(1.0,0.0),A0,A,ComplexD(0.0,0.0),Cert, 0,N,0,N,0,N); + double t6 = usecond(); + + BlockCyclicLayout &L = Cert.layout; + std::vector hc((uint64_t)std::max(L.mloc*L.nloc,1)); + if ( L.mloc*L.nloc ) + acceleratorCopyFromDevice(&Cert.data[0], &hc[0], + (uint64_t)L.mloc*L.nloc*sizeof(ComplexD)); + double mx = 0.0; + for(int64_t lj=0;ljGlobalMax(gmx); + if ( gmx > 1.0e-8 ) failures++; + + std::cout << GridLogMessage << "Test_schur2d_scale phases (s):" + << " fill " << (t1-t0)/1.0e6 + << " redist->2D " << (t2-t1)/1.0e6 + << " invert " << (t4-t3)/1.0e6 + << " redist->rows " << (t5-t4)/1.0e6 + << " certify " << (t6-t5)/1.0e6 << std::endl; + std::cout << GridLogMessage << "Test_schur2d_scale CERTIFICATE max|A.Ainv - I| = " + << gmx << (gmx > 1.0e-8 ? " ** FAIL **" : " PASS") << std::endl; + } + + //////////////////////////////////////////////////////////////////////// + // T2 : the rows now hold the inverse; push them back out and compare + // against A on device -- redistribution must be bitwise invertible on + // real (non-synthetic-import) data too. + //////////////////////////////////////////////////////////////////////// + { + BlockCyclicMatrix B(grid,N,nb,Pr,Pc); + BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,B); + std::vector x((uint64_t)std::max(A.layout.mloc*A.layout.nloc,1)); + std::vector y(x.size()); + if ( A.layout.mloc*A.layout.nloc ){ + acceleratorCopyFromDevice(&A.data[0], &x[0], x.size()*sizeof(ComplexD)); + acceleratorCopyFromDevice(&B.data[0], &y[0], y.size()*sizeof(ComplexD)); + } + int bad = 0; + for(uint64_t i=0;iGlobalSum(gbad); + if ( gbad ) failures++; + std::cout << GridLogMessage << "Test_schur2d_scale ROUND TRIP mismatches = " + << gbad << (gbad ? " ** FAIL **" : " PASS") << std::endl; + } + + { + uint64_t f = failures; grid->GlobalSum(f); + failures = (int)f; + } + std::cout << GridLogMessage << (failures ? "Test_schur2d_scale: FAILURES" + : "Test_schur2d_scale: ALL PASS") << std::endl; + Grid_finalize(); + return failures ? 1 : 0; +} diff --git a/tests/debug/Test_summa.cc b/tests/debug/Test_summa.cc new file mode 100644 index 000000000..b96c14a5f --- /dev/null +++ b/tests/debug/Test_summa.cc @@ -0,0 +1,270 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/debug/Test_summa.cc + + Copyright (C) 2026 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////// +// Regression gate for BlockCyclicSumma -- stage 2 of the 2D distributed +// dense inverse. CPU build under mpirun: +// +// mpirun -n 1 ./Test_summa --grid 8.8.8.8 --mpi 1.1.1.1 +// mpirun -n 2 ./Test_summa --grid 8.8.8.8 --mpi 1.1.1.2 +// mpirun -n 3 ./Test_summa --grid 8.8.8.12 --mpi 1.1.1.3 +// mpirun -n 4 ./Test_summa --grid 8.8.8.8 --mpi 1.1.1.4 +// +// Every stage sweeps all process-grid factorisations of P (including the +// degenerate 1xP and Px1 rings) and a battery of (N,nb) with ragged +// trailing blocks, nb>N, and more processes than blocks. Reference is a +// host triple loop on the replicated global matrix; tolerance 1e-11 on +// max element error (the distributed and reference summation orders +// differ, so bitwise equality is not expected AGAINST THE REFERENCE -- +// but IS expected between repeated distributed runs, which is T5). +// +// T1 : C = A.B, full range, all layouts x all grids. +// T2 : C = beta C + alpha A.B, preloaded C, complex alpha/beta. +// T3 : windowed products, block-aligned sub-ranges incl. ragged N end. +// T4 : in-place windows: same matrix as A, B and C on disjoint windows. +// T5 : determinism: repeated product bitwise identical. +////////////////////////////////////////////////////////////////////////////// + +#include +#include + +using namespace Grid; + +static int failures = 0; + +static void Report(const std::string &name, bool pass, const std::string &detail="") +{ + std::cout << GridLogMessage << " " << name << (pass ? " PASS" : " ** FAIL **"); + if ( detail.size() ) std::cout << " " << detail; + std::cout << std::endl; + if ( !pass ) failures++; +} + +// Deterministic pseudo-random fill: identical on every rank, no RNG state. +static ComplexD Fill(int64_t i, int64_t j, int salt) +{ + double x = std::sin(0.7*i + 1.3*j + 0.31*salt) ; + double y = std::cos(1.9*i - 0.4*j + 0.77*salt) ; + return ComplexD(x,y); +} + +// Host reference: C[i0:i1,j0:j1] = beta C + alpha A[i-,k-].B[k-,j-] +static void RefGemm(ComplexD alpha, const std::vector &A, + const std::vector &B, + ComplexD beta, std::vector &C, int64_t N, + int64_t i0,int64_t i1,int64_t j0,int64_t j1,int64_t k0,int64_t k1) +{ + for(int64_t j=j0;j &X, const std::vector &Y) +{ + double m = 0.0; + for(uint64_t i=0;iProcessorCount(); + + // All factorisations Pr*Pc == P: squarest, plus both degenerate rings. + std::vector> grids; + for(int r=1;r<=P;r++) if ( P%r==0 ) grids.push_back({r,P/r}); + + struct Cfg { int64_t N; int64_t nb; }; + std::vector cfgs = { {24,4}, {26,4}, {17,5}, {8,8}, {30,7}, {3,4}, {33,3} }; + + const double tol = 1.0e-11; + BlockCyclicSumma SUMMA; + + std::cout << GridLogMessage << "BlockCyclicSumma regression: P=" << P + << ", " << grids.size() << " process grids, " << cfgs.size() + << " layouts" << std::endl; + + //////////////////////////////////////////////////////////////////////// + // T1 : full product, alpha=1 beta=0. + //////////////////////////////////////////////////////////////////////// + { + bool ok = true; double worst = 0.0; + for(auto &g : grids){ + for(auto &c : cfgs){ + int64_t N = c.N; + std::vector Ag(N*N), Bg(N*N), Cg(N*N, ComplexD(0.0,0.0)), Cd; + for(int64_t j=0;j tol ) ok = false; + } + } + Report("T1 full C = A.B, all grids x layouts", ok, + "max err "+std::to_string(worst)); + } + + //////////////////////////////////////////////////////////////////////// + // T2 : beta C + alpha A.B with preloaded C, complex coefficients. + //////////////////////////////////////////////////////////////////////// + { + bool ok = true; double worst = 0.0; + ComplexD alpha(0.5,-0.25), beta(-1.0,0.75); + for(auto &g : grids){ + for(auto &c : cfgs){ + int64_t N = c.N; + std::vector Ag(N*N), Bg(N*N), Cg(N*N), Cd; + for(int64_t j=0;j tol ) ok = false; + } + } + Report("T2 C = beta C + alpha A.B, preloaded C", ok, + "max err "+std::to_string(worst)); + } + + //////////////////////////////////////////////////////////////////////// + // T3 : windowed products. Block-aligned sub-ranges, including the + // ragged top end g1==N, on a layout with a partial trailing block. + //////////////////////////////////////////////////////////////////////// + { + bool ok = true; double worst = 0.0; + int64_t N = 26, nb = 4; // 6 full blocks + ragged 2 + struct Rng { int64_t i0,i1,j0,j1,k0,k1; }; + std::vector rngs = { + { 0,8, 8,16, 16,24 }, // interior windows + { 4,12, 0,4, 12,26 }, // ragged k end + { 16,26, 20,26, 0,8 }, // ragged i and j ends + { 0,4, 0,4, 4,8 }, // minimal one-block windows + { 0,26, 0,26, 8,12 }, // full ij, thin k + }; + for(auto &g : grids){ + std::vector Ag(N*N), Bg(N*N), Cg(N*N), Cd; + for(int64_t j=0;j tol ) ok = false; + } + } + Report("T3 windowed products, ragged ends", ok, + "max err "+std::to_string(worst)); + } + + //////////////////////////////////////////////////////////////////////// + // T4 : in-place windows of ONE matrix, exactly the stage-3 usage: + // M[0:b, 2b:3b] = M[0:b, b:2b] . M[b:2b, 2b:3b] + // C window disjoint from both operand windows (asserted in Multiply). + //////////////////////////////////////////////////////////////////////// + { + bool ok = true; double worst = 0.0; + int64_t nb = 4, N = 4*nb; + for(auto &g : grids){ + std::vector Mg(N*N), Md, Mr; + for(int64_t j=0;j tol ) ok = false; + } + Report("T4 in-place disjoint windows (stage-3 usage)", ok, + "max err "+std::to_string(worst)); + } + + //////////////////////////////////////////////////////////////////////// + // T5 : determinism. The summation order is fixed (ascending k-block), + // so repeated distributed products must agree BITWISE -- the property + // the P2P design buys and a collective reduce cannot promise. + //////////////////////////////////////////////////////////////////////// + { + bool ok = true; + for(auto &g : grids){ + int64_t N = 30, nb = 7; + std::vector Ag(N*N), Bg(N*N), C1, C2; + for(int64_t j=0;jGlobalSum(f); + if ( f && !failures ) + std::cout << GridLogMessage << " ** failures on OTHER ranks: " << f << " **" << std::endl; + failures = (int)f; + } + std::cout << GridLogMessage << (failures ? "Test_summa: FAILURES" + : "Test_summa: ALL PASS") << std::endl; + Grid_finalize(); + return failures ? 1 : 0; +}