SparseDiffTools.jl
SparseDiffTools.jl copied to clipboard
`BandedMatrix` sparsity computes only diagonal elements
When using a BandedMatrix
as the sparsity pattern for a Jacobian, only the diagonal elements are computed:
using BandedMatrices, SparseArrays, SparseDiffTools, LinearAlgebra
A = spdiagm(0 => rand(5), -1 => rand(4), 1 => rand(4))
rand_ode(du, u) = mul!(du, A, u)
u0 = ones(size(A, 1))
jac_proto = similar(A)
banded_jac_proto = BandedMatrix(A)
J = forwarddiff_color_jacobian!(similar(A), rand_ode, u0, sparsity = jac_proto)
banded_J = forwarddiff_color_jacobian!(similar(A), rand_ode, u0, sparsity = banded_jac_proto)
display(A.-J)
display(A.-banded_J)
display(matrix_colors(jac_proto))
display(matrix_colors(banded_jac_proto))
5×5 SparseMatrixCSC{Float64, Int64} with 0 stored entries:
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
5×5 SparseMatrixCSC{Float64, Int64} with 8 stored entries:
⋅ 0.145873 ⋅ ⋅ ⋅
0.104158 ⋅ 0.992912 ⋅ ⋅
⋅ 0.438004 ⋅ 0.821486 ⋅
⋅ ⋅ 0.828545 ⋅ 0.594859
⋅ ⋅ ⋅ 0.522268 ⋅
5-element Vector{Int64}:
1
2
3
1
2
5-element Vector{Int64}:
1
2
3
1
2
These issues do not occur if using ArrayInterfaceBandedMatrices
is included.
There is a big notice in the README about this kind of thing: https://github.com/JuliaDiff/SparseDiffTools.jl#note-about-sparse-differentiation-of-gpuarrays-bandedmatrices-and-blockbandedmatrices , though I am surprised that this case hits some kind of bad fallback for the matrix decompression instead of just erroring.
This is fixed by ArrayInterface v7