From b5b430f570e1755e63d5e87aa60b4d8dd63921d6 Mon Sep 17 00:00:00 2001
From: StuartBenjamin <stubenj9@mit.edu>
Date: Wed, 7 Oct 2026 09:35:19 -0700
Subject: [PATCH 1/6] ForceFreeStates.Galerkin - FEATURE - Add the banded
 Cholesky solve for gal_solver = "cholesky"

gal_solver = "cholesky" called LinearAlgebra.LAPACK.pbtrf!/pbtrs!, which do not exist, so the
option died with an UndefVarError. Call LAPACK zpbtrf('L')/zpbtrs('L') directly through
libblastrampoline, unscaled, on the lower band the assembly already builds for Cholesky. Unlike
the reference implementation, a non-positive-definite matrix (zpbtrf info > 0) raises an error
naming gal_solver = "LU" instead of continuing with a partial factor.

Add examples/LAR_gal_cholesky_example (LAR analytic equilibrium, positive-definite Galerkin
matrix), regression case lar_gal_cholesky, test/runtests_galerkin_cholesky.jl, and a solver note
in docs/src/galerkin.md.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 docs/src/galerkin.md                          |  6 ++
 examples/LAR_gal_cholesky_example/gpec.toml   | 64 +++++++++++++++
 .../cases/lar_gal_cholesky.toml               | 78 +++++++++++++++++++
 src/ForceFreeStates/Galerkin/GalerkinSolve.jl | 28 ++++++-
 .../Galerkin/GalerkinStructs.jl               |  2 +-
 test/runtests.jl                              |  1 +
 test/runtests_galerkin_cholesky.jl            | 42 ++++++++++
 7 files changed, 218 insertions(+), 3 deletions(-)
 create mode 100644 examples/LAR_gal_cholesky_example/gpec.toml
 create mode 100644 regression-harness/cases/lar_gal_cholesky.toml
 create mode 100644 test/runtests_galerkin_cholesky.jl

diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md
index 04c4aba22..806fd9ef3 100644
--- a/docs/src/galerkin.md
+++ b/docs/src/galerkin.md
@@ -20,6 +20,12 @@ control surface (when `vac_flag`) and produces no free-boundary energies or ODE
 Setting `gal_match_flag` additionally matches the inner layer, giving a driven ξ solution that
 `PerturbedEquilibrium` consumes in place of a forward solution.
 
+`gal_solver` picks the banded solver. `"LU"` (the default) handles any assembled matrix and is
+required for the coil-response columns (`gal_rpec_flag`) and matching. `"cholesky"` factors the
+lower band as Hermitian positive definite (LAPACK `zpbtrf`/`zpbtrs`) and stops with an error when
+the matrix is not positive definite (as on the DIII-D-like examples);
+`examples/LAR_gal_cholesky_example` is a case where it applies.
+
 The implementation lives in `src/ForceFreeStates/Galerkin/`:
 
 | File | Content |
diff --git a/examples/LAR_gal_cholesky_example/gpec.toml b/examples/LAR_gal_cholesky_example/gpec.toml
new file mode 100644
index 000000000..45282bcea
--- /dev/null
+++ b/examples/LAR_gal_cholesky_example/gpec.toml
@@ -0,0 +1,64 @@
+# LAR Cholesky Δ′: large-aspect-ratio analytic TJ equilibrium, epsilon = a/R0 = 0.2, n = 1.
+#
+# Outer-region Galerkin Δ′ with the banded Cholesky solver (gal_solver = "cholesky"). On this
+# low-beta circular equilibrium the assembled Galerkin matrix is Hermitian positive definite, so
+# Cholesky applies; it reproduces the LU Δ′ and agrees with the Riccati Δ′ to about 0.5% on both
+# rational surfaces (q = 2, 3). Cholesky refuses a matrix that is not positive definite; use gal_solver = "LU" there.
+# Cholesky cannot carry the coil-response columns, so gal_rpec_flag/gal_match_flag stay off.
+# This run terminates after the stability stage.
+
+[Equilibrium]
+eq_type   = "tj_analytic"      # Type of the input 2D equilibrium (analytic large-aspect-ratio model)
+jac_type  = "hamada"           # Coordinate system (hamada, pest, boozer, equal_arc, park, custom)
+grid_type = "ldp"              # Radial grid packing type (ldp = linear-derivative packing toward rationals)
+psilow    = 0.01               # Lower limit of normalized flux coordinate
+psihigh   = 0.995              # Upper limit of normalized flux coordinate
+mpsi      = 128                # Number of radial grid points (0 = auto-compute from psi_accuracy)
+mtheta    = 512                # Number of poloidal grid points
+
+[TJ_ANALYTIC_INPUT]
+lar_r0 = 2.0                   # Major radius R0 [m]
+lar_a  = 0.4                   # Minor radius a [m] (epsilon = a/R0 = 0.2)
+qc     = 1.5                   # On-axis safety factor
+qa     = 3.6                   # Edge safety factor
+pc     = 0.01                  # Normalized on-axis pressure
+mu     = 2.0                   # Pressure peaking exponent: p2 = pc*(1-r^2)^mu
+B0     = 12.0                  # On-axis toroidal field [T]
+ma     = 128                   # Radial grid points
+mtau   = 128                   # Poloidal grid points
+
+[Wall]
+shape = "conformal"            # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
+a     = 20                     # Distance from plasma (conformal) or shape parameter
+
+[ForceFreeStates]
+vac_flag             = true    # Compute plasma, vacuum, and total energies for free-boundary modes
+
+qlow        = 1.02             # Integration initiated at q determined by min(q0, qlow)
+qhigh       = 3.6              # Integration terminated at q limit determined by min(qa, qhigh)
+sing_start  = 0                # Start integration at the sing_start'th rational from the axis (psilow)
+nn_low      = 1                # Smallest toroidal mode number to include
+nn_high     = 1                # Largest toroidal mode number to include
+delta_mlow  = 8                # Expands lower bound of Fourier harmonics
+delta_mhigh = 8                # Expands upper bound of Fourier harmonics
+mthvac      = 960              # Number of points used in splines over poloidal angle at plasma-vacuum interface
+
+integrator            = "galerkin" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match)
+set_psilim_via_dmlim  = false  # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim)
+force_termination     = true   # Terminate after the stability stage (no perturbed-equilibrium section here)
+write_outputs_to_HDF5 = true   # Write stability outputs to HDF5
+
+gal_solver             = "cholesky" # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
+gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
+gal_nq                 = 6     # Gauss-Lobatto quadrature order per element
+gal_pfac               = 0.03  # Grid packing ratio toward singular surfaces (smaller = tighter packing)
+gal_dx0                = 5e-4  # Resonant-element integration truncation distance from the rational (times 1/|n q'|)
+gal_dx1                = 1e-3  # Resonant-element size (times 1/|n q'|)
+gal_dx2                = 1e-3  # Extension-element size (times 1/|n q'|)
+gal_cutoff             = 10    # Number of elements carrying the large (driving) solution as a source term
+gal_tol                = 1e-10 # Resonant-cell adaptive-quadrature (QuadGK) tolerance
+gal_gnstep             = 5000  # Maximum resonant-integration steps
+gal_dx1dx2_flag        = true  # Enable the special dx1/dx2 element sizing for resonant/extension cells
+gal_sing_order         = 6     # Base power-series order for the Galerkin singular asymptotics
+gal_sing_order_ceiling = true  # Auto-raise the series order by ceil(2*Re(alpha)) on high-Mercier-index surfaces
+gal_rpec_flag          = false # Append mpert coil-response columns (unit edge sources) to the solve, giving the delta_coil block
diff --git a/regression-harness/cases/lar_gal_cholesky.toml b/regression-harness/cases/lar_gal_cholesky.toml
new file mode 100644
index 000000000..b993bca59
--- /dev/null
+++ b/regression-harness/cases/lar_gal_cholesky.toml
@@ -0,0 +1,78 @@
+# Regression case: LAR analytic equilibrium, n=1, outer-region Galerkin Δ′ with the banded
+# Cholesky solver (gal_solver = "cholesky", zpbtrf/zpbtrs). The Galerkin matrix is positive
+# definite here, so this is the Cholesky path's end-to-end track: resonant-surface inventory,
+# the PEST-3 Δ matrix (diagonal and every element), the raw outer Δ′ norm, and the asymptotics.
+# Each [quantities.*] block names an HDF5 path in the run output, how to extract it,
+# and the noise floor below which a difference is treated as zero.
+[case]
+name = "lar_gal_cholesky"
+description = "LAR analytic equilibrium, n=1, Galerkin Δ′ via the banded Cholesky solver"
+example_dir = "examples/LAR_gal_cholesky_example"
+
+# Number and location of the resonant surfaces in the gal domain
+[quantities.gal_msing]
+h5path = "ForceFreeStates/Solutions/GalerkinIntegration/rational_count"
+type = "int_scalar"
+extract = "value"
+label = "gal # singular surfaces"
+noise_threshold = 0
+order = 10
+
+[quantities.gal_sing_q]
+h5path = "ForceFreeStates/Solutions/GalerkinIntegration/rational_q"
+type = "real_vector"
+extract = "all_real"
+label = "gal singular q values"
+noise_threshold = 1e-8
+order = 11
+
+# PEST-3 Δ matching matrix — per-surface diagonal (the tearing Δ′) and every element
+[quantities.gal_pest3_delta_diag]
+h5path = "SingularSurfaces/Delta_prime_matrix"
+type = "complex_matrix"
+extract = "diagonal_complex"
+label = "gal PEST3 Δ diagonal"
+noise_threshold = 1e-6
+order = 20
+
+[quantities.gal_pest3_delta_all]
+h5path = "SingularSurfaces/Delta_prime_matrix"
+type = "complex_matrix"
+extract = "all_complex"
+label = "gal PEST3 Δ (all)"
+noise_threshold = 1e-6
+order = 21
+
+# Raw outer Δ′ matrix (2·msing × 2·msing) — Frobenius norm catches any element drift
+[quantities.gal_delta_norm]
+h5path = "SingularSurfaces/Delta_prime_raw"
+type = "complex_matrix"
+extract = "norm"
+label = "||gal Δ′ matrix||"
+noise_threshold = 1e-6
+order = 22
+
+# Mercier index and resonant exponents per surface
+[quantities.gal_di]
+h5path = "ForceFreeStates/Solutions/GalerkinIntegration/D_I"
+type = "real_vector"
+extract = "all_real"
+label = "gal D_I per surface"
+noise_threshold = 1e-8
+order = 30
+
+[quantities.gal_alpha]
+h5path = "ForceFreeStates/Solutions/GalerkinIntegration/alpha"
+type = "complex_vector"
+extract = "all_complex"
+label = "gal α per surface"
+noise_threshold = 1e-8
+order = 31
+
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999
diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
index 963c107f0..99d734548 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
@@ -150,8 +150,8 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines
         LinearAlgebra.LAPACK.gbtrs!('N', ws.kl, ws.ku, ws.ndim, ab, ipiv, ws.sol)
     else
         ctrl.verbose && @info "Galerkin Cholesky banded factorization + solve"
