Fix handling of singular Jacobians in default solver

Fixes #635

Author: ChrisRackauckas

src/poly_algs.jl

@@ -34,14 +34,16 @@ function FastShortcutNonlinearPolyalg(
 ) where {T}
     start_index = 1
     common_kwargs = (; concrete_jac, linsolve, autodiff, vjp_autodiff, jvp_autodiff)
+    common_kwargs_nocj = (; linsolve, autodiff, vjp_autodiff, jvp_autodiff)
     if must_use_jacobian isa Val{true}
         if T <: Complex
             algs = (NewtonRaphson(; common_kwargs...),)
         else
             algs = (
                 NewtonRaphson(; common_kwargs...),
                 TrustRegion(; common_kwargs...),
-                TrustRegion(; common_kwargs..., radius_update_scheme = RUS.Bastin)
+                TrustRegion(; common_kwargs..., radius_update_scheme = RUS.Bastin),
+                LevenbergMarquardt(; common_kwargs_nocj...)
             )
         end
     else
@@ -61,7 +63,8 @@ function FastShortcutNonlinearPolyalg(
                     SimpleKlement(),
                     NewtonRaphson(; common_kwargs...),
                     TrustRegion(; common_kwargs...),
-                    TrustRegion(; common_kwargs..., radius_update_scheme = RUS.Bastin)
+                    TrustRegion(; common_kwargs..., radius_update_scheme = RUS.Bastin),
+                    LevenbergMarquardt(; common_kwargs_nocj...)
                 )
             end
         else
@@ -79,7 +82,8 @@ function FastShortcutNonlinearPolyalg(
                     Klement(; linsolve, autodiff),
                     NewtonRaphson(; common_kwargs...),
                     TrustRegion(; common_kwargs...),
-                    TrustRegion(; common_kwargs..., radius_update_scheme = RUS.Bastin)
+                    TrustRegion(; common_kwargs..., radius_update_scheme = RUS.Bastin),
+                    LevenbergMarquardt(; common_kwargs_nocj...)
                 )
             end
         end

test/core_tests.jl

@@ -435,4 +435,11 @@ end
     sol = solve(prob)
     @test SciMLBase.successful_retcode(sol)
     @test sol.retcode == ReturnCode.StalledSuccess
+end
+
+@testitem "Default Algorithm Singular Handling" tags=[:core] begin
+    f(u, p) = [u[1]^2 - 2u[1] + 1, sum(u)]
+    prob = NonlinearProblem(f, [1.0, 1.0])
+    sol = solve(prob)
+    @test SciMLBase.successful_retcode(sol)
 end