bench: solver-improvement test wave verdicts

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).

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")