From 7bdf88e251d288e36088c055591ff7d2b0ffa8ba Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 24 Sep 2024 17:11:33 -0400 Subject: [PATCH 1/3] use LinearSolve precs setup --- docs/src/native/solvers.md | 4 -- docs/src/tutorials/large_systems.md | 46 ++++++---------- src/algorithms/gauss_newton.jl | 6 +-- src/algorithms/klement.jl | 6 +-- src/algorithms/levenberg_marquardt.jl | 5 +- src/algorithms/pseudo_transient.jl | 6 +-- src/algorithms/raphson.jl | 6 +-- src/algorithms/trust_region.jl | 6 +-- src/default.jl | 76 +++++++++++++-------------- src/descent/damped_newton.jl | 4 +- src/descent/dogleg.jl | 10 ++-- src/descent/newton.jl | 5 +- src/descent/steepest.jl | 4 +- src/internal/linear_solve.jl | 55 +++---------------- test/core/rootfind_tests.jl | 27 ++++------ 15 files changed, 94 insertions(+), 172 deletions(-) diff --git a/docs/src/native/solvers.md b/docs/src/native/solvers.md index d2c0fc6e5..e28476b60 100644 --- a/docs/src/native/solvers.md +++ b/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 diff --git a/docs/src/tutorials/large_systems.md b/docs/src/tutorials/large_systems.md index 368535287..128533630 100644 --- a/docs/src/tutorials/large_systems.md +++ b/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 ``` diff --git a/src/algorithms/gauss_newton.jl b/src/algorithms/gauss_newton.jl index 54b7193fc..ce14724d7 100644 --- a/src/algorithms/gauss_newton.jl +++ b/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 diff --git a/src/algorithms/klement.jl b/src/algorithms/klement.jl index 66daec1b7..d0777f868 100644 --- a/src/algorithms/klement.jl +++ b/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 diff --git a/src/algorithms/levenberg_marquardt.jl b/src/algorithms/levenberg_marquardt.jl index 66554ee2d..493064598 100644 --- a/src/algorithms/levenberg_marquardt.jl +++ b/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)) diff --git a/src/algorithms/pseudo_transient.jl b/src/algorithms/pseudo_transient.jl index fda39a3b9..26df2650d 100644 --- a/src/algorithms/pseudo_transient.jl +++ b/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) diff --git a/src/algorithms/raphson.jl b/src/algorithms/raphson.jl index c6847f54c..67e17bc46 100644 --- a/src/algorithms/raphson.jl +++ b/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 diff --git a/src/algorithms/trust_region.jl b/src/algorithms/trust_region.jl index d68e2d9ec..49eb41783 100644 --- a/src/algorithms/trust_region.jl +++ b/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 diff --git a/src/default.jl b/src/default.jl index 8a9aac6af..5a2ccb8ad 100644 --- a/src/default.jl +++ b/src/default.jl @@ -338,7 +338,7 @@ end """ RobustMultiNewton(::Type{T} = Float64; concrete_jac = nothing, linsolve = nothing, - precs = DEFAULT_PRECS, autodiff = nothing) + autodiff = nothing) A polyalgorithm focused on robustness. It uses a mixture of Newton methods with different globalizing techniques (trust region updates, line searches, etc.) in order to find a @@ -354,20 +354,20 @@ or more precision / more stable linear solver choice is required). are compatible with the problem type. Defaults to `Float64`. """ function RobustMultiNewton(::Type{T} = Float64; concrete_jac = nothing, linsolve = nothing, - precs = DEFAULT_PRECS, autodiff = nothing) where {T} + autodiff = nothing) where {T} if __is_complex(T) # Let's atleast have something here for complex numbers - algs = (NewtonRaphson(; concrete_jac, linsolve, precs, autodiff),) + algs = (NewtonRaphson(; concrete_jac, linsolve, autodiff),) else - algs = (TrustRegion(; concrete_jac, linsolve, precs, autodiff), - TrustRegion(; concrete_jac, linsolve, precs, autodiff, + algs = (TrustRegion(; concrete_jac, linsolve, autodiff), + TrustRegion(; concrete_jac, linsolve, autodiff, radius_update_scheme = RadiusUpdateSchemes.Bastin), - NewtonRaphson(; concrete_jac, linsolve, precs, autodiff), - NewtonRaphson(; concrete_jac, linsolve, precs, + NewtonRaphson(; concrete_jac, linsolve, autodiff), + NewtonRaphson(; concrete_jac, linsolve, linesearch = LineSearchesJL(; method = BackTracking()), autodiff), - TrustRegion(; concrete_jac, linsolve, precs, + TrustRegion(; concrete_jac, linsolve, radius_update_scheme = RadiusUpdateSchemes.NLsolve, autodiff), - TrustRegion(; concrete_jac, linsolve, precs, + TrustRegion(; concrete_jac, linsolve, radius_update_scheme = RadiusUpdateSchemes.Fan, autodiff)) end return NonlinearSolvePolyAlgorithm(algs, Val(:NLS)) @@ -375,7 +375,7 @@ end """ FastShortcutNonlinearPolyalg(::Type{T} = Float64; concrete_jac = nothing, - linsolve = nothing, precs = DEFAULT_PRECS, must_use_jacobian::Val = Val(false), + linsolve = nothing, must_use_jacobian::Val = Val(false), prefer_simplenonlinearsolve::Val{SA} = Val(false), autodiff = nothing, u0_len::Union{Int, Nothing} = nothing) where {T} @@ -395,19 +395,19 @@ for more performance and then tries more robust techniques if the faster ones fa """ function FastShortcutNonlinearPolyalg( ::Type{T} = Float64; concrete_jac = nothing, linsolve = nothing, - precs = DEFAULT_PRECS, must_use_jacobian::Val{JAC} = Val(false), + must_use_jacobian::Val{JAC} = Val(false), prefer_simplenonlinearsolve::Val{SA} = Val(false), u0_len::Union{Int, Nothing} = nothing, autodiff = nothing) where {T, JAC, SA} start_index = 1 if JAC if __is_complex(T) - algs = (NewtonRaphson(; concrete_jac, linsolve, precs, autodiff),) + algs = (NewtonRaphson(; concrete_jac, linsolve, autodiff),) else - algs = (NewtonRaphson(; concrete_jac, linsolve, precs, autodiff), - NewtonRaphson(; concrete_jac, linsolve, precs, + algs = (NewtonRaphson(; concrete_jac, linsolve, autodiff), + NewtonRaphson(; concrete_jac, linsolve, linesearch = LineSearchesJL(; method = BackTracking()), autodiff), - TrustRegion(; concrete_jac, linsolve, precs, autodiff), - TrustRegion(; concrete_jac, linsolve, precs, + TrustRegion(; concrete_jac, linsolve, autodiff), + TrustRegion(; concrete_jac, linsolve, radius_update_scheme = RadiusUpdateSchemes.Bastin, autodiff)) end else @@ -418,35 +418,35 @@ function FastShortcutNonlinearPolyalg( algs = (SimpleBroyden(), Broyden(; init_jacobian = Val(:true_jacobian), autodiff), SimpleKlement(), - NewtonRaphson(; concrete_jac, linsolve, precs, autodiff)) + NewtonRaphson(; concrete_jac, linsolve, autodiff)) else start_index = u0_len !== nothing ? (u0_len ≤ 25 ? 4 : 1) : 1 algs = (SimpleBroyden(), Broyden(; init_jacobian = Val(:true_jacobian), autodiff), SimpleKlement(), - NewtonRaphson(; concrete_jac, linsolve, precs, autodiff), - NewtonRaphson(; concrete_jac, linsolve, precs, + NewtonRaphson(; concrete_jac, linsolve, autodiff), + NewtonRaphson(; concrete_jac, linsolve, linesearch = LineSearchesJL(; method = BackTracking()), autodiff), - TrustRegion(; concrete_jac, linsolve, precs, + TrustRegion(; concrete_jac, linsolve, radius_update_scheme = RadiusUpdateSchemes.Bastin, autodiff)) end else if __is_complex(T) algs = ( Broyden(), Broyden(; init_jacobian = Val(:true_jacobian), autodiff), - Klement(; linsolve, precs, autodiff), - NewtonRaphson(; concrete_jac, linsolve, precs, autodiff)) + Klement(; linsolve, autodiff), + NewtonRaphson(; concrete_jac, linsolve, autodiff)) else # TODO: This number requires a bit rigorous testing start_index = u0_len !== nothing ? (u0_len ≤ 25 ? 4 : 1) : 1 algs = (Broyden(; autodiff), Broyden(; init_jacobian = Val(:true_jacobian), autodiff), - Klement(; linsolve, precs, autodiff), - NewtonRaphson(; concrete_jac, linsolve, precs, autodiff), - NewtonRaphson(; concrete_jac, linsolve, precs, + Klement(; linsolve, autodiff), + NewtonRaphson(; concrete_jac, linsolve, autodiff), + NewtonRaphson(; concrete_jac, linsolve, linesearch = LineSearchesJL(; method = BackTracking()), autodiff), - TrustRegion(; concrete_jac, linsolve, precs, autodiff), - TrustRegion(; concrete_jac, linsolve, precs, + TrustRegion(; concrete_jac, linsolve, autodiff), + TrustRegion(; concrete_jac, linsolve, radius_update_scheme = RadiusUpdateSchemes.Bastin, autodiff)) end end @@ -456,7 +456,7 @@ end """ FastShortcutNLLSPolyalg(::Type{T} = Float64; concrete_jac = nothing, linsolve = nothing, - precs = DEFAULT_PRECS, autodiff = nothing, kwargs...) + autodiff = nothing, kwargs...) A polyalgorithm focused on balancing speed and robustness. It first tries less robust methods for more performance and then tries more robust techniques if the faster ones fail. @@ -468,23 +468,23 @@ for more performance and then tries more robust techniques if the faster ones fa """ function FastShortcutNLLSPolyalg( ::Type{T} = Float64; concrete_jac = nothing, linsolve = nothing, - precs = DEFAULT_PRECS, autodiff = nothing, kwargs...) where {T} + autodiff = nothing, kwargs...) where {T} if __is_complex(T) - algs = (GaussNewton(; concrete_jac, linsolve, precs, autodiff, kwargs...), + algs = (GaussNewton(; concrete_jac, linsolve, autodiff, kwargs...), LevenbergMarquardt(; - linsolve, precs, autodiff, disable_geodesic = Val(true), kwargs...), - LevenbergMarquardt(; linsolve, precs, autodiff, kwargs...)) + linsolve, autodiff, disable_geodesic = Val(true), kwargs...), + LevenbergMarquardt(; linsolve, autodiff, kwargs...)) else - algs = (GaussNewton(; concrete_jac, linsolve, precs, autodiff, kwargs...), + algs = (GaussNewton(; concrete_jac, linsolve, autodiff, kwargs...), LevenbergMarquardt(; - linsolve, precs, disable_geodesic = Val(true), autodiff, kwargs...), - TrustRegion(; concrete_jac, linsolve, precs, autodiff, kwargs...), - GaussNewton(; concrete_jac, linsolve, precs, + linsolve, disable_geodesic = Val(true), autodiff, kwargs...), + TrustRegion(; concrete_jac, linsolve, autodiff, kwargs...), + GaussNewton(; concrete_jac, linsolve, linesearch = LineSearchesJL(; method = BackTracking()), autodiff, kwargs...), - TrustRegion(; concrete_jac, linsolve, precs, + TrustRegion(; concrete_jac, linsolve, radius_update_scheme = RadiusUpdateSchemes.Bastin, autodiff, kwargs...), - LevenbergMarquardt(; linsolve, precs, autodiff, kwargs...)) + LevenbergMarquardt(; linsolve, autodiff, kwargs...)) end return NonlinearSolvePolyAlgorithm(algs, Val(:NLLS)) end diff --git a/src/descent/damped_newton.jl b/src/descent/damped_newton.jl index 77b204b13..2d852560d 100644 --- a/src/descent/damped_newton.jl +++ b/src/descent/damped_newton.jl @@ -1,5 +1,5 @@ """ - DampedNewtonDescent(; linsolve = nothing, precs = DEFAULT_PRECS, initial_damping, + DampedNewtonDescent(; linsolve = nothing, initial_damping, damping_fn) A Newton descent algorithm with damping. The damping factor is computed using the @@ -19,7 +19,6 @@ The damping factor returned must be a non-negative number. """ @kwdef @concrete struct DampedNewtonDescent <: AbstractDescentAlgorithm linsolve = nothing - precs = DEFAULT_PRECS initial_damping damping_fn end @@ -27,7 +26,6 @@ end function Base.show(io::IO, d::DampedNewtonDescent) modifiers = String[] d.linsolve !== nothing && push!(modifiers, "linsolve = $(d.linsolve)") - d.precs !== DEFAULT_PRECS && push!(modifiers, "precs = $(d.precs)") push!(modifiers, "initial_damping = $(d.initial_damping)") push!(modifiers, "damping_fn = $(d.damping_fn)") print(io, "DampedNewtonDescent($(join(modifiers, ", ")))") diff --git a/src/descent/dogleg.jl b/src/descent/dogleg.jl index 5d0bb1c7c..c77535593 100644 --- a/src/descent/dogleg.jl +++ b/src/descent/dogleg.jl @@ -1,5 +1,5 @@ """ - Dogleg(; linsolve = nothing, precs = DEFAULT_PRECS) + Dogleg(; linsolve = nothing) Switch between Newton's method and the steepest descent method depending on the size of the trust region. The trust region is specified via keyword argument `trust_region` to @@ -20,17 +20,17 @@ end supports_trust_region(::Dogleg) = true get_linear_solver(alg::Dogleg) = get_linear_solver(alg.newton_descent) -function Dogleg(; linsolve = nothing, precs = DEFAULT_PRECS, damping = False, +function Dogleg(; linsolve = nothing, damping = False, damping_fn = missing, initial_damping = missing, kwargs...) if damping === False - return Dogleg(NewtonDescent(; linsolve, precs), SteepestDescent(; linsolve, precs)) + return Dogleg(NewtonDescent(; linsolve), SteepestDescent(; linsolve)) end if damping_fn === missing || initial_damping === missing throw(ArgumentError("`damping_fn` and `initial_damping` must be supplied if \ `damping = Val(true)`.")) end - return Dogleg(DampedNewtonDescent(; linsolve, precs, damping_fn, initial_damping), - SteepestDescent(; linsolve, precs)) + return Dogleg(DampedNewtonDescent(; linsolve, damping_fn, initial_damping), + SteepestDescent(; linsolve)) end @concrete mutable struct DoglegCache{pre_inverted, normalform} <: AbstractDescentCache diff --git a/src/descent/newton.jl b/src/descent/newton.jl index 2fe7abf9a..c02db35f1 100644 --- a/src/descent/newton.jl +++ b/src/descent/newton.jl @@ -1,5 +1,5 @@ """ - NewtonDescent(; linsolve = nothing, precs = DEFAULT_PRECS) + NewtonDescent(; linsolve = nothing) Compute the descent direction as ``J δu = -fu``. For non-square Jacobian problems, this is commonly referred to as the Gauss-Newton Descent. @@ -8,13 +8,10 @@ See also [`Dogleg`](@ref), [`SteepestDescent`](@ref), [`DampedNewtonDescent`](@r """ @kwdef @concrete struct NewtonDescent <: AbstractDescentAlgorithm linsolve = nothing - precs = DEFAULT_PRECS end function Base.show(io::IO, d::NewtonDescent) modifiers = String[] - d.linsolve !== nothing && push!(modifiers, "linsolve = $(d.linsolve)") - d.precs !== DEFAULT_PRECS && push!(modifiers, "precs = $(d.precs)") print(io, "NewtonDescent($(join(modifiers, ", ")))") end diff --git a/src/descent/steepest.jl b/src/descent/steepest.jl index cc2f4d128..7f63ac408 100644 --- a/src/descent/steepest.jl +++ b/src/descent/steepest.jl @@ -1,5 +1,5 @@ """ - SteepestDescent(; linsolve = nothing, precs = DEFAULT_PRECS) + SteepestDescent(; linsolve = nothing) Compute the descent direction as ``δu = -Jᵀfu``. The linear solver and preconditioner are only used if `J` is provided in the inverted form. @@ -8,13 +8,11 @@ See also [`Dogleg`](@ref), [`NewtonDescent`](@ref), [`DampedNewtonDescent`](@ref """ @kwdef @concrete struct SteepestDescent <: AbstractDescentAlgorithm linsolve = nothing - precs = DEFAULT_PRECS end function Base.show(io::IO, d::SteepestDescent) modifiers = String[] d.linsolve !== nothing && push!(modifiers, "linsolve = $(d.linsolve)") - d.precs !== DEFAULT_PRECS && push!(modifiers, "precs = $(d.precs)") print(io, "SteepestDescent($(join(modifiers, ", ")))") end diff --git a/src/internal/linear_solve.jl b/src/internal/linear_solve.jl index 27e6a780f..9decc5507 100644 --- a/src/internal/linear_solve.jl +++ b/src/internal/linear_solve.jl @@ -19,7 +19,7 @@ handled: ```julia (cache::LinearSolverCache)(; A = nothing, b = nothing, linu = nothing, du = nothing, p = nothing, - weight = nothing, cachedata = nothing, reuse_A_if_factorization = false, kwargs...) + cachedata = nothing, reuse_A_if_factorization = false, kwargs...) ``` Returns the solution of the system `u` and stores the updated cache in `cache.lincache`. @@ -49,7 +49,6 @@ not mutated, we do this by copying over `A` to a preconstructed cache. additional_lincache::Any A b - precs stats::NLStats end @@ -75,24 +74,15 @@ function LinearSolverCache(alg, linsolve, A, b, u; stats, kwargs...) (linsolve === nothing && A isa SMatrix) || (A isa Diagonal) || (linsolve isa typeof(\)) - return LinearSolverCache(nothing, nothing, nothing, A, b, nothing, stats) + return LinearSolverCache(nothing, nothing, nothing, A, b, stats) end @bb u_ = copy(u_fixed) linprob = LinearProblem(A, b; u0 = u_, kwargs...) - weight = __init_ones(u_fixed) - if __hasfield(alg, Val(:precs)) - precs = alg.precs - Pl_, Pr_ = precs(A, nothing, u, ntuple(Returns(nothing), 6)...) - else - precs, Pl_, Pr_ = nothing, nothing, nothing - end - Pl, Pr = __wrapprecs(Pl_, Pr_, weight) - # Unalias here, we will later use these as caches - lincache = init(linprob, linsolve; alias_A = false, alias_b = false, Pl, Pr) + lincache = init(linprob, linsolve; alias_A = false, alias_b = false) - return LinearSolverCache(lincache, linsolve, nothing, nothing, nothing, precs, stats) + return LinearSolverCache(lincache, linsolve, nothing, nothing, nothing, stats) end @kwdef @concrete struct LinearSolveResult @@ -120,35 +110,13 @@ end # Use LinearSolve.jl function (cache::LinearSolverCache)(; A = nothing, b = nothing, linu = nothing, du = nothing, - p = nothing, weight = nothing, cachedata = nothing, + p = nothing, cachedata = nothing, reuse_A_if_factorization = false, verbose = true, kwargs...) cache.stats.nsolve += 1 __update_A!(cache, A, reuse_A_if_factorization) b !== nothing && (cache.lincache.b = b) linu !== nothing && __set_lincache_u!(cache, linu) - - Plprev = cache.lincache.Pl isa ComposePreconditioner ? cache.lincache.Pl.outer : - cache.lincache.Pl - Prprev = cache.lincache.Pr isa ComposePreconditioner ? cache.lincache.Pr.outer : - cache.lincache.Pr - - if cache.precs === nothing - _Pl, _Pr = nothing, nothing - else - _Pl, _Pr = cache.precs(cache.lincache.A, du, linu, p, nothing, - A !== nothing, Plprev, Prprev, cachedata) - end - - if (_Pl !== nothing || _Pr !== nothing) - _weight = weight === nothing ? - (cache.lincache.Pr isa Diagonal ? cache.lincache.Pr.diag : - cache.lincache.Pr.inner.diag) : weight - Pl, Pr = __wrapprecs(_Pl, _Pr, _weight) - cache.lincache.Pl = Pl - cache.lincache.Pr = Pr - end - linres = solve!(cache.lincache) cache.lincache = linres.cache # Unfortunately LinearSolve.jl doesn't have the most uniform ReturnCode handling @@ -169,7 +137,7 @@ function (cache::LinearSolverCache)(; linprob = LinearProblem(A, b; u0 = linres.u) cache.additional_lincache = init( linprob, QRFactorization(ColumnNorm()); alias_u0 = false, - alias_A = false, alias_b = false, cache.lincache.Pl, cache.lincache.Pr) + alias_A = false, alias_b = false) else cache.additional_lincache.A = A cache.additional_lincache.b = b @@ -242,17 +210,6 @@ function __set_lincache_A(lincache, new_A) end end -@inline function __wrapprecs(_Pl, _Pr, weight) - Pl = _Pl !== nothing ? - ComposePreconditioner(InvPreconditioner(Diagonal(_vec(weight))), _Pl) : - InvPreconditioner(Diagonal(_vec(weight))) - - Pr = _Pr !== nothing ? ComposePreconditioner(Diagonal(_vec(weight)), _Pr) : - Diagonal(_vec(weight)) - - return Pl, Pr -end - @inline __needs_square_A(_, ::Number) = false @inline __needs_square_A(::Nothing, ::Number) = false @inline __needs_square_A(::Nothing, _) = false diff --git a/test/core/rootfind_tests.jl b/test/core/rootfind_tests.jl index bca02a44b..518adbec5 100644 --- a/test/core/rootfind_tests.jl +++ b/test/core/rootfind_tests.jl @@ -70,26 +70,23 @@ end NewtonRaphson(), abstol = 1e-9) @test (@ballocated solve!($cache)) < 200 end - - precs = [(u0) -> NonlinearSolve.DEFAULT_PRECS, - u0 -> ((args...) -> (Diagonal(rand!(similar(u0))), nothing))] - - @testset "[IIP] u0: $(typeof(u0)) precs: $(_nameof(prec)) linsolve: $(_nameof(linsolve))" for u0 in ([ + + precs = (A, p=nothing) -> (Diagonal(randn!(similar(A, size(A,1)))), nothing) + @testset "[IIP] u0: $(typeof(u0)) linsolve: $(_nameof(linsolve))" for u0 in ([ 1.0, 1.0],), - prec in precs, - linsolve in (nothing, KrylovJL_GMRES(), \) + linsolve in (nothing, KrylovJL(), KrylovJL(;precs), \) ad isa AutoZygote && continue if prec === :Random prec = (args...) -> (Diagonal(randn!(similar(u0))), nothing) end - solver = NewtonRaphson(; linsolve, precs = prec(u0), linesearch) + solver = NewtonRaphson(; linsolve, linesearch) sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver) @test SciMLBase.successful_retcode(sol) @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0), - NewtonRaphson(; linsolve, precs = prec(u0)), abstol = 1e-9) + NewtonRaphson(; linsolve), abstol = 1e-9) @test (@ballocated solve!($cache)) ≤ 64 end end @@ -412,18 +409,14 @@ end @test (@ballocated solve!($cache)) < 200 end - precs = [NonlinearSolve.DEFAULT_PRECS, :Random] + precs = (A, p=nothing) -> (Diagonal(randn!(similar(A, size(A,1)))), nothing) - @testset "[IIP] u0: $(typeof(u0)) precs: $(_nameof(prec)) linsolve: $(_nameof(linsolve))" for u0 in ([ + @testset "[IIP] u0: $(typeof(u0)) linsolve: $(_nameof(linsolve))" for u0 in ([ 1.0, 1.0],), - prec in precs, - linsolve in (nothing, KrylovJL_GMRES(), \) + linsolve in (nothing, KrylovJL(), KrylovJL(;precs), \) ad isa AutoZygote && continue - if prec === :Random - prec = (args...) -> (Diagonal(randn!(similar(u0))), nothing) - end - solver = PseudoTransient(; alpha_initial = 10.0, linsolve, precs = prec) + solver = PseudoTransient(; alpha_initial = 10.0, linsolve) sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver) @test SciMLBase.successful_retcode(sol) @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) From 9463461eda59895e9bb4360766da86c77b2e5cd4 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 24 Sep 2024 19:18:54 -0400 Subject: [PATCH 2/3] don't import preconditioner --- src/NonlinearSolve.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 2e420a761..982fbcd08 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -31,8 +31,8 @@ using LinearAlgebra: LinearAlgebra, ColumnNorm, Diagonal, I, LowerTriangular, Sy UpperTriangular, axpy!, cond, diag, diagind, dot, issuccess, istril, istriu, lu, mul!, norm, pinv, tril!, triu! using LineSearches: LineSearches -using LinearSolve: LinearSolve, LUFactorization, QRFactorization, ComposePreconditioner, - InvPreconditioner, needs_concrete_A, AbstractFactorization, +using LinearSolve: LinearSolve, LUFactorization, QRFactorization, + needs_concrete_A, AbstractFactorization, DefaultAlgorithmChoice, DefaultLinearSolver using MaybeInplace: @bb using Printf: @printf From c4e6defc3c1db681ba15c47246b42008afb3db64 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Wed, 25 Sep 2024 10:37:34 -0400 Subject: [PATCH 3/3] typo --- test/core/rootfind_tests.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/test/core/rootfind_tests.jl b/test/core/rootfind_tests.jl index 518adbec5..da06109f3 100644 --- a/test/core/rootfind_tests.jl +++ b/test/core/rootfind_tests.jl @@ -70,16 +70,13 @@ end NewtonRaphson(), abstol = 1e-9) @test (@ballocated solve!($cache)) < 200 end - + precs = (A, p=nothing) -> (Diagonal(randn!(similar(A, size(A,1)))), nothing) @testset "[IIP] u0: $(typeof(u0)) linsolve: $(_nameof(linsolve))" for u0 in ([ 1.0, 1.0],), linsolve in (nothing, KrylovJL(), KrylovJL(;precs), \) ad isa AutoZygote && continue - if prec === :Random - prec = (args...) -> (Diagonal(randn!(similar(u0))), nothing) - end solver = NewtonRaphson(; linsolve, linesearch) sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver) @test SciMLBase.successful_retcode(sol)