From 965b52b882b9fb5e924b695653f20351a2913962 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sun, 5 Jul 2026 17:12:45 -0400 Subject: [PATCH] exp: Pade 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 new GESVFactorization (v4.2) 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 the 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 also reproduces the gesv! results bit-for-bit, keeps the same per-call allocation bytes, and restores the historical SingularException throw on a singular denominator (the default algorithm's QR fallback resolved it silently instead). Requires LinearSolve 4.2 for GESVFactorization. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 --- Project.toml | 2 +- src/exp_baseexp.jl | 23 +++++++++++------------ 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/Project.toml b/Project.toml index 37b0bfe..a5aa5b6 100644 --- a/Project.toml +++ b/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.2" PrecompileTools = "1" Printf = "1.10" Random = "1" diff --git a/src/exp_baseexp.jl b/src/exp_baseexp.jl index 889247e..e2fa045 100644 --- a/src/exp_baseexp.jl +++ b/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 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