Use approximate sparsity detection by default

Author: avik-pal

Project.toml

@@ -1,7 +1,7 @@
 name = "NonlinearSolve"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 authors = ["SciML"]
-version = "3.0.1"
+version = "3.0.2"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -35,6 +35,7 @@ FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce"
 LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
 MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [extensions]
@@ -43,6 +44,7 @@ NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
 NonlinearSolveLeastSquaresOptimExt = "LeastSquaresOptim"
 NonlinearSolveMINPACKExt = "MINPACK"
 NonlinearSolveNLsolveExt = "NLsolve"
+NonlinearSolveSymbolicsExt = "Symbolics"
 NonlinearSolveZygoteExt = "Zygote"
 
 [compat]

docs/src/tutorials/large_systems.md

@@ -60,12 +60,14 @@ broadcast). Use `dx=dy=1/32`.
 The resulting `NonlinearProblem` definition is:
 
 ```@example ill_conditioned_nlprob
-using NonlinearSolve, LinearAlgebra, SparseArrays, LinearSolve, Symbolics
+using NonlinearSolve, LinearAlgebra, SparseArrays, LinearSolve
 
 const N = 32
 const xyd_brusselator = range(0, stop = 1, length = N)
+
 brusselator_f(x, y) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * 5.0
 limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a
+
 function brusselator_2d_loop(du, u, p)
     A, B, alpha, dx = p
     alpha = alpha / dx^2
@@ -96,6 +98,7 @@ function init_brusselator_2d(xyd)
     end
     u
 end
+
 u0 = init_brusselator_2d(xyd_brusselator)
 prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop, u0, p)
 ```
@@ -120,6 +123,26 @@ include:
   - BandedMatrix ([BandedMatrices.jl](https://github.com/JuliaLinearAlgebra/BandedMatrices.jl))
   - BlockBandedMatrix ([BlockBandedMatrices.jl](https://github.com/JuliaLinearAlgebra/BlockBandedMatrices.jl))
 
+## Approximate Sparsity Detection & Sparse Jacobians
+
+In the next section, we will discuss how to declare a sparse Jacobian and how to use
+[Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl), to compute exact sparse
+jacobians. For this to work it is important to not load `Symbolics.jl`. This is triggered
+if you pass in a sparse autodiff type such as `AutoSparseForwardDiff()`.
+
+```@example ill_conditioned_nlprob
+using BenchmarkTools # for @btime
+
+@btime solve(prob_brusselator_2d, NewtonRaphson());
+@btime solve(prob_brusselator_2d, NewtonRaphson(; autodiff = AutoSparseForwardDiff()));
+@btime solve(prob_brusselator_2d,
+    NewtonRaphson(; autodiff = AutoSparseForwardDiff(),
+        linsolve = KLUFactorization()));
+@btime solve(prob_brusselator_2d,
+    NewtonRaphson(; autodiff = AutoSparseForwardDiff(),
+        linsolve = KrylovJL_GMRES()));
+```
+
 ## Declaring a Sparse Jacobian with Automatic Sparsity Detection
 
 Jacobian sparsity is declared by the `jac_prototype` argument in the `NonlinearFunction`.
@@ -156,7 +179,6 @@ prob_brusselator_2d_sparse = NonlinearProblem(ff, u0, p)
 Now let's see how the version with sparsity compares to the version without:
 
 ```@example ill_conditioned_nlprob
-using BenchmarkTools # for @btime
 @btime solve(prob_brusselator_2d, NewtonRaphson());
 @btime solve(prob_brusselator_2d_sparse, NewtonRaphson());
 @btime solve(prob_brusselator_2d_sparse, NewtonRaphson(linsolve = KLUFactorization()));
@@ -258,10 +280,8 @@ function algebraicmultigrid2(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
     if newW === nothing || newW
         A = convert(AbstractMatrix, W)
         Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A,
-            presmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
-                1))),
-            postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
-                1)))))
+            presmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1))),
+            postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1)))))
     else
         Pl = Plprev
     end

ext/NonlinearSolveSymbolicsExt.jl

@@ -0,0 +1,7 @@
+module NonlinearSolveSymbolicsExt
+
+import NonlinearSolve, Symbolics
+
+NonlinearSolve.is_extension_loaded(::Val{:Symbolics}) = true
+
+end

src/jacobian.jl

@@ -12,7 +12,15 @@ sparsity_detection_alg(_, _) = NoSparsityDetection()
 function sparsity_detection_alg(f, ad::AbstractSparseADType)
     if f.sparsity === nothing
         if f.jac_prototype === nothing
-            return SymbolicsSparsityDetection()
+            if is_extension_loaded(Val(:Symbolics))
+                return SymbolicsSparsityDetection()
+            else
+                @warn "Symbolics.jl is not loaded and sparse AD mode $(ad) is being used. \
+                       Using approximate sparsity detection using ForwardDiff. This can \
+                       potentially fail or generate incorrect sparsity pattern for \
+                       complicated problems. Use with caution." maxlog=1
+                return ApproximateJacobianSparsity()
+            end
         else
             jac_prototype = f.jac_prototype
         end