diff --git a/Project.toml b/Project.toml
index 809ad81d1a..6769d1500d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -121,8 +121,8 @@ OrdinaryDiffEqQPRK = "1"
 OrdinaryDiffEqRKN = "1"
 OrdinaryDiffEqRosenbrock = "1"
 OrdinaryDiffEqSDIRK = "1"
-OrdinaryDiffEqStabilizedIRK = "1"
 OrdinaryDiffEqSSPRK = "1"
+OrdinaryDiffEqStabilizedIRK = "1"
 OrdinaryDiffEqStabilizedRK = "1"
 OrdinaryDiffEqSymplecticRK = "1"
 OrdinaryDiffEqTsit5 = "1"
diff --git a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl
index 5a2981219b..04ccc94195 100644
--- a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl
+++ b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl
@@ -7,8 +7,7 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!,
                            OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
                            OrdinaryDiffEqNewtonAdaptiveAlgorithm,
                            OrdinaryDiffEqNewtonAlgorithm,
-                           AbstractController, DEFAULT_PRECS,
-                           CompiledFloats, uses_uprev,
+                           AbstractController, CompiledFloats, uses_uprev,
                            alg_cache, _vec, _reshape, @cache,
                            isfsal, full_cache,
                            constvalue, isadaptive, error_constant,
diff --git a/lib/OrdinaryDiffEqBDF/src/algorithms.jl b/lib/OrdinaryDiffEqBDF/src/algorithms.jl
index 494bbcfa8a..0ab1955d5b 100644
--- a/lib/OrdinaryDiffEqBDF/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqBDF/src/algorithms.jl
@@ -10,7 +10,6 @@ function BDF_docstring(description::String,
         concrete_jac = nothing,
         diff_type = Val{:forward},
         linsolve = nothing,
-        precs = DEFAULT_PRECS,
         """ * "\n" * extra_keyword_default
 
     keyword_default_description = """
@@ -34,40 +33,6 @@ function BDF_docstring(description::String,
       For example, to use [KLU.jl](https://github.com/JuliaSparse/KLU.jl), specify
       `$name(linsolve = KLUFactorization()`).
        When `nothing` is passed, uses `DefaultLinearSolver`.
-    - `precs`: Any [LinearSolve.jl-compatible preconditioner](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/)
-      can be used as a left or right preconditioner.
-      Preconditioners are specified by the `Pl,Pr = precs(W,du,u,p,t,newW,Plprev,Prprev,solverdata)`
-      function where the arguments are defined as:
-        - `W`: the current Jacobian of the nonlinear system. Specified as either
-            ``I - \\gamma J`` or ``I/\\gamma - J`` depending on the algorithm. This will
-            commonly be a `WOperator` type defined by OrdinaryDiffEq.jl. It is a lazy
-            representation of the operator. Users can construct the W-matrix on demand
-            by calling `convert(AbstractMatrix,W)` to receive an `AbstractMatrix` matching
-            the `jac_prototype`.
-        - `du`: the current ODE derivative
-        - `u`: the current ODE state
-        - `p`: the ODE parameters
-        - `t`: the current ODE time
-        - `newW`: a `Bool` which specifies whether the `W` matrix has been updated since
-            the last call to `precs`. It is recommended that this is checked to only
-            update the preconditioner when `newW == true`.
-        - `Plprev`: the previous `Pl`.
-        - `Prprev`: the previous `Pr`.
-        - `solverdata`: Optional extra data the solvers can give to the `precs` function.
-            Solver-dependent and subject to change.
-      The return is a tuple `(Pl,Pr)` of the LinearSolve.jl-compatible preconditioners.
-      To specify one-sided preconditioning, simply return `nothing` for the preconditioner
-      which is not used. Additionally, `precs` must supply the dispatch:
-      ```julia
-      Pl, Pr = precs(W, du, u, p, t, ::Nothing, ::Nothing, ::Nothing, solverdata)
-      ```
-      which is used in the solver setup phase to construct the integrator
-      type with the preconditioners `(Pl,Pr)`.
-      The default is `precs=DEFAULT_PRECS` where the default preconditioner function
-      is defined as:
-      ```julia
-      DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
-      ```
         """ * "/n" * extra_keyword_description
     generic_solver_docstring(
         description, name, "Multistep Method.", references,
@@ -101,11 +66,10 @@ end
     controller = :Standard,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -115,14 +79,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
 
@@ -157,11 +121,10 @@ like `KenCarp4`, but instead using a multistep BDF approach",
     ark = false,
     order,
     """)
-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
@@ -171,14 +134,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,
@@ -189,15 +151,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,
@@ -278,11 +239,10 @@ Optional parameter kappa defaults to Shampine's accuracy-optimal -0.1850.",
     controller = :Standard,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -291,15 +251,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,
@@ -333,11 +292,10 @@ end
     controller = :Standard,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -346,15 +304,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,
@@ -393,12 +350,11 @@ Utilizes Shampine's accuracy-optimal kappa values as defaults (has a keyword arg
     controller = :Standard,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -410,16 +366,16 @@ 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
 
@@ -446,22 +402,20 @@ TruncatedStacktraces.@truncate_stacktrace QNDF
     nlsolve = NLNewton(),
     extrapolant = :constant,
     """)
-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
 
@@ -493,12 +447,11 @@ Utilizes Shampine's accuracy-optimal kappa values as defaults (has a keyword arg
     step_limiter! = trivial_limiter!,
     max_order::Val{MO} = Val{5}(),
     """)
-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
@@ -509,14 +462,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
 
@@ -614,23 +567,22 @@ It uses an apriori error estimator for adaptivity based on a finite differencing
     extrapolant = :constant,
     controller = :Standard,
     """)
-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
 
 @doc BDF_docstring("Fully implicit implementation of BDF2.",
@@ -653,36 +605,34 @@ end
     extrapolant = :constant,
     controller = :Standard,
     """)
-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
 
 #=
-struct DBDF{CS,AD,F,F2,P,FDT,ST,CJ} <: DAEAlgorithm{CS,AD,FDT,ST,CJ}
+struct DBDF{CS,AD,F,F2,FDT,ST,CJ} <: DAEAlgorithm{CS,AD,FDT,ST,CJ}
   linsolve::F
   nlsolve::F2
-  precs::P
   extrapolant::Symbol
 end
 
 DBDF(;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) =
