refactor: move stuff around a bit

Author: avik-pal

Project.toml

@@ -9,31 +9,25 @@ ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LineSearch = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-MaybeInplace = "bb5d69b7-63fc-4a16-80bd-7e42200c7bdb"
 NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
 NonlinearSolveQuasiNewton = "9a2c21bd-3a47-402d-9113-8faf9a0ee114"
 NonlinearSolveSpectralMethods = "26075421-4e9a-44e1-8bd1-420ed7ad02b2"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
-RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-SciMLJacobianOperators = "19f34311-ddf3-4b8b-af20-060888a46c0e"
-SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [weakdeps]
 FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce"
@@ -71,7 +65,6 @@ CUDA = "5.5"
 CommonSolve = "0.2.4"
 ConcreteStructs = "0.2.3"
 DiffEqBase = "6.155.3"
-DifferentiationInterface = "0.6.16"
 Enzyme = "0.13.11"
 ExplicitImports = "1.5"
 FastClosures = "0.3.2"
@@ -84,11 +77,10 @@ InteractiveUtils = "<0.0.1, 1"
 LeastSquaresOptim = "0.8.5"
 LineSearch = "0.1.4"
 LineSearches = "7.3"
-LinearAlgebra = "1.10"
-LinearSolve = "2.35"
+LinearAlgebra = "1.11.0"
+LinearSolve = "2.36.1"
 MINPACK = "1.2"
 MPI = "0.20.22"
-MaybeInplace = "0.1.4"
 NLSolvers = "0.5"
 NLsolve = "4.5"
 NaNMath = "1"
@@ -101,12 +93,9 @@ PrecompileTools = "1.2"
 Preferences = "1.4"
 Random = "1.10"
 ReTestItems = "1.24"
-RecursiveArrayTools = "3.27"
 Reexport = "1.2"
 SIAMFANLEquations = "1.0.1"
 SciMLBase = "2.54.0"
-SciMLJacobianOperators = "0.1"
-SciMLOperators = "0.3.10"
 SimpleNonlinearSolve = "2"
 SparseArrays = "1.10"
 SparseConnectivityTracer = "0.6.5"
@@ -118,7 +107,6 @@ StaticArraysCore = "1.4"
 Sundials = "4.23.1"
 SymbolicIndexingInterface = "0.3.31"
 Test = "1.10"
-TimerOutputs = "0.5.23"
 Zygote = "0.6.69"
 julia = "1.10"
 

lib/NonlinearSolveBase/src/abstract_types.jl

@@ -18,7 +18,7 @@ function reinit_self!(x::Any; kwargs...)
     return
 end
 
-function reinit_self!(stats::NLStats)
+function reinit!(stats::NLStats)
     stats.nf = 0
     stats.nsteps = 0
     stats.nfactors = 0
@@ -257,8 +257,8 @@ Abstract Type for all NonlinearSolveBase Caches.
   - `SciMLBase.set_u!(cache, u)`: set the current state.
   - `SciMLBase.reinit!(cache, u0; kwargs...)`: reinitialize the cache with the initial state
     `u0` and any additional keyword arguments.
-  - `SciMLBase.step!(cache; kwargs...)`: See [`SciMLBase.step!`](@ref) for more details.
   - `SciMLBase.isinplace(cache)`: whether or not the solver is inplace.
+  - `CommonSolve.step!(cache; kwargs...)`: See [`CommonSolve.step!`](@ref) for more details.
 
 Additionally implements `SymbolicIndexingInterface` interface Functions.
 
@@ -577,7 +577,8 @@ macro internal_caches(cType, internal_cache_names...)
             return
         end
         function NonlinearSolveBase.InternalAPI.reinit!(
-                cache::$(cType), args...; kwargs...)
+                cache::$(cType), args...; kwargs...
+        )
             $(reinit_caches...)
             $(InternalAPI.reinit_self!)(cache, args...; kwargs...)
             return

