From 1e1d56941ecfd9293311c3d8cda83bf72688afb4 Mon Sep 17 00:00:00 2001
From: oscarddssmith <oscar.smith@juliacomputing.com>
Date: Wed, 7 Aug 2024 17:40:08 -0400
Subject: [PATCH] simplify dolinsolve

---
 lib/OrdinaryDiffEqBDF/src/algorithms.jl       |  91 +++----
 .../src/OrdinaryDiffEqCore.jl                 |   2 +-
 .../src/OrdinaryDiffEqDifferentiation.jl      |   2 +-
 .../src/linsolve_utils.jl                     |  47 ++--
 .../src/algorithms.jl                         |  41 ++-
 .../src/extrapolation_caches.jl               |  32 +--
 lib/OrdinaryDiffEqFIRK/src/algorithms.jl      |  24 +-
 lib/OrdinaryDiffEqFIRK/src/firk_caches.jl     |  24 +-
 .../src/algorithms.jl                         |  17 +-
 lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl |  12 +-
 lib/OrdinaryDiffEqPDIRK/src/algorithms.jl     |   9 +-
 .../src/algorithms.jl                         |  30 +--
 .../src/generic_rosenbrock.jl                 |  30 +--
 .../src/rosenbrock_caches.jl                  |  97 ++-----
 .../src/rosenbrock_perform_step.jl            |  77 +-----
 lib/OrdinaryDiffEqSDIRK/src/algorithms.jl     | 251 ++++++++----------
 .../src/algorithms.jl                         |   9 +-
 test/interface/default_solver_tests.jl        |   4 +-
 test/interface/linear_nonlinear_tests.jl      |  69 ++---
 test/interface/linear_solver_test.jl          |   4 +-
 test/interface/preconditioners.jl             |  79 +++---
 21 files changed, 348 insertions(+), 603 deletions(-)

diff --git a/lib/OrdinaryDiffEqBDF/src/algorithms.jl b/lib/OrdinaryDiffEqBDF/src/algorithms.jl
index 1a4c960c7b..73327af5d4 100644
--- a/lib/OrdinaryDiffEqBDF/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqBDF/src/algorithms.jl
@@ -6,11 +6,10 @@ an Adaptive BDF2 Formula and Comparison with The MATLAB Ode15s. Procedia Compute
 ABDF2: Multistep Method
 An adaptive order 2 L-stable fixed leading coefficient multistep BDF method.
 """
-struct ABDF2{CS, AD, F, F2, P, FDT, ST, CJ, K, T, StepLimiter} <:
+struct ABDF2{CS, AD, F, F2, FDT, ST, CJ, K, T, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     κ::K
     tol::T
     smooth_est::Bool
@@ -20,14 +19,14 @@ struct ABDF2{CS, AD, F, F2, P, FDT, ST, CJ, K, T, StepLimiter} <:
 end
 function ABDF2(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        κ = nothing, tol = nothing, linsolve = nothing, precs = DEFAULT_PRECS,
+        κ = nothing, tol = nothing, linsolve = nothing,
         nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :Standard, step_limiter! = trivial_limiter!)
     ABDF2{
         _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
-        typeof(κ), typeof(tol), typeof(step_limiter!)}(linsolve, nlsolve, precs, κ, tol,
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(κ), typeof(tol), typeof(step_limiter!)}(linsolve, nlsolve, κ, tol,
         smooth_est, extrapolant, controller, step_limiter!)
 end
 
@@ -36,11 +35,10 @@ Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. Implicit-Explicit Methods fo
 Dependent Partial Differential Equations. 1995 Society for Industrial and Applied Mathematics
 Journal on Numerical Analysis, 32(3), pp 797-823, 1995. doi: https://doi.org/10.1137/0732037
 """
-struct SBDF{CS, AD, F, F2, P, FDT, ST, CJ, K, T} <:
+struct SBDF{CS, AD, F, F2, FDT, ST, CJ, K, T} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     κ::K
     tol::T
     extrapolant::Symbol
@@ -50,14 +48,13 @@ end
 
 function SBDF(order; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing,
+        linsolve = nothing, nlsolve = NLNewton(), κ = nothing,
         tol = nothing,
         extrapolant = :linear, ark = false)
     SBDF{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
         typeof(κ), typeof(tol)}(linsolve,
         nlsolve,
-        precs,
         κ,
         tol,
         extrapolant,
@@ -68,15 +65,14 @@ end
 # All keyword form needed for remake
 function SBDF(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing,
+        linsolve = nothing, nlsolve = NLNewton(), κ = nothing,
         tol = nothing,
         extrapolant = :linear,
         order, ark = false)
     SBDF{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
         typeof(κ), typeof(tol)}(linsolve,
         nlsolve,
-        precs,
         κ,
         tol,
         extrapolant,
@@ -136,11 +132,10 @@ Optional parameter kappa defaults to Shampine's accuracy-optimal -0.1850.
 
 See also `QNDF`.
 """
-struct QNDF1{CS, AD, F, F2, P, FDT, ST, CJ, κType, StepLimiter} <:
+struct QNDF1{CS, AD, F, F2, FDT, ST, CJ, κType, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     kappa::κType
     controller::Symbol
@@ -149,15 +144,14 @@ end
 
 function QNDF1(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear, kappa = -37//200,
         controller = :Standard, step_limiter! = trivial_limiter!)
     QNDF1{
         _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
         typeof(kappa), typeof(step_limiter!)}(linsolve,
         nlsolve,
-        precs,
         extrapolant,
         kappa,
         controller,
@@ -170,11 +164,10 @@ An adaptive order 2 quasi-constant timestep L-stable numerical differentiation f
 
 See also `QNDF`.
 """
-struct QNDF2{CS, AD, F, F2, P, FDT, ST, CJ, κType, StepLimiter} <:
+struct QNDF2{CS, AD, F, F2, FDT, ST, CJ, κType, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     kappa::κType
     controller::Symbol
@@ -183,15 +176,14 @@ end
 
 function QNDF2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear, kappa = -1 // 9,
         controller = :Standard, step_limiter! = trivial_limiter!)
     QNDF2{
         _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
         typeof(kappa), typeof(step_limiter!)}(linsolve,
         nlsolve,
-        precs,
         extrapolant,
         kappa,
         controller,
@@ -214,12 +206,11 @@ year={1997},
 publisher={SIAM}
 }
 """
-struct QNDF{MO, CS, AD, F, F2, P, FDT, ST, CJ, K, T, κType, StepLimiter} <:
+struct QNDF{MO, CS, AD, F, F2, FDT, ST, CJ, K, T, κType, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     max_order::Val{MO}
     linsolve::F
     nlsolve::F2
-    precs::P
     κ::K
     tol::T
     extrapolant::Symbol
@@ -231,15 +222,15 @@ end
 function QNDF(; max_order::Val{MO} = Val{5}(), chunk_size = Val{0}(),
         autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing,
+        linsolve = nothing, nlsolve = NLNewton(), κ = nothing,
         tol = nothing,
         extrapolant = :linear, kappa = (-37//200, -1//9, -823//10000, -83//2000, 0//1),
         controller = :Standard, step_limiter! = trivial_limiter!) where {MO}
     QNDF{MO, _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac),
         typeof(κ), typeof(tol), typeof(kappa), typeof(step_limiter!)}(
-        max_order, linsolve, nlsolve, precs, κ, tol,
+        max_order, linsolve, nlsolve, κ, tol,
         extrapolant, kappa, controller, step_limiter!)
 end
 
@@ -250,22 +241,20 @@ MEBDF2: Multistep Method
 The second order Modified Extended BDF method, which has improved stability properties over the standard BDF.
 Fixed timestep only.
 """
-struct MEBDF2{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct MEBDF2{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
 end
 function MEBDF2(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :constant)
     MEBDF2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac)}(linsolve,
         nlsolve,
-        precs,
         extrapolant)
 end
 
@@ -282,12 +271,11 @@ year={2002},
 publisher={Walter de Gruyter GmbH \\& Co. KG}
 }
 """
-struct FBDF{MO, CS, AD, F, F2, P, FDT, ST, CJ, K, T, StepLimiter} <:
+struct FBDF{MO, CS, AD, F, F2, FDT, ST, CJ, K, T, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     max_order::Val{MO}
     linsolve::F
     nlsolve::F2
-    precs::P
     κ::K
     tol::T
     extrapolant::Symbol
@@ -298,14 +286,14 @@ end
 function FBDF(; max_order::Val{MO} = Val{5}(), chunk_size = Val{0}(),
         autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing,
+        linsolve = nothing, nlsolve = NLNewton(), κ = nothing,
         tol = nothing,
         extrapolant = :linear, controller = :Standard, step_limiter! = trivial_limiter!) where {MO}
     FBDF{MO, _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac),
         typeof(κ), typeof(tol), typeof(step_limiter!)}(
-        max_order, linsolve, nlsolve, precs, κ, tol, extrapolant,
+        max_order, linsolve, nlsolve, κ, tol, extrapolant,
         controller, step_limiter!)
 end
 
@@ -389,41 +377,39 @@ See also `SBDF`, `IMEXEuler`.
 """
 IMEXEulerARK(; kwargs...) = SBDF(1; ark = true, kwargs...)
 
-struct DImplicitEuler{CS, AD, F, F2, P, FDT, ST, CJ} <: DAEAlgorithm{CS, AD, FDT, ST, CJ}
+struct DImplicitEuler{CS, AD, F, F2, FDT, ST, CJ} <: DAEAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     controller::Symbol
 end
 function DImplicitEuler(;
         chunk_size = Val{0}(), autodiff = true, standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :constant,
         controller = :Standard)
     DImplicitEuler{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac)}(linsolve,
-        nlsolve, precs, extrapolant, controller)
+        nlsolve, extrapolant, controller)
 end
 
-struct DABDF2{CS, AD, F, F2, P, FDT, ST, CJ} <: DAEAlgorithm{CS, AD, FDT, ST, CJ}
+struct DABDF2{CS, AD, F, F2, FDT, ST, CJ} <: DAEAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     controller::Symbol
 end
 function DABDF2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :constant,
         controller = :Standard)
     DABDF2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac)}(linsolve,
-        nlsolve, precs, extrapolant, controller)
+        nlsolve, extrapolant, controller)
 end
 
 #=