-     DBDF{_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)
+     linsolve=nothing,nlsolve=NLNewton(),extrapolant=:linear) =
+     DBDF{_unwrap_val(chunk_size),_unwrap_val(autodiff),typeof(linsolve),typeof(nlsolve),diff_type,_unwrap_val(standardtag),_unwrap_val(concrete_jac)}(
+     linsolve,nlsolve,extrapolant)
 =#
 
 @doc BDF_docstring("Fully implicit implementation of FBDF based on Shampine's",
@@ -709,11 +659,10 @@ DBDF(;chunk_size=Val{0}(),autodiff=Val{true}(), standardtag = Val{true}(), concr
     controller = :Standard,
     max_order::Val{MO} = Val{5}(),
     """)
-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
@@ -722,13 +671,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 bc36e5c7a0..929571c11f 100644
--- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
+++ b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
@@ -105,10 +105,10 @@ end
     Divergence = -2
 end
 const TryAgain = SlowConvergence
-
-DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
 isdiscretecache(cache) = false
 
+# unused. Delete this once StocasticDiffEq doesn't use it
+DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
 include("doc_utils.jl")
 include("misc_utils.jl")
 
diff --git a/lib/OrdinaryDiffEqCore/src/doc_utils.jl b/lib/OrdinaryDiffEqCore/src/doc_utils.jl
index 75eb00dc7d..a75221fcc8 100644
--- a/lib/OrdinaryDiffEqCore/src/doc_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/doc_utils.jl
@@ -87,7 +87,6 @@ function differentiation_rk_docstring(description::String,
     concrete_jac = nothing,
     diff_type = Val{:forward},
     linsolve = nothing,
-    precs = DEFAULT_PRECS,
     """ * extra_keyword_default
 
     keyword_default_description = """
@@ -111,40 +110,7 @@ function differentiation_rk_docstring(description::String,
       For example, to use [KLU.jl](https://github.com/JuliaSparse/KLU.jl), specify
       `$name(linsolve = KLUFactorization()`).
        When `nothing` is passed, uses `DefaultLinearSolver`.
-    - `precs`: Any [LinearSolve.jl-compatible preconditioner](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/)
-      can be used as a left or right preconditioner.
-      Preconditioners are specified by the `Pl,Pr = precs(W,du,u,p,t,newW,Plprev,Prprev,solverdata)`
-      function where the arguments are defined as:
-        - `W`: the current Jacobian of the nonlinear system. Specified as either
-            ``I - \\gamma J`` or ``I/\\gamma - J`` depending on the algorithm. This will
-            commonly be a `WOperator` type defined by OrdinaryDiffEq.jl. It is a lazy
-            representation of the operator. Users can construct the W-matrix on demand
-            by calling `convert(AbstractMatrix,W)` to receive an `AbstractMatrix` matching
-            the `jac_prototype`.
-        - `du`: the current ODE derivative
-        - `u`: the current ODE state
-        - `p`: the ODE parameters
-        - `t`: the current ODE time
-        - `newW`: a `Bool` which specifies whether the `W` matrix has been updated since
-            the last call to `precs`. It is recommended that this is checked to only
-            update the preconditioner when `newW == true`.
-        - `Plprev`: the previous `Pl`.
-        - `Prprev`: the previous `Pr`.
-        - `solverdata`: Optional extra data the solvers can give to the `precs` function.
-            Solver-dependent and subject to change.
-      The return is a tuple `(Pl,Pr)` of the LinearSolve.jl-compatible preconditioners.
-      To specify one-sided preconditioning, simply return `nothing` for the preconditioner
-      which is not used. Additionally, `precs` must supply the dispatch:
-      ```julia
-      Pl, Pr = precs(W, du, u, p, t, ::Nothing, ::Nothing, ::Nothing, solverdata)
-      ```
-      which is used in the solver setup phase to construct the integrator
-      type with the preconditioners `(Pl,Pr)`.
-      The default is `precs=DEFAULT_PRECS` where the default preconditioner function
-      is defined as:
-      ```julia
-      DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
-      ```
+
     """ * extra_keyword_description
 
     generic_solver_docstring(
diff --git a/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl b/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl
index 4191dbf4d8..a33ff38269 100644
--- a/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl
+++ b/lib/OrdinaryDiffEqDefault/test/default_solver_tests.jl
@@ -49,7 +49,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]
@@ -75,7 +75,7 @@ for n in (100, 600)
         vcat([1.0, 0.0, 0.0], ones(n)), (0.0, 100.0), (0.04, 3e7, 1e4))
     global 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/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/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
index e351ab54c3..b9e3695b69 100644
--- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
+++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
@@ -883,7 +883,7 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits},
         # Thus setup JacVec and a concrete J, using sparsity when possible
         _f = islin ? (isode ? f.f : f.f1.f) : f
         J = if f.jac_prototype === nothing
-            ArrayInterface.undefmatrix(u)
+            ArrayInterface.zeromatrix(u)
         else
             deepcopy(f.jac_prototype)
         end
@@ -907,7 +907,7 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits},
                 f.jac(uprev, p, t)
             end
         elseif f.jac_prototype === nothing
-            ArrayInterface.undefmatrix(u)
+            ArrayInterface.zeromatrix(u)
         else
             deepcopy(f.jac_prototype)
         end
@@ -1003,4 +1003,4 @@ function resize_J_W!(cache, integrator, i)
     end
 
     nothing
-end
\ No newline at end of file
+end
diff --git a/lib/OrdinaryDiffEqDifferentiation/src/linsolve_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/linsolve_utils.jl
index aa48161e1b..82482a530a 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,15 @@ 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 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/OrdinaryDiffEqExtrapolation.jl b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl
index 29bbc6ca8b..6b5726a71a 100644
--- a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl
+++ b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl
@@ -13,7 +13,7 @@ import OrdinaryDiffEqCore: alg_order, alg_maximum_order, get_current_adaptive_or
                            OrdinaryDiffEqAdaptiveAlgorithm,
                            OrdinaryDiffEqAdaptiveImplicitAlgorithm,
                            alg_cache, CompiledFloats, @threaded, stepsize_controller!,
-                           DEFAULT_PRECS, full_cache,
+                           full_cache,
                            constvalue, PolyesterThreads, Sequential, BaseThreads,
                            _digest_beta1_beta2, timedepentdtmin, _unwrap_val,
                            _reshape, _vec, get_fsalfirstlast, generic_solver_docstring,
diff --git a/lib/OrdinaryDiffEqExtrapolation/src/algorithms.jl b/lib/OrdinaryDiffEqExtrapolation/src/algorithms.jl
index 205e8e1f3d..266287c234 100644
--- a/lib/OrdinaryDiffEqExtrapolation/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqExtrapolation/src/algorithms.jl
@@ -55,10 +55,9 @@ Similar to Hairer's SEULEX.",
     thread = OrdinaryDiffEq.False(),
     sequence = :harmonic
     """)
-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
@@ -69,7 +68,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 &&
@@ -100,11 +98,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
 
 @doc generic_solver_docstring("Midpoint extrapolation using Barycentric coordinates.",
@@ -195,10 +192,9 @@ end
     thread = OrdinaryDiffEq.False(),
     sequence = :harmonic,
     """)
