diff --git a/Project.toml b/Project.toml
index 2110b26c1..e2fbe83cc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -58,7 +58,7 @@ ForwardDiffChainRules = "0.3.0"
 GraphRecipes = "0.5.13"
 Integrals = "4.5.0, 5.0.0"
 LinearAlgebra = "1.10.0"
-LinearSolve = "3.19.0"
+LinearSolve = "3.83.0"
 MatterPower = "0.5.0"
 ModelingToolkit = "11.0.0"
 ModelingToolkitBase = "1.0.0"
diff --git a/docs/src/benchmarks.md b/docs/src/benchmarks.md
index d0ed2a8f9..e18d80006 100644
--- a/docs/src/benchmarks.md
+++ b/docs/src/benchmarks.md
@@ -106,7 +106,7 @@ The points on each curve correspond to a sequence of tolerances.
 ```@example bench
 # TODO: test different nlsolve # hide
 # TODO: add AdaptiveRadau/RadauIIA5 when they support sparse J: https://github.com/SciML/OrdinaryDiffEq.jl/issues/2892 # hide
-linsolve = KLUFactorization()
+linsolve = PureKLUFactorization()
 ptalgs = [algtype(; linsolve) for algtype in [TRBDF2, KenCarp4, KenCarp47, KenCarp5, Kvaerno5, Rodas4P, Rodas5P, Rodas6P, QNDF, FBDF]]
 ptprobgen = SymBoltz.setuppt(prob.pt, bgsol)
 refalg = Rodas5P(; linsolve)
@@ -185,15 +185,15 @@ prob_jac = prob # CosmologyProblem(M, pars; jac = true, sparse = true)
 prob_nojac = CosmologyProblem(M, pars; jac = false, sparse = true)
 
 bgopts = (alg = Rodas5P(linsolve = RFLUFactorization(),),)
-ptopts = (alg = Rodas5P(linsolve = KLUFactorization(),), save_everystep = false) # generate function for J symbolically
+ptopts = (alg = Rodas5P(linsolve = PureKLUFactorization(),), save_everystep = false) # generate function for J symbolically
 bench["symbolic"] = @benchmarkable $solve($prob_jac, $ks; bgopts = $bgopts, ptopts = $ptopts) samples=5 seconds=30
 
 bgopts = (alg = Rodas5P(linsolve = RFLUFactorization(), autodiff = SymBoltz.AutoForwardDiff()),)
-ptopts = (alg = Rodas5P(linsolve = KLUFactorization(), autodiff = SymBoltz.AutoForwardDiff()), save_everystep = false) # compute J with forward-mode AD
+ptopts = (alg = Rodas5P(linsolve = PureKLUFactorization(), autodiff = SymBoltz.AutoForwardDiff()), save_everystep = false) # compute J with forward-mode AD
 bench["forward diff"] = @benchmarkable $solve($prob_nojac, $ks; bgopts = $bgopts, ptopts = $ptopts) samples=5 seconds=30
 
 bgopts = (alg = Rodas5P(linsolve = RFLUFactorization(), autodiff = SymBoltz.AutoForwardDiff()),) # fails with finite diff background J
-ptopts = (alg = Rodas5P(linsolve = KLUFactorization(), autodiff = SymBoltz.AutoFiniteDiff()), save_everystep = false) # compute J with finite differences
+ptopts = (alg = Rodas5P(linsolve = PureKLUFactorization(), autodiff = SymBoltz.AutoFiniteDiff()), save_everystep = false) # compute J with finite differences
 bench["finite diff"] = @benchmarkable $solve($prob_nojac, $ks; bgopts = $bgopts, ptopts = $ptopts) samples=5 seconds=30
 
 results = run(bench; verbose = true)
@@ -240,12 +240,12 @@ This plot compares the time to solve several perturbation $k$-modes with differe
 ks = 10 .^ range(-2, 4, length=75)
 ptopts1 = (alg = Rodas5P(linsolve = LUFactorization()), save_everystep = false)
 ptopts2 = (alg = Rodas5P(linsolve = RFLUFactorization()), save_everystep = false)
-ptopts3 = (alg = Rodas5P(linsolve = KLUFactorization()), save_everystep = false)
+ptopts3 = (alg = Rodas5P(linsolve = PureKLUFactorization()), save_everystep = false)
 ptopts4 = (alg = Rodas5P(linsolve = UMFPACKFactorization()), save_everystep = false)
-ts1 = [minimum(@elapsed solve(prob, ks; ptopts = ptopts1, verbose = true) for i in 1:3) for prob in probs_dense]
-ts2 = [minimum(@elapsed solve(prob, ks; ptopts = ptopts2, verbose = true) for i in 1:3) for prob in probs_dense]
-ts3 = [minimum(@elapsed solve(prob, ks; ptopts = ptopts3, verbose = true) for i in 1:3) for prob in probs_sparse]
-ts4 = [minimum(@elapsed solve(prob, ks; ptopts = ptopts4, verbose = true) for i in 1:3) for prob in probs_sparse]
+ts1 = [minimum(@elapsed solve(prob, ks; ptopts = ptopts1) for i in 1:3) for prob in probs_dense]
+ts2 = [minimum(@elapsed solve(prob, ks; ptopts = ptopts2) for i in 1:3) for prob in probs_dense]
+ts3 = [minimum(@elapsed solve(prob, ks; ptopts = ptopts3) for i in 1:3) for prob in probs_sparse]
+ts4 = [minimum(@elapsed solve(prob, ks; ptopts = ptopts4) for i in 1:3) for prob in probs_sparse]
 
 p1 = plot(ylabel = "time / s", xticks = (lmaxs, ""), ylims = (0.0, ceil(max(maximum(ts1), maximum(ts2), maximum(ts3), maximum(ts4)))))
 marker = :circle
@@ -273,7 +273,7 @@ plot!(p2, lmaxs, speedups4; marker, label = nothing)
 plot(p1, p2; size = (800, 600), layout = grid(2, 1, heights=(3//4, 1//4)))
 ```
 
