refactor: make nothing the default preconditioner

Author: avik-pal

lib/NonlinearSolveBase/Project.toml

@@ -23,6 +23,7 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLJacobianOperators = "19f34311-ddf3-4b8b-af20-060888a46c0e"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [weakdeps]
@@ -67,6 +68,7 @@ SciMLOperators = "0.3.10"
 SparseArrays = "1.10"
 SparseMatrixColorings = "0.4.8"
 StaticArraysCore = "1.4"
+SymbolicIndexingInterface = "0.3.31"
 Test = "1.10"
 TimerOutputs = "0.5.23"
 julia = "1.10"

lib/NonlinearSolveBase/src/abstract_types.jl

@@ -3,6 +3,7 @@ module InternalAPI
 function init end
 function solve! end
 function reinit! end
+function step! end
 
 end
 
@@ -134,6 +135,8 @@ for fname in (:preinverted_jacobian, :normal_form)
     end
 end
 
+
+
 """
     AbstractDampingFunction
 
@@ -145,7 +148,7 @@ Abstract Type for Damping Functions in DampedNewton.
 InternalAPI.init(
     prob::AbstractNonlinearProblem, f::AbstractDampingFunction, initial_damping,
     J, fu, u, args...;
-    internalnorm::F = L2_NORM, kwargs...
+    internalnorm = L2_NORM, kwargs...
 )::AbstractDampingFunctionCache
 ```
 

src/NonlinearSolve.jl

@@ -25,16 +25,15 @@ using NonlinearSolveBase: NonlinearSolveBase,
                           AbstractSafeBestNonlinearTerminationMode,
                           select_forward_mode_autodiff, select_reverse_mode_autodiff,
                           select_jacobian_autodiff,
-                          construct_linear_solver, construct_jacobian_cache,
+                          construct_jacobian_cache,
                           DescentResult,
                           SteepestDescent, NewtonDescent, DampedNewtonDescent, Dogleg,
                           GeodesicAcceleration,
                           reset_timer!, @static_timeit
 
 # XXX: Remove
 import NonlinearSolveBase: InternalAPI, concrete_jac, supports_line_search,
-                           supports_trust_region, set_du!, last_step_accepted,
-                           get_linear_solver,
+                           supports_trust_region, last_step_accepted, get_linear_solver,
                            AbstractDampingFunction, AbstractDampingFunctionCache,
                            requires_normal_form_jacobian, requires_normal_form_rhs,
                            returns_norm_form_damping, get_timer_output
@@ -59,8 +58,7 @@ using ADTypes: ADTypes, AbstractADType, AutoFiniteDiff, AutoForwardDiff,
 using DifferentiationInterface: DifferentiationInterface
 using FiniteDiff: FiniteDiff
 using ForwardDiff: ForwardDiff, Dual
-using SciMLJacobianOperators: AbstractJacobianOperator, VecJacOperator,
-                              JacVecOperator, StatefulJacobianOperator
+using SciMLJacobianOperators: VecJacOperator, JacVecOperator, StatefulJacobianOperator
 
 ## Sparse AD Support
 using SparseArrays: AbstractSparseMatrix, SparseMatrixCSC

src/algorithms/gauss_newton.jl

@@ -1,13 +1,13 @@
 """