-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
@@ -207,8 +203,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:
@@ -243,9 +238,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
@@ -342,10 +337,9 @@ end
     thread = OrdinaryDiffEq.False(),
     sequence = :harmonic,
     """)
-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
@@ -356,7 +350,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)
@@ -392,11 +386,10 @@ 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
 
 @doc differentiation_rk_docstring("Euler extrapolation using Barycentric coordinates,
@@ -420,10 +413,9 @@ end
     sequence = :harmonic,
     sequence_factor = 2,
     """)
-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
@@ -436,7 +428,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,
@@ -473,10 +465,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 30c8fc16c9..168318e928 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/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
index af3e4a7dd2..bc6bf9e12a 100644
--- a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
@@ -10,7 +10,7 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!,
                            constvalue, _unwrap_val,
                            differentiation_rk_docstring, trivial_limiter!,
                            _ode_interpolant!, _ode_addsteps!, AbstractController,
-                           qmax_default, alg_adaptive_order, DEFAULT_PRECS,
+                           qmax_default, alg_adaptive_order,
                            stepsize_controller!, step_accept_controller!,
                            step_reject_controller!,
                            PredictiveController, alg_can_repeat_jac, NewtonAlgorithm,
diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
index 844f50ea96..ab9e2baab7 100644
--- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
@@ -1,4 +1,3 @@
-
 hairer1999stiff = """@article{hairer1999stiff,
     title={Stiff differential equations solved by Radau methods},
     author={Hairer, Ernst and Wanner, Gerhard},
@@ -26,10 +25,9 @@ Similar to Hairer's SEULEX.",
     references = hairer1999stiff,
     extra_keyword_description = extra_keyword_description,
     extra_keyword_default = extra_keyword_default)
-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
@@ -42,16 +40,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,
@@ -69,10 +66,9 @@ Similar to Hairer's SEULEX.",
     references = hairer1999stiff,
     extra_keyword_description = extra_keyword_description,
     extra_keyword_default = extra_keyword_default)
-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
@@ -86,16 +82,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,
         κ,
@@ -114,10 +109,9 @@ Similar to Hairer's SEULEX.",
     references = hairer1999stiff,
     extra_keyword_description = extra_keyword_description,
     extra_keyword_default = extra_keyword_default)
-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
@@ -131,16 +125,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 47fdad5374..ec178ee765 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/OrdinaryDiffEqIMEXMultistep.jl b/lib/OrdinaryDiffEqIMEXMultistep/src/OrdinaryDiffEqIMEXMultistep.jl
index 567e6b5eae..6771de9751 100644
--- a/lib/OrdinaryDiffEqIMEXMultistep/src/OrdinaryDiffEqIMEXMultistep.jl
+++ b/lib/OrdinaryDiffEqIMEXMultistep/src/OrdinaryDiffEqIMEXMultistep.jl
@@ -1,8 +1,7 @@
 module OrdinaryDiffEqIMEXMultistep
 
 import OrdinaryDiffEqCore: alg_order, issplit, OrdinaryDiffEqNewtonAlgorithm, _unwrap_val,
-                           DEFAULT_PRECS, OrdinaryDiffEqConstantCache,
-                           OrdinaryDiffEqMutableCache,
+                           OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache,
                            @cache, alg_cache, initialize!, perform_step!, @unpack,
                            full_cache, get_fsalfirstlast,
                            generic_solver_docstring
diff --git a/lib/OrdinaryDiffEqIMEXMultistep/src/algorithms.jl b/lib/OrdinaryDiffEqIMEXMultistep/src/algorithms.jl
index 69e0d7ab57..75272c874b 100644
--- a/lib/OrdinaryDiffEqIMEXMultistep/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqIMEXMultistep/src/algorithms.jl
@@ -20,24 +20,21 @@
     pages={647--659},
     year={2010},
     publisher={Wiley Online Library}}", "", "")
-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
 
@@ -60,22 +57,20 @@ end
     pages={263--276},
     year={2015},
     publisher={Elsevier}}", "", "")
-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/newton.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl
index c1c8f67db7..20c92ae5aa 100644
--- a/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl
+++ b/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl
@@ -220,15 +220,9 @@ end
         reltol = eps(eltype(dz))
     end
 
-    if is_always_new(nlsolver) || (iter == 1 && new_W)
-        linres = dolinsolve(integrator, linsolve; A = W, b = _vec(b), linu = _vec(dz),
-            reltol = reltol)
-    else
-        linres = dolinsolve(
-            integrator, linsolve; A = nothing, b = _vec(b), linu = _vec(dz),
-            reltol = reltol)
-    end
-
+    make_new_W = is_always_new(nlsolver) || (iter == 1 && new_W)
+    linres = dolinsolve(integrator, linsolve; A = make_new_W ? W : nothing, b = _vec(b),
+        linu = dz, reltol)
     cache.linsolve = linres.cache
 
     if DiffEqBase.has_stats(integrator)
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/OrdinaryDiffEqPDIRK.jl b/lib/OrdinaryDiffEqPDIRK/src/OrdinaryDiffEqPDIRK.jl
index 18a4e56205..5248348b0b 100644
--- a/lib/OrdinaryDiffEqPDIRK/src/OrdinaryDiffEqPDIRK.jl
+++ b/lib/OrdinaryDiffEqPDIRK/src/OrdinaryDiffEqPDIRK.jl
@@ -3,7 +3,7 @@ module OrdinaryDiffEqPDIRK
 import OrdinaryDiffEqCore: isfsal, alg_order, _unwrap_val,
                            OrdinaryDiffEqNewtonAlgorithm, OrdinaryDiffEqConstantCache,
                            OrdinaryDiffEqMutableCache, constvalue, alg_cache,
-                           uses_uprev, @unpack, unwrap_alg, @cache, DEFAULT_PRECS,
+                           uses_uprev, @unpack, unwrap_alg, @cache,
                            @threaded, initialize!, perform_step!, isthreaded,
                            full_cache, get_fsalfirstlast, differentiation_rk_docstring
 import StaticArrays: SVector
diff --git a/lib/OrdinaryDiffEqPDIRK/src/algorithms.jl b/lib/OrdinaryDiffEqPDIRK/src/algorithms.jl
index 0e93763cfd..fa781ce0e8 100644
--- a/lib/OrdinaryDiffEqPDIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqPDIRK/src/algorithms.jl
@@ -21,20 +21,19 @@
     extrapolant = :constant,
     thread = OrdinaryDiffEq.True(),
     """)