@@ -440,11 +426,10 @@ DBDF(;chunk_size=Val{0}(),autodiff=Val{true}(), standardtag = Val{true}(), concr
      linsolve,nlsolve,precs,extrapolant)
 =#
 
-struct DFBDF{MO, CS, AD, F, F2, P, FDT, ST, CJ, K, T} <: DAEAlgorithm{CS, AD, FDT, ST, CJ}
+struct DFBDF{MO, CS, AD, F, F2, FDT, ST, CJ, K, T} <: DAEAlgorithm{CS, AD, FDT, ST, CJ}
     max_order::Val{MO}
     linsolve::F
     nlsolve::F2
-    precs::P
     κ::K
     tol::T
     extrapolant::Symbol
@@ -453,13 +438,13 @@ end
 function DFBDF(; max_order::Val{MO} = Val{5}(), chunk_size = Val{0}(),
         autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing,
+        linsolve = nothing, nlsolve = NLNewton(), κ = nothing,
         tol = nothing,
         extrapolant = :linear, controller = :Standard) where {MO}
     DFBDF{MO, _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac),
-        typeof(κ), typeof(tol)}(max_order, linsolve, nlsolve, precs, κ, tol, extrapolant,
+        typeof(κ), typeof(tol)}(max_order, linsolve, nlsolve, κ, tol, extrapolant,
         controller)
 end
 
diff --git a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
index 6cbb4ec20a..d948c0816e 100644
--- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
+++ b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
@@ -107,7 +107,7 @@ end
 end
 const TryAgain = SlowConvergence
 
-DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
+DEFAULT_PRECS(W, p) = nothing, nothing
 isdiscretecache(cache) = false
 
 include("doc_utils.jl")
diff --git a/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl b/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
index bc1e42339d..bc87f456b6 100644
--- a/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
+++ b/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
@@ -27,7 +27,7 @@ import StaticArrays: SArray, MVector, SVector, @SVector, StaticArray, MMatrix, S
 using DiffEqBase: TimeGradientWrapper,
                   UJacobianWrapper, TimeDerivativeWrapper,
                   UDerivativeWrapper
-using SciMLBase: AbstractSciMLOperator
+using SciMLBase: AbstractSciMLOperator, DEIntegrator
 import OrdinaryDiffEqCore
 using OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, OrdinaryDiffEqAdaptiveImplicitAlgorithm,
                           DAEAlgorithm,
diff --git a/lib/OrdinaryDiffEqDifferentiation/src/linsolve_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/linsolve_utils.jl
index aa48161e1b..443ebca5fa 100644
--- a/lib/OrdinaryDiffEqDifferentiation/src/linsolve_utils.jl
+++ b/lib/OrdinaryDiffEqDifferentiation/src/linsolve_utils.jl
@@ -3,27 +3,20 @@ issuccess_W(W::Number) = !iszero(W)
 issuccess_W(::Any) = true
 
 function dolinsolve(integrator, linsolve; A = nothing, linu = nothing, b = nothing,
-        du = nothing, u = nothing, p = nothing, t = nothing,
-        weight = nothing, solverdata = nothing,
         reltol = integrator === nothing ? nothing : integrator.opts.reltol)
-    A !== nothing && (linsolve.A = A)
     b !== nothing && (linsolve.b = b)
     linu !== nothing && (linsolve.u = linu)
 
-    Plprev = linsolve.Pl isa LinearSolve.ComposePreconditioner ? linsolve.Pl.outer :
-             linsolve.Pl
-    Prprev = linsolve.Pr isa LinearSolve.ComposePreconditioner ? linsolve.Pr.outer :
-             linsolve.Pr
-
     _alg = unwrap_alg(integrator, true)
-
-    _Pl, _Pr = _alg.precs(linsolve.A, du, u, p, t, A !== nothing, Plprev, Prprev,
-        solverdata)
-    if (_Pl !== nothing || _Pr !== nothing)
-        __Pl = _Pl === nothing ? SciMLOperators.IdentityOperator(length(integrator.u)) : _Pl
-        __Pr = _Pr === nothing ? SciMLOperators.IdentityOperator(length(integrator.u)) : _Pr
-        linsolve.Pl = __Pl
-        linsolve.Pr = __Pr
+    if !isnothing(A)
+        if integrator isa DEIntegrator
+            (;u, p, t) = integrator
+            du = hasproperty(integrator, :du) ? integrator.du : nothing
+            p = (du, u, p, t)
+            reinit!(linsolve; A, p)
+        else
+            reinit!(linsolve; A)
+        end
     end
 
     linres = solve!(linsolve; reltol)
@@ -44,16 +37,18 @@ function dolinsolve(integrator, linsolve; A = nothing, linu = nothing, b = nothi
     return linres
 end
 
-function wrapprecs(_Pl::Nothing, _Pr::Nothing, weight, u)
-    Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight)))
-    Pr = Diagonal(_vec(weight))
-    Pl, Pr
-end
-
-function wrapprecs(_Pl, _Pr, weight, u)
-    Pl = _Pl === nothing ? SciMLOperators.IdentityOperator(length(u)) : _Pl
-    Pr = _Pr === nothing ? SciMLOperators.IdentityOperator(length(u)) : _Pr
-    Pl, Pr
+function wrapprecs(linsolver, W, weight)
+    if isnothing(linsolver)
+        linsolver = LinearSolve.defaultalg(W, weight, LinearSolve.OperatorAssumptions(true))
+    end
+    if hasproperty(linsolver, :precs) && isnothing(linsolver.precs)
+        Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight)))
+        Pr = Diagonal(_vec(weight))
+        precs = Returns((Pl, Pr))
+        return remake(linsolver; precs)
+    else
+        return linsolver
+    end
 end
 
 Base.resize!(p::LinearSolve.LinearCache, i) = p
diff --git a/lib/OrdinaryDiffEqExtrapolation/src/algorithms.jl b/lib/OrdinaryDiffEqExtrapolation/src/algorithms.jl
index 7b84e36bdc..d0780b7542 100644
--- a/lib/OrdinaryDiffEqExtrapolation/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqExtrapolation/src/algorithms.jl
@@ -21,10 +21,9 @@ ImplicitEulerExtrapolation: Parallelized Implicit Extrapolation Method
 Extrapolation of implicit Euler method with Romberg sequence.
 Similar to Hairer's SEULEX.
 """
-struct ImplicitEulerExtrapolation{CS, AD, F, P, FDT, ST, CJ, TO} <:
+struct ImplicitEulerExtrapolation{CS, AD, F, FDT, ST, CJ, TO} <:
        OrdinaryDiffEqImplicitExtrapolationAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
-    precs::P
     max_order::Int
     min_order::Int
     init_order::Int
@@ -35,7 +34,6 @@ end
 function ImplicitEulerExtrapolation(; chunk_size = Val{0}(), autodiff = true,
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward}, linsolve = nothing,
-        precs = DEFAULT_PRECS,
         max_order = 12, min_order = 3, init_order = 5,
         threading = false, sequence = :harmonic)
     linsolve = (linsolve === nothing &&
@@ -66,11 +64,10 @@ Initial order: " * lpad(init_order, 2, " ") * " --> " * lpad(init_order, 2, " ")
         sequence = :harmonic
     end
     ImplicitEulerExtrapolation{_unwrap_val(chunk_size), _unwrap_val(autodiff),
-        typeof(linsolve), typeof(precs), diff_type,
+        typeof(linsolve), diff_type,
         _unwrap_val(standardtag), _unwrap_val(concrete_jac),
-        typeof(threading)}(linsolve, precs, max_order, min_order,
-        init_order,
-        threading, sequence)
+        typeof(threading)}(linsolve, max_order, min_order,
+        init_order, threading, sequence)
 end
 """
 ExtrapolationMidpointDeuflhard: Parallelized Explicit Extrapolation Method
@@ -129,10 +126,9 @@ end
 ImplicitDeuflhardExtrapolation: Parallelized Implicit Extrapolation Method
 Midpoint extrapolation using Barycentric coordinates
 """
-struct ImplicitDeuflhardExtrapolation{CS, AD, F, P, FDT, ST, CJ, TO} <:
+struct ImplicitDeuflhardExtrapolation{CS, AD, F, FDT, ST, CJ, TO} <:
        OrdinaryDiffEqImplicitExtrapolationAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
-    precs::P
     min_order::Int # Minimal extrapolation order
     init_order::Int # Initial extrapolation order
     max_order::Int # Maximal extrapolation order
@@ -141,8 +137,7 @@ struct ImplicitDeuflhardExtrapolation{CS, AD, F, P, FDT, ST, CJ, TO} <:
 end
 function ImplicitDeuflhardExtrapolation(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
-        linsolve = nothing, precs = DEFAULT_PRECS,
-        diff_type = Val{:forward},
+        linsolve = nothing, diff_type = Val{:forward},
         min_order = 1, init_order = 5, max_order = 10,
         sequence = :harmonic, threading = false)
     # Enforce 1 <=  min_order <= init_order <= max_order:
@@ -177,9 +172,9 @@ Initial order: " * lpad(init_order, 2, " ") * " --> " * lpad(init_order, 2, " ")
 
     # Initialize algorithm
     ImplicitDeuflhardExtrapolation{_unwrap_val(chunk_size), _unwrap_val(autodiff),
-        typeof(linsolve), typeof(precs), diff_type,
+        typeof(linsolve), diff_type,
         _unwrap_val(standardtag), _unwrap_val(concrete_jac),
-        typeof(threading)}(linsolve, precs, min_order,
+        typeof(threading)}(linsolve, min_order,
         init_order, max_order,
         sequence, threading)
 end
@@ -242,10 +237,9 @@ end
 ImplicitHairerWannerExtrapolation: Parallelized Implicit Extrapolation Method
 Midpoint extrapolation using Barycentric coordinates, following Hairer's SODEX in the adaptivity behavior.
 """
