Fix aliasing issue

SciML · Nov 21, 2023 · 29d06b7 · 29d06b7
1 parent b79228d
commit 29d06b7
Show file tree

Hide file tree

Showing 3 changed files with 24 additions and 9 deletions.
diff --git a/src/levenberg.jl b/src/levenberg.jl
@@ -366,9 +366,10 @@ function perform_step!(cache::LevenbergMarquardtCache{false, fastls}) where {fas
         if linsolve === nothing
             cache.v = -cache.mat_tmp \ (J' * fu1)
         else
-            linres = dolinsolve(alg.precs, linsolve; A = -__maybe_symmetric(cache.mat_tmp),
+            linres = dolinsolve(alg.precs, linsolve; A = __maybe_symmetric(cache.mat_tmp),
                 b = _vec(J' * _vec(fu1)), linu = _vec(cache.v), p, reltol = cache.abstol)
             cache.linsolve = linres.cache
+            cache.v .*= -1
         end
     end
 
@@ -384,9 +385,11 @@ function perform_step!(cache::LevenbergMarquardtCache{false, fastls}) where {fas
         if linsolve === nothing
             cache.a = -cache.mat_tmp \ _vec(J' * rhs_term)
         else
-            linres = dolinsolve(alg.precs, linsolve; b = _mutable(_vec(J' * rhs_term)),
-                linu = _vec(cache.a), p, reltol = cache.abstol)
+            linres = dolinsolve(alg.precs, linsolve; A = __maybe_symmetric(cache.mat_tmp),
+                b = _mutable(_vec(J' * rhs_term)), linu = _vec(cache.a), p,
+                reltol = cache.abstol, reuse_A_if_factorization = true)
             cache.linsolve = linres.cache
+            cache.a .*= -1
         end
     end
     cache.stats.nsolve += 1

diff --git a/src/utils.jl b/src/utils.jl
@@ -82,16 +82,28 @@ end
 DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, cachedata) = nothing, nothing
 
 function dolinsolve(precs::P, linsolve; A = nothing, linu = nothing, b = nothing,
-        du = nothing, u = nothing, p = nothing, t = nothing, weight = nothing,
-        cachedata = nothing, reltol = nothing) where {P}
-    A !== nothing && (linsolve.A = A)
+        du = nothing, p = nothing, weight = nothing, cachedata = nothing, reltol = nothing,
+        reuse_A_if_factorization = false) where {P}
+    # Some Algorithms would reuse factorization but it causes the cache to not reset in
+    # certain cases
+    if A !== nothing
+        alg = linsolve.alg
+        if (alg isa LinearSolve.AbstractFactorization) ||
+           (alg isa LinearSolve.DefaultLinearSolver && !(alg ==
+              LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.KrylovJL_GMRES)))
+            # Factorization Algorithm
+            !reuse_A_if_factorization && (linsolve.A = A)
+        else
+            linsolve.A = A
+        end
+    end
     b !== nothing && (linsolve.b = b)
     linu !== nothing && (linsolve.u = linu)
 
     Plprev = linsolve.Pl isa ComposePreconditioner ? linsolve.Pl.outer : linsolve.Pl
     Prprev = linsolve.Pr isa ComposePreconditioner ? linsolve.Pr.outer : linsolve.Pr
 
-    _Pl, _Pr = precs(linsolve.A, du, u, p, nothing, A !== nothing, Plprev, Prprev,
+    _Pl, _Pr = precs(linsolve.A, du, linu, p, nothing, A !== nothing, Plprev, Prprev,
         cachedata)
     if (_Pl !== nothing || _Pr !== nothing)
         _weight = weight === nothing ?

diff --git a/test/infeasible.jl b/test/infeasible.jl
@@ -29,8 +29,8 @@ function f1(u, p)
     v_x = 8.550491684548064e-12 + u[1]
     v_y = 6631.60076191005 + u[2]
     v_z = 3600.665431405663 + u[3]
-    r = @SVector [x, y, z]
-    v = @SVector [v_x, v_y, v_z]
+    r = [x, y, z]
+    v = [v_x, v_y, v_z]
     h = cross(r, v)
     ev = cross(v, h) / μ - r / norm(r)
     i = acos(h[3] / norm(h))