exp: Padé denominator solve through GESVFactorization

Benchmarks of the four candidate solve paths on identical per-cell
matrices (n = 5..250, scales 0.3/4.0, best-of-N with preallocated caches,
Julia 1.12, OpenBLAS single-threaded) showed LinearSolve's
GESVFactorization (v4.3) is the only variant at parity with the raw
LAPACK.gesv! this code historically used, in every cell: -1.8% to +2.3%
end-to-end, versus up to +5% for a pinned LUFactorization and up to +4.5%
for the default algorithm choice (whose safety fallback adds an O(n^2)
backup copy per refactorization). GESVFactorization reproduces the gesv!
results bit-for-bit and keeps the same per-call allocation bytes; on a
singular denominator it returns a Failure retcode, which this code turns
back into the historical SingularException (the default algorithm's QR
fallback resolved it silently instead).

Requires LinearSolve 4.3 for GESVFactorization.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

Author: ChrisRackauckas

Project.toml

@@ -30,7 +30,7 @@ ForwardDiff = "0.10.13, 1"
 GPUArraysCore = "0.1, 0.2"
 GenericSchur = "0.5.3"
 LinearAlgebra = "1.10"
-LinearSolve = "4.1"
+LinearSolve = "4.3"
 PrecompileTools = "1"
 Printf = "1.10"
 Random = "1"

src/exp_baseexp.jl

@@ -1,4 +1,4 @@
-using LinearSolve: LinearSolve, LinearProblem
+using LinearSolve: LinearSolve, LinearProblem, GESVFactorization
 
 """
     ExpMethodHigham2005Base()
@@ -16,15 +16,16 @@ function alloc_mem(
     U = Matrix{T}(undef, n, n)
     V = Matrix{T}(undef, n, n)
     temp = Matrix{T}(undef, n, n)
-    # Cached LinearSolve workspace for the Padé denominator solves, using
-    # LinearSolve's default algorithm choice (size-dependent: e.g. generic LU
-    # for tiny matrices, LAPACK/RecursiveFactorization above). The buffers are
-    # workspace-owned, so aliasing lets the factorization refactorize in place
-    # on every `solve!` instead of copying the n×n matrix.
+    # Cached LinearSolve workspace for the Padé denominator solves.
+    # GESVFactorization (LinearSolve 4.3) mirrors the raw LAPACK.gesv! this code historically
+    # used (single combined factorize-and-solve call, bit-identical results)
+    # and benchmarks at parity with it from n = 5 through n = 250. The buffers
+    # are workspace-owned, so aliasing lets the factorization overwrite Abuf in
+    # place on every `solve!` instead of copying the n×n matrix.
     Abuf = Matrix{T}(undef, n, n)
     Bbuf = Matrix{T}(undef, n, n)
     linsolve = LinearSolve.init(
-        LinearProblem(Abuf, Bbuf);
+        LinearProblem(Abuf, Bbuf), GESVFactorization();
         alias = LinearSolve.LinearAliasSpecifier(alias_A = true, alias_b = true)
     )
     return (A2, P, U, V, temp, linsolve)
@@ -37,13 +38,11 @@ function _pade_linsolve!(
     ) where {T}
     Abuf = linsolve.A
     copyto!(Abuf, temp)
-    linsolve.A = Abuf # flag the refactorization; lu! overwrites Abuf
+    linsolve.A = Abuf # flag the refactorization; gesv! overwrites Abuf
     copyto!(linsolve.b, X)
     sol = LinearSolve.solve!(linsolve)
-    # The Pade denominator is nonsingular by construction for the branch norm
-    # bounds, so this only triggers if even the default algorithm's safety
-    # fallback failed; surface it like LAPACK.gesv! did rather than silently
-    # propagating garbage.
+    # LAPACK.gesv! threw on a singular denominator; keep that contract instead
+    # of silently propagating garbage from a failed factorization.
     if !LinearSolve.SciMLBase.successful_retcode(sol.retcode)
         throw(LinearAlgebra.SingularException(0))
     end