Initialize sparse Jacobian caches with stored values#3833
Conversation
| J = similar(f.jac_prototype) | ||
| W = similar(J) | ||
| fill!(J, zero(eltype(J))) | ||
| fill!(W, zero(eltype(W))) |
There was a problem hiding this comment.
Should use zero instead of similar to start it. though actually, if it's zero then dropzeros is wrong... so maybe something like fill with ones? Benchmark and see if zeroing has adverse behavior when mixed with a sprase matrix and the default dropzeros of LinearSolve
|
@ChrisRackauckas If anything, should it not initialize with a nonzero number like ones? Initializing with zeros sounds like it would trigger the Alternatively, would it be better if OrdinaryDiffEq just sets the nonstructural zeros operator assumption to e.g. |
|
There's a lot of cases we found where the dropzeros is actually a major performance improvement, and also a numerical improvement. So we want to have it on and scan for it. But we should try to not do it on many problems, yes. The issue is that it hits init, so we should probably just init to ones. |
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
d67e68b to
8a4f275
Compare
This PR should be ignored until reviewed by @ChrisRackauckas.
Summary
nzvalentries for sparse matricesJandWkeep their sparse pattern and stored values immediately afterinitFixes #3830.
Local verification
~/.juliaup/bin/julia +1.12 --project=. -e 'using OrdinaryDiffEqRosenbrock, SparseArrays, LinearAlgebra; function f!(du,u,p,t); du .= u; end; jacproto = sparse(1.0I, 4, 4); prob = ODEProblem(ODEFunction(f!; jac_prototype=jacproto), ones(4), (0.0, 1.0)); integ = init(prob, Rodas5P()); @show nnz(integ.cache.J) nnz(integ.cache.W); @show nonzeros(integ.cache.J); @show nonzeros(integ.cache.W); @assert nnz(integ.cache.J) == 4; @assert nnz(integ.cache.W) == 4; @assert all(==(1), nonzeros(integ.cache.J)); @assert all(==(1), nonzeros(integ.cache.W))'ODEDIFFEQ_TEST_GROUP=Sparse ~/.juliaup/bin/julia +1.12 --project=lib/OrdinaryDiffEqDifferentiation -e 'using Pkg; Pkg.test(; test_args=String[])'~/.juliaup/bin/julia +1.12 -e 'using Runic; exit(Runic.main(ARGS))' -- --check lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl lib/OrdinaryDiffEqDifferentiation/test/sparse/nonfulldiagonal_sparse_tests.jlNote: repository-wide Runic check was also run with
--check --diff --verbose .; it fails on unrelated files. A clean-master check reproduced that failure and bisected it to661f605637/ PR #3777.