-struct ImplicitHairerWannerExtrapolation{CS, AD, F, P, FDT, ST, CJ, TO} <:
+struct ImplicitHairerWannerExtrapolation{CS, AD, F, FDT, ST, CJ, TO} <:
        OrdinaryDiffEqImplicitExtrapolationAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
-    precs::P
     min_order::Int # Minimal extrapolation order
     init_order::Int # Initial extrapolation order
     max_order::Int # Maximal extrapolation order
@@ -256,7 +250,7 @@ end
 function ImplicitHairerWannerExtrapolation(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(),
         concrete_jac = nothing,
-        linsolve = nothing, precs = DEFAULT_PRECS,
+        linsolve = nothing,
         diff_type = Val{:forward},
         min_order = 2, init_order = 5, max_order = 10,
         sequence = :harmonic, threading = false)
@@ -292,21 +286,19 @@ Initial order: " * lpad(init_order, 2, " ") * " --> " * lpad(init_order, 2, " ")
 
     # Initialize algorithm
     ImplicitHairerWannerExtrapolation{_unwrap_val(chunk_size), _unwrap_val(autodiff),
-        typeof(linsolve), typeof(precs), diff_type,
+        typeof(linsolve), diff_type,
         _unwrap_val(standardtag), _unwrap_val(concrete_jac),
-        typeof(threading)}(linsolve, precs, min_order,
-        init_order,
-        max_order, sequence, threading)
+        typeof(threading)}(linsolve, min_order,
+        init_order, max_order, sequence, threading)
 end
 
 """
 ImplicitEulerBarycentricExtrapolation: Parallelized Implicit Extrapolation Method
 Euler extrapolation using Barycentric coordinates, following Hairer's SODEX in the adaptivity behavior.
 """
-struct ImplicitEulerBarycentricExtrapolation{CS, AD, F, P, FDT, ST, CJ, TO} <:
+struct ImplicitEulerBarycentricExtrapolation{CS, AD, F, FDT, ST, CJ, TO} <:
        OrdinaryDiffEqImplicitExtrapolationAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
-    precs::P
     min_order::Int # Minimal extrapolation order
     init_order::Int # Initial extrapolation order
     max_order::Int # Maximal extrapolation order