lib/NonlinearSolveBase/src/descent/damped_newton.jl

@@ -42,7 +42,7 @@ supports_trust_region(::DampedNewtonDescent) = true
     mode <: Union{Val{:normal_form}, Val{:least_squares}, Val{:simple}}
 end
 
-NonlinearSolveBase.@internal_caches DampedNewtonDescentCache :lincache :damping_fn_cache
+@internal_caches DampedNewtonDescentCache :lincache :damping_fn_cache
 
 function InternalAPI.init(
         prob::AbstractNonlinearProblem, alg::DampedNewtonDescent, J, fu, u; stats,

lib/NonlinearSolveBase/src/descent/dogleg.jl

@@ -44,7 +44,7 @@ end
     normal_form <: Union{Val{false}, Val{true}}
 end
 
-NonlinearSolveBase.@internal_caches DoglegCache :newton_cache :cauchy_cache
+@internal_caches DoglegCache :newton_cache :cauchy_cache
 
 function InternalAPI.init(
         prob::AbstractNonlinearProblem, alg::Dogleg, J, fu, u;
@@ -111,15 +111,15 @@ function InternalAPI.solve!(
 
         l_grad = cache.internalnorm(δu_cauchy)
         @bb cache.δu_cache_mul = JᵀJ × vec(δu_cauchy)
-        δuJᵀJδu = Utils.dot(cache.δu_cache_mul, cache.δu_cache_mul)
+        δuJᵀJδu = Utils.safe_dot(cache.δu_cache_mul, cache.δu_cache_mul)
     else
         δu_cauchy = InternalAPI.solve!(
             cache.cauchy_cache, J, fu, u, idx; skip_solve, kwargs...
         ).δu
         J_ = preinverted_jacobian(cache) ? inv(J) : J
         l_grad = cache.internalnorm(δu_cauchy)
         @bb cache.Jᵀδu_cache = J_ × vec(δu_cauchy)
-        δuJᵀJδu = Utils.dot(cache.Jᵀδu_cache, cache.Jᵀδu_cache)
+        δuJᵀJδu = Utils.safe_dot(cache.Jᵀδu_cache, cache.Jᵀδu_cache)
     end
     d_cauchy = (l_grad^3) / δuJᵀJδu
 

lib/NonlinearSolveBase/src/descent/geodesic_acceleration.jl

@@ -54,7 +54,7 @@ function InternalAPI.reinit_self!(cache::GeodesicAccelerationCache; p = cache.p,
     cache.last_step_accepted = false
 end
 
-NonlinearSolveBase.@internal_caches GeodesicAccelerationCache :descent_cache
+@internal_caches GeodesicAccelerationCache :descent_cache
 
 function get_velocity(cache::GeodesicAccelerationCache)
     return SciMLBase.get_du(cache.descent_cache, Val(1))

lib/NonlinearSolveBase/src/descent/newton.jl

@@ -24,7 +24,7 @@ supports_line_search(::NewtonDescent) = true
     normal_form <: Union{Val{false}, Val{true}}
 end
 
-NonlinearSolveBase.@internal_caches NewtonDescentCache :lincache
+@internal_caches NewtonDescentCache :lincache
 
 function InternalAPI.init(
         prob::AbstractNonlinearProblem, alg::NewtonDescent, J, fu, u; stats,

lib/NonlinearSolveBase/src/descent/steepest.jl

@@ -21,7 +21,7 @@ supports_line_search(::SteepestDescent) = true
     preinverted_jacobian <: Union{Val{false}, Val{true}}
 end
 
-NonlinearSolveBase.@internal_caches SteepestDescentCache :lincache
+@internal_caches SteepestDescentCache :lincache
 
 function InternalAPI.init(
         prob::AbstractNonlinearProblem, alg::SteepestDescent, J, fu, u;

lib/NonlinearSolveBase/src/solve.jl

@@ -2,13 +2,13 @@ function SciMLBase.__solve(
         prob::AbstractNonlinearProblem, alg::AbstractNonlinearSolveAlgorithm, args...;
         kwargs...
 )
-    cache = SciMLBase.init(prob, alg, args...; kwargs...)
-    return SciMLBase.solve!(cache)
+    cache = SciMLBase.__init(prob, alg, args...; kwargs...)
+    return CommonSolve.solve!(cache)
 end
 
 function CommonSolve.solve!(cache::AbstractNonlinearSolveCache)
     while not_terminated(cache)
-        SciMLBase.step!(cache)
+        CommonSolve.step!(cache)
     end
 
     # The solver might have set a different `retcode`
@@ -47,7 +47,7 @@ Performs one step of the nonlinear solver.
     respectively. For algorithms that don't use jacobian information, this keyword is
     ignored with a one-time warning.
 """
-function SciMLBase.step!(cache::AbstractNonlinearSolveCache, args...; kwargs...)
+function CommonSolve.step!(cache::AbstractNonlinearSolveCache, args...; kwargs...)
     not_terminated(cache) || return
 
     has_time_limit(cache) && (time_start = time())

lib/NonlinearSolveBase/src/utils.jl

@@ -2,7 +2,7 @@ module Utils
 
 using ArrayInterface: ArrayInterface
 using FastClosures: @closure
-using LinearAlgebra: LinearAlgebra, Diagonal, Symmetric, norm, dot, cond, diagind, pinv, inv
+using LinearAlgebra: LinearAlgebra, Diagonal, Symmetric, norm, dot, cond, diagind, pinv
 using MaybeInplace: @bb
 using RecursiveArrayTools: AbstractVectorOfArray, ArrayPartition
 using SciMLOperators: AbstractSciMLOperator
@@ -117,6 +117,7 @@ is_default_value(::Any, ::Symbol, ::Nothing) = true
 is_default_value(::Any, ::Symbol, ::Missing) = true
 is_default_value(::Any, ::Symbol, val::Int) = val == typemax(typeof(val))
 is_default_value(::Any, ::Symbol, ::Any) = false
+is_default_value(::Any, ::Any, ::Any) = false
 
 maybe_symmetric(x) = Symmetric(x)
 maybe_symmetric(x::Number) = x

lib/NonlinearSolveFirstOrder/Project.toml

@@ -19,6 +19,7 @@ NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLJacobianOperators = "19f34311-ddf3-4b8b-af20-060888a46c0e"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 

lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl

@@ -22,6 +22,7 @@ using NonlinearSolveBase: NonlinearSolveBase, AbstractNonlinearSolveAlgorithm,
                           NewtonDescent, DampedNewtonDescent, GeodesicAcceleration,
                           Dogleg
 using SciMLBase: SciMLBase, AbstractNonlinearProblem, NLStats, ReturnCode, NonlinearFunction
+using SciMLJacobianOperators: VecJacOperator, JacVecOperator, StatefulJacobianOperator
 using Setfield: @set!
 using StaticArraysCore: SArray
 
@@ -41,9 +42,8 @@ include("solve.jl")
         @__DIR__, "..", "..", "..", "common", "nlls_problem_workloads.jl"
     ))
 
-    # XXX: TrustRegion
-    nlp_algs = [NewtonRaphson(), LevenbergMarquardt()]
-    nlls_algs = [GaussNewton(), LevenbergMarquardt()]
+    nlp_algs = [NewtonRaphson(), TrustRegion(), LevenbergMarquardt()]
+    nlls_algs = [GaussNewton(), TrustRegion(), LevenbergMarquardt()]
 
     @compile_workload begin
         for prob in nonlinear_problems, alg in nlp_algs
@@ -61,6 +61,8 @@ end
 export NewtonRaphson, PseudoTransient
 export GaussNewton, LevenbergMarquardt, TrustRegion
 
+export RadiusUpdateSchemes
+
 export GeneralizedFirstOrderAlgorithm
 
 end