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 <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: ChrisRackauckas

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