From 9873597ea73872f5799a0d837292b28bdcfb0da8 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 4 Jul 2026 22:42:23 -0400 Subject: [PATCH 1/5] exp: Pade denominator solve via LinearSolve default algorithm MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ExpMethodHigham2005Base solved (V - U) X = (V + U) with raw LAPACK.gesv! calls; route both sites (small-norm and squaring branches) through LinearSolve's default algorithm choice with the batched matrix RHS instead. This picks the best factorization per size/type (e.g. RecursiveFactorization for small matrices when loaded), and fails with a retcode rather than a LAPACKException on pathological input. Verified exp ≈ Base.exp for n = 5, 40, 120 at small and large norms (both code paths) against LinearSolve 4.0.0. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh --- src/exp_baseexp.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/exp_baseexp.jl b/src/exp_baseexp.jl index 6069d93..9afa212 100644 --- a/src/exp_baseexp.jl +++ b/src/exp_baseexp.jl @@ -80,7 +80,9 @@ function exponential!( U, temp = temp, U # equivalent to U = A * U @. X = V + U @. temp = V - U - LAPACK.gesv!(temp, X) + # Pade denominator solve through LinearSolve's default algorithm choice + # (matrix RHS, LinearSolve >= 4); replaces the raw LAPACK.gesv! call. + X .= LinearSolve.solve(LinearProblem(temp, X)).u else s = log2(nA / 5.4) # power of 2 later reversed by squaring si = 0 # always defined so the s > 0 squaring loop is type-stable @@ -109,7 +111,9 @@ function exponential!( U, temp = temp, U # equivalent to U = A * U @. X = V + U @. temp = V - U - LAPACK.gesv!(temp, X) + # Pade denominator solve through LinearSolve's default algorithm choice + # (matrix RHS, LinearSolve >= 4); replaces the raw LAPACK.gesv! call. + X .= LinearSolve.solve(LinearProblem(temp, X)).u if s > 0 # squaring to reverse dividing by power of 2 for t in 1:si From a7f54961c1344c1bc83f6215fdc27b652342e29c Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sun, 5 Jul 2026 02:25:28 -0400 Subject: [PATCH 2/5] exp: cache the Pade denominator LinearSolve workspace in alloc_mem The per-call solve(LinearProblem(temp, X)) added init, matrix copies, and a fresh factorization allocation on every exponential! call, which measured slower than the raw LAPACK.gesv! it replaced for n <= 60. Build a LinearSolve.init workspace once in alloc_mem with workspace-owned n-by-n buffers and alias = LinearAliasSpecifier( alias_A = true, alias_b = true), so each call just refills the buffers and solve! refactorizes in place (lu!) with no O(n^2) copy (LinearSolve >= 4.1). Both Pade branches (small-norm and scaling-and-squaring) route through the shared helper, which restores gesv!'s SingularException contract on a failed factorization. Best-of-N per-call times (Julia 1.12, OpenBLAS, this machine), old gesv! vs cached LinearSolve: ~1.3-1.8x slower for n <= 30 (fixed solve!/property overhead of a few microseconds), parity at n = 60, and 1.2-1.5x faster for n >= 120. Allocations on cache reuse stay O(n) (e.g. 1.4 KB at n = 60). Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 --- Project.toml | 2 +- src/exp_baseexp.jl | 42 ++++++++++++++++++++++++++++++++++-------- 2 files changed, 35 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 881cb95..37b0bfe 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" +LinearSolve = "4.1" PrecompileTools = "1" Printf = "1.10" Random = "1" diff --git a/src/exp_baseexp.jl b/src/exp_baseexp.jl index 9afa212..45e9331 100644 --- a/src/exp_baseexp.jl +++ b/src/exp_baseexp.jl @@ -1,3 +1,5 @@ +using LinearSolve: LinearSolve, LinearProblem, LUFactorization + """ ExpMethodHigham2005Base() @@ -14,7 +16,35 @@ function alloc_mem( U = Matrix{T}(undef, n, n) V = Matrix{T}(undef, n, n) temp = Matrix{T}(undef, n, n) - return (A2, P, U, V, temp) + # Cached LinearSolve workspace for the Padé denominator solves. The buffers + # are workspace-owned, so aliasing lets LUFactorization refactorize in + # place (`lu!`) 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), LUFactorization(); + alias = LinearSolve.LinearAliasSpecifier(alias_A = true, alias_b = true) + ) + return (A2, P, U, V, temp, linsolve) +end + +# X .= temp \ X through the cached LinearSolve workspace: refill the +# workspace-owned buffers, mark A as replaced, and solve in place. +function _pade_linsolve!( + X::StridedMatrix{T}, temp::StridedMatrix{T}, linsolve + ) where {T} + Abuf = linsolve.A + copyto!(Abuf, temp) + linsolve.A = Abuf # flag the refactorization; lu! overwrites Abuf + copyto!(linsolve.b, X) + sol = LinearSolve.solve!(linsolve) + # 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 + copyto!(X, sol.u) + return X end ## Destructive matrix exponential using algorithm from Higham, 2008, @@ -32,7 +62,7 @@ function exponential!( # return copytri!(parent(exp(Hermitian(A))), 'U', true) # end - A2, P, U, V, temp = cache + A2, P, U, V, temp, linsolve = cache fill!(P, zero(T)) fill!(@diagview(P), one(T)) # P = Inn @@ -80,9 +110,7 @@ function exponential!( U, temp = temp, U # equivalent to U = A * U @. X = V + U @. temp = V - U - # Pade denominator solve through LinearSolve's default algorithm choice - # (matrix RHS, LinearSolve >= 4); replaces the raw LAPACK.gesv! call. - X .= LinearSolve.solve(LinearProblem(temp, X)).u + _pade_linsolve!(X, temp, linsolve) else s = log2(nA / 5.4) # power of 2 later reversed by squaring si = 0 # always defined so the s > 0 squaring loop is type-stable @@ -111,9 +139,7 @@ function exponential!( U, temp = temp, U # equivalent to U = A * U @. X = V + U @. temp = V - U - # Pade denominator solve through LinearSolve's default algorithm choice - # (matrix RHS, LinearSolve >= 4); replaces the raw LAPACK.gesv! call. - X .= LinearSolve.solve(LinearProblem(temp, X)).u + _pade_linsolve!(X, temp, linsolve) if s > 0 # squaring to reverse dividing by power of 2 for t in 1:si From 719ea4cc7a50579d33523446459e3bb283c2ecdf Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sun, 5 Jul 2026 05:22:38 -0400 Subject: [PATCH 3/5] Retrigger CI: LinearSolve 4.1.0 registered Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 From 849d3069bb1c37d72953f5b8bd45de1e79856a98 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sun, 5 Jul 2026 05:22:38 -0400 Subject: [PATCH 4/5] Retrigger CI: LinearSolve 4.1.0 registered Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 From 70b6b93484f37aa73844c07a4e73fa535c3ce8e5 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sun, 5 Jul 2026 05:23:29 -0400 Subject: [PATCH 5/5] exp: Pade denominator workspace uses LinearSolve's default algorithm Drop the pinned LUFactorization in alloc_mem and let LinearSolve.init size-select the algorithm (GenericLUFactorization for n <= 10, tuned LAPACK/RecursiveFactorization choices above). The tiny-matrix generic LU refactorizes in place with a cached pivot vector, making the warm solve! allocation-free and at parity with (n = 5) or faster than (n = 10) the raw LAPACK.gesv! this path originally used. Best-of-N end-to-end exponential! times against the old gesv! code (Julia 1.12, OpenBLAS single-threaded, this machine) are within ~3% at all of n = 5, 10, 30, 60, 120, 250 for scales 0.3 and 4.0, with fewer per-call allocation bytes at n <= 10. One behavioral note: an exactly singular Pade denominator is now resolved by the default algorithm's QR safety fallback instead of throwing, matching generic \ semantics. The denominator is nonsingular by construction for the branch norm bounds, so the retcode guard only fires if the fallback itself fails. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 --- src/exp_baseexp.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/exp_baseexp.jl b/src/exp_baseexp.jl index 45e9331..889247e 100644 --- a/src/exp_baseexp.jl +++ b/src/exp_baseexp.jl @@ -1,4 +1,4 @@ -using LinearSolve: LinearSolve, LinearProblem, LUFactorization +using LinearSolve: LinearSolve, LinearProblem """ ExpMethodHigham2005Base() @@ -16,13 +16,15 @@ 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. The buffers - # are workspace-owned, so aliasing lets LUFactorization refactorize in - # place (`lu!`) on every `solve!` instead of copying the n×n matrix. + # 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. Abuf = Matrix{T}(undef, n, n) Bbuf = Matrix{T}(undef, n, n) linsolve = LinearSolve.init( - LinearProblem(Abuf, Bbuf), LUFactorization(); + LinearProblem(Abuf, Bbuf); alias = LinearSolve.LinearAliasSpecifier(alias_A = true, alias_b = true) ) return (A2, P, U, V, temp, linsolve) @@ -38,8 +40,10 @@ function _pade_linsolve!( linsolve.A = Abuf # flag the refactorization; lu! overwrites Abuf copyto!(linsolve.b, X) sol = LinearSolve.solve!(linsolve) - # LAPACK.gesv! threw on a singular denominator; keep that contract instead - # of silently propagating garbage from a failed factorization. + # 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. if !LinearSolve.SciMLBase.successful_retcode(sol.retcode) throw(LinearAlgebra.SingularException(0)) end