-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/OrdinaryDiffEqRosenbrock.jl b/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl
index ed5502cedd..491f61e869 100644
--- a/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl
+++ b/lib/OrdinaryDiffEqRosenbrock/src/OrdinaryDiffEqRosenbrock.jl
@@ -1,7 +1,7 @@
 module OrdinaryDiffEqRosenbrock
 
 import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, isWmethod, isfsal, _unwrap_val,
-                           DEFAULT_PRECS, OrdinaryDiffEqRosenbrockAlgorithm, @cache,
+                           OrdinaryDiffEqRosenbrockAlgorithm, @cache,
                            alg_cache, initialize!, @unpack,
                            calculate_residuals!, OrdinaryDiffEqMutableCache,
                            OrdinaryDiffEqConstantCache, _ode_interpolant, _ode_interpolant!,
@@ -52,7 +52,6 @@ function rosenbrock_wanner_docstring(description::String,
     concrete_jac = nothing,
     diff_type = Val{:central},
     linsolve = nothing,
-    precs = DEFAULT_PRECS,
     """ * extra_keyword_default
 
     keyword_default_description = """
@@ -76,41 +75,7 @@ function rosenbrock_wanner_docstring(description::String,
       For example, to use [KLU.jl](https://github.com/JuliaSparse/KLU.jl), specify
       `$name(linsolve = KLUFactorization()`).
        When `nothing` is passed, uses `DefaultLinearSolver`.
-    - `precs`: Any [LinearSolve.jl-compatible preconditioner](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/)
-      can be used as a left or right preconditioner.
-      Preconditioners are specified by the `Pl,Pr = precs(W,du,u,p,t,newW,Plprev,Prprev,solverdata)`
-      function where the arguments are defined as:
-        - `W`: the current Jacobian of the nonlinear system. Specified as either
-            ``I - \\gamma J`` or ``I/\\gamma - J`` depending on the algorithm. This will
-            commonly be a `WOperator` type defined by OrdinaryDiffEq.jl. It is a lazy
-            representation of the operator. Users can construct the W-matrix on demand
-            by calling `convert(AbstractMatrix,W)` to receive an `AbstractMatrix` matching
-            the `jac_prototype`.
-        - `du`: the current ODE derivative
-        - `u`: the current ODE state
-        - `p`: the ODE parameters
-        - `t`: the current ODE time
-        - `newW`: a `Bool` which specifies whether the `W` matrix has been updated since
-            the last call to `precs`. It is recommended that this is checked to only
-            update the preconditioner when `newW == true`.
-        - `Plprev`: the previous `Pl`.
-        - `Prprev`: the previous `Pr`.
-        - `solverdata`: Optional extra data the solvers can give to the `precs` function.
-            Solver-dependent and subject to change.
-      The return is a tuple `(Pl,Pr)` of the LinearSolve.jl-compatible preconditioners.
-      To specify one-sided preconditioning, simply return `nothing` for the preconditioner
-      which is not used. Additionally, `precs` must supply the dispatch:
-      ```julia
-      Pl, Pr = precs(W, du, u, p, t, ::Nothing, ::Nothing, ::Nothing, solverdata)
-      ```
-      which is used in the solver setup phase to construct the integrator
-      type with the preconditioners `(Pl,Pr)`.
-      The default is `precs=DEFAULT_PRECS` where the default preconditioner function
-      is defined as:
-      ```julia
-      DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
-      ```
-    """ * extra_keyword_description
+     """ * extra_keyword_description
 
     if with_step_limiter
         keyword_default *= "step_limiter! = OrdinaryDiffEq.trivial_limiter!,\n"
@@ -150,41 +115,7 @@ function rosenbrock_docstring(description::String,
       For example, to use [KLU.jl](https://github.com/JuliaSparse/KLU.jl), specify
       `$name(linsolve = KLUFactorization()`).
        When `nothing` is passed, uses `DefaultLinearSolver`.
-    - `precs`: Any [LinearSolve.jl-compatible preconditioner](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/)
-      can be used as a left or right preconditioner.
-      Preconditioners are specified by the `Pl,Pr = precs(W,du,u,p,t,newW,Plprev,Prprev,solverdata)`
-      function where the arguments are defined as:
-        - `W`: the current Jacobian of the nonlinear system. Specified as either
-            ``I - \\gamma J`` or ``I/\\gamma - J`` depending on the algorithm. This will
-            commonly be a `WOperator` type defined by OrdinaryDiffEq.jl. It is a lazy
-            representation of the operator. Users can construct the W-matrix on demand
-            by calling `convert(AbstractMatrix,W)` to receive an `AbstractMatrix` matching
-            the `jac_prototype`.
-        - `du`: the current ODE derivative
-        - `u`: the current ODE state
-        - `p`: the ODE parameters
-        - `t`: the current ODE time
-        - `newW`: a `Bool` which specifies whether the `W` matrix has been updated since
-            the last call to `precs`. It is recommended that this is checked to only
-            update the preconditioner when `newW == true`.
-        - `Plprev`: the previous `Pl`.
-        - `Prprev`: the previous `Pr`.
-        - `solverdata`: Optional extra data the solvers can give to the `precs` function.
-            Solver-dependent and subject to change.
-      The return is a tuple `(Pl,Pr)` of the LinearSolve.jl-compatible preconditioners.
-      To specify one-sided preconditioning, simply return `nothing` for the preconditioner
-      which is not used. Additionally, `precs` must supply the dispatch:
-      ```julia
-      Pl, Pr = precs(W, du, u, p, t, ::Nothing, ::Nothing, ::Nothing, solverdata)
-      ```
-      which is used in the solver setup phase to construct the integrator
-      type with the preconditioners `(Pl,Pr)`.
-      The default is `precs=DEFAULT_PRECS` where the default preconditioner function
-      is defined as:
-      ```julia
-      DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
-      ```
-    """ * extra_keyword_default
+     """ * extra_keyword_default
 
     keyword_default_description = """
     - `chunk_size`: TBD
@@ -193,7 +124,6 @@ function rosenbrock_docstring(description::String,
     - `concrete_jac`: function of the form `jac!(J, u, p, t)`
     - `diff_type`: TBD
     - `linsolve`: custom solver for the inner linear systems
-    - `precs`: custom preconditioner for the inner linear solver
     """ * extra_keyword_description
 
     if with_step_limiter
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 a67b4a9548..cda6b4e944 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)))
diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl
index 4b6841a1bc..bfce56d6b3 100644
--- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl
+++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl
@@ -137,12 +137,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)
@@ -182,13 +178,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))
@@ -331,12 +323,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)
@@ -417,12 +405,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)
@@ -610,12 +594,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)
@@ -654,13 +634,9 @@ 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,
-        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,
+         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)
     Rodas3PCache(u, uprev, dense1, dense2, dense3, du, du1, du2, k1, k2, k3, k4, k5,
@@ -761,16 +737,10 @@ function alg_cache(alg::Union{Rodas4, Rodas42, Rodas4P, Rodas4P2}, u, rate_proto
     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,
+         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)
 
@@ -868,12 +838,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)
@@ -932,13 +898,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 51c938ee39..03a8c4c443 100644
--- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl
+++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl
@@ -48,21 +48,10 @@ 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
-
-    vecu = _vec(linres.u)
+    linres = dolinsolve(integrator, cache.linsolve; A = repeat_step ? nothing : W, b = _vec(linsolve_tmp))
     veck₁ = _vec(k₁)
 