@@ -319,7 +311,7 @@ function ImplicitEulerBarycentricExtrapolation(; chunk_size = Val{0}(),
         autodiff = Val{true}(),
         standardtag = Val{true}(),
         concrete_jac = nothing,
-        linsolve = nothing, precs = DEFAULT_PRECS,
+        linsolve = nothing,
         diff_type = Val{:forward},
         min_order = 3, init_order = 5,
         max_order = 12, sequence = :harmonic,
@@ -356,10 +348,9 @@ Initial order: " * lpad(init_order, 2, " ") * " --> " * lpad(init_order, 2, " ")
 
     # Initialize algorithm
     ImplicitEulerBarycentricExtrapolation{_unwrap_val(chunk_size), _unwrap_val(autodiff),
-        typeof(linsolve), typeof(precs), diff_type,
+        typeof(linsolve), diff_type,
         _unwrap_val(standardtag),
         _unwrap_val(concrete_jac), typeof(threading)}(linsolve,
-        precs,
         min_order,
         init_order,
         max_order,
diff --git a/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_caches.jl b/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_caches.jl
index 7ede4e5b1a..88385fc466 100644
--- a/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_caches.jl
+++ b/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_caches.jl
@@ -263,18 +263,14 @@ function alg_cache(alg::ImplicitEulerExtrapolation, u, rate_prototype,
         linsolve_tmps[i] = zero(rate_prototype)
     end
 
-    linprob = LinearProblem(W[1], _vec(linsolve_tmps[1]); u0 = _vec(k_tmps[1]))
+    linprob = LinearProblem(W[1], _vec(linsolve_tmps[1]), (nothing, u, p, t); u0 = _vec(k_tmps[1]))
     linsolve1 = init(linprob, alg.linsolve, alias_A = true, alias_b = true)
-    #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-    #Pr = Diagonal(_vec(weight)))
 
     linsolve = Array{typeof(linsolve1), 1}(undef, Threads.nthreads())
     linsolve[1] = linsolve1
     for i in 2:Threads.nthreads()
-        linprob = LinearProblem(W[i], _vec(linsolve_tmps[i]); u0 = _vec(k_tmps[i]))
+        linprob = LinearProblem(W[i], _vec(linsolve_tmps[i]), (nothing, u, p, t); u0 = _vec(k_tmps[i]))
         linsolve[i] = init(linprob, alg.linsolve, alias_A = true, alias_b = true)
-        #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-        #Pr = Diagonal(_vec(weight)))
     end
 
     res = uEltypeNoUnits.(zero(u))
@@ -1150,18 +1146,14 @@ function alg_cache(alg::ImplicitDeuflhardExtrapolation, u, rate_prototype,
         linsolve_tmps[i] = zero(rate_prototype)
     end
 
-    linprob = LinearProblem(W[1], _vec(linsolve_tmps[1]); u0 = _vec(k_tmps[1]))
+    linprob = LinearProblem(W[1], _vec(linsolve_tmps[1]), (nothing, u, p, t); u0 = _vec(k_tmps[1]))
     linsolve1 = init(linprob, alg.linsolve, alias_A = true, alias_b = true)
-    #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-    #Pr = Diagonal(_vec(weight)))
 
     linsolve = Array{typeof(linsolve1), 1}(undef, Threads.nthreads())
     linsolve[1] = linsolve1
     for i in 2:Threads.nthreads()
-        linprob = LinearProblem(W[i], _vec(linsolve_tmps[i]); u0 = _vec(k_tmps[i]))
+        linprob = LinearProblem(W[i], _vec(linsolve_tmps[i]), (nothing, u, p, t); u0 = _vec(k_tmps[i]))
         linsolve[i] = init(linprob, alg.linsolve, alias_A = true, alias_b = true)
-        #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-        #Pr = Diagonal(_vec(weight)))
     end
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, du1, du2)
@@ -1478,18 +1470,14 @@ function alg_cache(alg::ImplicitHairerWannerExtrapolation, u, rate_prototype,
         linsolve_tmps[i] = zero(rate_prototype)
     end
 
-    linprob = LinearProblem(W[1], _vec(linsolve_tmps[1]); u0 = _vec(k_tmps[1]))
+    linprob = LinearProblem(W[1], _vec(linsolve_tmps[1]), (nothing, u, p, t); u0 = _vec(k_tmps[1]))
     linsolve1 = init(linprob, alg.linsolve, alias_A = true, alias_b = true)
-    #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-    #Pr = Diagonal(_vec(weight)))
 
     linsolve = Array{typeof(linsolve1), 1}(undef, Threads.nthreads())
     linsolve[1] = linsolve1
     for i in 2:Threads.nthreads()
-        linprob = LinearProblem(W[i], _vec(linsolve_tmps[i]); u0 = _vec(k_tmps[i]))
+        linprob = LinearProblem(W[i], _vec(linsolve_tmps[i]), (nothing, u, p, t); u0 = _vec(k_tmps[i]))
         linsolve[i] = init(linprob, alg.linsolve, alias_A = true, alias_b = true)
-        #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-        #Pr = Diagonal(_vec(weight)))
     end
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, du1, du2)
@@ -1674,18 +1662,14 @@ function alg_cache(alg::ImplicitEulerBarycentricExtrapolation, u, rate_prototype
         linsolve_tmps[i] = zero(rate_prototype)
     end
 
-    linprob = LinearProblem(W[1], _vec(linsolve_tmps[1]); u0 = _vec(k_tmps[1]))
+    linprob = LinearProblem(W[1], _vec(linsolve_tmps[1]), (nothing, u, p, t); u0 = _vec(k_tmps[1]))
     linsolve1 = init(linprob, alg.linsolve, alias_A = true, alias_b = true)
-    #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-    #Pr = Diagonal(_vec(weight)))
 
     linsolve = Array{typeof(linsolve1), 1}(undef, Threads.nthreads())
     linsolve[1] = linsolve1
     for i in 2:Threads.nthreads()
-        linprob = LinearProblem(W[i], _vec(linsolve_tmps[i]); u0 = _vec(k_tmps[i]))
+        linprob = LinearProblem(W[i], _vec(linsolve_tmps[i]), (nothing, u, p, t); u0 = _vec(k_tmps[i]))
         linsolve[i] = init(linprob, alg.linsolve, alias_A = true, alias_b = true)
-        #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-        #Pr = Diagonal(_vec(weight)))
     end
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, du1, du2)
diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
index d95430f5a5..e9ba0fc116 100644
--- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
@@ -15,10 +15,9 @@ publisher={Elsevier}
 RadauIIA3: Fully-Implicit Runge-Kutta Method
 An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
 """
-struct RadauIIA3{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
+struct RadauIIA3{CS, AD, F, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
-    precs::P
     extrapolant::Symbol
     κ::Tol
     maxiters::Int
@@ -31,16 +30,15 @@ end
 function RadauIIA3(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS,
+        linsolve = nothing,
         extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
         new_W_γdt_cutoff = 1 // 5,
         controller = :Predictive, κ = nothing, maxiters = 10,
         step_limiter! = trivial_limiter!)
     RadauIIA3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
         typeof(κ), typeof(fast_convergence_cutoff),
         typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
-        precs,
         extrapolant,
         κ,
         maxiters,
@@ -65,10 +63,9 @@ publisher={Elsevier}
 RadauIIA5: Fully-Implicit Runge-Kutta Method
 An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
 """
-struct RadauIIA5{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
+struct RadauIIA5{CS, AD, F, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     κ::Tol
@@ -82,16 +79,15 @@ end
 function RadauIIA5(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS,
+        linsolve = nothing,
         extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
         new_W_γdt_cutoff = 1 // 5,
         controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
         step_limiter! = trivial_limiter!)
     RadauIIA5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
         typeof(κ), typeof(fast_convergence_cutoff),
         typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
-        precs,
         smooth_est,
         extrapolant,
         κ,
@@ -117,10 +113,9 @@ publisher={Elsevier}
 RadauIIA9: Fully-Implicit Runge-Kutta Method
 An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
 """
-struct RadauIIA9{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
+struct RadauIIA9{CS, AD, F, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     κ::Tol
@@ -134,16 +129,15 @@ end
 function RadauIIA9(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS,
+        linsolve = nothing,
         extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
         new_W_γdt_cutoff = 1 // 5,
         controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
         step_limiter! = trivial_limiter!)
     RadauIIA9{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
         typeof(κ), typeof(fast_convergence_cutoff),
         typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
-        precs,
         smooth_est,
         extrapolant,
         κ,
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
index 516775e085..98293fe81d 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
@@ -107,11 +107,9 @@ function alg_cache(alg::RadauIIA3, u, rate_prototype, ::Type{uEltypeNoUnits},
     recursivefill!(atmp, false)
     jac_config = jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, dw12)
 
-    linprob = LinearProblem(W1, _vec(cubuff); u0 = _vec(dw12))
+    linprob = LinearProblem(W1, _vec(cubuff), (nothing,u,p,t); u0 = _vec(dw12))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
-    #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-    #Pr = Diagonal(_vec(weight)))
 
     rtol = reltol isa Number ? reltol : zero(reltol)
     atol = reltol isa Number ? reltol : zero(reltol)
@@ -252,16 +250,12 @@ function alg_cache(alg::RadauIIA5, u, rate_prototype, ::Type{uEltypeNoUnits},
     recursivefill!(atmp, false)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, dw1)
 
-    linprob = LinearProblem(W1, _vec(ubuff); u0 = _vec(dw1))
+    linprob = LinearProblem(W1, _vec(ubuff), (nothing,u,p,t); u0 = _vec(dw1))
     linsolve1 = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
-    #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-    #Pr = Diagonal(_vec(weight)))
-    linprob = LinearProblem(W2, _vec(cubuff); u0 = _vec(dw23))
+    linprob = LinearProblem(W2, _vec(cubuff), (nothing,u,p,t); u0 = _vec(dw23))
     linsolve2 = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
-    #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-    #Pr = Diagonal(_vec(weight)))
 
     rtol = reltol isa Number ? reltol : zero(reltol)
     atol = reltol isa Number ? reltol : zero(reltol)
@@ -441,21 +435,15 @@ function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits},
     recursivefill!(atmp, false)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, dw1)
 
-    linprob = LinearProblem(W1, _vec(ubuff); u0 = _vec(dw1))
+    linprob = LinearProblem(W1, _vec(ubuff), (nothing,u,p,t); u0 = _vec(dw1))
     linsolve1 = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
-    #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-    #Pr = Diagonal(_vec(weight)))
-    linprob = LinearProblem(W2, _vec(cubuff1); u0 = _vec(dw23))
+    linprob = LinearProblem(W2, _vec(cubuff1), (nothing,u,p,t); u0 = _vec(dw23))
     linsolve2 = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
-    #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-    #Pr = Diagonal(_vec(weight)))
-    linprob = LinearProblem(W3, _vec(cubuff2); u0 = _vec(dw45))
+    linprob = LinearProblem(W3, _vec(cubuff2), (nothing,u,p,t); u0 = _vec(dw45))
     linsolve3 = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
-    #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
-    #Pr = Diagonal(_vec(weight)))
 
     rtol = reltol isa Number ? reltol : zero(reltol)
     atol = reltol isa Number ? reltol : zero(reltol)
diff --git a/lib/OrdinaryDiffEqIMEXMultistep/src/algorithms.jl b/lib/OrdinaryDiffEqIMEXMultistep/src/algorithms.jl
index 2aca6a5bd6..1af29b8689 100644
--- a/lib/OrdinaryDiffEqIMEXMultistep/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqIMEXMultistep/src/algorithms.jl
@@ -1,42 +1,37 @@
 # IMEX Multistep methods
 
-struct CNAB2{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct CNAB2{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
 end
 
 function CNAB2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
-        extrapolant = :linear)
+        linsolve = nothing, nlsolve = NLNewton(), extrapolant = :linear)
     CNAB2{
         _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(
         linsolve,
         nlsolve,
-        precs,
         extrapolant)
 end
 
-struct CNLF2{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct CNLF2{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
 end
 function CNLF2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear)
     CNLF2{
         _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(
         linsolve,
         nlsolve,
-        precs,
         extrapolant)
 end
diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl
index 3bed6a453a..bf11eb0d49 100644
--- a/lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl
+++ b/lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl
@@ -191,14 +191,10 @@ function build_nlsolver(
             end
             jac_config = build_jac_config(alg, nf, uf, du1, uprev, u, ztmp, dz)
         end
-        linprob = LinearProblem(W, _vec(k); u0 = _vec(dz))
-        Pl, Pr = wrapprecs(
-            alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-                nothing)...,
-            weight, dz)
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-            Pl = Pl, Pr = Pr,
-            assumptions = LinearSolve.OperatorAssumptions(true))
+        linprob = LinearProblem(W, _vec(k), (isdae ? du1 : nothing,u,p,t); u0 = _vec(dz))
+        linsolve = init(linprob,
+                        wrapprecs(alg.linsolve, W, weight),
+                        (isdae ? du1 : nothing,u,p,t); alias_A = true, alias_b = true)
 
         tType = typeof(t)
         invγdt = inv(oneunit(t) * one(uTolType))
diff --git a/lib/OrdinaryDiffEqPDIRK/src/algorithms.jl b/lib/OrdinaryDiffEqPDIRK/src/algorithms.jl
index 264feb64a0..b0813d30ec 100644
--- a/lib/OrdinaryDiffEqPDIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqPDIRK/src/algorithms.jl
@@ -5,20 +5,19 @@
 PDIRK44: Parallel Diagonally Implicit Runge-Kutta Method
 A 2 processor 4th order diagonally non-adaptive implicit method.
 """
-struct PDIRK44{CS, AD, F, F2, P, FDT, ST, CJ, TO} <:
+struct PDIRK44{CS, AD, F, F2, FDT, ST, CJ, TO} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     threading::TO
 end
 function PDIRK44(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :constant, threading = true)
     PDIRK44{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac), typeof(threading)}(linsolve, nlsolve, precs,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac), typeof(threading)}(linsolve, nlsolve,
         extrapolant, threading)
 end
diff --git a/lib/OrdinaryDiffEqRosenbrock/src/algorithms.jl b/lib/OrdinaryDiffEqRosenbrock/src/algorithms.jl
index 03cbfae9d4..59a493b904 100644
--- a/lib/OrdinaryDiffEqRosenbrock/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqRosenbrock/src/algorithms.jl
@@ -104,22 +104,21 @@ for Alg in [
     :Rodas5Pe,
     :Rodas5Pr]
     @eval begin
-        struct $Alg{CS, AD, F, P, FDT, ST, CJ, StepLimiter, StageLimiter} <:
+        struct $Alg{CS, AD, F, FDT, ST, CJ, StepLimiter, StageLimiter} <:
                OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
             linsolve::F
-            precs::P
             step_limiter!::StepLimiter
             stage_limiter!::StageLimiter
         end
         function $Alg(; chunk_size = Val{0}(), autodiff = Val{true}(),
                 standardtag = Val{true}(), concrete_jac = nothing,
                 diff_type = Val{:forward}, linsolve = nothing,
-                precs = DEFAULT_PRECS, step_limiter! = trivial_limiter!,
+                step_limiter! = trivial_limiter!,
                 stage_limiter! = trivial_limiter!)
             $Alg{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-                typeof(precs), diff_type, _unwrap_val(standardtag),
+                diff_type, _unwrap_val(standardtag),
                 _unwrap_val(concrete_jac), typeof(step_limiter!),
-                typeof(stage_limiter!)}(linsolve, precs, step_limiter!,
+                typeof(stage_limiter!)}(linsolve, step_limiter!,
                 stage_limiter!)
         end
     end
@@ -147,20 +146,17 @@ end
     references = """
     https://doi.org/10.1016/j.cam.2009.09.017
     """)
-struct RosenbrockW6S4OS{CS, AD, F, P, FDT, ST, CJ} <:
+struct RosenbrockW6S4OS{CS, AD, F, FDT, ST, CJ} <:
        OrdinaryDiffEqRosenbrockAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
-    precs::P
 end
 function RosenbrockW6S4OS(; chunk_size = Val{0}(), autodiff = true,
         standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:central},
-        linsolve = nothing,
-        precs = DEFAULT_PRECS)
+        linsolve = nothing)
     RosenbrockW6S4OS{_unwrap_val(chunk_size),
-        _unwrap_val(autodiff), typeof(linsolve), typeof(precs), diff_type,
-        _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(linsolve,
-        precs)
+        _unwrap_val(autodiff), typeof(linsolve), diff_type,
+        _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(linsolve)
 end
 
 for Alg in [
@@ -185,18 +181,16 @@ for Alg in [
     :GRK4A,
     :Ros4LStab]
     @eval begin
-        struct $Alg{CS, AD, F, P, FDT, ST, CJ} <:
+        struct $Alg{CS, AD, F, FDT, ST, CJ} <:
                OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
             linsolve::F
-            precs::P
         end
         function $Alg(; chunk_size = Val{0}(), autodiff = Val{true}(),
                 standardtag = Val{true}(), concrete_jac = nothing,
-                diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS)
+                diff_type = Val{:forward}, linsolve = nothing)
             $Alg{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-                typeof(precs), diff_type, _unwrap_val(standardtag),
-                _unwrap_val(concrete_jac)}(linsolve,
-                precs)
+                diff_type, _unwrap_val(standardtag),
+                _unwrap_val(concrete_jac)}(linsolve)
         end
     end
 end
