Add tests

Author: avik-pal

src/broyden.jl

@@ -36,9 +36,6 @@ An implementation of `Broyden` with resetting and line search.
       + `Val(:diagonal)`: Only update the diagonal of the Jacobian. This algorithm may be
         useful for specific problems, but whether it will work may depend strongly on the
         problem.
-      + `Val(:exciting_mixing)`: Update the diagonal of the Jacobian. This is based off
-        SciPy's implementation of Broyden's method. This algorithm may be useful for
-        specific problems, but whether it will work may depend strongly on the problem.
 """
 @concrete struct GeneralBroyden{IJ, UR, CJ, AD} <: AbstractNewtonAlgorithm{CJ, AD}
     ad::AD
@@ -65,7 +62,7 @@ function GeneralBroyden(; max_resets = 100, linesearch = nothing, reset_toleranc
         init_jacobian::Val = Val(:identity), autodiff = nothing, alpha = nothing,
         update_rule = Val(:good_broyden))
     UR = _unwrap_val(update_rule)
-    @assert UR ∈ (:good_broyden, :bad_broyden, :diagonal, :exciting_mixing)
+    @assert UR ∈ (:good_broyden, :bad_broyden, :diagonal)
     IJ = _unwrap_val(init_jacobian)
     @assert IJ ∈ (:identity, :true_jacobian)
     linesearch = linesearch isa LineSearch ? linesearch : LineSearch(; method = linesearch)
@@ -167,7 +164,12 @@ function perform_step!(cache::GeneralBroydenCache{iip, IJ, UR}) where {iip, IJ,
     T = eltype(cache.u)
 
     if IJ === :true_jacobian && cache.stats.nsteps == 0
-        cache.J⁻¹ = __safe_inv(jacobian!!(cache.J⁻¹, cache)) # This allocates
+        if __isdiag(cache.J⁻¹) && cache.J⁻¹_cache !== nothing
+            cache.J⁻¹_cache = __safe_inv(jacobian!!(cache.J⁻¹_cache, cache))
+            cache.J⁻¹ = __get_diagonal!!(cache.J⁻¹, cache.J⁻¹_cache)
+        else
+            cache.J⁻¹ = __safe_inv(jacobian!!(cache.J⁻¹, cache))
+        end
     end
 
     if __isdiag(cache.J⁻¹)
@@ -195,7 +197,7 @@ function perform_step!(cache::GeneralBroydenCache{iip, IJ, UR}) where {iip, IJ,
             return nothing
         end
         if IJ === :true_jacobian
-            if __isdiag(cache.J⁻¹)
+            if __isdiag(cache.J⁻¹) && cache.J⁻¹_cache !== nothing
                 cache.J⁻¹_cache = __safe_inv(jacobian!!(cache.J⁻¹_cache, cache))
                 cache.J⁻¹ = __get_diagonal!!(cache.J⁻¹, cache.J⁻¹_cache)
             else

test/23_test_problems.jl

@@ -94,19 +94,25 @@ end
         GeneralBroyden(; init_jacobian = Val(:true_jacobian)),
         GeneralBroyden(; update_rule = Val(:bad_broyden)),
         GeneralBroyden(; init_jacobian = Val(:true_jacobian),
-            update_rule = Val(:bad_broyden)))
+            update_rule = Val(:bad_broyden)),
+        GeneralBroyden(; update_rule = Val(:diagonal)),
+        GeneralBroyden(; init_jacobian = Val(:true_jacobian), update_rule = Val(:diagonal)))
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
     broken_tests[alg_ops[1]] = [1, 5, 6, 11]
     broken_tests[alg_ops[2]] = [1, 5, 6, 8, 11, 18]
     broken_tests[alg_ops[3]] = [1, 4, 5, 6, 9, 11]
     broken_tests[alg_ops[4]] = [1, 5, 6, 8, 11]
+    broken_tests[alg_ops[5]] = [1, 3, 4, 5, 6, 8, 9, 11, 12, 15, 21]
+    broken_tests[alg_ops[6]] = [3, 4, 5, 6, 8, 9, 11, 12, 21]
 
     skip_tests = Dict(alg => Int[] for alg in alg_ops)
     skip_tests[alg_ops[1]] = [2, 22]
     skip_tests[alg_ops[2]] = [2, 22]
     skip_tests[alg_ops[3]] = [2, 22]
     skip_tests[alg_ops[4]] = [2, 22]
+    skip_tests[alg_ops[5]] = [2, 22]
+    skip_tests[alg_ops[6]] = [2, 22]
 
     test_on_library(problems, dicts, alg_ops, broken_tests; skip_tests)
 end

test/basictests.jl

@@ -692,42 +692,50 @@ end
 # --- GeneralBroyden tests ---
 
 @testset "GeneralBroyden" begin
-    function benchmark_nlsolve_oop(f, u0, p = 2.0; linesearch = nothing)
+    function benchmark_nlsolve_oop(f, u0, p = 2.0; linesearch = nothing,
+        init_jacobian = Val(:identity), update_rule = Val(:good_broyden))
         prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, GeneralBroyden(; linesearch), abstol = 1e-9)
+        return solve(prob, GeneralBroyden(; linesearch, init_jacobian, update_rule);
+            abstol = 1e-9)
     end
 
-    function benchmark_nlsolve_iip(f, u0, p = 2.0; linesearch = nothing)
+    function benchmark_nlsolve_iip(f, u0, p = 2.0; linesearch = nothing,
+            init_jacobian = Val(:identity), update_rule = Val(:good_broyden))
         prob = NonlinearProblem{true}(f, u0, p)
-        return solve(prob, GeneralBroyden(; linesearch), abstol = 1e-9)
+        return solve(prob, GeneralBroyden(; linesearch, init_jacobian, update_rule);
+            abstol = 1e-9)
     end
 
-    @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in (Static(),
+    @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad)) Init Jacobian: $(init_jacobian) Update Rule: $(update_rule)" for lsmethod in (Static(),
             StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente(),
             LiFukushimaLineSearch()),
-        ad in (AutoFiniteDiff(), AutoZygote())
+        ad in (AutoFiniteDiff(), AutoZygote()),
+        init_jacobian in (Val(:identity), Val(:true_jacobian)),
+        update_rule in (Val(:good_broyden), Val(:bad_broyden), Val(:diagonal))
 
         linesearch = LineSearch(; method = lsmethod, autodiff = ad)
         u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
 
         @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
-            sol = benchmark_nlsolve_oop(quadratic_f, u0; linesearch)
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; linesearch, update_rule,
+                init_jacobian)
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
             cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0),
-                GeneralBroyden(; linesearch), abstol = 1e-9)
+                GeneralBroyden(; linesearch, update_rule, init_jacobian), abstol = 1e-9)
             @test (@ballocated solve!($cache)) < 200
         end
 
         @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
             ad isa AutoZygote && continue
-            sol = benchmark_nlsolve_iip(quadratic_f!, u0; linesearch)
+            sol = benchmark_nlsolve_iip(quadratic_f!, u0; linesearch, update_rule,
+                init_jacobian)
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
             cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0),
-                GeneralBroyden(; linesearch), abstol = 1e-9)
+                GeneralBroyden(; linesearch, update_rule, init_jacobian), abstol = 1e-9)
             @test (@ballocated solve!($cache)) ≤ 64
         end
     end