-        LinearAlgebra.LAPACK.pbtrf!('L', ws.kl, ws.mat)
-        LinearAlgebra.LAPACK.pbtrs!('L', ws.kl, ws.mat, ws.sol)
+        gal_zpbtrf!(ws.mat, ws.kl)
+        gal_zpbtrs!(ws.mat, ws.kl, ws.sol)
     end
 
     # --- extract Δ′ from the small resonant coefficients (gal.f) ---
@@ -340,3 +340,27 @@ function annotate_galerkin!(out_h5)
     Utilities.HDF5Annotations.annotate!(out_h5, GALERKIN_H5_ANNOTATIONS)
     return nothing
 end
+
+# Lower-band Hermitian Cholesky, unscaled, as rdcon gal.f: zpbtrf('L') then zpbtrs('L').
+# LinearAlgebra.LAPACK has no pbtrf!/pbtrs! wrapper, so call LAPACK directly.
+function gal_zpbtrf!(AB::Matrix{ComplexF64}, kd::Int)
+    info = Ref{LinearAlgebra.BlasInt}()
+    ccall((LinearAlgebra.BLAS.@blasfunc(zpbtrf_), LinearAlgebra.libblastrampoline), Cvoid,
+        (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{ComplexF64}, Ref{LinearAlgebra.BlasInt},
+            Ref{LinearAlgebra.BlasInt}, Clong),
+        'L', size(AB, 2), kd, AB, max(1, stride(AB, 2)), info, 1)
+    # gal.f ignores info; a non-positive-definite matrix would leave a partial factor and a wrong Δ′.
+    info[] > 0 && error("gal_solver=\"cholesky\": Galerkin matrix is not positive definite (zpbtrf info=$(info[])); use gal_solver=\"LU\"")
+    LAPACK.chklapackerror(info[])
+    return AB
+end
+
+function gal_zpbtrs!(AB::Matrix{ComplexF64}, kd::Int, B::AbstractVecOrMat{ComplexF64})
+    info = Ref{LinearAlgebra.BlasInt}()
+    ccall((LinearAlgebra.BLAS.@blasfunc(zpbtrs_), LinearAlgebra.libblastrampoline), Cvoid,
+        (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{ComplexF64},
+            Ref{LinearAlgebra.BlasInt}, Ptr{ComplexF64}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Clong),
+        'L', size(AB, 2), kd, size(B, 2), AB, max(1, stride(AB, 2)), B, max(1, stride(B, 2)), info, 1)
+    LAPACK.chklapackerror(info[])
+    return B
+end
diff --git a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl
index f4e79deb0..866c60f2c 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl
@@ -87,7 +87,7 @@ Global workspace for the Galerkin solve. Port of Fortran `gal_type` (gal.f).
 `kl = ku = mpert*(np+1)`. The banded `mat` storage layout depends on `solver` (gal.f):
 
   - `"LU"`       → `ldab = 2*kl + ku + 1`, full band (LAPACK `gbtrf!`/`gbtrs!`)