diff --git a/lib/OrdinaryDiffEqRosenbrock/src/generic_rosenbrock.jl b/lib/OrdinaryDiffEqRosenbrock/src/generic_rosenbrock.jl
index 5600b1bc27..0e0b54a238 100644
--- a/lib/OrdinaryDiffEqRosenbrock/src/generic_rosenbrock.jl
+++ b/lib/OrdinaryDiffEqRosenbrock/src/generic_rosenbrock.jl
@@ -247,7 +247,7 @@ function gen_algcache(cacheexpr::Expr,constcachename::Symbol,algname::Symbol,tab
             tf = TimeGradientWrapper(f,uprev,p)
             uf = UJacobianWrapper(f,t,p)
             linsolve_tmp = zero(rate_prototype)
-            linprob = LinearProblem(W,_vec(linsolve_tmp); u0=_vec(tmp))
+            linprob = LinearProblem(W,_vec(linsolve_tmp), (nothing, u, p, t); u0=_vec(tmp))
             linsolve = init(linprob,alg.linsolve,alias_A=true,alias_b=true,
                             Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
                             Pr = Diagonal(_vec(weight)))
@@ -995,7 +995,7 @@ references = """
 """,
 "Rodas3",
 references = """
-- Steinebach G. Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks 
+- Steinebach G. Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks
   within the Julia Differential Equations package.
   In: BIT Numerical Mathematics, 63(2), 2023
 """,
@@ -1056,7 +1056,7 @@ lower if not corrected).
 """,
 "Rodas4P",
 references = """
-- Steinebach G., Rodas23W / Rodas32P - a Rosenbrock-type method for DAEs with additional error estimate 
+- Steinebach G., Rodas23W / Rodas32P - a Rosenbrock-type method for DAEs with additional error estimate
   for dense output and Julia implementation,
   In progress.
 """,
@@ -1070,7 +1070,7 @@ of Roadas4P and in case of inexact Jacobians a second order W method.
 """,
 "Rodas4P2",
 references = """
-- Steinebach G., Rodas23W / Rodas32P - a Rosenbrock-type method for DAEs with additional error estimate 
+- Steinebach G., Rodas23W / Rodas32P - a Rosenbrock-type method for DAEs with additional error estimate
   for dense output and Julia implementation,
   In progress.
 """,
@@ -1095,7 +1095,7 @@ Has improved stability in the adaptive time stepping embedding.
 """,
 "Rodas5P",
 references = """
-- Steinebach G. Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks 
+- Steinebach G. Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks
   within the Julia Differential Equations package.
   In: BIT Numerical Mathematics, 63(2), 2023
 """,
@@ -1225,10 +1225,10 @@ references = """
     ROS2PRTableau()
 
 2nd order stiffly accurate Rosenbrock method with 3 internal stages with (Rinf=0).
