useW-sparsity

Author: vyudu

src/systems/diffeqs/abstractodesystem.jl

@@ -294,9 +294,13 @@ function jacobian_dae_sparsity(sys::AbstractODESystem)
 end
 
 function W_sparsity(sys::AbstractODESystem) 
-    jac_sparsity = jacobian_sparsity(sys)
+    if has_tearing_state(sys)
+        jac_sparsity = jacobian_sparsity(sys)
+    else
+        jac = calculate_jacobian(sys; sparse = true)
+        jac_sparsity = sparse((!iszero).(jac))
+    end
     (n, n) = size(jac_sparsity)
-
     M = calculate_massmatrix(sys)
     M_sparsity = M isa UniformScaling ? sparse(I(n)) : sparse((!iszero).(M))
     jac_sparsity .|| M_sparsity
@@ -408,14 +412,8 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem,
 
     if sparse
         uElType = u0 === nothing ? Float64 : eltype(u0)
-        if jac
-            jac_prototype = similar(calculate_jacobian(sys; sparse), uElType)
-        else
-            jac_prototype = similar(jacobian_sparsity(sys), uElType)
-        end
         W_prototype = similar(W_sparsity(sys), uElType)
     else
-        jac_prototype = nothing
         W_prototype = nothing
     end
 
@@ -427,9 +425,8 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem,
         tgrad = _tgrad === nothing ? nothing : _tgrad,
         mass_matrix = _M,
         jac_prototype = W_prototype,
-        W_prototype = W_prototype,
         observed = observedfun,
-        sparsity = sparsity ? jacobian_sparsity(sys) : nothing,
+        sparsity = sparsity ? W_sparsity(sys) : nothing,
         analytic = analytic,
         initialization_data)
 end

src/systems/diffeqs/sdesystem.jl

@@ -645,10 +645,8 @@ function DiffEqBase.SDEFunction{iip, specialize}(sys::SDESystem, dvs = unknowns(
     M = calculate_massmatrix(sys)
     if sparse
         uElType = u0 === nothing ? Float64 : eltype(u0)
-        jac_prototype = similar(calculate_jacobian(sys; sparse), uElType)
-        W_prototype = similar(jac_prototype .+ M, uElType)
+        W_prototype = similar(W_sparsity(sys), uElType)
     else
-        jac_prototype = nothing
         W_prototype = nothing
     end
 
@@ -664,7 +662,7 @@ function DiffEqBase.SDEFunction{iip, specialize}(sys::SDESystem, dvs = unknowns(
         mass_matrix = _M,
         jac_prototype = W_prototype,
         observed = observedfun,
-        sparsity = sparsity ? jacobian_sparsity(sys) : nothing,
+        sparsity = sparsity ? W_sparsity(sys) : nothing,
         analytic = analytic,
         Wfact = _Wfact === nothing ? nothing : _Wfact,
         Wfact_t = _Wfact_t === nothing ? nothing : _Wfact_t,