Fix TR scalars

SciML · Oct 27, 2023 · 2d4a88d · 2d4a88d
1 parent 9c50e5f
commit 2d4a88d
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 7 deletions.
diff --git a/src/default.jl b/src/default.jl
@@ -127,8 +127,7 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip},
     alg::FastShortcutNonlinearPolyalg, args...; kwargs...) where {uType, iip}
     @unpack adkwargs, linsolve, precs = alg
 
-    algs = (
-        GeneralKlement(; linsolve, precs),
+    algs = (GeneralKlement(; linsolve, precs),
         GeneralBroyden(),
         NewtonRaphson(; linsolve, precs, adkwargs...),
         NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...),

diff --git a/src/trustRegion.jl b/src/trustRegion.jl
@@ -388,11 +388,16 @@ function perform_step!(cache::TrustRegionCache{false})
         cache.g = _restructure(fu, J' * _vec(fu))
         cache.stats.njacs += 1
 
-        # do not use A = cache.H, b = _vec(cache.g) since it is equivalent
-        # to  A = cache.J, b = _vec(fu) as long as the Jacobian is non-singular
-        linres = dolinsolve(cache.alg.precs, cache.linsolve, A = cache.J, b = -_vec(fu),
-            linu = _vec(cache.u_gauss_newton), p = p, reltol = cache.abstol)
-        cache.linsolve = linres.cache
+        if cache.linsolve === nothing
+            # Scalar
+            cache.u_gauss_newton = cache.H \ cache.g
+        else
+            # do not use A = cache.H, b = _vec(cache.g) since it is equivalent
+            # to  A = cache.J, b = _vec(fu) as long as the Jacobian is non-singular
+            linres = dolinsolve(cache.alg.precs, cache.linsolve, A = cache.J, b = -_vec(fu),
+                linu = _vec(cache.u_gauss_newton), p = p, reltol = cache.abstol)
+            cache.linsolve = linres.cache
+        end
     end
 
     # Compute the Newton step.