-Except for models with a very small perturbation system, it is a good idea to generate the sparse Jacobian and use the sparse `KLUFactorization` linear solver.
+Except for models with a very small perturbation system, it is a good idea to generate the sparse Jacobian and use the sparse `PureKLUFactorization` linear solver.
 
 ```@setup
 # TODO: tune Krylov with verbose = 1, ILU, ..., atol, rtol # hide
diff --git a/src/solve.jl b/src/solve.jl
index 54c61ddc6..12e55f1c8 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -312,7 +312,7 @@ issparse(prob::ODEProblem) = issparse(prob.f.jac_prototype)
 
 function bgalg(prob::ODEProblem; stiff = true)
     if issparse(prob)
-        linsolve = KLUFactorization()
+        linsolve = PureKLUFactorization()
     else
         linsolve = RFLUFactorization(throwerror = true)
     end
@@ -326,7 +326,7 @@ bgalg(prob::CosmologyProblem; kwargs...) = bgalg(prob.bg; kwargs...)
 
 function ptalg(prob::ODEProblem; accuracy = 2)
     if issparse(prob)
-        linsolve = KLUFactorization()
+        linsolve = PureKLUFactorization()
     else
         linsolve = RFLUFactorization(throwerror = true)
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index 84a785823..a1f9f1c9d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -593,14 +593,14 @@ end
 
     # with ΛCDM model
     k = [1e0, 1e1, 1e2, 1e3]
-    sol = solve(prob_sparse, k; bgopts = (alg = SymBoltz.Rodas4P(linsolve = SymBoltz.LUFactorization()),), ptopts = (alg = SymBoltz.KenCarp4(linsolve = SymBoltz.KLUFactorization()),))
+    sol = solve(prob_sparse, k; bgopts = (alg = SymBoltz.Rodas4P(linsolve = SymBoltz.LUFactorization()),), ptopts = (alg = SymBoltz.KenCarp4(linsolve = SymBoltz.PureKLUFactorization()),))
     @test issuccess(sol)
 
     M2 = RMΛ()
     pars2 = Dict(M2.m.Ω₀ => 0.3, M2.r.Ω₀ => 1e-5, M2.g.h => NaN, M2.r.T₀ => NaN)
     prob2 = CosmologyProblem(M2, pars2; jac = true, sparse = true, bgopts = (sparse = true,)) # demand sparse background
-    bgopts = (alg = SymBoltz.Rodas4P(linsolve = SymBoltz.KLUFactorization()),)
-    ptopts = (alg = SymBoltz.KenCarp4(linsolve = SymBoltz.KLUFactorization()),)
+    bgopts = (alg = SymBoltz.Rodas4P(linsolve = SymBoltz.PureKLUFactorization()),)
+    ptopts = (alg = SymBoltz.KenCarp4(linsolve = SymBoltz.PureKLUFactorization()),)
     sol = solve(prob2, k; bgopts, ptopts)
     @test issuccess(sol)
 end
@@ -612,7 +612,7 @@ end
     @test issuccess(solve(prob_sparse, 1.0))
     @test issuccess(solve(prob_dense, 1.0; bgopts = (alg = SymBoltz.Rodas5P(),), ptopts = (alg = SymBoltz.Rodas5P(),))) # should automatically find compatible linsolves
     @test issuccess(solve(prob_sparse, 1.0; bgopts = (alg = SymBoltz.Rodas5P(),), ptopts = (alg = SymBoltz.Rodas5P(),)))
-    @test_throws "dense Jacobian must be solved with dense" solve(prob_dense; bgopts = (alg = SymBoltz.Rodas5P(linsolve = SymBoltz.KLUFactorization()),)) # has dense background
+    @test_throws "dense Jacobian must be solved with dense" solve(prob_dense; bgopts = (alg = SymBoltz.Rodas5P(linsolve = SymBoltz.PureKLUFactorization()),)) # has dense background
     @test_throws "sparse Jacobian must be solved with sparse" solve(prob_sparse, 1.0; ptopts = (alg = SymBoltz.Rodas5P(linsolve = SymBoltz.RFLUFactorization()),)) # has sparse perturbations
     @test issuccess(solve(prob_dense, 1.0; bgopts = (alg = SymBoltz.bgalg(prob_dense),), ptopts = (alg = SymBoltz.ptalg(prob_dense; accuracy = 0),)))
     @test issuccess(solve(prob_dense, 1.0; bgopts = (alg = SymBoltz.bgalg(prob_dense),), ptopts = (alg = SymBoltz.ptalg(prob_dense; accuracy = 1),)))