From c13371651c1bb5f39d9cb7856164c870424054d2 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sun, 5 Jul 2026 17:12:45 -0400 Subject: [PATCH] =?UTF-8?q?exp:=20Pad=C3=A9=20denominator=20solve=20via=20?= =?UTF-8?q?LinearSolve=20default=20algorithm?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Route ExpMethodHigham2005Base's Padé denominator solve (V - U) X = (V + U) through a cached LinearSolve workspace using the *default* algorithm choice with alias_A/alias_b, replacing the raw LAPACK.gesv! calls. The default size-selects the fastest available LU — notably RecursiveFactorization when loaded (the common case in SciML/exponential- integrator stacks), whose solve is ~2-2.7x faster than LAPACK LU at small/medium n. alias_A lets it refactorize in place (no O(n^2) copy); per-call allocation stays small and O(1)-ish. On a singular denominator the solve returns a Failure retcode, which this code turns back into the SingularException exp! has always thrown. An earlier revision pinned GESVFactorization for raw-gesv! parity; with alias_A the default is allocation-competitive and faster (RFLU), needs no 4.3-only feature, and auto-benefits from future LinearSolve improvements, so it wins here. Requires LinearSolve 4.1 (alias_A in-place refactorization). Verified vs Base.exp to ~2e-15; exp-methods 40/40 and Phi 284/284 on LinearSolve 4.1. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh --- src/exp_baseexp.jl | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/src/exp_baseexp.jl b/src/exp_baseexp.jl index 889247e..34d9745 100644 --- a/src/exp_baseexp.jl +++ b/src/exp_baseexp.jl @@ -16,11 +16,12 @@ 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, using the + # default algorithm choice: it size-selects the fastest available LU + # (RecursiveFactorization when loaded — markedly faster than LAPACK at + # small/medium n). The buffers are workspace-owned, so aliasing lets the + # factorization overwrite Abuf in place on every `solve!` (no n×n copy), + # keeping each call allocation-free. Abuf = Matrix{T}(undef, n, n) Bbuf = Matrix{T}(undef, n, n) linsolve = LinearSolve.init( @@ -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 refactorization; alias_A lets it overwrite 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. + # exp! historically threw on a singular denominator (via LAPACK.gesv!); + # keep that contract rather than propagating a failed factorization. if !LinearSolve.SciMLBase.successful_retcode(sol.retcode) throw(LinearAlgebra.SingularException(0)) end