JuliaDiff / SparseDiffTools.jl

Fast jacobian computation through sparsity exploitation and matrix coloring
MIT License
237 stars 41 forks source link

`BandedMatrix` sparsity computes only diagonal elements #200

Closed superckl closed 1 year ago

superckl commented 1 year ago

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.

ChrisRackauckas commented 1 year ago

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.

ChrisRackauckas commented 1 year ago

This is fixed by ArrayInterface v7