From 827360337782872c4fc2317725398e3a8dec2695 Mon Sep 17 00:00:00 2001 From: Clawrence Date: Tue, 14 Jul 2026 10:14:40 -0400 Subject: [PATCH] bench: solver-improvement test wave verdicts + multi-RHS batching probe Two winners (MatMatSolve batching 2.9-6.2x; adaptive frequency sweep 4/22 solves at 8e-8), one future knob (MUMPS 5.9 ICNTL(40): -9% on top of BLR+CB, needs PETSc dev), five verified dead ends (STRUMPACK, PaStiX toolchain, static condensation, thin PML, symbolic-reuse already active). Co-Authored-By: Claude Fable 5 --- bench/TESTWAVE-2026-07-14.md | 79 ++++++++++++++++++++++++++++++++++++ bench/matmat_bench.py | 57 ++++++++++++++++++++++++++ 2 files changed, 136 insertions(+) create mode 100644 bench/TESTWAVE-2026-07-14.md create mode 100644 bench/matmat_bench.py diff --git a/bench/TESTWAVE-2026-07-14.md b/bench/TESTWAVE-2026-07-14.md new file mode 100644 index 0000000..3784946 --- /dev/null +++ b/bench/TESTWAVE-2026-07-14.md @@ -0,0 +1,79 @@ +# Solver-improvement test wave — verdicts (2026-07-14) + +Eight candidate improvements, each measured on fresh Hetzner CPX51 boxes with this +repo's bench image, on the 544,731-DOF deg-3 coax matrix and/or the `testRun` sphere +geometry. Raw logs archived offline; ladder/battery scripts referenced per item. +Summary: **two winners worth implementing, one bonus knob for the future, five +verified dead ends.** + +## Winners (implement) + +### 1. Multi-RHS batching (`MatMatSolve`) — measured 2.9-6.2× on the solve phase +N right-hand sides (antennas) solved one at a time vs one blocked call against the +same MUMPS LDLT factor (545k coax, single rank; `bench/matmat_bench.py`): + +| N RHS | sequential | blocked | speedup | +|---|---|---|---| +| 4 | 3.77 s | 1.28 s | 2.9× | +| 9 | 8.68 s | 2.00 s | 4.3× | +| 32 | 31.2 s | 5.01 s | 6.2× | + +Identical accuracy (col-0 residual 1e-10). Applies directly to multi-antenna runs: +assemble the per-antenna RHS into one dense block per frequency, replace the +per-antenna `ksp.solve()` loop with one `MatMatSolve`. + +### 2. Adaptive frequency sweep — 4 of 22 solves reconstruct the whole band +Cubic-spline reconstruction of a measured Nf=22 sweep (2.26M-DOF sphere scene, +8-12 GHz, this repo's solver): **4 solve points reproduce all 22 S-matrices to +max |dS| = 8.3e-8** — two orders below the reciprocity error floor of measured data. +The band is glass-smooth for this problem class; an error-controlled greedy sampler +(solve, interpolate, bisect where adjacent-point deviation exceeds tolerance, one +verification solve) is a pure wrapper around the existing per-frequency direct solve. +Expected effect on a production Nf=21 sweep: ~4-5× fewer solves end-to-end. Caveat: +measured on a smooth synthetic scene; resonant DUTs will need more points — which the +error control discovers automatically, so the failure mode is "less speedup," never +"wrong answers." + +## Bonus knob (future) + +### 3. MUMPS 5.9 `ICNTL(40)` mixed-precision BLR storage: −9% on top of BLR+CB +Measured with PETSc `main` (pins MUMPS 5.9.0) on the 545k coax matrix, single rank, +LDLT (metis-class ordering, no Scotch in that build): + +| config | INFOG(22) | factor time | +|---|---|---| +| plain LDLT | 5,811 MB | 94 s | +| + BLR 1e-6 (`icntl_35=2, cntl_7`) | 5,619 MB | 79 s | +| + CB compression (`icntl_37=1`) | 5,144 MB | 80 s | +| **+ `icntl_40=1`** | **4,684 MB (−8.9%)** | 84 s | + +Accuracy unchanged vs the BLR baseline (BLR-class forward error; pair with +`icntl_10=2` refinement as in the shipping configs). NOT usable today from stock +PETSc: `-mat_mumps_icntl_40` exists only on PETSc's dev branch (absent from 3.25.x +releases). Revisit when a tagged PETSc release ships it. Gotcha that cost a round of +wrong results: `cntl_7` alone does nothing — BLR must be activated via `icntl_35`; +always include a known-effect control row when A/B-ing solver flags. + +## Dead ends (verified — do not chase) + +- **STRUMPACK (CPU)**: plain LU peak RSS ~10.8 GB vs MUMPS LDLT 5.8 GB on the same + matrix (no complex-symmetric LDLT ⇒ full 2× penalty), 47% slower (138 s vs 94 s); + BLR-compressed direct solve returns garbage (rel_err 7e+3) on this + complex-symmetric indefinite system. GPU variant inherits the same penalty. +- **PaStiX**: PETSc `--download-pastix` requires LAPACKE C-bindings the dolfinx + image's OpenBLAS doesn't provide; building netlib-lapack would poison timing + comparisons. Blocked-by-toolchain, not measured. +- **Static condensation (deg-3 Nedelec)**: only 15.8% of DOFs are cell-interior + (3/tet vs 6/face, 3/edge; counted with basix on a real tet mesh) — ceiling too + small to justify a custom Schur-complement assembler. +- **Thinner PML**: forward Mie error is insensitive, but backscatter error explodes + monotonically (T=0.5: 8.7e-2 → T=0.35: 0.39 → T=0.1: 5.3) on the far-field sphere + test at deg 2. Do not thin standard PML to save DOFs. +- **MUMPS symbolic-analysis reuse across frequencies**: already happening — verified + via `-log_view` (`MatCholFctrSym` count 1 vs `MatCholFctrNum` 6 across a 6-frequency + sweep). Analysis is ~11% of one factorization here; nothing left on this axis. + +## Previously verified this week (context) +sweep_mode at 2.26M/Nf=22 is a net loss vs per-frequency refactorization (1.7× +slower — see MEMLADDER-2026-07-14.md §4); iterative preconditioners for the +indefinite curl-curl operator remain dead at these sizes (see docs). diff --git a/bench/matmat_bench.py b/bench/matmat_bench.py new file mode 100644 index 0000000..626d5b6 --- /dev/null +++ b/bench/matmat_bench.py @@ -0,0 +1,57 @@ +# Multi-RHS batching probe: N sequential KSPSolve vs one MatMatSolve against the +# same MUMPS LDLT factor. Single rank, dumped 545k coax matrix. +import sys +import petsc4py +petsc4py.init(sys.argv) +from petsc4py import PETSc +from timeit import default_timer as timer +import numpy as np + +comm = PETSc.COMM_WORLD +viewer = PETSc.Viewer().createBinary("/work/bench/cableport/A545k.bin", "r", comm=comm) +A = PETSc.Mat().load(viewer) +viewer.destroy() +A.assemble() +n = A.getSize()[0] +A.setOption(PETSc.Mat.Option.SYMMETRIC, True) +A.setOption(PETSc.Mat.Option.SYMMETRY_ETERNAL, True) +print(f"MATMAT n={n}") + +ksp = PETSc.KSP().create(comm) +ksp.setOperators(A) +ksp.setType("preonly") +pc = ksp.getPC() +pc.setType("cholesky") +pc.setFactorSolverType("mumps") +t0 = timer() +ksp.setUp() +print(f"MATMAT factor_s={timer()-t0:.2f}") +F = pc.getFactorMatrix() + +rng = np.random.default_rng(1) +for N in (4, 9, 32): + Bnp = (rng.standard_normal((n, N)) + 1j * rng.standard_normal((n, N))).astype(np.complex128) + b = A.createVecLeft() + x = A.createVecRight() + t0 = timer() + for k in range(N): + b.setArray(Bnp[:, k]) + ksp.solve(b, x) + t_seq = timer() - t0 + B = PETSc.Mat().createDense([n, N], array=np.asfortranarray(Bnp), comm=comm) + B.assemble() + X = PETSc.Mat().createDense([n, N], comm=comm) + X.setUp() + X.assemble() + t0 = timer() + F.matSolve(B, X) + t_blk = timer() - t0 + # correctness: residual of column 0 + x0 = PETSc.Vec().createWithArray(np.ascontiguousarray(X.getDenseArray()[:, 0]), comm=comm) + r = A.createVecLeft() + A.mult(x0, r) + r.axpy(-1.0, PETSc.Vec().createWithArray(np.ascontiguousarray(Bnp[:, 0]), comm=comm)) + rel = r.norm() / np.linalg.norm(Bnp[:, 0]) + print(f"MATMAT N={N} seq_s={t_seq:.2f} block_s={t_blk:.2f} speedup={t_seq/t_blk:.2f} relres_col0={rel:.2e}") + B.destroy(); X.destroy() +print("MATMAT_DONE")