matt-pharr · Oct 27, 2026 · Oct 7, 2026 · Oct 9, 2026 · Oct 9, 2026 · Oct 9, 2026
diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md
@@ -7,6 +7,11 @@ a single global banded system.  Cells adjacent to each rational surface ("resona
 "extension" cells) carry the large/small power-law asymptotic solutions analytically, so the
 singular behavior is built into the basis rather than resolved numerically.
 
+The ideal coefficient is reconstructed as ``F = (QL)(QL)^\dagger`` from the interpolated
+Cholesky factor ``L``, with ``Q = \operatorname{diag}(m-nq)``. This uses the same coefficient
+representation as the Euler–Lagrange integration and singular-surface expansions, including
+between spline knots, and preserves Hermitian positive semidefiniteness of ``F``.
+
 The solve produces the inter-surface Δ′ matrix and the PEST-3 matching blocks
 (``A'``, ``B'``, ``\Gamma'``, ``\Delta'``), written to the HDF5 output under the same
 `SingularSurfaces/` paths the Riccati Δ′ uses (`Delta_prime_matrix`, `Delta_prime_raw`,
@@ -21,6 +26,15 @@ Setting `gal_match_flag` additionally matches the inner layer, giving a driven
 composite inner-region profiles are a scripting diagnostic: they need `Galerkin(; cut_solution=true)`
 followed by `solve(MatchProblem(...))`, and a deck run does not write them.
 
+`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 free block's off-diagonal entries must be Hermitian within
+`1e-10 * sqrt(abs(A_ii)) * sqrt(abs(A_jj))`; the diagonal's imaginary part is deliberately
+discarded. Unequal bandwidths, off-diagonal asymmetry, or failure of positive definiteness
+trigger a reported fallback to Jacobi-scaled LU. `"LU"` skips the Cholesky attempt;
+`"cholesky"` stops with an error instead of falling back. The factorization actually used is
+recorded as `ForceFreeStates/Solutions/GalerkinIntegration/solver` in `gpec.h5`.
+
 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
@@ -67,7 +67,7 @@ rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rati
 
 # 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
@@ -68,7 +68,7 @@ rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rati
 
 # 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
@@ -67,7 +67,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
@@ -68,7 +68,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
@@ -148,7 +148,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 `true`, guarding the cancellation in the PEST3 combination. Set `false` to use Float64 after checking the complex Δ' matrix against Double64 at the chosen matching distance.
 """
 @kwdef struct ForceFreeStatesControl
     verbose::Bool = true
@@ -181,10 +181,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 = true
 
     # --- 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/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl
@@ -22,8 +22,7 @@ the Galerkin solver's variant of the same matrix.
   - `F_spline_lower` - reduced F̄ as its lower-triangular Cholesky factor L, F̄ = L·Lᴴ. The ideal
     Euler-Lagrange kernel only ever solves against F̄, so the factor is stored instead of F̄ itself.
   - `F_spline_prim` - primitive F before the Schur complement (for kinetic).
-  - `F_spline_gal` - reduced Hermitian F̄, un-factored, for the Galerkin solver, which needs F̄
-    directly rather than through a solve.
+  - `F_spline_gal` - independently interpolated reduced Hermitian F̄, retained for coefficient diagnostics.
   - `K_spline`, `G_spline` - reduced K̄ and Ḡ of the Euler-Lagrange system
     d/dψ(F̄·ξ′ + K̄·ξ) = K̄†·ξ′ + Ḡ·ξ (Eqs A5-A7).
   - `J_spline` - Jacobian Fourier band (2·mpert−1 conjugate-symmetric coefficients per surface), used

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
 #
@@ -63,7 +63,8 @@ end
 
 Evaluate the `mpert×mpert` matrices `F = Q F̄ Qᴴ`, `K = Q K̄`, `G = Ḡ` at flux `x` with safety factor
 `q`, where `Q = diag(singfac)` and `singfac = m - n q` (direct). Port of `gal_get_fkg` (gal.f).
-Uses the un-factored reduced `mats.ideal.F_spline_gal` (F̄), `mats.ideal.K_spline` (K̄), `mats.ideal.G_spline` (Ḡ). F/K/G are the
+Reconstructs F̄ = LLᴴ from `mats.ideal.F_spline_lower`, matching the Euler–Lagrange kernel and singular
+expansions, and uses `mats.ideal.K_spline` (K̄), `mats.ideal.G_spline` (Ḡ). F/K/G are the
 ideal-MHD Euler–Lagrange coefficient matrices of the outer-region weak form (Glasser 2016, PoP 23, 112506).
 """
 function gal_get_fkg(mats::MatrixSplines, intr::ForceFreeStatesInternal, x::Float64, q::Float64;
@@ -75,7 +76,9 @@ function gal_get_fkg(mats::MatrixSplines, intr::ForceFreeStatesInternal, x::Floa
     F = Fbuf === nothing ? Matrix{ComplexF64}(undef, N, N) : Fbuf
     K = Kbuf === nothing ? Matrix{ComplexF64}(undef, N, N) : Kbuf
     G = Gbuf === nothing ? Matrix{ComplexF64}(undef, N, N) : Gbuf
-    mats.ideal.F_spline_gal(vec(F), x; hint=mats._hint)
+    # Reuse G as factor scratch before evaluating Ḡ.
+    mats.ideal.F_spline_lower(vec(G), x; hint=mats._hint)
+    mul!(F, LowerTriangular(G), adjoint(LowerTriangular(G)))
     mats.ideal.K_spline(vec(K), x; hint=mats._hint)
     mats.ideal.G_spline(vec(G), x; hint=mats._hint)
 
@@ -119,11 +122,12 @@ function gal_gauss_quad!(cell::GalCell, mats::MatrixSplines, profiles, intr::For
         pb = swap_edge ? (pbt[1], pbt[2], pbt[4], pbt[3]) : pbt
         qb = swap_edge ? (qbt[1], qbt[2], qbt[4], qbt[3]) : qbt
         for ip in 0:np, ipert in 1:N, jp in 0:np, jpert in 1:N
-            cell.mat[ipert, jpert, ip+1, jp+1] += w * (
-                F[ipert, jpert] * qb[ip+1] * qb[jp+1] +
-                K[ipert, jpert] * qb[ip+1] * pb[jp+1] +
-                conj(K[jpert, ipert]) * pb[ip+1] * qb[jp+1] +
-                G[ipert, jpert] * pb[ip+1] * pb[jp+1])
+            cell.mat[ipert, jpert, ip+1, jp+1] +=
+                w * (
+                    F[ipert, jpert] * qb[ip+1] * qb[jp+1] +
+                    K[ipert, jpert] * qb[ip+1] * pb[jp+1] +
+                    conj(K[jpert, ipert]) * pb[ip+1] * qb[jp+1] +
+                    G[ipert, jpert] * pb[ip+1] * pb[jp+1])
         end
     end
     return cell
@@ -345,8 +349,8 @@ function gal_resonant!(cell::GalCell, ising::Int, mats::MatrixSplines, profiles,
     if verbose
         @info "  resonant jsing=$jsing side=$(cell.extra) qerr=$(qerr) res1=$(raw[1]) res2=$(raw[2])"
     end
-    hbig = reshape(@view(raw[3:2+N*(np+1)]), N, np + 1)
-    hsmall = reshape(@view(raw[3+N*(np+1):end]), N, np + 1)
+    hbig = reshape(@view(raw[3:(2+N*(np+1))]), N, np + 1)
+    hsmall = reshape(@view(raw[(3+N*(np+1)):end]), N, np + 1)
 
     # Apply the right-side negation and the gal_make_arrays signs (gal.f gal_make_arrays)
     cell.erhs = -sgn_u * raw[1]
@@ -359,15 +363,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 +379,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
@@ -449,17 +441,15 @@ including its full three-way edge branch (rpec count `ncoil = ws.nsol - 2*msing`
   - free boundary (`wv_edge !== nothing`): add the vacuum block `wvac·psio²` at the edge value DOF
     (gal.f).
   - fixed boundary (else): identity edge (gal.f).
-`wv_edge` is `nothing` for the rpec and fixed cases.
+    `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,20 +459,18 @@ 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
         for ipert in 1:mpert
-            ws.rhs[cell.map[ipert, 4], 2 * msing + ipert] = 1
+            ws.rhs[cell.map[ipert, 4], 2*msing+ipert] = 1
         end
     elseif wv_edge !== nothing
         # free boundary: vacuum block wvac·psio² (gal.f)
         @views cell.mat[:, :, 4, 4] .+= wv_edge
     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