Safe condition number

Author: avik-pal

src/broyden.jl

@@ -50,8 +50,8 @@ end
 
 function __alg_print_modifiers(alg::GeneralBroyden{IJ, UR}) where {IJ, UR}
     modifiers = String[]
-    IJ !== :identity && push!(modifiers, "init_jacobian = :$(IJ)")
-    UR !== :good_broyden && push!(modifiers, "update_rule = :$(UR)")
+    IJ !== :identity && push!(modifiers, "init_jacobian = Val(:$(IJ))")
+    UR !== :good_broyden && push!(modifiers, "update_rule = Val(:$(UR))")
     alg.alpha !== nothing && push!(modifiers, "alpha = $(alg.alpha)")
     return modifiers
 end
@@ -141,8 +141,8 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::GeneralBroyd
 
     abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fu, u,
         termination_condition)
-    trace = init_nonlinearsolve_trace(alg, u, fu, J⁻¹, du; uses_jac_inverse = Val(true),
-        kwargs...)
+    trace = init_nonlinearsolve_trace(alg, u, fu, ApplyArray(__zero, J⁻¹), du;
+        uses_jac_inverse = Val(true), kwargs...)
 
     return GeneralBroydenCache{iip, IJ, UR}(f, alg, u, u_cache, du, fu, fu_cache, dfu, p,
         uf, J⁻¹, J⁻¹dfu, inv_alpha, alg.alpha, false, 0, alg.max_resets, maxiters,

src/klement.jl

@@ -49,7 +49,7 @@ end
 
 function __alg_print_modifiers(alg::GeneralKlement{IJ}) where {IJ}
     modifiers = String[]
-    IJ !== :identity && push!(modifiers, "init_jacobian = :$(IJ)")
+    IJ !== :identity && push!(modifiers, "init_jacobian = Val(:$(IJ))")
     alg.alpha !== nothing && push!(modifiers, "alpha = $(alg.alpha)")
     return modifiers
 end
@@ -231,7 +231,7 @@ function perform_step!(cache::GeneralKlementCache{iip, IJ}) where {iip, IJ}
     else
         # Klement Updates to the Full Jacobian don't work for most problems, we should
         # probably be using the Broyden Update Rule here
-        @bb @. cache.J_cache = cache.J' ^ 2
+        @bb @. cache.J_cache = cache.J'^2
         @bb @. cache.Jdu = cache.du^2
         @bb cache.Jdu_cache = cache.J_cache × vec(cache.Jdu)
         @bb cache.Jdu = cache.J × vec(cache.du)

src/trace.jl

@@ -174,7 +174,7 @@ function init_nonlinearsolve_trace(alg, ::Val{show_trace},
         print("\nAlgorithm: ")
         Base.printstyled(alg, "\n\n"; color = :green, bold = true)
     end
-    J_ = uses_jac_inverse ? (trace_level isa TraceMinimal ? J : inv(J)) : J
+    J_ = uses_jac_inverse ? (trace_level isa TraceMinimal ? J : __safe_inv(J)) : J
     history = __init_trace_history(Val{show_trace}(), trace_level, Val{store_trace}(), u,
         fu, J_, δu)
     return NonlinearSolveTrace{show_trace, store_trace}(history, trace_level)
@@ -231,7 +231,7 @@ function update_trace!(cache::AbstractNonlinearSolveCache, α = true)
                 nothing, cache.du, α)
         else
             update_trace!(trace, cache.stats.nsteps + 1, get_u(cache), get_fu(cache),
-                ApplyArray(inv, J_inv), cache.du, α)
+                ApplyArray(__safe_inv, J_inv), cache.du, α)
         end
     else
         update_trace!(trace, cache.stats.nsteps + 1, get_u(cache), get_fu(cache), J,

src/utils.jl

@@ -528,6 +528,7 @@ end
 @inline __diag(x::Number) = x
 
 # Safe Inverse: Try to use `inv` but if lu fails use `pinv`
+__safe_inv(A::AbstractMatrix) = pinv(A)
 function __safe_inv(A::StridedMatrix{T}) where {T}
     LinearAlgebra.checksquare(A)
     if istriu(A)