use LinearSolve precs setup

docs/src/native/solvers.md

@@ -18,10 +18,6 @@ documentation.
  uses the LinearSolve.jl default algorithm choice. For more information on available
  algorithm choices, see the
  [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/).
- - `precs`: the choice of preconditioners for the linear solver. Defaults to using no
- preconditioners. For more information on specifying preconditioners for LinearSolve
- algorithms, consult the
- [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/).
  - `linesearch`: the line search algorithm to use. Defaults to [`NoLineSearch()`](@ref),
  which means that no line search is performed. Algorithms from
  [`LineSearches.jl`](https://github.com/JuliaNLSolvers/LineSearches.jl/) must be

docs/src/tutorials/large_systems.md

@@ -221,35 +221,27 @@ choices, see the
 
 Any [LinearSolve.jl-compatible preconditioner](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/)
 can be used as a preconditioner in the linear solver interface. To define preconditioners,
-one must define a `precs` function in compatible with nonlinear solvers which returns the
+one must define a `precs` function in compatible with linear solver which returns the
 left and right preconditioners, matrices which approximate the inverse of `W = I - gamma*J`
 used in the solution of the ODE. An example of this with using
 [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl) is as follows:
 
 ```julia
-# FIXME: On 1.10+ this is broken. Skipping this for now.
 using IncompleteLU
 
-function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
- if newW === nothing || newW
- Pl = ilu(W, τ = 50.0)
- else
- Pl = Plprev
- end
+function incompletelu(W, p=nothing)
+ Pl = ilu(W, τ = 50.0)
  Pl, nothing
 end
 
 @btime solve(prob_brusselator_2d_sparse,
- NewtonRaphson(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true));
-nothing # hide
+ NewtonRaphson(linsolve = KrylovJL_GMRES(precs = incompletelu), concrete_jac = true));
+nothing
 ```
 
 Notice a few things about this preconditioner. This preconditioner uses the sparse Jacobian,
 and thus we set `concrete_jac = true` to tell the algorithm to generate the Jacobian
-(otherwise, a Jacobian-free algorithm is used with GMRES by default). Then `newW = true`
-whenever a new `W` matrix is computed, and `newW = nothing` during the startup phase of the
-solver. Thus, we do a check `newW === nothing || newW` and when true, it's only at these
-points when we update the preconditioner, otherwise we just pass on the previous version.
+(otherwise, a Jacobian-free algorithm is used with GMRES by default).
 We use `convert(AbstractMatrix,W)` to get the concrete `W` matrix (matching `jac_prototype`,
 thus `SpraseMatrixCSC`) which we can use in the preconditioner's definition. Then we use
 `IncompleteLU.ilu` on that sparse matrix to generate the preconditioner. We return
@@ -265,39 +257,31 @@ which is more automatic. The setup is very similar to before:
 ```@example ill_conditioned_nlprob
 using AlgebraicMultigrid
 
-function algebraicmultigrid(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
- if newW === nothing || newW
- Pl = aspreconditioner(ruge_stuben(convert(AbstractMatrix, W)))
- else
- Pl = Plprev
- end
+function algebraicmultigrid(W, p=nothing)
+ Pl = aspreconditioner(ruge_stuben(convert(AbstractMatrix, W)))
  Pl, nothing
 end
 
 @btime solve(prob_brusselator_2d_sparse,
  NewtonRaphson(
- linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid, concrete_jac = true));
+ linsolve = KrylovJL_GMRES(precs = algebraicmultigrid), concrete_jac = true));
 nothing # hide
 ```
 
 or with a Jacobi smoother:
 
 ```@example ill_conditioned_nlprob
-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
+function algebraicmultigrid2(W, p=nothing)
+ 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, nothing
 end
 
 @btime solve(prob_brusselator_2d_sparse,
  NewtonRaphson(
- linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid2, concrete_jac = true));
+ linsolve = KrylovJL_GMRES(precs = algebraicmultigrid2), concrete_jac = true));
 nothing # hide
 ```
 

src/algorithms/gauss_newton.jl

@@ -1,14 +1,14 @@
 """
  GaussNewton(; concrete_jac = nothing, linsolve = nothing, linesearch = NoLineSearch(),
- precs = DEFAULT_PRECS, adkwargs...)
+ adkwargs...)
 
 An advanced GaussNewton implementation with support for efficient handling of sparse
 matrices via colored automatic differentiation and preconditioned linear solvers. Designed
 for large-scale and numerically-difficult nonlinear least squares problems.
 """
-function GaussNewton(; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS,
+function GaussNewton(; concrete_jac = nothing, linsolve = nothing,
  linesearch = NoLineSearch(), vjp_autodiff = nothing, autodiff = nothing)
- descent = NewtonDescent(; linsolve, precs)
+ descent = NewtonDescent(; linsolve)
  return GeneralizedFirstOrderAlgorithm(; concrete_jac, name = :GaussNewton, descent,
  jacobian_ad = autodiff, reverse_ad = vjp_autodiff, linesearch)
 end

src/algorithms/klement.jl

@@ -1,6 +1,6 @@
 """
  Klement(; max_resets = 100, linsolve = NoLineSearch(), linesearch = nothing,
- precs = DEFAULT_PRECS, alpha = nothing, init_jacobian::Val = Val(:identity),
+ alpha = nothing, init_jacobian::Val = Val(:identity),
  autodiff = nothing)
 
 An implementation of `Klement` [klement2014using](@citep) with line search, preconditioning
@@ -25,7 +25,7 @@ over this.
  differentiable problems.
 """
 function Klement(; max_resets::Int = 100, linsolve = nothing, alpha = nothing,
- linesearch = NoLineSearch(), precs = DEFAULT_PRECS,
+ linesearch = NoLineSearch(),
  autodiff = nothing, init_jacobian::Val{IJ} = Val(:identity)) where {IJ}
  if !(linesearch isa AbstractNonlinearSolveLineSearchAlgorithm)
  Base.depwarn(
@@ -48,7 +48,7 @@ function Klement(; max_resets::Int = 100, linsolve = nothing, alpha = nothing,
  CJ = IJ === :true_jacobian || IJ === :true_jacobian_diagonal
 
  return ApproximateJacobianSolveAlgorithm{CJ, :Klement}(;
- linesearch, descent = NewtonDescent(; linsolve, precs),
+ linesearch, descent = NewtonDescent(; linsolve),
  update_rule = KlementUpdateRule(),
  reinit_rule = IllConditionedJacobianReset(), max_resets, initialization)
 end

src/algorithms/levenberg_marquardt.jl

@@ -1,6 +1,6 @@
 """
  LevenbergMarquardt(; linsolve = nothing,
- precs = DEFAULT_PRECS, damping_initial::Real = 1.0, α_geodesic::Real = 0.75,
+ damping_initial::Real = 1.0, α_geodesic::Real = 0.75,
  damping_increase_factor::Real = 2.0, damping_decrease_factor::Real = 3.0,
  finite_diff_step_geodesic = 0.1, b_uphill::Real = 1.0, autodiff = nothing,
  min_damping_D::Real = 1e-8, disable_geodesic = Val(false))
@@ -31,7 +31,7 @@ For the remaining arguments, see [`GeodesicAcceleration`](@ref) and
 [`NonlinearSolve.LevenbergMarquardtTrustRegion`](@ref) documentations.
 """
 function LevenbergMarquardt(;
- concrete_jac = missing, linsolve = nothing, precs = DEFAULT_PRECS,
+ concrete_jac = missing, linsolve = nothing,
  damping_initial::Real = 1.0, α_geodesic::Real = 0.75,
  damping_increase_factor::Real = 2.0, damping_decrease_factor::Real = 3.0,
  finite_diff_step_geodesic = 0.1, b_uphill::Real = 1.0,
@@ -45,7 +45,6 @@ function LevenbergMarquardt(;
  end
 
  descent = DampedNewtonDescent(; linsolve,
- precs,
  initial_damping = damping_initial,
  damping_fn = LevenbergMarquardtDampingFunction(
  damping_increase_factor, damping_decrease_factor, min_damping_D))

src/algorithms/pseudo_transient.jl

@@ -1,7 +1,7 @@
 """
  PseudoTransient(; concrete_jac = nothing, linsolve = nothing,
  linesearch::AbstractNonlinearSolveLineSearchAlgorithm = NoLineSearch(),
- precs = DEFAULT_PRECS, autodiff = nothing)
+ autodiff = nothing)
 
 An implementation of PseudoTransient Method [coffey2003pseudotransient](@cite) that is used
 to solve steady state problems in an accelerated manner. It uses an adaptive time-stepping
@@ -17,8 +17,8 @@ This implementation specifically uses "switched evolution relaxation"
 """
 function PseudoTransient(; concrete_jac = nothing, linsolve = nothing,
  linesearch::AbstractNonlinearSolveLineSearchAlgorithm = NoLineSearch(),
- precs = DEFAULT_PRECS, autodiff = nothing, alpha_initial = 1e-3)
- descent = DampedNewtonDescent(; linsolve, precs, initial_damping = alpha_initial,
+ autodiff = nothing, alpha_initial = 1e-3)
+ descent = DampedNewtonDescent(; linsolve, initial_damping = alpha_initial,
  damping_fn = SwitchedEvolutionRelaxation())
  return GeneralizedFirstOrderAlgorithm(;
  concrete_jac, name = :PseudoTransient, linesearch, descent, jacobian_ad = autodiff)

src/algorithms/raphson.jl

@@ -1,14 +1,14 @@
 """
  NewtonRaphson(; concrete_jac = nothing, linsolve = nothing, linesearch = NoLineSearch(),
- precs = DEFAULT_PRECS, autodiff = nothing)
+ autodiff = nothing)
 
 An advanced NewtonRaphson implementation with support for efficient handling of sparse
 matrices via colored automatic differentiation and preconditioned linear solvers. Designed
 for large-scale and numerically-difficult nonlinear systems.
 """
 function NewtonRaphson(; concrete_jac = nothing, linsolve = nothing,
- linesearch = NoLineSearch(), precs = DEFAULT_PRECS, autodiff = nothing)
- descent = NewtonDescent(; linsolve, precs)
+ linesearch = NoLineSearch(), autodiff = nothing)
+ descent = NewtonDescent(; linsolve)
  return GeneralizedFirstOrderAlgorithm(;
  concrete_jac, name = :NewtonRaphson, linesearch, descent, jacobian_ad = autodiff)
 end

src/algorithms/trust_region.jl

@@ -1,5 +1,5 @@
 """
- TrustRegion(; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS,
+ TrustRegion(; concrete_jac = nothing, linsolve = nothing,
  radius_update_scheme = RadiusUpdateSchemes.Simple, max_trust_radius::Real = 0 // 1,
  initial_trust_radius::Real = 0 // 1, step_threshold::Real = 1 // 10000,
  shrink_threshold::Real = 1 // 4, expand_threshold::Real = 3 // 4,
@@ -19,13 +19,13 @@ for large-scale and numerically-difficult nonlinear systems.
 For the remaining arguments, see [`NonlinearSolve.GenericTrustRegionScheme`](@ref)
 documentation.
 """
-function TrustRegion(; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS,
+function TrustRegion(; concrete_jac = nothing, linsolve = nothing,
  radius_update_scheme = RadiusUpdateSchemes.Simple, max_trust_radius::Real = 0 // 1,
  initial_trust_radius::Real = 0 // 1, step_threshold::Real = 1 // 10000,
  shrink_threshold::Real = 1 // 4, expand_threshold::Real = 3 // 4,
  shrink_factor::Real = 1 // 4, expand_factor::Real = 2 // 1,
  max_shrink_times::Int = 32, autodiff = nothing, vjp_autodiff = nothing)
- descent = Dogleg(; linsolve, precs)
+ descent = Dogleg(; linsolve)
  if autodiff !== nothing && ADTypes.mode(autodiff) isa ADTypes.ForwardMode
  forward_ad = autodiff
  else