ForceFreeStates.Galerkin - FEATURE! - Solve the Galerkin system by banded Cholesky, with LU fallback

[StuartBenjamin]

Release note

• Audience: users
• Numerical impact: the Galerkin Δ′ solve now uses Cholesky wherever the matrix is positive definite. On gal_resistive_diiid the PEST3 Δ diagonal moves by 0.40 % and the Δ′ matrix norm by 2e-5 (relative) (harness @ cd290c6)
• Migration: to reproduce the previous results, set gal_solver = "LU".

The Galerkin Δ′ solve now factors the matrix by banded Cholesky and falls back to LU when it is not positive definite (gal_solver = "auto", the new default, as @matt-pharr suggested). Before, gal_solver = "cholesky" crashed (UndefVarError: pbtrf!).

Regression report

regress --cases gal_resistive_diiid,gal_resistive_pe --refs 6fc89d021,cd290c659 (develop merge-base vs this head; julia 1.11.7, pinned manifest, 2 threads). gal_resistive_pe is N/A on both refs, as on develop.

Case: gal_resistive_diiid — DIII-D-like, n=1, RDCON outer-region Galerkin Δ′ with rpec coil columns (delta_coil)
================================================================

Regression Report: gal_resistive_diiid
=====================================================================================
Ref 1: 6fc89d021  @ 6fc89d021 (2026-10-09)
       env: julia 1.11.7, x86_64-linux-gnu, manifest fdfdedfa (pinned), 2 threads/2 BLAS, jit build
Ref 2: cd290c659  @ cd290c659 (2026-10-13)
       env: julia 1.11.7, x86_64-linux-gnu, manifest fdfdedfa (pinned), 2 threads/2 BLAS, jit build
-------------------------------------------------------------------------------------
Quantity                 6fc89d021     cd290c659     Diff               Status       
-------------------------------------------------------------------------------------
gal # singular surfaces  4             4             0.0e+00            OK           
gal singular q values    [4 elem]      [4 elem]      0.0e+00            OK           
gal PEST3 Δ diagonal     [4 elem]      [4 elem]      6.140e+01 (0.40%)  ** CHANGED **
||gal Δ′ matrix||        1.243825e+04  1.243849e+04  2.367e-01 (0.00%)  ** CHANGED **
||gal Δ_coil block||     2.498141e+02  2.498168e+02  2.686e-03 (0.00%)  ** CHANGED **
gal D_I per surface      [4 elem]      [4 elem]      0.0e+00            OK           
gal α per surface        [4 elem]      [4 elem]      0.0e+00            OK           
||gal inner-layer Δ||    1.474630e+06  1.474630e+06  0.0e+00            OK           
||gal match cout||       3.061595e-02  3.061584e-02  1.2e-07            OK           
gal match residual       1.901643e-16  2.612202e-16  7.1e-17            OK           
Runtime (s)              226.0s        185.4s                           --           
=====================================================================================
Summary: 3 changed, 7 unchanged

Notes for reviewers