-    @.. broadcast=false veck₁=-vecu
+    @.. broadcast=false veck₁=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false u=uprev + dto2 * k₁
@@ -79,10 +68,9 @@ end
     @.. broadcast=false linsolve_tmp=f₁ - tmp
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
-    vecu = _vec(linres.u)
     veck2 = _vec(k₂)
 
-    @.. broadcast=false veck2=-vecu
+    @.. broadcast=false veck2=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false k₂+=k₁
@@ -105,9 +93,8 @@ end
         end
 
         linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
-        vecu = _vec(linres.u)
         veck3 = _vec(k₃)
-        @.. broadcast=false veck3=-vecu
+        @.. broadcast=false veck3=-linres.u
 
         integrator.stats.nsolve += 1
 
@@ -160,21 +147,10 @@ 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₁)
-
-    @.. broadcast=false veck₁=-vecu
+    @.. broadcast=false veck₁=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false u=uprev + dto2 * k₁
@@ -191,10 +167,8 @@ end
     @.. broadcast=false linsolve_tmp=f₁ - tmp
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
-    vecu = _vec(linres.u)
     veck2 = _vec(k₂)
-
-    @.. broadcast=false veck2=-vecu
+    @.. broadcast=false veck2=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false k₂+=k₁
@@ -213,10 +187,9 @@ end
     end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
-    vecu = _vec(linres.u)
     veck3 = _vec(k₃)
 
-    @.. broadcast=false veck3=-vecu
+    @.. broadcast=false veck3=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false u=uprev + dto6 * (k₁ + 4k₂ + k₃)
@@ -517,21 +490,9 @@ 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
-
-    vecu = _vec(linres.u)
+    linres = dolinsolve(integrator, cache.linsolve; A = repeat_step ? nothing : W, b = _vec(linsolve_tmp))
     veck1 = _vec(k1)
-
-    @.. broadcast=false veck1=-vecu
+    @.. broadcast=false veck1=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false u=uprev + a21 * k1
@@ -548,10 +509,9 @@ end
     end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
-    vecu = _vec(linres.u)
     veck2 = _vec(k2)
 
-    @.. broadcast=false veck2=-vecu
+    @.. broadcast=false veck2=-linres.u
 
     integrator.stats.nsolve += 1
 
@@ -569,10 +529,9 @@ end
     end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
-    vecu = _vec(linres.u)
     veck3 = _vec(k3)
 
-    @.. broadcast=false veck3=-vecu
+    @.. broadcast=false veck3=-linres.u
 
     integrator.stats.nsolve += 1
 
@@ -712,21 +671,10 @@ 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
-
-    vecu = _vec(linres.u)
+    linres = dolinsolve(integrator, cache.linsolve; A = repeat_step ? nothing : W, b = _vec(linsolve_tmp))
     veck1 = _vec(k1)
 
-    @.. broadcast=false veck1=-vecu
+    @.. broadcast=false veck1=-linres.u
     integrator.stats.nsolve += 1
 
     #=
@@ -747,7 +695,7 @@ end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
     veck2 = _vec(k2)
-    @.. broadcast=false veck2=-vecu
+    @.. broadcast=false veck2=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false u=uprev + a31 * k1 + a32 * k2
@@ -765,7 +713,7 @@ end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
     veck3 = _vec(k3)
-    @.. broadcast=false veck3=-vecu
+    @.. broadcast=false veck3=-linres.u
     integrator.stats.nsolve += 1
     @.. broadcast=false u=uprev + a41 * k1 + a42 * k2 + a43 * k3
     stage_limiter!(u, integrator, p, t + dt)
@@ -783,7 +731,7 @@ end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
     veck4 = _vec(k4)
-    @.. broadcast=false veck4=-vecu
+    @.. broadcast=false veck4=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false u=uprev + b1 * k1 + b2 * k2 + b3 * k3 + b4 * k4
@@ -1020,16 +968,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
 
@@ -1320,16 +1259,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))
 
     @.. $(_vec(ks[1]))=-linres.u
     integrator.stats.nsolve += 1
@@ -1665,21 +1595,11 @@ 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)
 
-    @.. broadcast=false veck1=-vecu
+    @.. broadcast=false veck1=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false u=uprev + a21 * k1
@@ -1697,7 +1617,7 @@ end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
     veck2 = _vec(k2)
-    @.. broadcast=false veck2=-vecu
+    @.. broadcast=false veck2=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false u=uprev + a31 * k1 + a32 * k2
@@ -1715,7 +1635,7 @@ end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
     veck3 = _vec(k3)
-    @.. broadcast=false veck3=-vecu
+    @.. broadcast=false veck3=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false u=uprev + a41 * k1 + a42 * k2 + a43 * k3
@@ -1734,7 +1654,7 @@ end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
     veck4 = _vec(k4)
-    @.. broadcast=false veck4=-vecu
+    @.. broadcast=false veck4=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false u=uprev + a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4
@@ -1753,7 +1673,7 @@ end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
     veck5 = _vec(k5)
-    @.. broadcast=false veck5=-vecu
+    @.. broadcast=false veck5=-linres.u
     integrator.stats.nsolve += 1
 
     @.. broadcast=false u=uprev + a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5
@@ -1773,7 +1693,7 @@ end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
     veck6 = _vec(k6)
-    @.. broadcast=false veck6=-vecu
+    @.. broadcast=false veck6=-linres.u
     integrator.stats.nsolve += 1
 
     u .+= k6
@@ -1793,7 +1713,7 @@ end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
     veck7 = _vec(k7)
-    @.. broadcast=false veck7=-vecu
+    @.. broadcast=false veck7=-linres.u
     integrator.stats.nsolve += 1
 
     u .+= k7
@@ -1812,7 +1732,7 @@ end
 
     linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
     veck8 = _vec(k8)
-    @.. broadcast=false veck8=-vecu
+    @.. broadcast=false veck8=-linres.u
     integrator.stats.nsolve += 1
 
     du .= k8
