[StuartBenjamin]

STRIDE solves the Δ′ boundary value problem as a dense matrix of size nMat = (2 + 4 msing) mpert. It factors it with ZGETRF and keeps two dense copies (M and MInv). Time grows as nMat³ and memory as nMat², so at high n with many rational surfaces this one step dominates STRIDE. On the DIII-D resistive example at n = 8 it takes 632 s of the 750 s run and most of the 5.1 GiB peak. The matrix is block-banded, though: each surface's block couples only to its neighbours.

This PR solves it as a banded matrix. Each singular surface's rows now end with that surface's two big-solution rows, instead of all big-solution rows sitting at the end of the matrix. That makes the matrix banded with kl = ku = 3 mpert − 1. It is factored and solved with ZGBTRF / ZGBCON / ZGBTRS, and the residual check uses the sparse entries, so no dense matrix is formed. The change is in stride/ode.F only. The printed condition number now uses the norm of the matrix itself, as ZGBCON expects (develop took the norm after factorization).

Δ′ is identical to the bit with develop in every run below (reference LAPACK). GPEC examples, same machine, 2 threads, develop → this PR:

example	n	nMat	Δ′ step [s]	STRIDE wall [s]	peak memory [GiB]
solovev_resistive	4	720	0.18 → 0.03	6.3 → 6.9	0.10 → 0.09
DIIID_resistive	2	1 520	1.7 → 0.09	7.5 → 6.7	0.18 → 0.12
DIIID_resistive	4	4 056	25 → 0.45	46 → 21	0.77 → 0.31
DIIID_resistive	6	7 296	150 → 1.4	203 → 54	2.15 → 0.62
DIIID_resistive	8	11 704	632 → 4.0	750 → 123	5.1 → 1.1

We also ran 28 reconstructed DIII-D equilibria at n = 1–4 (112 runs) on our GPEC branch with and without this change. Δ′ was bit-identical in all 112. At n = 4 the Δ′ step went from 215 s median (947 s max) to 2.4 s (5.8 s), and the largest STRIDE peak memory went from 7.3 to 1.8 GiB. The attached plot shows both sets against nMat. The zip holds the script that runs the examples (run_demo.sh DEVELOP_DIR PR_DIR OUT), the plotting script and the two tables.

stride_banded_dprime_demo.zip

Review check, on this PR merged onto current develop (with #325's axis-interval split, #328 and the EQUIL fixes), against current develop, DIII-D resistive example: delta_prime.out is byte-identical at n = 1, 2 and 4 on both Linux (glibc 2.39, reference LAPACK, Ubuntu 24.04 container) and macOS (Apple Accelerate).

platform	n	peak memory	STRIDE wall [s]
Linux	2	254 → 197 MB	2.7 → 1.7
Linux	4	940 → 467 MB	23.2 → 5.1
macOS	2	189 → 154 MB (footprint)	1.8 → 1.5
macOS	4	1050 → 358 MB (footprint)	9.5 → 4.5

On the branch's original base, which predates #325, n = 4 was not bit-identical (median 7e-5 relative on the Δ′ diagonal, up to 5% on the innermost surface), because there all 38 driven solves failed the big-solution residual check and the result was dominated by roundoff. With #325 only 2 fail, at about 1e-6, and the banded and dense solves agree to the bit.

🤖 Generated with Claude Code

[matt-pharr]

Ran this myself and can confirm the positive claims of this PR. We mostly have been moving away from banded matrices, but that only applies to the regions where banded matrices are an approximation; they are exact here and can save quite a bit of runtime for higher $n$ problems. Approve.