-  - `"cholesky"` → `ldab = kl + 1`, lower band only (LAPACK `pbtrf!`/`pbtrs!`, `'L'`)
+  - `"cholesky"` → `ldab = kl + 1`, lower band only (LAPACK `zpbtrf`/`zpbtrs`, `'L'`, unscaled as gal.f)
 """
 mutable struct GalWorkspace
     solver::String
diff --git a/test/runtests.jl b/test/runtests.jl
index 389370209..e1b697bf8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -32,6 +32,7 @@ else
     include("./runtests_parallel_integration.jl")
     include("./runtests_result_struct.jl")
     include("./runtests_solve_api.jl")
+    include("./runtests_galerkin_cholesky.jl")
     include("./runtests_dominant_coupling.jl")
     include("./runtests_error_fields.jl")
     include("./runtests_tolerance_toml.jl")
diff --git a/test/runtests_galerkin_cholesky.jl b/test/runtests_galerkin_cholesky.jl
new file mode 100644
index 000000000..1caedbd84
--- /dev/null
+++ b/test/runtests_galerkin_cholesky.jl
@@ -0,0 +1,42 @@
+using Test
+using LinearAlgebra
+
+# gal_solver = "cholesky": unscaled lower-band zpbtrf/zpbtrs (rdcon gal.f) on GalWorkspace's
+# `ldab = kl + 1` lower band storage.
+@testset "Galerkin banded Cholesky" begin
+    FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
+
+    n, kl = 24, 3
+    A = zeros(ComplexF64, n, n)
+    for j in 1:n, i in max(1, j - kl):min(n, j + kl)
+        A[i, j] = i == j ? 4 + sin(j) : 0.3 * cis(0.3i + 0.7j)
+    end
+    A = (A + A') / 2  # Hermitian, diagonally dominant → positive definite
+
+    function lowerband(A)
+        m = size(A, 2)
+        ab = zeros(ComplexF64, kl + 1, m)
+        for j in 1:m, i in j:min(m, j + kl)
+            ab[1+i-j, j] = A[i, j]
+        end
+        return ab
+    end
+
+    x = ComplexF64[cis(0.1j) * (1 + j / n) for j in 1:n]
+    B = A * hcat(x, -3x)
+    ab = lowerband(A)
+    FFS.gal_zpbtrf!(ab, kl)
+    FFS.gal_zpbtrs!(ab, kl, B)
+    @test B[:, 1] ≈ x rtol = 1e-13
+    @test B[:, 2] ≈ -3x rtol = 1e-13
+
+    # Lower factor matches the dense Cholesky factor
+    L = cholesky(Hermitian(A, :L)).L
+    @test all(ab[1+i-j, j] ≈ L[i, j] for j in 1:n for i in j:min(n, j + kl))
+
+    @testset "indefinite matrix is refused" begin
+        Aind = copy(A)
+        Aind[5, 5] = -1
+        @test_throws ErrorException FFS.gal_zpbtrf!(lowerband(Aind), kl)
+    end
+end

From fc5eb1b7956bfcb686dc2170baf9baf10b8c4574 Mon Sep 17 00:00:00 2001
From: Matthew Pharr <m.pharr@protonmail.com>
Date: Fri, 9 Oct 2026 11:04:19 +0200
Subject: [PATCH 2/6] ForceFreeStates.Galerkin - MINOR - Cleanup of
 comments/example cases

---
 examples/LAR_gal_cholesky_example/gpec.toml   | 64 ---------------
 .../cases/lar_gal_cholesky.toml               | 78 -------------------
 .../Galerkin/GalerkinStructs.jl               |  2 +-
 3 files changed, 1 insertion(+), 143 deletions(-)
 delete mode 100644 examples/LAR_gal_cholesky_example/gpec.toml
 delete mode 100644 regression-harness/cases/lar_gal_cholesky.toml

diff --git a/examples/LAR_gal_cholesky_example/gpec.toml b/examples/LAR_gal_cholesky_example/gpec.toml
deleted file mode 100644
index 45282bcea..000000000
--- a/examples/LAR_gal_cholesky_example/gpec.toml
+++ /dev/null
@@ -1,64 +0,0 @@
-# LAR Cholesky Δ′: large-aspect-ratio analytic TJ equilibrium, epsilon = a/R0 = 0.2, n = 1.
-#
-# Outer-region Galerkin Δ′ with the banded Cholesky solver (gal_solver = "cholesky"). On this
-# low-beta circular equilibrium the assembled Galerkin matrix is Hermitian positive definite, so
-# Cholesky applies; it reproduces the LU Δ′ and agrees with the Riccati Δ′ to about 0.5% on both
-# rational surfaces (q = 2, 3). Cholesky refuses a matrix that is not positive definite; use gal_solver = "LU" there.
-# Cholesky cannot carry the coil-response columns, so gal_rpec_flag/gal_match_flag stay off.
-# This run terminates after the stability stage.
-
-[Equilibrium]
-eq_type   = "tj_analytic"      # Type of the input 2D equilibrium (analytic large-aspect-ratio model)
-jac_type  = "hamada"           # Coordinate system (hamada, pest, boozer, equal_arc, park, custom)
-grid_type = "ldp"              # Radial grid packing type (ldp = linear-derivative packing toward rationals)
-psilow    = 0.01               # Lower limit of normalized flux coordinate
-psihigh   = 0.995              # Upper limit of normalized flux coordinate
-mpsi      = 128                # Number of radial grid points (0 = auto-compute from psi_accuracy)
-mtheta    = 512                # Number of poloidal grid points
-
-[TJ_ANALYTIC_INPUT]
-lar_r0 = 2.0                   # Major radius R0 [m]
-lar_a  = 0.4                   # Minor radius a [m] (epsilon = a/R0 = 0.2)
-qc     = 1.5                   # On-axis safety factor
-qa     = 3.6                   # Edge safety factor
-pc     = 0.01                  # Normalized on-axis pressure
-mu     = 2.0                   # Pressure peaking exponent: p2 = pc*(1-r^2)^mu
-B0     = 12.0                  # On-axis toroidal field [T]
-ma     = 128                   # Radial grid points
-mtau   = 128                   # Poloidal grid points
-
-[Wall]
-shape = "conformal"            # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
-a     = 20                     # Distance from plasma (conformal) or shape parameter
-
-[ForceFreeStates]
-vac_flag             = true    # Compute plasma, vacuum, and total energies for free-boundary modes
-
-qlow        = 1.02             # Integration initiated at q determined by min(q0, qlow)
-qhigh       = 3.6              # Integration terminated at q limit determined by min(qa, qhigh)
-sing_start  = 0                # Start integration at the sing_start'th rational from the axis (psilow)
-nn_low      = 1                # Smallest toroidal mode number to include
-nn_high     = 1                # Largest toroidal mode number to include
-delta_mlow  = 8                # Expands lower bound of Fourier harmonics
-delta_mhigh = 8                # Expands upper bound of Fourier harmonics
-mthvac      = 960              # Number of points used in splines over poloidal angle at plasma-vacuum interface
-
-integrator            = "galerkin" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match)
-set_psilim_via_dmlim  = false  # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim)
-force_termination     = true   # Terminate after the stability stage (no perturbed-equilibrium section here)
-write_outputs_to_HDF5 = true   # Write stability outputs to HDF5
-
-gal_solver             = "cholesky" # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
-gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
-gal_nq                 = 6     # Gauss-Lobatto quadrature order per element
-gal_pfac               = 0.03  # Grid packing ratio toward singular surfaces (smaller = tighter packing)
-gal_dx0                = 5e-4  # Resonant-element integration truncation distance from the rational (times 1/|n q'|)
-gal_dx1                = 1e-3  # Resonant-element size (times 1/|n q'|)
-gal_dx2                = 1e-3  # Extension-element size (times 1/|n q'|)
-gal_cutoff             = 10    # Number of elements carrying the large (driving) solution as a source term
-gal_tol                = 1e-10 # Resonant-cell adaptive-quadrature (QuadGK) tolerance
-gal_gnstep             = 5000  # Maximum resonant-integration steps
-gal_dx1dx2_flag        = true  # Enable the special dx1/dx2 element sizing for resonant/extension cells
-gal_sing_order         = 6     # Base power-series order for the Galerkin singular asymptotics
-gal_sing_order_ceiling = true  # Auto-raise the series order by ceil(2*Re(alpha)) on high-Mercier-index surfaces
-gal_rpec_flag          = false # Append mpert coil-response columns (unit edge sources) to the solve, giving the delta_coil block
diff --git a/regression-harness/cases/lar_gal_cholesky.toml b/regression-harness/cases/lar_gal_cholesky.toml
deleted file mode 100644
index b993bca59..000000000
--- a/regression-harness/cases/lar_gal_cholesky.toml
+++ /dev/null
@@ -1,78 +0,0 @@
-# Regression case: LAR analytic equilibrium, n=1, outer-region Galerkin Δ′ with the banded
-# Cholesky solver (gal_solver = "cholesky", zpbtrf/zpbtrs). The Galerkin matrix is positive
-# definite here, so this is the Cholesky path's end-to-end track: resonant-surface inventory,
-# the PEST-3 Δ matrix (diagonal and every element), the raw outer Δ′ norm, and the asymptotics.
-# Each [quantities.*] block names an HDF5 path in the run output, how to extract it,
-# and the noise floor below which a difference is treated as zero.
-[case]
-name = "lar_gal_cholesky"
-description = "LAR analytic equilibrium, n=1, Galerkin Δ′ via the banded Cholesky solver"
-example_dir = "examples/LAR_gal_cholesky_example"
-
-# Number and location of the resonant surfaces in the gal domain
-[quantities.gal_msing]
-h5path = "ForceFreeStates/Solutions/GalerkinIntegration/rational_count"
-type = "int_scalar"
-extract = "value"
-label = "gal # singular surfaces"
-noise_threshold = 0
-order = 10
-
-[quantities.gal_sing_q]
-h5path = "ForceFreeStates/Solutions/GalerkinIntegration/rational_q"
-type = "real_vector"
-extract = "all_real"
-label = "gal singular q values"
-noise_threshold = 1e-8
-order = 11
-
-# PEST-3 Δ matching matrix — per-surface diagonal (the tearing Δ′) and every element
-[quantities.gal_pest3_delta_diag]
-h5path = "SingularSurfaces/Delta_prime_matrix"
-type = "complex_matrix"
-extract = "diagonal_complex"
-label = "gal PEST3 Δ diagonal"
-noise_threshold = 1e-6
-order = 20
-
-[quantities.gal_pest3_delta_all]
-h5path = "SingularSurfaces/Delta_prime_matrix"
-type = "complex_matrix"
-extract = "all_complex"
-label = "gal PEST3 Δ (all)"
-noise_threshold = 1e-6
-order = 21
-
-# Raw outer Δ′ matrix (2·msing × 2·msing) — Frobenius norm catches any element drift
-[quantities.gal_delta_norm]
-h5path = "SingularSurfaces/Delta_prime_raw"
-type = "complex_matrix"
-extract = "norm"
-label = "||gal Δ′ matrix||"
-noise_threshold = 1e-6
-order = 22
-
-# Mercier index and resonant exponents per surface
-[quantities.gal_di]
-h5path = "ForceFreeStates/Solutions/GalerkinIntegration/D_I"
-type = "real_vector"
-extract = "all_real"
-label = "gal D_I per surface"
-noise_threshold = 1e-8
-order = 30
-
-[quantities.gal_alpha]
-h5path = "ForceFreeStates/Solutions/GalerkinIntegration/alpha"
-type = "complex_vector"
-extract = "all_complex"
-label = "gal α per surface"
-noise_threshold = 1e-8
-order = 31
-
-[quantities.runtime]
-h5path = ""
-type = "runtime"
-extract = "value"
-label = "Runtime (s)"
-noise_threshold = 0.0
-order = 999
diff --git a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl
index 866c60f2c..fd85df1e5 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl
@@ -87,7 +87,7 @@ Global workspace for the Galerkin solve. Port of Fortran `gal_type` (gal.f).
 `kl = ku = mpert*(np+1)`. The banded `mat` storage layout depends on `solver` (gal.f):
 
   - `"LU"`       → `ldab = 2*kl + ku + 1`, full band (LAPACK `gbtrf!`/`gbtrs!`)
-  - `"cholesky"` → `ldab = kl + 1`, lower band only (LAPACK `zpbtrf`/`zpbtrs`, `'L'`, unscaled as gal.f)
+  - `"cholesky"` → `ldab = kl + 1`, lower band only (LAPACK `zpbtrf`/`zpbtrs`, `'L'`)
 """
 mutable struct GalWorkspace
     solver::String

From 12456572bec5a0046bbb4b06e9da5dd6adf60ad3 Mon Sep 17 00:00:00 2001
From: Matthew Pharr <m.pharr@protonmail.com>
Date: Fri, 9 Oct 2026 12:34:42 +0200
Subject: [PATCH 3/6] ForceFreeStates.galerkin - FEATURE - make coil RHS
 positive definite to enable Cholesky solver for ncoil > 0

---
 src/ForceFreeStates/Galerkin/GalerkinAssembly.jl | 6 ++++++
 src/ForceFreeStates/Galerkin/GalerkinSolve.jl    | 4 ----
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl
index 1ad84cb8b..1fbb94b78 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl
@@ -468,6 +468,12 @@ function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing,
     cell = ws.intvl[msing+1].cells[ws.nx]
     if ncoil > 0
         # rpec: identity edge + unit coil sources (gal.f)