diff --git a/lib/OrdinaryDiffEqSDIRK/src/OrdinaryDiffEqSDIRK.jl b/lib/OrdinaryDiffEqSDIRK/src/OrdinaryDiffEqSDIRK.jl
index f9fc3bbd9e..201ee7c1db 100644
--- a/lib/OrdinaryDiffEqSDIRK/src/OrdinaryDiffEqSDIRK.jl
+++ b/lib/OrdinaryDiffEqSDIRK/src/OrdinaryDiffEqSDIRK.jl
@@ -7,7 +7,6 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!,
                            OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
                            OrdinaryDiffEqNewtonAdaptiveAlgorithm,
                            OrdinaryDiffEqNewtonAlgorithm,
-                           DEFAULT_PRECS,
                            OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev,
                            alg_cache, _vec, _reshape, @cache, isfsal, full_cache,
                            constvalue, _unwrap_val, _ode_interpolant,
diff --git a/lib/OrdinaryDiffEqSDIRK/src/algorithms.jl b/lib/OrdinaryDiffEqSDIRK/src/algorithms.jl
index 226f92e57f..bcd31c410d 100644
--- a/lib/OrdinaryDiffEqSDIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqSDIRK/src/algorithms.jl
@@ -10,7 +10,6 @@ function SDIRK_docstring(description::String,
         concrete_jac = nothing,
         diff_type = Val{:forward},
         linsolve = nothing,
-        precs = DEFAULT_PRECS,
         nlsolve = NLNewton(),
         """ * extra_keyword_default
 
@@ -35,41 +34,7 @@ function SDIRK_docstring(description::String,
       For example, to use [KLU.jl](https://github.com/JuliaSparse/KLU.jl), specify
       `$name(linsolve = KLUFactorization()`).
        When `nothing` is passed, uses `DefaultLinearSolver`.
-    - `precs`: Any [LinearSolve.jl-compatible preconditioner](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/)
-      can be used as a left or right preconditioner.
-      Preconditioners are specified by the `Pl,Pr = precs(W,du,u,p,t,newW,Plprev,Prprev,solverdata)`
-      function where the arguments are defined as:
-        - `W`: the current Jacobian of the nonlinear system. Specified as either
-            ``I - \\gamma J`` or ``I/\\gamma - J`` depending on the algorithm. This will
-            commonly be a `WOperator` type defined by OrdinaryDiffEq.jl. It is a lazy
-            representation of the operator. Users can construct the W-matrix on demand
-            by calling `convert(AbstractMatrix,W)` to receive an `AbstractMatrix` matching
-            the `jac_prototype`.
-        - `du`: the current ODE derivative
-        - `u`: the current ODE state
-        - `p`: the ODE parameters
-        - `t`: the current ODE time
-        - `newW`: a `Bool` which specifies whether the `W` matrix has been updated since
-            the last call to `precs`. It is recommended that this is checked to only
-            update the preconditioner when `newW == true`.
-        - `Plprev`: the previous `Pl`.
-        - `Prprev`: the previous `Pr`.
-        - `solverdata`: Optional extra data the solvers can give to the `precs` function.
-            Solver-dependent and subject to change.
-      The return is a tuple `(Pl,Pr)` of the LinearSolve.jl-compatible preconditioners.
-      To specify one-sided preconditioning, simply return `nothing` for the preconditioner
-      which is not used. Additionally, `precs` must supply the dispatch:
-      ```julia
-      Pl, Pr = precs(W, du, u, p, t, ::Nothing, ::Nothing, ::Nothing, solverdata)
-      ```
-      which is used in the solver setup phase to construct the integrator
-      type with the preconditioners `(Pl,Pr)`.
-      The default is `precs=DEFAULT_PRECS` where the default preconditioner function
-      is defined as:
-      ```julia
-      DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
-      ```
-    - `nlsolve`: TBD
+     - `nlsolve`: TBD
         """ * extra_keyword_description
 
     generic_solver_docstring(
@@ -98,11 +63,10 @@ end
     controller = :PI,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -111,13 +75,13 @@ 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
 
 @doc SDIRK_docstring("A second order A-stable symplectic and symmetric implicit solver.
@@ -137,11 +101,10 @@ end
     extrapolant = :linear,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -149,13 +112,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
@@ -178,11 +140,10 @@ estimate (the SPICE approximation strategy).""",
     controller = :PI,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -191,14 +152,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!)
@@ -229,11 +189,10 @@ end
     controller = :PI,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -242,12 +201,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
 
@@ -276,11 +235,10 @@ TruncatedStacktraces.@truncate_stacktrace TRBDF2
     controller = :PI,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -289,13 +247,13 @@ 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
@@ -318,11 +276,10 @@ end
     controller = :PI,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -331,14 +288,13 @@ 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!)
@@ -370,11 +326,10 @@ but where reaction equations are sufficient stiff to justify implicit integratio
     extrapolant = :constant,
     controller = :PI,
     """)
-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
@@ -383,12 +338,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
 
@@ -415,11 +370,10 @@ end
     controller = :PI,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -428,12 +382,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
 
@@ -457,11 +411,10 @@ end
     controller = :PI,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -470,12 +423,12 @@ 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
 
@@ -496,23 +449,21 @@ end
     extra_keyword_default = """
     extrapolant = :linear,
     """)
-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
 
@@ -539,11 +490,10 @@ end
     controller = :PI,
     embedding = 3,
     """)
-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
@@ -551,15 +501,14 @@ 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,
@@ -583,23 +532,21 @@ end
     extra_keyword_default = """
     extrapolant = :linear,
     """)
-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
 
@@ -620,24 +567,22 @@ end
     extra_keyword_default = """
     extrapolant = :linear,
     """)
-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
 
@@ -658,24 +603,22 @@ end
     extra_keyword_default = """
     extrapolant = :linear,
     """)
-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
 
@@ -696,24 +639,22 @@ end
     extra_keyword_default = """
     extrapolant = :linear,
     """)
-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
 
@@ -734,24 +675,22 @@ end
     extra_keyword_default = """
     extrapolant = :linear,
     """)
-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
 
@@ -770,11 +709,10 @@ end
     extrapolant = :linear,
     controller = :PI,
     """)
-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
@@ -782,12 +720,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
 
@@ -806,11 +744,10 @@ end
     extrapolant = :linear,
     controller = :PI,
     """)
-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
@@ -818,12 +755,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
 
@@ -850,11 +787,10 @@ end
     controller = :PI,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -863,12 +799,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
 
@@ -895,11 +831,10 @@ end
     controller = :PI,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -908,12 +843,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
 
@@ -937,11 +872,10 @@ end
     controller = :PI,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -950,12 +884,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
 
@@ -982,11 +916,10 @@ TruncatedStacktraces.@truncate_stacktrace KenCarp4
     extrapolant = :linear,
     controller = :PI,
     """)
-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
@@ -994,12 +927,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
 