• gal_solver: "auto" tries Cholesky, then LU; "LU" skips Cholesky; "cholesky" errors instead of falling back. The key stays because the TPS wrapper passes it.
• The matrix is always assembled with the full band, so the LU fallback sees the same matrix as before; the Cholesky-only assembly branches are gone. The Dirichlet rows (axis, fixed edge, coil edge) are identity rows with their columns kept; the Cholesky path moves those columns to the right-hand side, which leaves a Hermitian system.
• Cholesky uses only the lower band and the real part of the diagonal. The resonant diagonal entries carry a small imaginary error (up to 5.6e-5 relative on gal_resistive_diiid, 1.6e-7 with ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces #517), which LU keeps and Cholesky drops: that is the 0.40 % change above.
• zpbtrf's info is checked. Fortran RDCON's gal.f ignores it and solves with a partial factor (fixed in RDCON - BUGFIX - Galerkin factorization failures: Cholesky falls back to LU, a failed LU stops PrincetonUniversity/GPEC#301).
• Unit test runtests_galerkin_cholesky.jl: a banded system with Dirichlet rows and a nonzero prescribed value is solved by Cholesky and matches a dense solve; made indefinite, it falls back to LU with the same answer, and "cholesky" errors.

---

🤖 Generated with Claude Code

[matt-pharr]

Couple cleanup items and one main question, being whether we should do away with the flag altogether and just make every galerkin solve try to use cholesky first and then use LU as a fallback. See the relevant comment.

[d-burg]

#517 was merged into this branch, so I'm reviewing its exponent change here, where it will reach develop.

The exponent change is correct, and it is the paper's definition rather than a departure from it, so that's good. Glasser (2016) defines the powers by det[M₀ − αI] = 0 (Eq. 46). Of α = ±√(−D_I) (Eqs. 48–49) the paper says it could not be proved analytically, only that the roots converge to it as Fourier components are added, and it proposes that agreement as a Fourier-convergence check. So sqrt(−det) was Eq. 46 under an unstated tr = 0 assumption.

From claude: regarding @matt-pharr's question (more right, or only more self-consistent): as far as I can derive it, the resonant block has the form [−κ+σ/2 φ; γ conj(κ)−σ/2] with φ, γ real, so tr = −2i·Im κ. That is a genuine term of the truncated operator, not an assembly slip. Matching to a basis with any other powers makes Δ′ depend on the matching distance, so this is more right for the truncated system. It does hide one thing: the old imaginary Δ′ was the only visible symptom of a non-converged resonant block, so I'd like the trace itself reported. Log Im tr / 2 and α per surface (ideally also next to alpha in the HDF5 output), and warn when |tr| is not small against α. The debug print in Asymptotics.jl (lines 125–126, outside this diff) still prints di next to an α that is no longer √(−di); something like:

        block_det = m0mat[1, 1]*m0mat[2, 2] - m0mat[2, 1]*m0mat[1, 2]
        block_tr = m0mat[1, 1] + m0mat[2, 2]
        msg *= @sprintf("  det= %+.12e, tr= %+.12e %+.12ei\n", real(block_det), real(block_tr), imag(block_tr))
        msg *= @sprintf("  alpha= %+.12e %+.12ei\n", real(alpha[1]), imag(alpha[1]))

Requests before this goes to develop:

1. Release note. This PR now carries #517's results-moving change and the removal of extended_precision_bvp, but the release note covers only gal_solver. Suggested Migration: "To reproduce the previous Galerkin results, set gal_solver = "LU". Remove extended_precision_bvp from [ForceFreeStates] (the TOML key is ignored with a warning; passing it as a keyword to solve is now an error); the Δ′ BVP always solves in Float64. Anything pinned to imaginary Δ′ or to Galerkin alpha must be re-baselined." The Numerical-impact line needs the exponent change too.
2. CI is red at test/runtests_eulerlagrange.jl:706: the opposed-powers testset still checks bitwise a branch that was removed (suggestion inline).
3. Harness at the current head. The stamp here (cd290c659) predates the #517 merge, and #517's own report was against a merge-base older than the axis-start change (#514). The near-axis surface is the one most exposed to that.
4. Wording: "Hermitian with a real diagonal" rather than "Δ′ must be real", as Matt suggested on #517.
5. Solovev. Solovev is up-down symmetric, so its blocks should be traceless and the exponent change should not move it, yet #517 reported Im Δ′ going from −1.8e9 to −27 there. What is Im tr at its q = 2 surface? If it is round-off level, that change needs another explanation; if not, the matrices have an imaginary part on a symmetric equilibrium.
6. Near-axis surface at n = 4. A residual Im Δ′ of −25 (about 1 % of the diagonal) means the series and the integrated solution disagree by about 1 % there; I'd treat that as the error bar. A cheap diagnostic: run with JULIA_DEBUG=GeneralizedPerturbedEquilibrium, compare Im tr, α and dpsi^(−2α) at q = 1.25 with the other surfaces, then a three-point singfac_min scan (3e-4, 1e-4, 3e-5) to see whether Δ′₁₁ plateaus.
7. Small: using DoubleFloats in src/ForceFreeStates/ForceFreeStates.jl (line 22) is now dead; GGJ has its own using.

[d-burg]

The opposed-powers branch this testset pinned is gone, so its bitwise check now fails CI (line 706). The second half of the testset already checks the general path exactly, so I'd drop the first half rather than relax it to ≈.

Suggested change

	@testset "column scaling stays /dpsi^α when the powers are opposed" begin
	α = ComplexF64(1.5)
	dpsi = 1.3e-4
	big = ones(ComplexF64, 2, 2)
	small = ones(ComplexF64, 2, 2)
	FFS._apply_column_powers!(big, small, -α, α, dpsi)
	pfac = dpsi^α
	big_ref = ones(ComplexF64, 2, 2)
	small_ref = ones(ComplexF64, 2, 2)
	big_ref ./= pfac
	small_ref .*= pfac
	@test big == big_ref
	@test small == small_ref
	@testset "column scaling multiplies each column by dpsi^power" begin
	dpsi = 1.3e-4

[d-burg]

With the traceless branch removed this is the only path, so a real block with complex roots throws DomainError and one with real roots returns Float64. Converting keeps it total and type-stable without changing any complex result.

Suggested change

	α = sqrt(tr * tr / 4 - det)
	α = sqrt(ComplexF64(tr * tr / 4 - det))

[d-burg]

When tr²/4 − det is negative (a Mercier-unstable surface), it is real by construction, so the sign of the root comes from the sign of a zero or round-off imaginary part. Left and right can then label opposite solutions "big"; a standalone copy of this helper does that for a real block and its mirror. The old shared α made the pairing consistent by construction, and it is now gone entirely. Could α be canonicalised the same way on both sides (Re α > 0, or Im α > 0 when Re α is zero to round-off), or the pair computed once per surface and passed to both sides?

[d-burg]

These tests pin the 2×2 algebra but not the claim. The assembled block always has a purely imaginary trace, which none of the inputs here has, and the unit-test fixtures are up-down symmetric, so CI never exercises the new path. I'd add a case with the real block structure (before the column-scaling testset):

    @testset "Hermitian-structured block: shared imaginary shift, real eigenvectors, same powers on both sides" begin
        # Shape the assembly produces: [−κ+σ/2 φ; γ conj(κ)−σ/2] with φ, γ real and κ → −κ on the left side.
        κ, φ, γ = 0.3 + 0.02im, 0.7, 2.5
        m_right = ComplexF64[-κ+0.5 φ; γ conj(κ)-0.5]
        m_left = ComplexF64[κ-0.5 φ; γ -conj(κ)+0.5]
        α_right, big_right, small_right = FFS.resonant_block_exponents(m_right, 1.0)
        α_left, big_left, small_left = FFS.resonant_block_exponents(m_left, -1.0)
        @test α_right ≈ sqrt((real(κ) - 0.5)^2 + φ * γ)
        @test α_left ≈ α_right
        @test imag(big_right) ≈ -imag(κ)
        @test imag(small_right) ≈ -imag(κ)
        @test big_left ≈ big_right
        @test small_left ≈ small_right
        # ξ′ entries of both leading eigenvectors are real, so the big/small cross-Wronskian 2α/φ is real.
        @test abs(imag(-(m_right[1, 1] - big_right) / φ)) < 1e-14
        @test abs(imag(-(m_right[1, 1] - small_right) / φ)) < 1e-14
        @test -(m_right[1, 1] - small_right) / φ + (m_right[1, 1] - big_right) / φ ≈ 2 * α_right / φ
    end

Ideally also one integration-level assertion on an asymmetric fixture: left and right powers agree, and the Riccati Δ′ diagonal is real to a stated tolerance.

[d-burg]

I don't think the evidence supports deleting extended_precision_bvp before this reaches develop, only changing its default. It covers only singfac_min = 1e-4; #517's own 2 % miss at singfac_min = 1e-6 (Solovev q₀ = 1.85) is a case where Float64 is measurably wrong, and the code can no longer detect or repair it. A singfac_min scan is also exactly how one would verify the new exponents. The exponent fix and the precision of this LU solve are separate error sources. Could we keep the option as opt-in (default false, which keeps the whole speedup), or replace it with one step of iterative refinement (Float64 LU, residual in Double64) plus a warning when the correction to the read-off entries is large? What do you think?