+        # Cholesky also zeroes the edge column, so first lift it into the coil RHS: rhs_I -= A_IB·e_k.
+        if chol
+            for k in 1:mpert, ip in 1:3, ipert in 1:mpert
+                ws.rhs[cell.map[ipert, ip], 2 * msing + k] -= cell.mat[ipert, k, ip, 4]
+            end
+        end
         cell.mat[:, :, 4, :] .= 0
         chol && (cell.mat[:, :, :, 4] .= 0)
         for idx in 1:mpert
diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
index 8a77510ab..b904920af 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
@@ -100,11 +100,7 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines
     ldab = ctrl.gal_solver == "LU" ? 2kl + ku + 1 : kl + 1
     # rpec_flag (RDCON, gal.f): append mpert coil-response columns. Each is a unit source at
     # the edge value DOF in one poloidal mode; the recorded plasma response is the coil block of Δ_gw.
-    # Cholesky's lower-only edge zeroing can't represent the rpec identity edge, so require LU.
     ncoil = ctrl.gal_rpec_flag ? mpert : 0
-    if ncoil > 0 && ctrl.gal_solver != "LU"
-        error("galerkin_solve: gal_rpec_flag=true requires gal_solver=\"LU\" (coil edge BC needs the full-band path)")
-    end
     nsol = 2 * msing + ncoil
     intvl = [GalInterval(zeros(Float64, nx + 1), zeros(Float64, nx + 1), [GalCell(mpert) for _ in 1:nx])
              for _ in 0:msing]

From e01619e157d2954ee7727143b606aafcb0dbf0df Mon Sep 17 00:00:00 2001
From: StuartBenjamin <stubenj9@mit.edu>
Date: Mon, 12 Oct 2026 16:43:15 -0700
Subject: [PATCH 4/6] ForceFreeStates.Galerkin - FEATURE - Try the banded
 Cholesky solve first and fall back to LU

gal_solver defaults to "auto", following @matt-pharr's review recommendation
on #513: the full-band Galerkin matrix is solved by banded Hermitian Cholesky
when it is Hermitian positive definite (Dirichlet rows' columns moved to the
right-hand side), else by the Jacobi-scaled LU. "LU" skips the attempt;
"cholesky" errors instead of falling back. The matrix is always assembled
full band, so the lower-band assembly branches are gone. The test now checks
the solver choice and the solution against a dense solve.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 docs/src/galerkin.md                          | 10 +--
 .../gpec.toml                                 |  2 +-
 .../gpec.toml                                 |  2 +-
 examples/LAR_ideal_match_test/gpec.toml       |  2 +-
 examples/LAR_resistive_match_test/gpec.toml   |  2 +-
 src/ForceFreeStates/CoreTypes.jl              |  2 +-
 .../Galerkin/GalerkinAssembly.jl              | 34 ++------
 src/ForceFreeStates/Galerkin/GalerkinSolve.jl | 78 +++++++++++++++----
 .../Galerkin/GalerkinStructs.jl               |  6 +-
 test/runtests_galerkin_cholesky.jl            | 72 +++++++++++------
 10 files changed, 126 insertions(+), 84 deletions(-)

diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md
index 806fd9ef3..b77bd25e1 100644
--- a/docs/src/galerkin.md
+++ b/docs/src/galerkin.md
@@ -20,11 +20,11 @@ control surface (when `vac_flag`) and produces no free-boundary energies or ODE
 Setting `gal_match_flag` additionally matches the inner layer, giving a driven ξ solution that
 `PerturbedEquilibrium` consumes in place of a forward solution.
 
-`gal_solver` picks the banded solver. `"LU"` (the default) handles any assembled matrix and is
-required for the coil-response columns (`gal_rpec_flag`) and matching. `"cholesky"` factors the
-lower band as Hermitian positive definite (LAPACK `zpbtrf`/`zpbtrs`) and stops with an error when
-the matrix is not positive definite (as on the DIII-D-like examples);
-`examples/LAR_gal_cholesky_example` is a case where it applies.
+`gal_solver` picks the banded solver. `"auto"` (the default) tries a Hermitian Cholesky
+factorization of the lower band (LAPACK `zpbtrf`/`zpbtrs`, with the Dirichlet rows' columns moved
+to the right-hand side) and falls back to Jacobi-scaled LU when the matrix is not Hermitian
+positive definite (as on the DIII-D-like examples). `"LU"` skips the Cholesky attempt;
+`"cholesky"` stops with an error instead of falling back.
 
 The implementation lives in `src/ForceFreeStates/Galerkin/`:
 
diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml
index 71d16afcb..890137170 100644
--- a/examples/DIIID-like_gal_resistive_example/gpec.toml
+++ b/examples/DIIID-like_gal_resistive_example/gpec.toml
@@ -57,7 +57,7 @@ dmlim               = 0.2      # Truncate integration at (last_rational_q + dmli
 
 # Outer-region singular Galerkin Δ′ solver with rpec coil columns.
 # This case exercises the gal solve AND the rpec coil/edge path (delta_coil) for regression tracking.
-gal_solver             = "LU"  # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
+gal_solver             = "auto"  # Banded solver: "auto" (Cholesky if Hermitian positive definite, else LU), "LU" or "cholesky"
 gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
 gal_nq                 = 6     # Gauss-Lobatto quadrature order per element
 gal_pfac               = 0.03  # Grid packing ratio toward singular surfaces (smaller = tighter packing)
diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
index 71e584905..4249bbe12 100644
--- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
+++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
@@ -57,7 +57,7 @@ dmlim               = 0.2      # Truncate integration at (last_rational_q + dmli
 
 # Outer-region singular Galerkin Δ′ solver with rpec coil columns.
 # This case exercises the gal solve AND the rpec coil/edge path (delta_coil) for regression tracking.
-gal_solver             = "LU"  # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
+gal_solver             = "auto"  # Banded solver: "auto" (Cholesky if Hermitian positive definite, else LU), "LU" or "cholesky"
 gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
 gal_nq                 = 6     # Gauss-Lobatto quadrature order per element
 gal_pfac               = 0.03  # Grid packing ratio toward singular surfaces (smaller = tighter packing)
diff --git a/examples/LAR_ideal_match_test/gpec.toml b/examples/LAR_ideal_match_test/gpec.toml
index 0706863b7..6b4a808cf 100644
--- a/examples/LAR_ideal_match_test/gpec.toml
+++ b/examples/LAR_ideal_match_test/gpec.toml
@@ -57,7 +57,7 @@ force_termination     = true   # Terminate after the stability stage (no perturb
 write_outputs_to_HDF5 = true   # Write stability outputs to HDF5
 
 # Outer-region singular Galerkin Delta-prime solver with rpec coil columns.
-gal_solver             = "LU"  # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
+gal_solver             = "LU"  # Banded solver: "auto" (Cholesky if Hermitian positive definite, else LU), "LU" or "cholesky"
 gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
 gal_nq                 = 6     # Gauss-Lobatto quadrature order per element
 gal_pfac               = 0.03  # Grid packing ratio toward singular surfaces (smaller = tighter packing)
diff --git a/examples/LAR_resistive_match_test/gpec.toml b/examples/LAR_resistive_match_test/gpec.toml
index 0571301c7..94106ee0d 100644
--- a/examples/LAR_resistive_match_test/gpec.toml
+++ b/examples/LAR_resistive_match_test/gpec.toml
@@ -58,7 +58,7 @@ force_termination     = true   # Terminate after the stability stage (no perturb
 write_outputs_to_HDF5 = true   # Write stability outputs to HDF5
 
 # Outer-region singular Galerkin Delta-prime solver with rpec coil columns.
-gal_solver             = "LU"  # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
+gal_solver             = "LU"  # Banded solver: "auto" (Cholesky if Hermitian positive definite, else LU), "LU" or "cholesky"
 gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
 gal_nq                 = 6     # Gauss-Lobatto quadrature order per element
 gal_pfac               = 0.03  # Grid packing ratio toward singular surfaces (smaller = tighter packing)
diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl
index 5544b0e60..5ebbe08d9 100644
--- a/src/ForceFreeStates/CoreTypes.jl
+++ b/src/ForceFreeStates/CoreTypes.jl
@@ -192,7 +192,7 @@ gpec.toml.
     extended_precision_bvp::Bool = true   # Promote Δ' BVP to Complex{Double64}; default on (Float64 drifts the imaginary Δ' by 2–5× on DIIID-class cases).
 
     # --- RDCON outer-region Galerkin Δ′ solver (gal_solve port); selected by integrator = "galerkin" ---
-    gal_solver::String = "LU"       # "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs)
+    gal_solver::String = "auto"     # "auto" (banded Cholesky if the matrix is Hermitian positive definite, else LU), "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs, error if not positive definite)
     gal_nx::Int = 256               # elements per interval between singular surfaces
     gal_nq::Int = 6                 # Gauss-Lobatto quadrature order per element
     gal_pfac::Float64 = 0.001       # grid packing ratio near singular surfaces
diff --git a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl
index 1fbb94b78..356c86958 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl
@@ -359,15 +359,13 @@ end
 """
     gal_assemble_mat!(ws, mpert)
 
-Scatter every cell's `mat`/`emat`/`ediag` into the banded global matrix `ws.mat`. Port of
-`gal_assemble_mat` (gal.f). Supports both `"LU"` (full band, `offset = kl+ku+1`) and
-`"cholesky"` (lower band only, `offset = 1`).
+Scatter every cell's `mat`/`emat`/`ediag` into the full-band global matrix `ws.mat`
+(`offset = kl+ku+1`). Port of `gal_assemble_mat` (gal.f).
 """
 function gal_assemble_mat!(ws::GalWorkspace, mpert::Int)
     np = ws.np
     msing = length(ws.intvl) - 1
