StuartBenjamin · Oct 7, 2026 · Oct 9, 2026 · Oct 9, 2026 · Oct 9, 2026 · Oct 12, 2026
diff --git 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. `"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; 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/`:
 
 | File | Content |

diff --git 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
@@ -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
@@ -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
@@ -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
@@ -153,7 +153,7 @@ gpec.toml.
   - `force_termination::Bool` - Terminate after force-free states (skip perturbed equilibrium calculations)
   - `integrator::String` - Which formalism integrates the Euler-Lagrange system. `"forward"` sweeps the plasma serially with Gaussian reduction and returns `u_store` / `du_store` / `xi_s_store` dense in the axis (EL) basis — the only convention PerturbedEquilibrium and FieldReconstruction consume correctly, and the only path that supports `kinetic_factor > 0`. `"riccati"` (default) runs the chunked fundamental-matrix propagator driver (Glasser 2018 Phys. Plasmas 25, 032507): chunks are integrated independently from identity initial conditions and assembled serially with Riccati-style crossings, which is the only way to obtain the singular-surface Δ' matrix for the tearing-mode solvers downstream, but leaves `u_store` as sparse chunk-endpoint Riccati states, so dense ξ profiles are unavailable. `"galerkin"` solves the same Euler-Lagrange system variationally instead of by radial ODE integration — the RDCON outer-region singular Galerkin method (Glasser, Wang & Park 2016 Phys. Plasmas 23, 112506), which discretizes the displacement on packed Hermite-cubic elements and solves one global banded system — producing the resistive Δ′ matrix and, when `gal_match_flag` is set, the RPEC inner-layer-matched ξ; it computes its own vacuum response and returns no free-boundary energies, and does not support `kinetic_factor > 0`. Requires `singfac_min != 0` for `"riccati"`.
   - `nchunks::Int` - Target number of Riccati integration chunks. `0` (the default) derives the count from problem structure alone: `max(2·msing + 3, 8·(msing + 1) + msing)`, enough sub-chunks per segment to keep the accumulated propagator products well-conditioned. An explicit value below `2·msing + 3` is clamped up with a warning. Chunk sizing never consults `Threads.nthreads()`, so Riccati outputs are identical whatever thread count `julia -t` provides; threads only change wall-clock.
-  - `extended_precision_bvp::Bool` - When `true` (default), promote the Δ' BVP linear system to `Complex{Double64}` (~31 digits) for the LU solve and PEST3 combination. Guards against catastrophic cancellation in the PEST3 four-term combination (dp_raw entries can be 10⁴–10⁵× larger than the result; the imaginary part of off-diagonal Δ' is particularly sensitive). Disabling (`false`) saves ~1.5–2× the BVP solve time but on DIIID-class equilibria the imaginary Δ' components can drift by factors of 2–5×; only disable for performance experiments on cases where Float64 has been validated against Double64.
+  - `extended_precision_bvp::Bool` - Solve the Δ' BVP and its PEST3 combination in `Complex{Double64}` (~31 digits). Default `false` (Float64, 1.5–5.4× faster on large matrices); set `true` to check Float64 at small `singfac_min`, where the four-term combination loses digits.
 """
 @kwdef struct ForceFreeStatesControl
     verbose::Bool = true
@@ -191,10 +191,10 @@ gpec.toml.
     force_termination::Bool = false
     integrator::String = "riccati"   # Default: unlocks SingularSurfaces/Delta_prime_matrix (STRIDE BVP Δ′ matrix) used by SLAYER/GGJ downstream. Use "forward" for dense ξ (PerturbedEquilibrium) or kinetic runs.
     nchunks::Int = 0                 # Riccati chunk-count target; 0 = auto (derived from msing alone, never from Threads.nthreads()).
-    extended_precision_bvp::Bool = true   # Promote Δ' BVP to Complex{Double64}; default on (Float64 drifts the imaginary Δ' by 2–5× on DIIID-class cases).
+    extended_precision_bvp::Bool = false   # Δ' BVP in Complex{Double64}; opt-in check of Float64 at small singfac_min.
 
     # --- 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 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/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
@@ -769,7 +769,7 @@ function cross_ideal_singular_surf!(
     # Compute direction-specific asymptotic power series for this singular surface
     singp = intr.sing[ising]
     sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=1.0)
-    sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0, alpha_override=sing_asymp_right.alpha)
+    sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0)
     dpsi = singp.psifac - odet.psifac # ψ_res - ψ (positive)
 
     # Get asymptotic coefficients before crossing (left side)

diff --git 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
 #
@@ -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
@@ -469,7 +455,6 @@ function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing,
     if ncoil > 0
         # rpec: identity edge + unit coil sources (gal.f)
         cell.mat[:, :, 4, :] .= 0
-        chol && (cell.mat[:, :, :, 4] .= 0)
         for idx in 1:mpert
             cell.mat[idx, idx, 4, 4] = 1
         end
@@ -482,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