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

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 |

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)

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)

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)

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)

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 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

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

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,14 +98,10 @@ 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.
-    # 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]
@@ -131,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
@@ -144,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"
-        LinearAlgebra.LAPACK.pbtrf!('L', ws.kl, ws.mat)
-        LinearAlgebra.LAPACK.pbtrs!('L', ws.kl, ws.mat, 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)
@@ -361,3 +352,71 @@ function annotate_galerkin!(out_h5)
     Utilities.HDF5Annotations.annotate!(out_h5, GALERKIN_H5_ANNOTATIONS)
     return nothing
 end
+
+"""
+    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 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 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. 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
+    n = size(mat, 2)
+    off = 2kl + 1
+    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]
+    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
+    for c in 1:n
+        fixed[c] || continue
+        for r in band(c)
+            r == c || fixed[r] || (@views sol[r, :] .-= a(r, c) .* sol[c, :])
+        end
+    end
+    gal_zpbtrs!(ab, kl, sol)
+    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)
+    LAPACK.chklapackerror(min(info[], 0))
+    return info[] == 0
+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

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 `pbtrf!`/`pbtrs!`, `'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

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 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

test/runtests.jl

@@ -33,6 +33,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")

test/runtests_galerkin_cholesky.jl

@@ -0,0 +1,38 @@
+using Test
+using LinearAlgebra
+
+# 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, 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  # diagonally dominant → positive definite
+    function system(A)
+        A = copy(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 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]
+
+    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