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 6069d93..889247e 100644
--- a/src/exp_baseexp.jl
+++ b/src/exp_baseexp.jl
@@ -1,3 +1,5 @@
+using LinearSolve: LinearSolve, LinearProblem
+
 """
     ExpMethodHigham2005Base()
 
@@ -14,7 +16,39 @@ 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, 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);
+        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)
+    # 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
+    copyto!(X, sol.u)
+    return X
 end
 
 ## Destructive matrix exponential using algorithm from Higham, 2008,
@@ -32,7 +66,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,7 +114,7 @@ function exponential!(
         U, temp = temp, U # equivalent to U = A * U
         @. X = V + U
         @. temp = V - U
-        LAPACK.gesv!(temp, X)
+        _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
@@ -109,7 +143,7 @@ function exponential!(
         U, temp = temp, U # equivalent to U = A * U
         @. X = V + U
         @. temp = V - U
-        LAPACK.gesv!(temp, X)
+        _pade_linsolve!(X, temp, linsolve)
 
         if s > 0            # squaring to reverse dividing by power of 2
             for t in 1:si