-    lu = ws.solver == "LU"
-    offset = lu ? ws.kl + ws.ku + 1 : 1
+    offset = ws.kl + ws.ku + 1
 
     for ising in 0:msing
         for cell in ws.intvl[ising+1].cells
@@ -377,24 +375,14 @@ function gal_assemble_mat!(ws::GalWorkspace, mpert::Int)
                 # nonresonant stiffness
                 for jp in 0:np, jpert in 1:mpert
                     j = cell.map[jpert, jp+1]
-                    if lu
-                        ws.mat[offset+i-j, j] += cell.mat[ipert, jpert, ip+1, jp+1]
-                    elseif j <= i
-                        ws.mat[offset+i-j, j] += cell.mat[ipert, jpert, ip+1, jp+1]
-                    end
+                    ws.mat[offset+i-j, j] += cell.mat[ipert, jpert, ip+1, jp+1]
                 end
                 # resonant coupling (res/ext only)
                 if res
                     j = cell.emap
                     e = cell.emat[ipert, ip+1]
-                    if lu
-                        ws.mat[offset+i-j, j] += e
-                        ws.mat[offset+j-i, i] += conj(e)
-                    elseif j < i
-                        ws.mat[offset+i-j, j] += e
-                    elseif j > i
-                        ws.mat[offset+j-i, i] += conj(e)
-                    end
+                    ws.mat[offset+i-j, j] += e
+                    ws.mat[offset+j-i, i] += conj(e)
                 end
             end
             if res
