Merge pull request #200 from avik-pal/ap/jac_prototype

Newton Raphson used to ignore jac-prototype
SciML · Aug 19, 2023 · e93fc58 · e93fc58
2 parents cb80ed2 + 503d579
commit e93fc58
Show file tree

Hide file tree

Showing 4 changed files with 13 additions and 3 deletions.
diff --git a/.gitignore b/.gitignore
@@ -23,3 +23,5 @@ docs/site/
 # environment.
 Manifest.toml
 docs/src/assets/Project.toml
+
+.vscode
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "NonlinearSolve"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 authors = ["SciML"]
-version = "1.9.0"
+version = "1.10.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"

diff --git a/src/raphson.jl b/src/raphson.jl
@@ -106,7 +106,11 @@ end
 
 function jacobian_caches(alg::NewtonRaphson, f, u, p, ::Val{true})
     uf = JacobianWrapper(f, p)
-    J = ArrayInterface.undefmatrix(u)
+    J = if f.jac_prototype === nothing
+        ArrayInterface.undefmatrix(u)
+    else
+        f.jac_prototype
+    end
 
     linprob = LinearProblem(J, _vec(zero(u)); u0 = _vec(zero(u)))
     weight = similar(u)

diff --git a/test/sparse.jl b/test/sparse.jl
@@ -41,11 +41,15 @@ sol = solve(prob_brusselator_2d, NewtonRaphson())
 du0 = copy(u0)
 jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p), du0,
     u0)
+jac_prototype = float.(jac_sparsity)
+fill!(jac_prototype, 0)
+@test all(iszero, jac_prototype)
 
-ff = NonlinearFunction(brusselator_2d_loop; jac_prototype = float.(jac_sparsity))
+ff = NonlinearFunction(brusselator_2d_loop; jac_prototype)
 prob_brusselator_2d = NonlinearProblem(ff, u0, p)
 sol = solve(prob_brusselator_2d, NewtonRaphson())
 @test norm(sol.resid) < 1e-8
+@test !all(iszero, jac_prototype)
 
 sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = false))
 @test norm(sol.resid) < 1e-6