diff --git a/tests/debug/Test_schur_dense_coarse.cc b/tests/debug/Test_schur_dense_coarse.cc new file mode 100644 index 000000000..a2168cab6 --- /dev/null +++ b/tests/debug/Test_schur_dense_coarse.cc @@ -0,0 +1,289 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: Test_schur_dense_coarse.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 */ + +// +// 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. +// +// 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. +// +// 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): +// mpirun -n 1 ./Test_schur_dense_coarse --grid 12.12.12.12 --mpi 1.1.1.1 +// mpirun -n 2 ./Test_schur_dense_coarse --grid 12.12.12.24 --mpi 1.1.1.2 +// mpirun -n 3 ./Test_schur_dense_coarse --grid 12.12.12.36 --mpi 1.1.1.3 +// mpirun -n 4 ./Test_schur_dense_coarse --grid 12.12.12.48 --mpi 1.1.1.4 +// +#include +#include +#include +#include + +using namespace std; +using namespace Grid; + +/////////////////////////////////////////////////////////////////////// +// MdagM + shift: Galerkin projection of a PD operator plus sigma I is +// safely invertible whatever the (random) subspace quality. +/////////////////////////////////////////////////////////////////////// +template +class ShiftedHermOpAdaptor : public LinearOperatorBase +{ + LinearOperatorBase &wrapped; + RealD shift; +public: + ShiftedHermOpAdaptor(LinearOperatorBase &wrapme, RealD s) + : wrapped(wrapme), shift(s) {}; + + void Op(const Field &in, Field &out) + { + wrapped.HermOp(in, out); + out = out + shift*in; + } + void AdjOp(const Field &in, Field &out) + { + Op(in, out); + } + void HermOp(const Field &in, Field &out) + { + Op(in, out); + } + void OpDiag(const Field &in, Field &out) { GRID_ASSERT(0); } + void OpDir (const Field &in, Field &out, int dir, int disp) { GRID_ASSERT(0); } + void OpDirAll(const Field &in, std::vector &out) { GRID_ASSERT(0); } + void HermOpAndNorm(const Field &in, Field &out, RealD &n1, RealD &n2) { GRID_ASSERT(0); } +}; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls = 4; + const int nbasis = 4; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + // Per-dimension blocking {4,4,3,3}: fine 12.12.12.12 -> coarse + // 3.3.4.4. Two constraints meet here (both MEASURED today): + // - coarse dims of 2 hit the probing pathology (health probe below) + // - GEN-simd lanes {1,1,2,2} must land on even coarse dims, so the + // odd production-like 3s go on the lane-free x,y axes (exactly the + // production [3,6,8,8] trick). + Coordinate blocks({4,4,3,3}); + Coordinate clatt = GridDefaultLatt(); + for(int d=0; d seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + std::vector cseeds({9,10,11,12}); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG CRNG(Coarse5d); CRNG.SeedFixedIntegers(cseeds); + + LatticeGaugeField Umu(UGrid); + SU::HotConfiguration(RNG4,Umu); + + RealD mass = 0.1; + RealD M5 = 1.8; + DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + + MdagMLinearOperator HermDefOp(Ddwf); + ShiftedHermOpAdaptor HOA(HermDefOp, 0.5); + + std::cout << GridLogMessage << "Building random aggregation space, nbasis " << nbasis << std::endl; + typedef Aggregation Subspace; + Subspace Aggregates(Coarse5d,FGrid,0); + Aggregates.CreateSubspaceRandom(RNG5); + + std::cout << GridLogMessage << "Coarsening shifted MdagM" << std::endl; + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); + LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse5d); + LittleDiracOp.CoarsenOperator(HOA,Aggregates); + + /////////////////////////////////////////////////////////////////////// + // Operator health probes (independent of DenseCoarseMatrix). + // + // MEASURED PATHOLOGY, banked 2026-08-14: on coarse dims of 2 (fine + // 8.8.8.8, block 4 -> coarse 2.2.2.2) the coarsened operator is + // rank 16/128 with 112 zero ROWS (output support = 2 of 16 sites) + // and Hermiticity violation 0.17 -- the probing construction breaks + // on the size-2 torus. The import certificate cannot see this + // (dense and M share _A). Out of scope here; coarse dims >= 3. + // + // Cheap any-size probes: output support + Hermiticity via inner + // products on random vectors. + /////////////////////////////////////////////////////////////////////// + { + CoarseVector px(Coarse5d); + CoarseVector py(Coarse5d); + CoarseVector Mx(Coarse5d); + CoarseVector My(Coarse5d); + random(CRNG, px); + random(CRNG, py); + LittleDiracOp.M(px, Mx); + LittleDiracOp.M(py, My); + ComplexD ip1 = innerProduct(py, Mx); // + ComplexD ip3 = innerProduct(px, My); // + RealD hermdev = abs(ip1 - conj(ip3)) / std::sqrt(norm2(Mx)*norm2(py)); + RealD support = norm2(Mx) / norm2(px); + std::cout << GridLogMessage << "Operator health: ||Mx||^2/||x||^2 = " << support + << " herm-dev " << hermdev << std::endl; + // Hermitian fine op => exactly Hermitian Galerkin coarse op. + // (A measured herm-dev of 1.3e-4 here was the CPU SIMT-lane + // CoarsenOperator bug -- fixed 2026-08-14, now 4e-15. A loud + // failure here means _A population is broken again.) + GRID_ASSERT( support > 1.0e-3 ); + GRID_ASSERT( hermdev < 1.0e-10 ); + } + + /////////////////////////////////////////////////////////////////////// + // Full-matrix conditioning probe at small N: dense columns by + // applying M to unit vectors, fp64 Eigen SVD. + /////////////////////////////////////////////////////////////////////// + { + int64_t Nprobe = Coarse5d->gSites() * nbasis; + if ( Nprobe <= 700 ) + { + Eigen::MatrixXcd eA(Nprobe, Nprobe); + CoarseVector e(Coarse5d); + CoarseVector Me(Coarse5d); + for(int64_t j=0; j_ndimension); + Lexicographic::CoorFromIndex(gcoor, gsite, Coarse5d->GlobalDimensions()); + typedef typename CoarseVector::vector_object::scalar_object csobj; + csobj s; + s = Zero(); + ((ComplexD *)&s)[b] = ComplexD(1.0,0.0); + pokeSite(s, e, gcoor); + LittleDiracOp.M(e, Me); + for(int64_t i=0; i_ndimension); + Lexicographic::CoorFromIndex(gci, gsi, Coarse5d->GlobalDimensions()); + csobj si; + peekSite(si, Me, gci); + eA(i,j) = ((ComplexD *)&si)[bi]; + } + } + Eigen::JacobiSVD svd(eA); + double smax = svd.singularValues()(0); + double smin = svd.singularValues()(Nprobe-1); + int64_t rank = 0; + for(int64_t i=0; i 1.0e-10*smax ) rank++; + } + double herm = (eA - eA.adjoint()).cwiseAbs().maxCoeff(); + int64_t zrows = 0; + int64_t zcols = 0; + for(int64_t i=0; i DenseCC; + DenseCC dcm(LittleDiracOp, Coarse5d); + + 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 + /////////////////////////////////////////////////////////////////////// + CoarseVector x(Coarse5d); + CoarseVector y(Coarse5d); + CoarseVector z(Coarse5d); + random(CRNG, x); + dcm(x, y); + LittleDiracOp.M(y, z); + z = z - x; + RealD rel = std::sqrt(norm2(z)/norm2(x)); + std::cout << GridLogMessage << "T6 round trip ||A Ainv x - x||/||x|| (random x) = " + << rel << std::endl; + GRID_ASSERT( rel < 1.0e-2 ); + + 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 new file mode 100644 index 000000000..3dd986c88 --- /dev/null +++ b/tests/debug/Test_schur_inverse.cc @@ -0,0 +1,763 @@ +/************************************************************************************* + + 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 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(); +}