@@ -452,14 +440,12 @@ including its full three-way edge branch (rpec count `ncoil = ws.nsol - 2*msing`
 `wv_edge` is `nothing` for the rpec and fixed cases.
 """
 function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing,Matrix{ComplexF64}})
-    chol = ws.solver == "cholesky"
     msing = length(ws.intvl) - 1
     ncoil = ws.nsol - 2 * msing
 
     # axis u(0): first cell of interval 0, Hermite value DOF ip=0 (Julia index 1)
     cell = ws.intvl[1].cells[1]
     cell.mat[:, :, 1, :] .= 0
-    chol && (cell.mat[:, :, :, 1] .= 0)
     for idx in 1:mpert
         cell.mat[idx, idx, 1, 1] = 1
     end
@@ -468,14 +454,7 @@ function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing,
     cell = ws.intvl[msing+1].cells[ws.nx]
     if ncoil > 0
         # rpec: identity edge + unit coil sources (gal.f)
-        # Cholesky also zeroes the edge column, so first lift it into the coil RHS: rhs_I -= A_IB·e_k.
-        if chol
-            for k in 1:mpert, ip in 1:3, ipert in 1:mpert
-                ws.rhs[cell.map[ipert, ip], 2 * msing + k] -= cell.mat[ipert, k, ip, 4]
-            end
-        end
         cell.mat[:, :, 4, :] .= 0
-        chol && (cell.mat[:, :, :, 4] .= 0)
         for idx in 1:mpert
             cell.mat[idx, idx, 4, 4] = 1
         end
@@ -488,7 +467,6 @@ function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing,
     else
         # fixed boundary: identity edge (gal.f)
         cell.mat[:, :, 4, :] .= 0
-        chol && (cell.mat[:, :, :, 4] .= 0)
         for idx in 1:mpert
             cell.mat[idx, idx, 4, 4] = 1
         end
diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
index b904920af..4211e9f3b 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
@@ -60,7 +60,8 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines
     intr::ForceFreeStatesInternal; wv=nothing)
 
     intr.npert == 1 || error("galerkin_solve: only single-n (npert == 1) is supported")
-    ctrl.gal_solver in ("LU", "cholesky") || error("galerkin_solve: gal_solver must be \"LU\" or \"cholesky\"")
+    ctrl.gal_solver in ("auto", "LU", "cholesky") ||
+        error("galerkin_solve: gal_solver must be \"auto\", \"LU\" or \"cholesky\"")
     nn = intr.nlow
     mpert = intr.numpert_total
     np = GAL_NP
@@ -97,7 +98,7 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines
     nx = ctrl.gal_nx
     kl = mpert * (np + 1)
     ku = kl
-    ldab = ctrl.gal_solver == "LU" ? 2kl + ku + 1 : kl + 1
+    ldab = 2kl + ku + 1
     # rpec_flag (RDCON, gal.f): append mpert coil-response columns. Each is a unit source at
     # the edge value DOF in one poloidal mode; the recorded plasma response is the coil block of Δ_gw.
     ncoil = ctrl.gal_rpec_flag ? mpert : 0
@@ -127,7 +128,7 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines
     gal_make_arrays!(ws, ctrl, equil, mats, intr, asymps, sings, nn, wv_edge)
 
     if ctrl.verbose
-        offdbg = ws.solver == "LU" ? ws.kl + ws.ku + 1 : 1
+        offdbg = ws.kl + ws.ku + 1
         @info "Galerkin assembly: ndim=$(ws.ndim) nsol=$nsol kl=$(ws.kl) ldab=$(ws.ldab) |rhs|=$(norm(ws.rhs))"
         for ising in 0:msing
             for cell in ws.intvl[ising+1].cells
@@ -140,14 +141,8 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines
 
     # --- banded solve (gal.f) ---
     ws.sol .= ws.rhs
-    if ws.solver == "LU"
-        ctrl.verbose && @info "Galerkin LU banded factorization + solve"
-        gal_scaled_lu_solve!(ws.mat, ws.sol, ws.kl, ws.ku)
-    else
-        ctrl.verbose && @info "Galerkin Cholesky banded factorization + solve"
-        gal_zpbtrf!(ws.mat, ws.kl)
-        gal_zpbtrs!(ws.mat, ws.kl, ws.sol)
-    end
+    ws.solver = gal_banded_solve!(ws.mat, ws.sol, ws.kl, ws.ku, ctrl.gal_solver)
+    ctrl.verbose && @info "Galerkin banded solve: $(ws.solver)"
 
     # --- extract Δ′ from the small resonant coefficients (gal.f) ---
     delta = zeros(ComplexF64, nsol, 2 * msing)
@@ -358,18 +353,67 @@ function annotate_galerkin!(out_h5)
     return nothing
 end
 
-# Lower-band Hermitian Cholesky, unscaled, as rdcon gal.f: zpbtrf('L') then zpbtrs('L').
-# LinearAlgebra.LAPACK has no pbtrf!/pbtrs! wrapper, so call LAPACK directly.
+"""
+    gal_banded_solve!(mat, sol, kl, ku, solver) -> String
+
+Solve the full-band system `mat` (LAPACK `gbtrf` storage) in place of `sol` and return the
+solver used. `"auto"` tries `gal_cholesky_solve!` and falls back to `gal_scaled_lu_solve!` when
+the matrix is not Hermitian positive definite; `"cholesky"` errors instead of falling back.
+"""
+function gal_banded_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{ComplexF64}, kl::Int, ku::Int, solver::String)
+    if solver != "LU"
+        gal_cholesky_solve!(mat, sol, kl, ku) && return "cholesky"
+        solver == "cholesky" && error("gal_solver=\"cholesky\": the Galerkin matrix is not Hermitian positive definite; use \"auto\" or \"LU\"")
+    end
+    gal_scaled_lu_solve!(mat, sol, kl, ku)
+    return "LU"
+end
+
+"""
+    gal_cholesky_solve!(mat, sol, kl, ku) -> Bool
+
+Banded Cholesky solve of the full-band system `mat`, leaving `mat` untouched. Rows that are
+identity rows (Dirichlet DOFs: axis, fixed or coil edge) have their columns moved to the
+right-hand side, which keeps the system Hermitian. Returns `false`, with `sol` unchanged, when
+the remaining matrix is not Hermitian (to 1e-12 of its largest entry) or not positive definite.
+"""
+function gal_cholesky_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{ComplexF64}, kl::Int, ku::Int)
+    n = size(mat, 2)
+    off = kl + ku + 1
+    a(i, j) = mat[off+i-j, j]
+    band(i) = max(1, i - kl):min(n, i + ku)
+    fixed = [all(j == i ? a(i, j) == 1 : iszero(a(i, j)) for j in band(i)) for i in 1:n]
+    tol = 1e-12 * maximum(abs, mat)
+    for j in 1:n, i in (j+1):min(n, j + kl)
+        fixed[i] || fixed[j] || abs(a(i, j) - conj(a(j, i))) <= tol || return false
+    end
+    ab = zeros(ComplexF64, kl + 1, n)
+    for j in 1:n, i in j:min(n, j + kl)
+        ab[1+i-j, j] = fixed[i] || fixed[j] ? (i == j ? one(ComplexF64) : zero(ComplexF64)) : a(i, j)
+    end
+    gal_zpbtrf!(ab, kl) || return false
+    x = copy(sol)
+    for c in 1:n
+        fixed[c] || continue
+        for r in band(c)
+            r == c || fixed[r] || (@views x[r, :] .-= a(r, c) .* x[c, :])
+        end
+    end
+    gal_zpbtrs!(ab, kl, x)
+    sol .= x
+    return true
+end
+
+# Lower-band Hermitian Cholesky as rdcon gal.f: zpbtrf('L') then zpbtrs('L'). LinearAlgebra.LAPACK
+# has no pbtrf!/pbtrs! wrapper, so call LAPACK directly. Returns false if not positive definite.
 function gal_zpbtrf!(AB::Matrix{ComplexF64}, kd::Int)
     info = Ref{LinearAlgebra.BlasInt}()
     ccall((LinearAlgebra.BLAS.@blasfunc(zpbtrf_), LinearAlgebra.libblastrampoline), Cvoid,
         (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{ComplexF64}, Ref{LinearAlgebra.BlasInt},
             Ref{LinearAlgebra.BlasInt}, Clong),
         'L', size(AB, 2), kd, AB, max(1, stride(AB, 2)), info, 1)
-    # gal.f ignores info; a non-positive-definite matrix would leave a partial factor and a wrong Δ′.
-    info[] > 0 && error("gal_solver=\"cholesky\": Galerkin matrix is not positive definite (zpbtrf info=$(info[])); use gal_solver=\"LU\"")
-    LAPACK.chklapackerror(info[])
-    return AB
+    LAPACK.chklapackerror(min(info[], 0))
+    return info[] == 0
 end
 
 function gal_zpbtrs!(AB::Matrix{ComplexF64}, kd::Int, B::AbstractVecOrMat{ComplexF64})
diff --git a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl
index fd85df1e5..a1de3187a 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl
@@ -84,10 +84,8 @@ end
 Global workspace for the Galerkin solve. Port of Fortran `gal_type` (gal.f).
 
 `intvl` has length `msing+1` (Fortran `intvl(0:msing)`, accessed as `intvl[ising+1]`).
-`kl = ku = mpert*(np+1)`. The banded `mat` storage layout depends on `solver` (gal.f):
-
-  - `"LU"`       → `ldab = 2*kl + ku + 1`, full band (LAPACK `gbtrf!`/`gbtrs!`)
-  - `"cholesky"` → `ldab = kl + 1`, lower band only (LAPACK `zpbtrf`/`zpbtrs`, `'L'`)
+`kl = ku = mpert*(np+1)`. `mat` is stored full band, `ldab = 2*kl + ku + 1` (LAPACK `gbtrf!`
+layout); `solver` records the factorization `gal_banded_solve!` used (`"cholesky"` or `"LU"`).
 """
 mutable struct GalWorkspace
     solver::String
diff --git a/test/runtests_galerkin_cholesky.jl b/test/runtests_galerkin_cholesky.jl
index 1caedbd84..f9c815b6e 100644
--- a/test/runtests_galerkin_cholesky.jl
+++ b/test/runtests_galerkin_cholesky.jl
@@ -1,42 +1,64 @@
 using Test
 using LinearAlgebra
 
-# gal_solver = "cholesky": unscaled lower-band zpbtrf/zpbtrs (rdcon gal.f) on GalWorkspace's
-# `ldab = kl + 1` lower band storage.
-@testset "Galerkin banded Cholesky" begin
+# gal_banded_solve! on full-band (gbtrf layout) systems built like the Galerkin matrix: Hermitian
+# with Dirichlet identity rows (row zeroed, column kept) and nonzero prescribed values.
+@testset "Galerkin banded solve: Cholesky first, LU fallback" begin
     FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
 
     n, kl = 24, 3
-    A = zeros(ComplexF64, n, n)
+    H = zeros(ComplexF64, n, n)
     for j in 1:n, i in max(1, j - kl):min(n, j + kl)
-        A[i, j] = i == j ? 4 + sin(j) : 0.3 * cis(0.3i + 0.7j)
+        H[i, j] = i == j ? 4 + sin(j) : 0.3 * cis(0.3i + 0.7j)
     end
-    A = (A + A') / 2  # Hermitian, diagonally dominant → positive definite
+    H = (H + H') / 2  # Hermitian, diagonally dominant → positive definite
+    fixed = [1, 2, n]
 
-    function lowerband(A)
-        m = size(A, 2)
-        ab = zeros(ComplexF64, kl + 1, m)
-        for j in 1:m, i in j:min(m, j + kl)
-            ab[1+i-j, j] = A[i, j]
+    # Dirichlet rows as gal_set_boundary! leaves them: identity row, column untouched
+    function dirichlet(A)
+        A = copy(A)
+        for i in fixed
+            A[i, :] .= 0
+            A[i, i] = 1
+        end
+        return A
+    end
+    function fullband(A)
+        ab = zeros(ComplexF64, 3kl + 1, n)
+        for j in 1:n, i in max(1, j - kl):min(n, j + kl)
+            ab[2kl+1+i-j, j] = A[i, j]
         end
         return ab
     end
+    rhs = ComplexF64[cis(0.1j) * (1 + j / n) for j in 1:n, _ in 1:2]
+    rhs[fixed, 2] .= [1.0, -0.5im, 2.0]  # prescribed values, as in the rpec coil columns
 
-    x = ComplexF64[cis(0.1j) * (1 + j / n) for j in 1:n]
-    B = A * hcat(x, -3x)
-    ab = lowerband(A)
-    FFS.gal_zpbtrf!(ab, kl)
-    FFS.gal_zpbtrs!(ab, kl, B)
-    @test B[:, 1] ≈ x rtol = 1e-13
-    @test B[:, 2] ≈ -3x rtol = 1e-13
+    function check(A, solver)
+        sol = copy(rhs)
+        used = FFS.gal_banded_solve!(fullband(A), sol, kl, kl, solver)
+        @test sol ≈ A \ rhs rtol = 1e-12
+        return used
+    end
 
-    # Lower factor matches the dense Cholesky factor
-    L = cholesky(Hermitian(A, :L)).L
-    @test all(ab[1+i-j, j] ≈ L[i, j] for j in 1:n for i in j:min(n, j + kl))
+    @test check(dirichlet(H), "auto") == "cholesky"
+    @test check(dirichlet(H), "cholesky") == "cholesky"
+    @test check(dirichlet(H), "LU") == "LU"
 
-    @testset "indefinite matrix is refused" begin
-        Aind = copy(A)
-        Aind[5, 5] = -1
-        @test_throws ErrorException FFS.gal_zpbtrf!(lowerband(Aind), kl)
+    indefinite = copy(H)
+    indefinite[5, 5] = -1
+    @test check(dirichlet(indefinite), "auto") == "LU"
+    @test_throws ErrorException FFS.gal_banded_solve!(fullband(dirichlet(indefinite)), copy(rhs), kl, kl, "cholesky")
+
+    nonhermitian = copy(H)
+    nonhermitian[6, 4] += 0.1
+    @test check(dirichlet(nonhermitian), "auto") == "LU"
+
+    # the lower-band factor equals the dense Cholesky factor
+    ab = zeros(ComplexF64, kl + 1, n)
+    for j in 1:n, i in j:min(n, j + kl)
+        ab[1+i-j, j] = H[i, j]
     end
+    @test FFS.gal_zpbtrf!(ab, kl)
+    L = cholesky(Hermitian(H, :L)).L
+    @test all(ab[1+i-j, j] ≈ L[i, j] for j in 1:n for i in j:min(n, j + kl))
 end

From b88037f4e037bde500d2f08a0d097799882a52fd Mon Sep 17 00:00:00 2001
From: StuartBenjamin <stubenj9@mit.edu>
Date: Mon, 12 Oct 2026 20:15:45 -0700
Subject: [PATCH 5/6] ForceFreeStates.Galerkin - BUGFIX - Test the Cholesky
 precondition entry by entry, diagonal included
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The Hermiticity check before the banded Cholesky skipped the diagonal and used one tolerance,
1e-12 of the largest entry. zpbtrf reads only the real part of the diagonal, so a complex diagonal
(the resonant elements' conj(u)·Lu) was dropped silently; and since the diagonal spans many decades,
the small rows passed any asymmetry. Each entry is now compared with 1e-10·√|aᵢᵢaⱼⱼ|, diagonal
included (the quadrature round-off is ~1e-11 of that scale). gal_cholesky_solve! requires kl == ku and solves in place of sol once zpbtrf succeeds.
Galerkin() now defaults to solver = "auto" like gal_solver, and its docstring lists "auto".

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 docs/src/galerkin.md                          |  2 +-
 .../Galerkin/GalerkinAssembly.jl              |  2 +-
 src/ForceFreeStates/Galerkin/GalerkinSolve.jl | 20 +++++++++----------
 src/ForceFreeStates/Integrators.jl            |  6 +++---
 test/runtests_galerkin_cholesky.jl            |  8 ++++++++
 5 files changed, 23 insertions(+), 15 deletions(-)

diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md
index b77bd25e1..c787323a3 100644
--- a/docs/src/galerkin.md
+++ b/docs/src/galerkin.md
@@ -23,7 +23,7 @@ Setting `gal_match_flag` additionally matches the inner layer, giving a driven 
 `gal_solver` picks the banded solver. `"auto"` (the default) tries a Hermitian Cholesky
 factorization of the lower band (LAPACK `zpbtrf`/`zpbtrs`, with the Dirichlet rows' columns moved
 to the right-hand side) and falls back to Jacobi-scaled LU when the matrix is not Hermitian
-positive definite (as on the DIII-D-like examples). `"LU"` skips the Cholesky attempt;
+positive definite. `"LU"` skips the Cholesky attempt;
 `"cholesky"` stops with an error instead of falling back.
 
 The implementation lives in `src/ForceFreeStates/Galerkin/`:
diff --git a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl
index 356c86958..593a15ef2 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl
@@ -6,7 +6,7 @@
 #   gal_gauss_quad    (gal.f) nonresonant stiffness via Gauss-Lobatto quadrature
 #   gal_extension     (gal.f) extension-cell matrix/rhs (big-solution driving + surface terms)
 #   gal_resonant      (gal.f)   resonant-cell integral (QuadGK in place of Fortran LSODE)
-#   gal_assemble_mat  (gal.f)   scatter cell blocks into the banded global matrix (LU/Cholesky)
+#   gal_assemble_mat  (gal.f)   scatter cell blocks into the full-band global matrix
 #   gal_assemble_rhs  (gal.f)   scatter cell rhs into the global RHS
 #   gal_set_boundary  (gal.f) axis identity + edge (vacuum or fixed) BCs
 #
diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
index 4211e9f3b..025ad4e8c 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
@@ -375,32 +375,32 @@ end
 Banded Cholesky solve of the full-band system `mat`, leaving `mat` untouched. Rows that are
 identity rows (Dirichlet DOFs: axis, fixed or coil edge) have their columns moved to the
 right-hand side, which keeps the system Hermitian. Returns `false`, with `sol` unchanged, when
-the remaining matrix is not Hermitian (to 1e-12 of its largest entry) or not positive definite.
+the remaining matrix is not Hermitian (each entry to 1e-10 of √|aᵢᵢaⱼⱼ|, diagonal included) or not
+positive definite.
 """
 function gal_cholesky_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{ComplexF64}, kl::Int, ku::Int)
+    kl == ku || return false
     n = size(mat, 2)
-    off = kl + ku + 1
+    off = 2kl + 1
     a(i, j) = mat[off+i-j, j]
-    band(i) = max(1, i - kl):min(n, i + ku)
+    band(i) = max(1, i - kl):min(n, i + kl)
     fixed = [all(j == i ? a(i, j) == 1 : iszero(a(i, j)) for j in band(i)) for i in 1:n]
-    tol = 1e-12 * maximum(abs, mat)
-    for j in 1:n, i in (j+1):min(n, j + kl)
-        fixed[i] || fixed[j] || abs(a(i, j) - conj(a(j, i))) <= tol || return false
+    # Scale-free test (the diagonal spans many decades); 1e-10 clears the quadrature round-off (~1e-11)
+    for j in 1:n, i in j:min(n, j + kl)
+        fixed[i] || fixed[j] || abs(a(i, j) - conj(a(j, i))) <= 1e-10 * sqrt(abs(a(i, i) * a(j, j))) || return false
     end
     ab = zeros(ComplexF64, kl + 1, n)
     for j in 1:n, i in j:min(n, j + kl)
         ab[1+i-j, j] = fixed[i] || fixed[j] ? (i == j ? one(ComplexF64) : zero(ComplexF64)) : a(i, j)
     end
     gal_zpbtrf!(ab, kl) || return false
-    x = copy(sol)
     for c in 1:n
         fixed[c] || continue
         for r in band(c)
-            r == c || fixed[r] || (@views x[r, :] .-= a(r, c) .* x[c, :])
+            r == c || fixed[r] || (@views sol[r, :] .-= a(r, c) .* sol[c, :])
         end
     end
-    gal_zpbtrs!(ab, kl, x)
-    sol .= x
+    gal_zpbtrs!(ab, kl, sol)
     return true
 end
 
diff --git a/src/ForceFreeStates/Integrators.jl b/src/ForceFreeStates/Integrators.jl
index 8b123b2d0..ac5406782 100644
--- a/src/ForceFreeStates/Integrators.jl
+++ b/src/ForceFreeStates/Integrators.jl
@@ -37,7 +37,7 @@ and never depends on the thread count. Maps onto `integrator = "riccati"`.
 end
 
 """
-    Galerkin(; solver="LU", nx=256, ...)
+    Galerkin(; solver="auto", nx=256, ...)
 
 RDCON outer-region singular Galerkin solver: solves the same Euler-Lagrange system
 variationally on a finite-element grid packed around the rational surfaces, producing Δ′
@@ -46,7 +46,7 @@ matching `gal_*` control key without the prefix.
 
 ## Fields
 
-  - `solver::String` - Banded linear solver, `"LU"` (zgbtrf/zgbtrs) or `"cholesky"` (zpbtrf/zpbtrs). Matching requires `"LU"`.
+  - `solver::String` - Banded linear solver: `"auto"` (Cholesky, zpbtrf/zpbtrs, when the matrix is Hermitian positive definite, else LU), `"LU"` (zgbtrf/zgbtrs) or `"cholesky"` (error if not positive definite).
   - `nx::Int` - Elements per interval between singular surfaces.
   - `nq::Int` - Gauss-Lobatto quadrature order per element.
   - `pfac::Float64` - Grid packing ratio near singular surfaces.
@@ -63,7 +63,7 @@ matching `gal_*` control key without the prefix.
   - `edge_onesided::Bool` - Pack the two end intervals one-sided toward their single rational end instead of the Fortran symmetric pack.
 """
 @kwdef struct Galerkin <: AbstractIntegrator
-    solver::String = "LU"
+    solver::String = "auto"
     nx::Int = 256
     nq::Int = 6
     pfac::Float64 = 0.001
diff --git a/test/runtests_galerkin_cholesky.jl b/test/runtests_galerkin_cholesky.jl
index f9c815b6e..4426185ff 100644
--- a/test/runtests_galerkin_cholesky.jl
+++ b/test/runtests_galerkin_cholesky.jl
@@ -52,6 +52,14 @@ using LinearAlgebra
     nonhermitian = copy(H)
     nonhermitian[6, 4] += 0.1
     @test check(dirichlet(nonhermitian), "auto") == "LU"
+    complexdiag = copy(H)
+    complexdiag[5, 5] += 0.5im
+    @test check(dirichlet(complexdiag), "auto") == "LU"
+    # a small block's asymmetry counts against its own scale, not the largest entry
+    D = Diagonal([j <= n ÷ 2 ? 1e6 : 1e-6 for j in 1:n])
+    graded = D * H * D
+    graded[n-2, n-3] *= 1.01
+    @test FFS.gal_banded_solve!(fullband(dirichlet(graded)), copy(rhs), kl, kl, "auto") == "LU"
 
     # the lower-band factor equals the dense Cholesky factor
     ab = zeros(ComplexF64, kl + 1, n)

From cd290c659b393c5bc1aa62b425b09d66d89dd2bb Mon Sep 17 00:00:00 2001
From: StuartBenjamin <stubenj9@mit.edu>
Date: Tue, 13 Oct 2026 00:46:48 -0700
Subject: [PATCH 6/6] ForceFreeStates.Galerkin - MINOR - Drop the Hermitian
 precondition; Cholesky falls back to LU only when not positive definite

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 docs/src/galerkin.md                          |  4 +-
 src/ForceFreeStates/CoreTypes.jl              |  2 +-
 src/ForceFreeStates/Galerkin/GalerkinSolve.jl | 13 +---
 src/ForceFreeStates/Integrators.jl            |  2 +-
 test/runtests_galerkin_cholesky.jl            | 76 +++++--------------
 5 files changed, 29 insertions(+), 68 deletions(-)

diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md
index c787323a3..344d643e0 100644
--- a/docs/src/galerkin.md
+++ b/docs/src/galerkin.md
@@ -22,8 +22,8 @@ Setting `gal_match_flag` additionally matches the inner layer, giving a driven 
 
 `gal_solver` picks the banded solver. `"auto"` (the default) tries a Hermitian Cholesky
 factorization of the lower band (LAPACK `zpbtrf`/`zpbtrs`, with the Dirichlet rows' columns moved
-to the right-hand side) and falls back to Jacobi-scaled LU when the matrix is not Hermitian
-positive definite. `"LU"` skips the Cholesky attempt;
+to the right-hand side; the diagonal's imaginary part is ignored) and falls back to Jacobi-scaled
+LU when the matrix is not positive definite. `"LU"` skips the Cholesky attempt;
 `"cholesky"` stops with an error instead of falling back.
 
 The implementation lives in `src/ForceFreeStates/Galerkin/`:
diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl
index 5ebbe08d9..9e0831265 100644
--- a/src/ForceFreeStates/CoreTypes.jl
+++ b/src/ForceFreeStates/CoreTypes.jl
@@ -192,7 +192,7 @@ gpec.toml.
     extended_precision_bvp::Bool = true   # Promote Δ' BVP to Complex{Double64}; default on (Float64 drifts the imaginary Δ' by 2–5× on DIIID-class cases).
 
     # --- RDCON outer-region Galerkin Δ′ solver (gal_solve port); selected by integrator = "galerkin" ---
-    gal_solver::String = "auto"     # "auto" (banded Cholesky if the matrix is Hermitian positive definite, else LU), "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs, error if not positive definite)
+    gal_solver::String = "auto"     # "auto" (banded Cholesky if the matrix is positive definite, else LU), "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs, error if not positive definite)
     gal_nx::Int = 256               # elements per interval between singular surfaces
     gal_nq::Int = 6                 # Gauss-Lobatto quadrature order per element
     gal_pfac::Float64 = 0.001       # grid packing ratio near singular surfaces
diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
index 025ad4e8c..6da63d2f1 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
@@ -358,12 +358,12 @@ end
 
 Solve the full-band system `mat` (LAPACK `gbtrf` storage) in place of `sol` and return the
 solver used. `"auto"` tries `gal_cholesky_solve!` and falls back to `gal_scaled_lu_solve!` when
-the matrix is not Hermitian positive definite; `"cholesky"` errors instead of falling back.
+the matrix is not positive definite; `"cholesky"` errors instead of falling back.
 """
 function gal_banded_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{ComplexF64}, kl::Int, ku::Int, solver::String)
     if solver != "LU"
         gal_cholesky_solve!(mat, sol, kl, ku) && return "cholesky"
-        solver == "cholesky" && error("gal_solver=\"cholesky\": the Galerkin matrix is not Hermitian positive definite; use \"auto\" or \"LU\"")
+        solver == "cholesky" && error("gal_solver=\"cholesky\": the Galerkin matrix is not positive definite; use \"auto\" or \"LU\"")
     end
     gal_scaled_lu_solve!(mat, sol, kl, ku)
     return "LU"
@@ -374,9 +374,8 @@ end
 
 Banded Cholesky solve of the full-band system `mat`, leaving `mat` untouched. Rows that are
 identity rows (Dirichlet DOFs: axis, fixed or coil edge) have their columns moved to the
-right-hand side, which keeps the system Hermitian. Returns `false`, with `sol` unchanged, when
-the remaining matrix is not Hermitian (each entry to 1e-10 of √|aᵢᵢaⱼⱼ|, diagonal included) or not
-positive definite.
+right-hand side, which keeps the system Hermitian. Only the lower band and the real part of the
+diagonal are used. Returns `false`, with `sol` unchanged, when the matrix is not positive definite.
 """
 function gal_cholesky_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{ComplexF64}, kl::Int, ku::Int)
     kl == ku || return false
@@ -385,10 +384,6 @@ function gal_cholesky_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{Comp
     a(i, j) = mat[off+i-j, j]
     band(i) = max(1, i - kl):min(n, i + kl)
     fixed = [all(j == i ? a(i, j) == 1 : iszero(a(i, j)) for j in band(i)) for i in 1:n]
-    # Scale-free test (the diagonal spans many decades); 1e-10 clears the quadrature round-off (~1e-11)
-    for j in 1:n, i in j:min(n, j + kl)
-        fixed[i] || fixed[j] || abs(a(i, j) - conj(a(j, i))) <= 1e-10 * sqrt(abs(a(i, i) * a(j, j))) || return false
-    end
     ab = zeros(ComplexF64, kl + 1, n)
     for j in 1:n, i in j:min(n, j + kl)
         ab[1+i-j, j] = fixed[i] || fixed[j] ? (i == j ? one(ComplexF64) : zero(ComplexF64)) : a(i, j)
diff --git a/src/ForceFreeStates/Integrators.jl b/src/ForceFreeStates/Integrators.jl
index ac5406782..05a79f81d 100644
--- a/src/ForceFreeStates/Integrators.jl
+++ b/src/ForceFreeStates/Integrators.jl
@@ -46,7 +46,7 @@ matching `gal_*` control key without the prefix.
 
 ## Fields
 
-  - `solver::String` - Banded linear solver: `"auto"` (Cholesky, zpbtrf/zpbtrs, when the matrix is Hermitian positive definite, else LU), `"LU"` (zgbtrf/zgbtrs) or `"cholesky"` (error if not positive definite).
+  - `solver::String` - Banded linear solver: `"auto"` (Cholesky, zpbtrf/zpbtrs, when the matrix is positive definite, else LU), `"LU"` (zgbtrf/zgbtrs) or `"cholesky"` (error if not positive definite).
   - `nx::Int` - Elements per interval between singular surfaces.
   - `nq::Int` - Gauss-Lobatto quadrature order per element.
   - `pfac::Float64` - Grid packing ratio near singular surfaces.
diff --git a/test/runtests_galerkin_cholesky.jl b/test/runtests_galerkin_cholesky.jl
index 4426185ff..01bc5deee 100644
--- a/test/runtests_galerkin_cholesky.jl
+++ b/test/runtests_galerkin_cholesky.jl
@@ -1,72 +1,38 @@
 using Test
 using LinearAlgebra
 
-# gal_banded_solve! on full-band (gbtrf layout) systems built like the Galerkin matrix: Hermitian
-# with Dirichlet identity rows (row zeroed, column kept) and nonzero prescribed values.
+# gal_banded_solve! on a full-band (gbtrf layout) system built like the Galerkin matrix: Hermitian, with
+# Dirichlet identity rows (row zeroed, column kept) and nonzero prescribed values, as in the rpec coil columns.
 @testset "Galerkin banded solve: Cholesky first, LU fallback" begin
     FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
-
-    n, kl = 24, 3
+    n, kl, fixed = 24, 3, [1, 2, 24]
     H = zeros(ComplexF64, n, n)
     for j in 1:n, i in max(1, j - kl):min(n, j + kl)
         H[i, j] = i == j ? 4 + sin(j) : 0.3 * cis(0.3i + 0.7j)
     end
-    H = (H + H') / 2  # Hermitian, diagonally dominant → positive definite
-    fixed = [1, 2, n]
-
-    # Dirichlet rows as gal_set_boundary! leaves them: identity row, column untouched
-    function dirichlet(A)
+    H = (H + H') / 2  # diagonally dominant → positive definite
+    function system(A)
         A = copy(A)
-        for i in fixed
-            A[i, :] .= 0
-            A[i, i] = 1
-        end
-        return A
-    end
-    function fullband(A)
+        A[fixed, :] .= 0
+        A[fixed, fixed] = I(length(fixed))
         ab = zeros(ComplexF64, 3kl + 1, n)
         for j in 1:n, i in max(1, j - kl):min(n, j + kl)
             ab[2kl+1+i-j, j] = A[i, j]
         end
-        return ab
+        return A, ab
     end
     rhs = ComplexF64[cis(0.1j) * (1 + j / n) for j in 1:n, _ in 1:2]
-    rhs[fixed, 2] .= [1.0, -0.5im, 2.0]  # prescribed values, as in the rpec coil columns
-
-    function check(A, solver)
-        sol = copy(rhs)
-        used = FFS.gal_banded_solve!(fullband(A), sol, kl, kl, solver)
-        @test sol ≈ A \ rhs rtol = 1e-12
-        return used
-    end
-
-    @test check(dirichlet(H), "auto") == "cholesky"
-    @test check(dirichlet(H), "cholesky") == "cholesky"
-    @test check(dirichlet(H), "LU") == "LU"
-
-    indefinite = copy(H)
-    indefinite[5, 5] = -1
-    @test check(dirichlet(indefinite), "auto") == "LU"
-    @test_throws ErrorException FFS.gal_banded_solve!(fullband(dirichlet(indefinite)), copy(rhs), kl, kl, "cholesky")
-
-    nonhermitian = copy(H)
-    nonhermitian[6, 4] += 0.1
-    @test check(dirichlet(nonhermitian), "auto") == "LU"
-    complexdiag = copy(H)
-    complexdiag[5, 5] += 0.5im
-    @test check(dirichlet(complexdiag), "auto") == "LU"
-    # a small block's asymmetry counts against its own scale, not the largest entry
-    D = Diagonal([j <= n ÷ 2 ? 1e6 : 1e-6 for j in 1:n])
-    graded = D * H * D
-    graded[n-2, n-3] *= 1.01
-    @test FFS.gal_banded_solve!(fullband(dirichlet(graded)), copy(rhs), kl, kl, "auto") == "LU"
-
-    # the lower-band factor equals the dense Cholesky factor
-    ab = zeros(ComplexF64, kl + 1, n)
-    for j in 1:n, i in j:min(n, j + kl)
-        ab[1+i-j, j] = H[i, j]
-    end
-    @test FFS.gal_zpbtrf!(ab, kl)
-    L = cholesky(Hermitian(H, :L)).L
-    @test all(ab[1+i-j, j] ≈ L[i, j] for j in 1:n for i in j:min(n, j + kl))
+    rhs[fixed, 2] .= [1.0, -0.5im, 2.0]
+
+    A, ab = system(H)
+    sol = copy(rhs)
+    @test FFS.gal_banded_solve!(ab, sol, kl, kl, "auto") == "cholesky"
+    @test sol ≈ A \ rhs rtol = 1e-12
+
+    H[5, 5] = -1
+    A, ab = system(H)
+    sol = copy(rhs)
+    @test FFS.gal_banded_solve!(copy(ab), sol, kl, kl, "auto") == "LU"
+    @test sol ≈ A \ rhs rtol = 1e-12
+    @test_throws ErrorException FFS.gal_banded_solve!(ab, copy(rhs), kl, kl, "cholesky")
 end