-For problems with medium stiffness the convergence behaviour is very poor and it is recommended to use 
+For problems with medium stiffness the convergence behaviour is very poor and it is recommended to use
 [`ROS2S`](@ref) instead.
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function ROS2PRTableau() # 2nd order
@@ -1248,7 +1248,7 @@ end
 @doc rosenbrock_docstring(
 """
 2nd order stiffly accurate Rosenbrock method with 3 internal stages with (Rinf=0).
-For problems with medium stiffness the convergence behaviour is very poor and it is recommended to use 
+For problems with medium stiffness the convergence behaviour is very poor and it is recommended to use
 [`ROS2S`](@ref) instead.
 """,
 "ROS2PR",
@@ -1267,7 +1267,7 @@ ROS2PR
 2nd order stiffly accurate Rosenbrock-Wanner W-method with 3 internal stages with B_PR consistent of order 2 with (Rinf=0).
 More Information at https://doi.org/10.24355/dbbs.084-201408121139-0
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function ROS2STableau() # 2nd order
@@ -1335,7 +1335,7 @@ references = """
 
 3nd order stiffly accurate Rosenbrock-Wanner method with 3 internal stages with B_PR consistent of order 3, which is strongly A-stable with Rinf~=-0.73.
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function ROS3PRTableau() # 3rd order
@@ -1371,7 +1371,7 @@ references = """
 3nd order stiffly accurate Rosenbrock method with 3 internal stages with B_PR consistent of order 3, which is strongly A-stable with Rinf~=-0.73
 Convergence with order 4 for the stiff case, but has a poor accuracy.
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function Scholz4_7Tableau() # 3rd order
@@ -1599,7 +1599,7 @@ references = """
 B_PR consistent of order 2 with Rinf=0.
 The order of convergence decreases if medium stiff problems are considered, but it has good results for very stiff cases.
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function ROS3PRLTableau() # 3rd order
@@ -1639,7 +1639,7 @@ references = """
 B_PR consistent of order 3.
 The order of convergence does NOT decreases if medium stiff problems are considered as it does for [`ROS3PRL`](@ref).
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function ROS3PRL2Tableau() # 3rd order
@@ -1677,7 +1677,7 @@ references = """
 
 4rd order L-stable Rosenbrock-Krylov method with 4 internal stages,
 with a 3rd order embedded method which is strongly A-stable with Rinf~=0.55. (when using exact Jacobians)
-Tranquilli, Paul and Sandu, Adrian (2014): Rosenbrock--Krylov Methods for Large Systems of Differential Equations 
+Tranquilli, Paul and Sandu, Adrian (2014): Rosenbrock--Krylov Methods for Large Systems of Differential Equations
 https://doi.org/10.1137/130923336
 """
 function ROK4aTableau() # 4rd order
@@ -1703,7 +1703,7 @@ with a 3rd order embedded method which is strongly A-stable with Rinf~=0.55. (wh
 """,
 "ROK4a",
 references = """
-- Tranquilli, Paul and Sandu, Adrian (2014): 
+- Tranquilli, Paul and Sandu, Adrian (2014):
   Rosenbrock--Krylov Methods for Large Systems of Differential Equations
   https://doi.org/10.1137/130923336
 """) ROK4a
diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl
index 4a266e77f3..214815c3fa 100644
--- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl
+++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl
@@ -99,12 +99,8 @@ function alg_cache(alg::Rosenbrock23, u, rate_prototype, ::Type{uEltypeNoUnits},
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
 
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
 
     grad_config = build_grad_config(alg, f, tf, du1, t)
@@ -144,13 +140,9 @@ function alg_cache(alg::Rosenbrock32, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
 
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2, Val(false))
@@ -293,12 +285,8 @@ function alg_cache(alg::ROS3P, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
@@ -379,12 +367,8 @@ function alg_cache(alg::Rodas3, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
@@ -572,12 +556,8 @@ function alg_cache(alg::Rodas23W, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
@@ -616,12 +596,8 @@ function alg_cache(alg::Rodas3P, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
@@ -739,12 +715,8 @@ function alg_cache(alg::Rodas4, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
@@ -799,12 +771,8 @@ function alg_cache(alg::Rodas42, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
@@ -859,12 +827,8 @@ function alg_cache(alg::Rodas4P, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
@@ -919,12 +883,8 @@ function alg_cache(alg::Rodas4P2, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
@@ -1036,12 +996,8 @@ function alg_cache(alg::Rodas5, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
@@ -1100,13 +1056,8 @@ function alg_cache(
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
-    Pl, Pr = wrapprecs(
-        alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
-            nothing)..., weight, tmp)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr,
-        assumptions = LinearSolve.OperatorAssumptions(true))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
+    linsolve = init(linprob, wrapprecs(alg.linsolve, W, weight), alias_A = true, alias_b = true)
     grad_config = build_grad_config(alg, f, tf, du1, t)
     jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
     Rosenbrock5Cache(u, uprev, dense1, dense2, dense3, du, du1, du2, k1, k2, k3, k4,
diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl
index a724abbc52..fac119f814 100644
--- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl
+++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl
@@ -48,16 +48,7 @@ end
         integrator.opts.abstol, integrator.opts.reltol,
         integrator.opts.internalnorm, t)
 
-    if repeat_step
-        linres = dolinsolve(
-            integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp),
-            du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = γ))
-    else
-        linres = dolinsolve(integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp),
-            du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = γ))
-    end
+    linres = dolinsolve(integrator, cache.linsolve; A = repeat_step ? nothing : W, b = _vec(linsolve_tmp))
 
     vecu = _vec(linres.u)
     veck₁ = _vec(k₁)
@@ -160,16 +151,7 @@ end
         integrator.opts.abstol, integrator.opts.reltol,
         integrator.opts.internalnorm, t)
 
-    if repeat_step
-        linres = dolinsolve(
-            integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp),
-            du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = γ))
-    else
-        linres = dolinsolve(integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp),
-            du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = γ))
-    end
+    linres = dolinsolve(integrator, cache.linsolve; A = repeat_step ? nothing : W, b = _vec(linsolve_tmp))
 
     vecu = _vec(linres.u)
     veck₁ = _vec(k₁)
@@ -517,16 +499,7 @@ end
         integrator.opts.abstol, integrator.opts.reltol,
         integrator.opts.internalnorm, t)
 
-    if repeat_step
-        linres = dolinsolve(
-            integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp),
-            du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = dtgamma))
-    else
-        linres = dolinsolve(integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp),
-            du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = dtgamma))
-    end
+    linres = dolinsolve(integrator, cache.linsolve; A = repeat_step ? nothing : W, b = _vec(linsolve_tmp))
 
     vecu = _vec(linres.u)
     veck1 = _vec(k1)
@@ -712,16 +685,7 @@ end
         integrator.opts.abstol, integrator.opts.reltol,
         integrator.opts.internalnorm, t)
 
-    if repeat_step
-        linres = dolinsolve(
-            integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp),
-            du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = dtgamma))
-    else
-        linres = dolinsolve(integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp),
-            du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = dtgamma))
-    end
+    linres = dolinsolve(integrator, cache.linsolve; A = repeat_step ? nothing : W, b = _vec(linsolve_tmp))
 
     vecu = _vec(linres.u)
     veck1 = _vec(k1)
@@ -1020,16 +984,7 @@ end
         integrator.opts.abstol, integrator.opts.reltol,
         integrator.opts.internalnorm, t)
 
-    if repeat_step
-        linres = dolinsolve(
-            integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp),
-            du = cache.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = dtgamma))
-    else
-        linres = dolinsolve(integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp),
-            du = cache.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = dtgamma))
-    end
+    linres = dolinsolve(integrator, cache.linsolve; A = repeat_step ? nothing : W, b = _vec(linsolve_tmp))
 
     @.. broadcast=false $(_vec(k1))=-linres.u
 
@@ -1383,16 +1338,7 @@ end
         integrator.opts.abstol, integrator.opts.reltol,
         integrator.opts.internalnorm, t)
 
-    if repeat_step
-        linres = dolinsolve(
-            integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp),
-            du = cache.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = dtgamma))
-    else
-        linres = dolinsolve(integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp),
-            du = cache.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = dtgamma))
-    end
+    linres = dolinsolve(integrator, cache.linsolve; A = repeat_step ? nothing : W, b = _vec(linsolve_tmp))
 
     @.. broadcast=false $(_vec(k1))=-linres.u
 
@@ -1786,16 +1732,7 @@ end
         integrator.opts.abstol, integrator.opts.reltol,
         integrator.opts.internalnorm, t)
 
-    if repeat_step
-        linres = dolinsolve(
-            integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp),
-            du = cache.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = dtgamma))
-    else
-        linres = dolinsolve(integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp),
-            du = cache.fsalfirst, u = u, p = p, t = t, weight = weight,
-            solverdata = (; gamma = dtgamma))
-    end
+    linres = dolinsolve(integrator, cache.linsolve; A = repeat_step ? nothing : W, b = _vec(linsolve_tmp))
 
     vecu = _vec(linres.u)
     veck1 = _vec(k1)
diff --git a/lib/OrdinaryDiffEqSDIRK/src/algorithms.jl b/lib/OrdinaryDiffEqSDIRK/src/algorithms.jl
index 01f796f650..8a105c7203 100644
--- a/lib/OrdinaryDiffEqSDIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqSDIRK/src/algorithms.jl
@@ -4,11 +4,10 @@ ImplicitEuler: SDIRK Method
 A 1st order implicit solver. A-B-L-stable. Adaptive timestepping through a divided differences estimate via memory.
 Strong-stability preserving (SSP).
 """
-struct ImplicitEuler{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct ImplicitEuler{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     controller::Symbol
     step_limiter!::StepLimiter
@@ -17,24 +16,23 @@ end
 function ImplicitEuler(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :constant,
         controller = :PI, step_limiter! = trivial_limiter!)
     ImplicitEuler{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve,
-        nlsolve, precs, extrapolant, controller, step_limiter!)
+        nlsolve, extrapolant, controller, step_limiter!)
 end
 """
 ImplicitMidpoint: SDIRK Method
 A second order A-stable symplectic and symmetric implicit solver.
 Good for highly stiff equations which need symplectic integration.
 """
-struct ImplicitMidpoint{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct ImplicitMidpoint{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     step_limiter!::StepLimiter
 end
@@ -42,13 +40,12 @@ end
 function ImplicitMidpoint(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear, step_limiter! = trivial_limiter!)
     ImplicitMidpoint{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve,
         nlsolve,
-        precs,
         extrapolant,
         step_limiter!)
 end
@@ -64,11 +61,10 @@ Also known as Crank-Nicolson when applied to PDEs. Adaptive timestepping via div
 differences approximation to the second derivative terms in the local truncation error
 estimate (the SPICE approximation strategy).
 """
-struct Trapezoid{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct Trapezoid{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     controller::Symbol
     step_limiter!::StepLimiter
@@ -77,14 +73,13 @@ end
 function Trapezoid(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear,
         controller = :PI, step_limiter! = trivial_limiter!)
     Trapezoid{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve,
         nlsolve,
-        precs,
         extrapolant,
         controller,
         step_limiter!)
@@ -106,11 +101,10 @@ TRBDF2: SDIRK Method
 A second order A-B-L-S-stable one-step ESDIRK method.
 Includes stiffness-robust error estimates for accurate adaptive timestepping, smoothed derivatives for highly stiff and oscillatory problems.
 """
-struct TRBDF2{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct TRBDF2{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -119,12 +113,12 @@ end
 
 function TRBDF2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI, step_limiter! = trivial_limiter!)
     TRBDF2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve,
         smooth_est, extrapolant, controller, step_limiter!)
 end
 
@@ -145,11 +139,10 @@ publisher={ACM}
 SDIRK2: SDIRK Method
 An A-B-L stable 2nd order SDIRK method
 """
-struct SDIRK2{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct SDIRK2{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -158,22 +151,21 @@ end
 
 function SDIRK2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI, step_limiter! = trivial_limiter!)
     SDIRK2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac), typeof(step_limiter!)}(
-        linsolve, nlsolve, precs, smooth_est, extrapolant,
+        linsolve, nlsolve, smooth_est, extrapolant,
         controller,
         step_limiter!)
 end
 
-struct SDIRK22{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct SDIRK22{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     controller::Symbol
     step_limiter!::StepLimiter
@@ -182,24 +174,22 @@ end
 function SDIRK22(;
         chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear,
         controller = :PI, step_limiter! = trivial_limiter!)
     Trapezoid{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve,
         nlsolve,
-        precs,
         extrapolant,
         controller,
         step_limiter!)
 end
 
-struct SSPSDIRK2{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct SSPSDIRK2{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} # Not adaptive
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -208,12 +198,12 @@ end
 function SSPSDIRK2(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :constant,
         controller = :PI)
     SSPSDIRK2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac)}(linsolve, nlsolve, smooth_est, extrapolant,
         controller)
 end
 
@@ -232,11 +222,10 @@ publisher={Springer}
 Kvaerno3: SDIRK Method
 An A-L stable stiffly-accurate 3rd order ESDIRK method
 """
-struct Kvaerno3{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct Kvaerno3{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -245,12 +234,12 @@ end
 function Kvaerno3(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI, step_limiter! = trivial_limiter!)
     Kvaerno3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve,
         smooth_est, extrapolant, controller, step_limiter!)
 end
 
@@ -265,11 +254,10 @@ publisher={National Aeronautics and Space Administration, Langley Research Cente
 KenCarp3: SDIRK Method
 An A-L stable stiffly-accurate 3rd order ESDIRK method with splitting
 """
-struct KenCarp3{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct KenCarp3{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -278,32 +266,30 @@ end
 function KenCarp3(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI, step_limiter! = trivial_limiter!)
     KenCarp3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve,
         smooth_est, extrapolant, controller, step_limiter!)
 end
 
-struct CFNLIRK3{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct CFNLIRK3{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
 end
 function CFNLIRK3(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear)
     CFNLIRK3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac)}(linsolve,
         nlsolve,
-        precs,
         extrapolant)
 end
 
@@ -322,11 +308,10 @@ publisher={ACM}
 Cash4: SDIRK Method
 An A-L stable 4th order SDIRK method
 """
-struct Cash4{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct Cash4{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     embedding::Int
@@ -334,122 +319,111 @@ struct Cash4{CS, AD, F, F2, P, FDT, ST, CJ} <:
 end
 function Cash4(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI, embedding = 3)
     Cash4{
         _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(
         linsolve,
         nlsolve,
-        precs,
         smooth_est,
         extrapolant,
         embedding,
         controller)
 end
 
-struct SFSDIRK4{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct SFSDIRK4{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
 end
 function SFSDIRK4(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear)
     SFSDIRK4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac)}(linsolve,
         nlsolve,
-        precs,
         extrapolant)
 end
 
-struct SFSDIRK5{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct SFSDIRK5{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
 end
 
 function SFSDIRK5(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear)
     SFSDIRK5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac)}(linsolve,
         nlsolve,
-        precs,
         extrapolant)
 end
 
-struct SFSDIRK6{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct SFSDIRK6{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
 end
 
 function SFSDIRK6(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear)
     SFSDIRK6{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac)}(linsolve,
         nlsolve,
-        precs,
         extrapolant)
 end
 
-struct SFSDIRK7{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct SFSDIRK7{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
 end
 
 function SFSDIRK7(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear)
     SFSDIRK7{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac)}(linsolve,
         nlsolve,
-        precs,
         extrapolant)
 end
 
-struct SFSDIRK8{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct SFSDIRK8{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
 end
 
 function SFSDIRK8(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear)
     SFSDIRK8{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
         _unwrap_val(concrete_jac)}(linsolve,
         nlsolve,
-        precs,
         extrapolant)
 end
 
@@ -461,11 +435,10 @@ Springer (1996)
 Hairer4: SDIRK Method
 An A-L stable 4th order SDIRK method
 """
-struct Hairer4{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct Hairer4{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -473,12 +446,12 @@ end
 function Hairer4(;
         chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI)
     Hairer4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac)}(linsolve, nlsolve, smooth_est, extrapolant,
         controller)
 end
 
@@ -490,11 +463,10 @@ Springer (1996)
 Hairer42: SDIRK Method
 An A-L stable 4th order SDIRK method
 """
-struct Hairer42{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct Hairer42{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -502,12 +474,12 @@ end
 function Hairer42(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI)
     Hairer42{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac)}(linsolve, nlsolve, smooth_est, extrapolant,
         controller)
 end
 
@@ -526,11 +498,10 @@ publisher={Springer}
 Kvaerno4: SDIRK Method
 An A-L stable stiffly-accurate 4th order ESDIRK method.
 """
-struct Kvaerno4{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct Kvaerno4{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -539,12 +510,12 @@ end
 function Kvaerno4(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI, step_limiter! = trivial_limiter!)
     Kvaerno4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve,
         smooth_est, extrapolant, controller, step_limiter!)
 end
 
@@ -563,11 +534,10 @@ publisher={Springer}
 Kvaerno5: SDIRK Method
 An A-L stable stiffly-accurate 5th order ESDIRK method
 """
-struct Kvaerno5{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct Kvaerno5{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -576,12 +546,12 @@ end
 function Kvaerno5(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI, step_limiter! = trivial_limiter!)
     Kvaerno5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve,
         smooth_est, extrapolant, controller, step_limiter!)
 end
 
@@ -596,11 +566,10 @@ publisher={National Aeronautics and Space Administration, Langley Research Cente
 KenCarp4: SDIRK Method
 An A-L stable stiffly-accurate 4th order ESDIRK method with splitting
 """
-struct KenCarp4{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct KenCarp4{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -609,12 +578,12 @@ end
 function KenCarp4(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI, step_limiter! = trivial_limiter!)
     KenCarp4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve,
         smooth_est, extrapolant, controller, step_limiter!)
 end
 
@@ -634,11 +603,10 @@ publisher={Elsevier}
 KenCarp47: SDIRK Method
 An A-L stable stiffly-accurate 4th order seven-stage ESDIRK method with splitting
 """
-struct KenCarp47{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct KenCarp47{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -646,12 +614,12 @@ end
 function KenCarp47(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI)
     KenCarp47{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac)}(linsolve, nlsolve, smooth_est, extrapolant,
         controller)
 end
 
@@ -666,11 +634,10 @@ publisher={National Aeronautics and Space Administration, Langley Research Cente
 KenCarp5: SDIRK Method
 An A-L stable stiffly-accurate 5th order ESDIRK method with splitting
 """
-struct KenCarp5{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+struct KenCarp5{CS, AD, F, F2, FDT, ST, CJ, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -679,12 +646,12 @@ end
 function KenCarp5(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI, step_limiter! = trivial_limiter!)
     KenCarp5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve,
         smooth_est, extrapolant, controller, step_limiter!)
 end
 """
@@ -701,11 +668,10 @@ publisher={Elsevier}
 KenCarp58: SDIRK Method
 An A-L stable stiffly-accurate 5th order eight-stage ESDIRK method with splitting
 """
-struct KenCarp58{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct KenCarp58{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
@@ -713,32 +679,31 @@ end
 function KenCarp58(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
         controller = :PI)
     KenCarp58{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac)}(linsolve, nlsolve, smooth_est, extrapolant,
         controller)
 end
 
 # `smooth_est` is not necessary, as the embedded method is also L-stable
-struct ESDIRK54I8L2SA{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct ESDIRK54I8L2SA{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     controller::Symbol
 end
 function ESDIRK54I8L2SA(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear, controller = :PI)
     ESDIRK54I8L2SA{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac)}(linsolve, nlsolve, extrapolant,
         controller)
 end
 
@@ -752,22 +717,21 @@ volume={146},
 pages={221-244}
 }
 """
-struct ESDIRK436L2SA2{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct ESDIRK436L2SA2{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     controller::Symbol
 end
 function ESDIRK436L2SA2(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear, controller = :PI)
     ESDIRK436L2SA2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac)}(linsolve, nlsolve, extrapolant,
         controller)
 end
 
@@ -781,22 +745,21 @@ volume={146},
 pages={221-244}
 }
 """
-struct ESDIRK437L2SA{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct ESDIRK437L2SA{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     controller::Symbol
 end
 function ESDIRK437L2SA(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear, controller = :PI)
     ESDIRK437L2SA{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac)}(linsolve, nlsolve, extrapolant,
         controller)
 end
 
@@ -810,22 +773,21 @@ volume={146},
 pages={221-244}
 }
 """
-struct ESDIRK547L2SA2{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct ESDIRK547L2SA2{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     controller::Symbol
 end
 function ESDIRK547L2SA2(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear, controller = :PI)
     ESDIRK547L2SA2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac)}(linsolve, nlsolve, extrapolant,
         controller)
 end
 
@@ -842,21 +804,20 @@ Currently has STABILITY ISSUES, causing it to fail the adaptive tests.
 Check issue https://github.com/SciML/OrdinaryDiffEq.jl/issues/1933 for more details.
 }
 """
-struct ESDIRK659L2SA{CS, AD, F, F2, P, FDT, ST, CJ} <:
+struct ESDIRK659L2SA{CS, AD, F, F2, FDT, ST, CJ} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     extrapolant::Symbol
     controller::Symbol
 end
 function ESDIRK659L2SA(; chunk_size = Val{0}(), autodiff = Val{true}(),
         standardtag = Val{true}(), concrete_jac = nothing,
         diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        linsolve = nothing, nlsolve = NLNewton(),
         extrapolant = :linear, controller = :PI)
     ESDIRK659L2SA{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant,
+        typeof(nlsolve), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac)}(linsolve, nlsolve, extrapolant,
         controller)
 end
diff --git a/lib/OrdinaryDiffEqStabilizedIRK/src/algorithms.jl b/lib/OrdinaryDiffEqStabilizedIRK/src/algorithms.jl
index 1a59689554..b23b397c9f 100644
--- a/lib/OrdinaryDiffEqStabilizedIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqStabilizedIRK/src/algorithms.jl
@@ -1,8 +1,7 @@
-struct IRKC{CS, AD, F, F2, P, FDT, ST, CJ, K, T, E} <:
+struct IRKC{CS, AD, F, F2, FDT, ST, CJ, K, T, E} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
-    precs::P
     κ::K
     tol::T
     extrapolant::Symbol
@@ -12,11 +11,11 @@ end
 
 function IRKC(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
         concrete_jac = nothing, diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing,
+        linsolve = nothing, nlsolve = NLNewton(), κ = nothing,
         tol = nothing,
         extrapolant = :linear, controller = :Standard, eigen_est = nothing)
     IRKC{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
-        typeof(κ), typeof(tol), typeof(eigen_est)}(linsolve, nlsolve, precs, κ, tol,
+        diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(κ), typeof(tol), typeof(eigen_est)}(linsolve, nlsolve, κ, tol,
         extrapolant, controller, eigen_est)
 end
diff --git a/test/interface/default_solver_tests.jl b/test/interface/default_solver_tests.jl
index 4cdef37e80..1e36a9b842 100644
--- a/test/interface/default_solver_tests.jl
+++ b/test/interface/default_solver_tests.jl
@@ -48,7 +48,7 @@ rosensol = solve(prob_rober, AutoTsit5(Rosenbrock23(autodiff = false)))
 sol = solve(prob_rober, reltol = 1e-7, abstol = 1e-7)
 rosensol = solve(
     prob_rober, AutoVern7(Rodas5P(autodiff = false)), reltol = 1e-7, abstol = 1e-7)
-# test that default has the same performance as AutoTsit5(Rosenbrock23()) (which we expect it to use for this).
+# test that default has the same performance as AutoTsit5(Rodas5P()) (which we expect it to use for this).
 @test sol.stats.naccept == rosensol.stats.naccept
 @test sol.stats.nf == rosensol.stats.nf
 @test unique(sol.alg_choice) == [2, 4]
@@ -74,7 +74,7 @@ for n in (100, 600)
         vcat([1.0, 0.0, 0.0], ones(n)), (0.0, 100.0), (0.04, 3e7, 1e4))
     sol = solve(prob_ex_rober)
     fsol = solve(prob_ex_rober, AutoTsit5(FBDF(; autodiff = false, linsolve)))
-    # test that default has the same performance as AutoTsit5(Rosenbrock23()) (which we expect it to use for this).
+    # test that default has the same performance as AutoTsit5(FBDF()) (which we expect it to use for this).
     @test sol.stats.naccept == fsol.stats.naccept
     @test sol.stats.nf == fsol.stats.nf
     @test unique(sol.alg_choice) == [1, stiffalg]
diff --git a/test/interface/linear_nonlinear_tests.jl b/test/interface/linear_nonlinear_tests.jl
index b948996346..e799ca2616 100644
--- a/test/interface/linear_nonlinear_tests.jl
+++ b/test/interface/linear_nonlinear_tests.jl
@@ -8,33 +8,22 @@ end
 u0 = rand(3)
 prob = ODEProblem(rn, u0, (0, 50.0))
 
-function precsl(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
-    if newW === nothing || newW
-        Pl = lu(convert(AbstractMatrix, W), check = false)
-    else
-        Pl = Plprev
-    end
-    Pl, nothing
+function precsl(W, p)
+    Pl = lu(convert(AbstractMatrix, W), check = false)
+    Pl, IdentityOperator(size(W, 1))
 end
 
-function precsr(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
-    if newW === nothing || newW
-        Pr = lu(convert(AbstractMatrix, W), check = false)
-    else
-        Pr = Prprev
-    end
-    nothing, Pr
+function precsr(W, p)
+    Pr = lu(convert(AbstractMatrix, W), check = false)
+    IdentityOperator(size(W, 1)), Pr
 end
 
-function precslr(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
-    if newW === nothing || newW
-        Pr = lu(convert(AbstractMatrix, W), check = false)
-    else
-        Pr = Prprev
-    end
+function precslr(W, p)
+    Pr = lu(convert(AbstractMatrix, W), check = false)
     Pr, Pr
 end
 
+
 sol = @test_nowarn solve(prob, TRBDF2(autodiff = false));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
@@ -45,16 +34,16 @@ solref = @test_nowarn solve(prob,
         smooth_est = false));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
-        precs = precsl, smooth_est = false, concrete_jac = true));
+    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(precs = precsl),
+        smooth_est = false, concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
-        precs = precsr, smooth_est = false, concrete_jac = true));
+    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(precs = precsr),
+           smooth_est = false, concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
-        precs = precslr, smooth_est = false, concrete_jac = true));
+    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(precs = precslr)
+        , smooth_est = false, concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
     QNDF(autodiff = false, linsolve = KrylovJL_GMRES(),
@@ -62,12 +51,12 @@ sol = @test_nowarn solve(prob,
 @test length(sol.t) < 25
 sol = @test_nowarn solve(prob,
     Rosenbrock23(autodiff = false,
-        linsolve = KrylovJL_GMRES(),
-        precs = precslr, concrete_jac = true));
+        linsolve = KrylovJL_GMRES(precs = precslr),
+        concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-    Rodas4(autodiff = false, linsolve = KrylovJL_GMRES(),
-        precs = precslr, concrete_jac = true));
+    Rodas4(autodiff = false, linsolve = KrylovJL_GMRES(precs = precslr),
+        concrete_jac = true));
 @test length(sol.t) < 20
 
 sol = @test_nowarn solve(prob, TRBDF2(autodiff = false));
@@ -79,26 +68,26 @@ sol = @test_nowarn solve(prob,
         smooth_est = false));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
-        precs = precsl, smooth_est = false, concrete_jac = true));
+    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(precs = precsl),
+        smooth_est = false, concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
-        precs = precsr, smooth_est = false, concrete_jac = true));
+    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(precs = precsr),
+        smooth_est = false, concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
-        precs = precslr, smooth_est = false, concrete_jac = true));
+    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(precs = precslr),
+        smooth_est = false, concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
     QNDF(autodiff = false, linsolve = KrylovJL_GMRES(),
         concrete_jac = true));
 @test length(sol.t) < 25
 sol = @test_nowarn solve(prob,
-    Rosenbrock23(autodiff = false, linsolve = KrylovJL_GMRES(),
-        precs = precslr, concrete_jac = true));
+    Rosenbrock23(autodiff = false, linsolve = KrylovJL_GMRES(precs = precslr),
+        concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-    Rodas4(autodiff = false, linsolve = KrylovJL_GMRES(),
-        precs = precslr, concrete_jac = true));
+    Rodas4(autodiff = false, linsolve = KrylovJL_GMRES(precs = precslr),
+        concrete_jac = true));
 @test length(sol.t) < 20
diff --git a/test/interface/linear_solver_test.jl b/test/interface/linear_solver_test.jl
index 12a449d603..516d5ab299 100644
--- a/test/interface/linear_solver_test.jl
+++ b/test/interface/linear_solver_test.jl
@@ -197,7 +197,7 @@ refsol = solve(probiip, FBDF(), abstol = 1e-12, reltol = 1e-12)
         @testset "$solname" for (solname, solver) in pairs(solvers)
             sol = solve(prob, solver, abstol = 1e-12, reltol = 1e-12, maxiters = 2e4)
             @test sol.retcode == ReturnCode.Success
-            @test isapprox(sol.u[end], refsol.u[end], rtol = 1e-8, atol = 1e-10)
+            @test isapprox(sol.u[end], refsol.u[end], rtol = 2e-8, atol = 1e-10)
         end
     end
 end
@@ -207,7 +207,7 @@ end
         @testset "$solname" for (solname, solver) in pairs(solvers)
             sol = solve(prob, solver, maxiters = 2e4)
             @test sol.retcode == ReturnCode.Success
-            @test isapprox(sol.u[end], refsol.u[end], rtol = 2e-3, atol = 1e-6)
+            @test isapprox(sol.u[end], refsol.u[end], rtol = 5e-3, atol = 1e-6)
         end
     end
 end
diff --git a/test/interface/preconditioners.jl b/test/interface/preconditioners.jl
index 53188d1a12..c41ada32e3 100644
--- a/test/interface/preconditioners.jl
+++ b/test/interface/preconditioners.jl
@@ -57,38 +57,25 @@ prob_ode_brusselator_2d_sparse = ODEProblem(
         jac_prototype = float.(jac)),
     u0, (0.0, 11.5), p)
 
-function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
-    if newW === nothing || newW
-        Pl = ilu(convert(AbstractMatrix, W), τ = 50.0)
-    else
-        Pl = Plprev
-    end
-    Pl, nothing
+function incompletelu(W, p)
+    Pl = ilu(convert(AbstractMatrix, W), τ = 50.0)
+    Pl, IdentityOperator(size(W, 1))
 end
 
-function algebraicmultigrid(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
-    if newW === nothing || newW
-        Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(convert(
-            AbstractMatrix,
-            W)))
-    else
-        Pl = Plprev
-    end
-    Pl, nothing
+function algebraicmultigrid(W, p)
+    A = convert(AbstractMatrix, W)
+    Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A))
+    Pl, IdentityOperator(size(W, 1))
 end
 
-function algebraicmultigrid2(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
-    if newW === nothing || newW
-        A = convert(AbstractMatrix, W)
-        Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A,
-            presmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
-                1))),
-            postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
-                1)))))
-    else
-        Pl = Plprev
-    end
-    Pl, nothing
+function algebraicmultigrid2(W, p)
+    A = convert(AbstractMatrix, W)
+    Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A,
+        presmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
+            1))),
+        postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
+            1)))))
+    Pl, IdentityOperator(size(W, 1))
 end
 
 iter[] = 0
@@ -97,17 +84,17 @@ sol1 = solve(prob_ode_brusselator_2d, KenCarp47(linsolve = KrylovJL_GMRES()),
 iter1 = iter[];
 iter[] = 0;
 sol2 = solve(prob_ode_brusselator_2d_sparse,
-    KenCarp47(linsolve = KrylovJL_GMRES(), precs = incompletelu,
+    KenCarp47(linsolve = KrylovJL_GMRES(precs = incompletelu),
         concrete_jac = true), save_everystep = false);
 iter2 = iter[];
 iter[] = 0;
 sol3 = solve(prob_ode_brusselator_2d_sparse,
-    KenCarp47(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid,
+    KenCarp47(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid),
         concrete_jac = true), save_everystep = false);
 iter3 = iter[];
 iter[] = 0;
 sol4 = solve(prob_ode_brusselator_2d_sparse,
-    KenCarp47(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid2,
+    KenCarp47(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid2),
         concrete_jac = true), save_everystep = false);
 iter4 = iter[];
 iter[] = 0;
@@ -126,17 +113,17 @@ sol1 = solve(prob_ode_brusselator_2d, Rosenbrock23(linsolve = KrylovJL_GMRES()),
 iter1 = iter[];
 iter[] = 0;
 sol2 = solve(prob_ode_brusselator_2d_sparse,
-    Rosenbrock23(linsolve = KrylovJL_GMRES(), precs = incompletelu,
+    Rosenbrock23(linsolve = KrylovJL_GMRES(precs = incompletelu),
         concrete_jac = true), save_everystep = false);
 iter2 = iter[];
 iter[] = 0;
 sol3 = solve(prob_ode_brusselator_2d_sparse,
-    Rosenbrock23(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid,
+    Rosenbrock23(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid),
         concrete_jac = true), save_everystep = false);
 iter3 = iter[];
 iter[] = 0;
 sol4 = solve(prob_ode_brusselator_2d_sparse,
-    Rosenbrock23(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid2,
+    Rosenbrock23(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid2),
         concrete_jac = true), save_everystep = false);
 iter4 = iter[];
 iter[] = 0;
@@ -155,17 +142,17 @@ sol1 = solve(prob_ode_brusselator_2d, Rodas4(linsolve = KrylovJL_GMRES()),
 iter1 = iter[];
 iter[] = 0;
 sol2 = solve(prob_ode_brusselator_2d_sparse,
-    Rodas4(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true),
+    Rodas4(linsolve = KrylovJL_GMRES(precs = incompletelu), concrete_jac = true),
     save_everystep = false);
 iter2 = iter[];
 iter[] = 0;
 sol3 = solve(prob_ode_brusselator_2d_sparse,
-    Rodas4(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid,
+    Rodas4(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid),
         concrete_jac = true), save_everystep = false);
 iter3 = iter[];
 iter[] = 0;
 sol4 = solve(prob_ode_brusselator_2d_sparse,
-    Rodas4(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid2,
+    Rodas4(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid2),
         concrete_jac = true), save_everystep = false);
 iter4 = iter[];
 iter[] = 0;
@@ -184,17 +171,17 @@ sol1 = solve(prob_ode_brusselator_2d, Rodas5(linsolve = KrylovJL_GMRES()),
 iter1 = iter[];
 iter[] = 0;
 sol2 = solve(prob_ode_brusselator_2d_sparse,
-    Rodas5(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true),
+    Rodas5(linsolve = KrylovJL_GMRES(precs = incompletelu), concrete_jac = true),
     save_everystep = false);
 iter2 = iter[];
 iter[] = 0;
 sol3 = solve(prob_ode_brusselator_2d_sparse,
-    Rodas5(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid,
+    Rodas5(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid),
         concrete_jac = true), save_everystep = false);
 iter3 = iter[];
 iter[] = 0;
 sol4 = solve(prob_ode_brusselator_2d_sparse,
-    Rodas5(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid2,
+    Rodas5(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid2),
         concrete_jac = true), save_everystep = false);
 iter4 = iter[];
 iter[] = 0;
@@ -213,17 +200,17 @@ sol1 = solve(prob_ode_brusselator_2d, TRBDF2(linsolve = KrylovJL_GMRES()),
 iter1 = iter[];
 iter[] = 0;
 sol2 = solve(prob_ode_brusselator_2d_sparse,
-    TRBDF2(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true),
+    TRBDF2(linsolve = KrylovJL_GMRES(precs = incompletelu), concrete_jac = true),
     save_everystep = false);
 iter2 = iter[];
 iter[] = 0;
 sol3 = solve(prob_ode_brusselator_2d_sparse,
-    TRBDF2(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid,
+    TRBDF2(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid),
         concrete_jac = true), save_everystep = false);
 iter3 = iter[];
 iter[] = 0;
 sol4 = solve(prob_ode_brusselator_2d_sparse,
-    TRBDF2(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid2,
+    TRBDF2(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid2),
         concrete_jac = true), save_everystep = false);
 iter4 = iter[];
 iter[] = 0;
@@ -242,17 +229,17 @@ sol1 = solve(prob_ode_brusselator_2d, TRBDF2(linsolve = KrylovJL_GMRES()),
 iter1 = iter[];
 iter[] = 0;
 sol2 = solve(prob_ode_brusselator_2d_sparse,
-    TRBDF2(linsolve = KrylovJL_GMRES(), precs = incompletelu,
+    TRBDF2(linsolve = KrylovJL_GMRES(precs = incompletelu),
         concrete_jac = true), save_everystep = false);
 iter2 = iter[];
 iter[] = 0;
 sol3 = solve(prob_ode_brusselator_2d_sparse,
-    TRBDF2(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid,
+    TRBDF2(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid),
         concrete_jac = true), save_everystep = false);
 iter3 = iter[];
 iter[] = 0;
 sol4 = solve(prob_ode_brusselator_2d_sparse,
-    TRBDF2(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid2,
+    TRBDF2(linsolve = KrylovJL_GMRES(precs = algebraicmultigrid2),
         concrete_jac = true), save_everystep = false);
 iter4 = iter[];
 iter[] = 0;