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