@@ -1023,11 +956,10 @@ end
     controller = :PI,
     step_limiter! = trivial_limiter!,
     """)
-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
@@ -1036,12 +968,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
 
@@ -1066,11 +998,10 @@ end
     extrapolant = :linear,
     controller = :PI,
     """)
-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
@@ -1078,12 +1009,12 @@ 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
 
@@ -1109,22 +1040,21 @@ but are still being fully evaluated in context.",
     extrapolant = :linear,
     controller = :PI,
     """)
-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
 
@@ -1149,22 +1079,21 @@ but are still being fully evaluated in context.",
     extrapolant = :linear,
     controller = :PI,
     """)
-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
 
@@ -1189,22 +1118,21 @@ but are still being fully evaluated in context.",
     extrapolant = :linear,
     controller = :PI,
     """)
-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
 
@@ -1229,22 +1157,21 @@ but are still being fully evaluated in context.",
     extrapolant = :linear,
     controller = :PI,
     """)
-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
 
@@ -1271,21 +1198,20 @@ Check issue https://github.com/SciML/OrdinaryDiffEq.jl/issues/1933 for more deta
     extrapolant = :linear,
     controller = :PI,
     """)
-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/OrdinaryDiffEqStabilizedIRK.jl b/lib/OrdinaryDiffEqStabilizedIRK/src/OrdinaryDiffEqStabilizedIRK.jl
index f16310de7e..6041a3aee8 100644
--- a/lib/OrdinaryDiffEqStabilizedIRK/src/OrdinaryDiffEqStabilizedIRK.jl
+++ b/lib/OrdinaryDiffEqStabilizedIRK/src/OrdinaryDiffEqStabilizedIRK.jl
@@ -9,7 +9,7 @@ import OrdinaryDiffEqCore: alg_order, alg_maximum_order,
                            OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
                            OrdinaryDiffEqAdaptiveAlgorithm,
                            OrdinaryDiffEqAdaptiveImplicitAlgorithm,
-                           alg_cache, _unwrap_val, DEFAULT_PRECS, @cache,
+                           alg_cache, _unwrap_val, @cache,
                            _reshape, _vec, full_cache, get_fsalfirstlast,
                            generic_solver_docstring
 
diff --git a/lib/OrdinaryDiffEqStabilizedIRK/src/algorithms.jl b/lib/OrdinaryDiffEqStabilizedIRK/src/algorithms.jl
index 85cf0825fa..73cfd9d9a1 100644
--- a/lib/OrdinaryDiffEqStabilizedIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqStabilizedIRK/src/algorithms.jl
@@ -7,11 +7,10 @@
     where `upper_bound` is an estimated upper bound on the spectral radius of the Jacobian matrix.
     If `eigen_est` is not provided, `upper_bound` will be estimated using the power iteration.",
     "eigen_est = nothing,")
-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
@@ -21,11 +20,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/lib/OrdinaryDiffEqStabilizedIRK/src/irkc_perform_step.jl b/lib/OrdinaryDiffEqStabilizedIRK/src/irkc_perform_step.jl
index da07d95c41..5696877a66 100644
--- a/lib/OrdinaryDiffEqStabilizedIRK/src/irkc_perform_step.jl
+++ b/lib/OrdinaryDiffEqStabilizedIRK/src/irkc_perform_step.jl
@@ -156,9 +156,8 @@ function perform_step!(integrator, cache::IRKCCache, repeat_step = false)
     # The the number of degree for Chebyshev polynomial
     #maxm = max(2,int(floor(sqrt(integrator.opts.internalnorm(integrator.opts.reltol,t)/(10 *eps(integrator.opts.internalnorm(uprev,t)))))))
     maxm = 50
-    mdeg = 1 + Int(floor(sqrt(1.54 * abs(dt) * integrator.eigen_est + 1)))
-    mdeg = (mdeg < minm) ? minm : mdeg
-    mdeg = (mdeg >= maxm) ? maxm : mdeg
+    mdeg = 1 + floor(Int, sqrt(1.54 * abs(dt) * integrator.eigen_est + 1))
+    mdeg = clamp(mdeg, minm, maxm)
 
     ω₀ = 1 + 2 / (13 * (mdeg^2))
     temp₁ = ω₀^2 - 1
diff --git a/lib/OrdinaryDiffEqStabilizedIRK/test/irkc_tests.jl b/lib/OrdinaryDiffEqStabilizedIRK/test/irkc_tests.jl
index 87a5ec3457..d6b92a7867 100644
--- a/lib/OrdinaryDiffEqStabilizedIRK/test/irkc_tests.jl
+++ b/lib/OrdinaryDiffEqStabilizedIRK/test/irkc_tests.jl
@@ -4,9 +4,8 @@ using OrdinaryDiffEqStabilizedIRK: maxeig!
 @testset "Power Iteration of Runge-Kutta-Chebyshev Tests" begin
     Random.seed!(123)
     eigen_est = (integrator) -> integrator.eigen_est = 1.5e2
-    for iip in [true, false], Alg in [IRKC]
+    @testset "iip=$iip, $Alg" for iip in [true, false], Alg in [IRKC]
         alg = Alg()
-        println(typeof(alg))
         A = randn(20, 20)
         B = randn(20, 20)
         test_f1 = !iip ? (u, p, t) -> A * u : (du, u, p, t) -> mul!(du, A, u)
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
index 503484ac33..5af395859c 100644
--- a/src/OrdinaryDiffEq.jl
+++ b/src/OrdinaryDiffEq.jl
@@ -35,10 +35,8 @@ import OrdinaryDiffEqCore: trivial_limiter!, CompositeAlgorithm, alg_order,
                            unwrap_alg, apply_step!, initialize_tstops, uses_uprev,
                            initialize_saveat,
                            isimplicit, initialize_d_discontinuities, isdtchangeable,
-                           _searchsortedfirst,
-                           _searchsortedlast,
-                           @unpack, ismultistep, DEFAULT_PRECS, isautoswitch,
-                           get_chunksize_int,
+                           _searchsortedfirst, _searchsortedlast,
+                           @unpack, ismultistep, isautoswitch, get_chunksize_int,
                            _unwrap_val, alg_autodiff, concrete_jac, alg_difftype,
                            standardtag,
                            alg_extrapolates, alg_maximum_order, alg_adaptive_order,
@@ -229,4 +227,8 @@ export AutoSwitch, AutoTsit5, AutoDP5,
 
 import OrdinaryDiffEqCore: IController, PIController, PIDController
 export IController, PIController, PIDController
+
+# unused, delete once StocasticDiffEq doesn't use this
+import OrdinaryDiffEqCore: DEFAULT_PRECS
+export DEFAULT_PRECS
 end # module
diff --git a/test/interface/linear_nonlinear_tests.jl b/test/interface/linear_nonlinear_tests.jl
index b948996346..537365a4e2 100644
--- a/test/interface/linear_nonlinear_tests.jl
+++ b/test/interface/linear_nonlinear_tests.jl
@@ -1,38 +1,43 @@
 using OrdinaryDiffEq, Test, Random, LinearAlgebra, LinearSolve
 Random.seed!(123)
 
-A = 0.01 * rand(3, 3)
+const A = 0.01 * rand(3, 3)
 rn = (du, u, p, t) -> begin
-    mul!(du, A, u)
+    du .= A * u
 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)
+    F = lu(convert(AbstractMatrix, W), check = false)
+    return F, 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)
+    F = lu(convert(AbstractMatrix, W), check = false)
+    IdentityOperator(size(W, 1)), F
+end
+
+function precslr(W, p)
+    F = lu(convert(AbstractMatrix, W), check = false)
+    F, F
 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
+@testset "precs" begin
+    @testset "$linsolve" for linsolve in (KrylovJL_GMRES(),
+                              KrylovJL_GMRES(precs = precsl),
+                              KrylovJL_GMRES(precs = precsr),
+                              KrylovJL_GMRES(precs = precslr))
+        sol = @test_nowarn solve(prob, TRBDF2(;linsolve,
+            smooth_est = false, concrete_jac = true), maxiters=20)
+        sol = @test_nowarn solve(prob, Rodas5P(autodiff = false;
+            linsolve, concrete_jac = true), maxiters=30)
+    end
+    @testset "$solver" for solver in (Rosenbrock23, FBDF, QNDF)
+        sol = @test_nowarn solve(prob, solver(
+            linsolve=KrylovJL_GMRES(precs = precslr), concrete_jac = true),
+            maxiters=30)
     end
-    Pr, Pr
 end
 
 sol = @test_nowarn solve(prob, TRBDF2(autodiff = false));
@@ -44,61 +49,7 @@ solref = @test_nowarn solve(prob,
     TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
         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));
-@test length(sol.t) < 20
-sol = @test_nowarn solve(prob,
-    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));
-@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));
-@test length(sol.t) < 20
-sol = @test_nowarn solve(prob,
-    Rodas4(autodiff = false, linsolve = KrylovJL_GMRES(),
-        precs = precslr, concrete_jac = true));
-@test length(sol.t) < 20
-
-sol = @test_nowarn solve(prob, TRBDF2(autodiff = false));
-@test length(sol.t) < 20
-sol = @test_nowarn solve(prob, TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES()));
-@test length(sol.t) < 20
-sol = @test_nowarn solve(prob,
-    TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
-        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));
-@test length(sol.t) < 20
-sol = @test_nowarn solve(prob,
-    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));
-@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));
-@test length(sol.t) < 20
-sol = @test_nowarn solve(prob,
-    Rodas4(autodiff = false, linsolve = KrylovJL_GMRES(),
-        precs = precslr, concrete_jac = true));
-@test length(sol.t) < 20
diff --git a/test/interface/linear_solver_split_ode_test.jl b/test/interface/linear_solver_split_ode_test.jl
index a959976899..83af64ff64 100644
--- a/test/interface/linear_solver_split_ode_test.jl
+++ b/test/interface/linear_solver_split_ode_test.jl
@@ -2,8 +2,6 @@ using Test
 using OrdinaryDiffEq
 using LinearAlgebra, LinearSolve
 
-import OrdinaryDiffEqCore.dolinsolve
-
 n = 8
 dt = 1 / 1000
 u0 = ones(n)
diff --git a/test/interface/linear_solver_test.jl b/test/interface/linear_solver_test.jl
index a615deb515..3471db21c7 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;
diff --git a/test/interface/sized_matrix_tests.jl b/test/interface/sized_matrix_tests.jl
index 79ba7b6e73..dd43f03ca2 100644
--- a/test/interface/sized_matrix_tests.jl
+++ b/test/interface/sized_matrix_tests.jl
@@ -7,17 +7,17 @@ end
 η0 = 1
 τ = 1
 α = 0.8
-p_giesekus = [η0, τ, α]
 
 ## Exercise code paths that are in-place but not `<:Array` but not AbstractVector
 σ0 = SizedMatrix{3, 3}([1.0 0.0 0.0; 0.0 2.0 0.0; 3.0 0.0 0.0])
 
-prob_giesekus = ODEProblem(dudt!, σ0, (0.0, 2.0), p_giesekus)
+prob_giesekus = ODEProblem(dudt!, σ0, (0.0, 2.0))
+prob_normal = ODEProblem(dudt!, Matrix(σ0), (0.0, 2.0))
+
 
 if VERSION >= v"1.9"
-    solve_giesekus = solve(prob_giesekus, Rodas4(), saveat = 0.2, abstol = 1e-14,
-        reltol = 1e-14)
-    for alg in [
+    @testset "giesekus" begin
+    @testset "$alg" for alg in [
         Rosenbrock23(),
         Rodas4(),
         Rodas4P(),
@@ -30,7 +30,13 @@ if VERSION >= v"1.9"
         Vern9(),
         DP5()
     ]
-        sol = solve(prob_giesekus, alg, saveat = 0.2, abstol = 1e-14, reltol = 1e-14)
-        @test Array(sol) ≈ Array(solve_giesekus)
+        sol = solve(prob_giesekus, alg, saveat = 0.2, abstol = 1e-10, reltol = 1e-10)
+        sol_normal = solve(prob_giesekus, alg, saveat = 0.2, abstol = 1e-10, reltol = 1e-10)
+        @test sol.retcode == sol_normal.retcode == ReturnCode.Success
+        @test Array(sol) ≈ Array(sol_normal)
+        # give ourselves some stats wiggle room
+        @test sol.stats.nreject <= 2*sol_normal.stats.nreject + 1
+        @test sol.stats.naccept <= 2*sol_normal.stats.naccept + 1
+    end
     end
 end
diff --git a/test/interface/utility_tests.jl b/test/interface/utility_tests.jl
index d0dbf2427b..f95e3aebd5 100644
--- a/test/interface/utility_tests.jl
+++ b/test/interface/utility_tests.jl
@@ -60,9 +60,7 @@ end
         jac_prototype = MatrixOperator(similar(A);
             update_func! = (J, u, p, t) -> J .= t .*
                                                 A))
-
-    for Alg in [ImplicitEuler, Rosenbrock23, Rodas5]
-        println(Alg)
+    @testset "$Alg" for Alg in [ImplicitEuler, Rosenbrock23, Rodas5]
         sol1 = solve(ODEProblem(fun1, u0, tspan), Alg(); adaptive = false, dt = 0.01)
         sol2 = solve(ODEProblem(fun2, u0, tspan), Alg(); adaptive = false, dt = 0.01)
         @test sol1(1.0)≈sol2(1.0) atol=1e-3