-    GaussNewton(; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS,
+    GaussNewton(; concrete_jac = nothing, linsolve = nothing, precs = nothing,
         linesearch = nothing, vjp_autodiff = nothing, autodiff = nothing,
         jvp_autodiff = nothing)
 
 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, precs = nothing,
         linesearch = nothing, vjp_autodiff = nothing, autodiff = nothing,
         jvp_autodiff = nothing)
     return GeneralizedFirstOrderAlgorithm{concrete_jac, :GaussNewton}(; linesearch,

src/algorithms/klement.jl

@@ -1,6 +1,6 @@
 """
     Klement(; max_resets = 100, linsolve = nothing, linesearch = nothing,
-        precs = DEFAULT_PRECS, alpha = nothing, init_jacobian::Val = Val(:identity),
+        precs = nothing, 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 = nothing, precs = DEFAULT_PRECS,
+        linesearch = nothing, precs = nothing,
         autodiff = nothing, init_jacobian::Val = Val(:identity))
     initialization = klement_init(init_jacobian, autodiff, alpha)
     CJ = init_jacobian isa Val{:true_jacobian} ||

src/algorithms/levenberg_marquardt.jl

@@ -1,6 +1,6 @@
 """
     LevenbergMarquardt(; linsolve = nothing,
-        precs = DEFAULT_PRECS, damping_initial::Real = 1.0, α_geodesic::Real = 0.75,
+        precs = 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, autodiff = nothing,
         min_damping_D::Real = 1e-8, disable_geodesic = Val(false), vjp_autodiff = nothing,
@@ -32,7 +32,7 @@ For the remaining arguments, see [`GeodesicAcceleration`](@ref) and
 [`NonlinearSolve.LevenbergMarquardtTrustRegion`](@ref) documentations.
 """
 function LevenbergMarquardt(;
-        linsolve = nothing, precs = DEFAULT_PRECS, damping_initial::Real = 1.0,
+        linsolve = nothing, precs = 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, min_damping_D::Real = 1e-8, disable_geodesic = False,

src/algorithms/pseudo_transient.jl

@@ -1,6 +1,6 @@
 """
     PseudoTransient(; concrete_jac = nothing, linsolve = nothing,
-        linesearch = NoLineSearch(), precs = DEFAULT_PRECS, autodiff = nothing,
+        linesearch = NoLineSearch(), precs = nothing, autodiff = nothing,
         jvp_autodiff = nothing, vjp_autodiff = nothing)
 
 An implementation of PseudoTransient Method [coffey2003pseudotransient](@cite) that is used
@@ -17,7 +17,7 @@ This implementation specifically uses "switched evolution relaxation"
 """
 function PseudoTransient(;
         concrete_jac = nothing, linsolve = nothing, linesearch = nothing,
-        precs = DEFAULT_PRECS, alpha_initial = 1e-3, autodiff = nothing,
+        precs = nothing, alpha_initial = 1e-3, autodiff = nothing,
         jvp_autodiff = nothing, vjp_autodiff = nothing)
     descent = DampedNewtonDescent(; linsolve, precs, initial_damping = alpha_initial,
         damping_fn = SwitchedEvolutionRelaxation())

src/algorithms/raphson.jl

@@ -1,14 +1,14 @@
 """
     NewtonRaphson(; concrete_jac = nothing, linsolve = nothing, linesearch = missing,
-        precs = DEFAULT_PRECS, autodiff = nothing, vjp_autodiff = nothing,
+        precs = nothing, autodiff = nothing, vjp_autodiff = nothing,
         jvp_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 = nothing,
-        precs = DEFAULT_PRECS, autodiff = nothing, vjp_autodiff = nothing,
+        precs = nothing, autodiff = nothing, vjp_autodiff = nothing,
         jvp_autodiff = nothing)
     return GeneralizedFirstOrderAlgorithm{concrete_jac, :NewtonRaphson}(; linesearch,
         descent = NewtonDescent(; linsolve, precs), autodiff, vjp_autodiff, jvp_autodiff)

src/algorithms/trust_region.jl

@@ -1,5 +1,5 @@
 """
-    TrustRegion(; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS,
+    TrustRegion(; concrete_jac = nothing, linsolve = nothing, precs = 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,
@@ -20,7 +20,7 @@ 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, precs = 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,

src/default.jl

@@ -338,7 +338,7 @@ end
 
 """
     RobustMultiNewton(::Type{T} = Float64; concrete_jac = nothing, linsolve = nothing,
-        precs = DEFAULT_PRECS, autodiff = nothing)
+        precs = 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,7 +354,7 @@ 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}
+        precs = nothing, autodiff = nothing) where {T}
     if __is_complex(T)
         # Let's atleast have something here for complex numbers
         algs = (NewtonRaphson(; concrete_jac, linsolve, precs, autodiff),)
@@ -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, precs = nothing, must_use_jacobian::Val = Val(false),
         prefer_simplenonlinearsolve::Val{SA} = Val(false), autodiff = nothing,
         u0_len::Union{Int, Nothing} = nothing) where {T}
 
@@ -395,7 +395,7 @@ 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),
+        precs = nothing, 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
@@ -457,7 +457,7 @@ end
 
 """
     FastShortcutNLLSPolyalg(::Type{T} = Float64; concrete_jac = nothing, linsolve = nothing,
-        precs = DEFAULT_PRECS, autodiff = nothing, kwargs...)
+        precs = nothing, 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.
@@ -469,7 +469,7 @@ 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}
+        precs = nothing, autodiff = nothing, kwargs...) where {T}
     if __is_complex(T)
         algs = (GaussNewton(; concrete_jac, linsolve, precs, autodiff, kwargs...),
             LevenbergMarquardt(;

src/utils.jl

@@ -1,6 +1,3 @@
-# Defaults
-@inline DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, cachedata) = nothing, nothing
-
 # Helper  Functions
 @generated function __getproperty(s::S, ::Val{X}) where {S, X}
     hasfield(S, X) && return :(s.$X)

test/core/rootfind_tests.jl

@@ -96,7 +96,7 @@ end
             @test (@ballocated solve!($cache)) < 200
         end
 
-        precs = [(u0) -> NonlinearSolve.DEFAULT_PRECS,
+        precs = [(u0) -> nothing,
             u0 -> ((args...) -> (Diagonal(rand!(similar(u0))), nothing))]
 
         @testset "[IIP] u0: $(typeof(u0)) precs: $(_nameof(prec)) linsolve: $(_nameof(linsolve))" for u0 in ([
@@ -437,7 +437,7 @@ end
             @test (@ballocated solve!($cache)) < 200
         end
 
-        precs = [NonlinearSolve.DEFAULT_PRECS, :Random]
+        precs = [nothing, :Random]
 
         @testset "[IIP] u0: $(typeof(u0)) precs: $(_nameof(prec)) linsolve: $(_nameof(linsolve))" for u0 in ([
                 1.0, 1.0],),