Skip to content

Commit 13aa479

Browse files
prepare_user_sparsity: unwrap operator mass matrices before findall
findall(!iszero, mm) and mm[idx] need a concrete matrix. A SciMLOperator mass matrix (e.g. MatrixOperator) supports neither, so a sparse jac_prototype + operator mass matrix crashed with MethodError: no method matching keys(::MatrixOperator...). Partially addresses #2929. convert(AbstractMatrix, mm) only reads the operator's currently cached array, not a live re-evaluation, so this is safe even for operators whose in-place update_coefficients! path is broken upstream (that part of #2929's repro is a separate SciMLOperators.jl bug, not fixable here).
1 parent 2111422 commit 13aa479

3 files changed

Lines changed: 58 additions & 3 deletions

File tree

lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -123,14 +123,20 @@ function prepare_user_sparsity(ad_alg, prob)
123123
@. @view(jac_prototype[idxs]) = 1
124124
end
125125
else
126-
idxs = findall(!iszero, prob.f.mass_matrix)
126+
# `findall`/`getindex` need a concrete matrix; a SciMLOperator mass
127+
# matrix (e.g. `MatrixOperator`) doesn't support either directly.
128+
# `convert` here only reads the operator's currently-cached array,
129+
# not a live re-evaluation (no `update_coefficients!` call).
130+
mm = prob.f.mass_matrix
131+
mm = mm isa AbstractSciMLOperator ? convert(AbstractMatrix, mm) : mm
132+
idxs = findall(!iszero, mm)
127133
for idx in idxs
128-
sparsity[idx] = prob.f.mass_matrix[idx]
134+
sparsity[idx] = mm[idx]
129135
end
130136

131137
if !isnothing(jac_prototype)
132138
for idx in idxs
133-
jac_prototype[idx] = prob.f.mass_matrix[idx]
139+
jac_prototype[idx] = mm[idx]
134140
end
135141
end
136142
end
Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
using OrdinaryDiffEqDifferentiation
2+
using SparseArrays
3+
using SciMLOperators
4+
using SciMLBase
5+
using ADTypes
6+
using LinearAlgebra
7+
using Test
8+
9+
f!(du, u, p, t) = (du .= u; nothing)
10+
ad_alg = AutoForwardDiff()
11+
12+
# `findall`/`getindex` on the mass matrix need a concrete matrix; a SciMLOperator
13+
# (e.g. `MatrixOperator`) supports neither directly, so it must be unwrapped first.
14+
@testset "MatrixOperator mass matrix" begin
15+
M = sparse([1.0 2.0; 0.0 1.0])
16+
mass_matrix = MatrixOperator(M; update_func = (A, u, p, t) -> M)
17+
jac_prototype = sparse([1.0 0.0; 0.0 1.0])
18+
sparsity = copy(jac_prototype)
19+
odef = ODEFunction(f!; mass_matrix, jac_prototype, sparsity)
20+
prob = ODEProblem(odef, ones(2), (0.0, 1.0))
21+
22+
OrdinaryDiffEqDifferentiation.prepare_user_sparsity(ad_alg, prob)
23+
@test Matrix(prob.f.sparsity) == Matrix(M)
24+
@test Matrix(prob.f.jac_prototype) == Matrix(M)
25+
end
26+
27+
@testset "UniformScaling mass matrix (unchanged path)" begin
28+
jac_prototype = sparse([1.0 0.0; 0.0 1.0])
29+
sparsity = copy(jac_prototype)
30+
odef = ODEFunction(f!; jac_prototype, sparsity)
31+
prob = ODEProblem(odef, ones(2), (0.0, 1.0))
32+
33+
OrdinaryDiffEqDifferentiation.prepare_user_sparsity(ad_alg, prob)
34+
@test Matrix(prob.f.sparsity) == Matrix(jac_prototype)
35+
@test Matrix(prob.f.jac_prototype) == Matrix(jac_prototype)
36+
end
37+
38+
@testset "Plain sparse mass matrix (unchanged path)" begin
39+
M = sparse([1.0 2.0; 0.0 1.0])
40+
jac_prototype = sparse([1.0 0.0; 0.0 1.0])
41+
sparsity = copy(jac_prototype)
42+
odef = ODEFunction(f!; mass_matrix = M, jac_prototype, sparsity)
43+
prob = ODEProblem(odef, ones(2), (0.0, 1.0))
44+
45+
OrdinaryDiffEqDifferentiation.prepare_user_sparsity(ad_alg, prob)
46+
@test Matrix(prob.f.sparsity) == Matrix(M)
47+
@test Matrix(prob.f.jac_prototype) == Matrix(M)
48+
end

lib/OrdinaryDiffEqDifferentiation/test/runtests.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ end
3030
# Run functional tests
3131
if TEST_GROUP ("QA", "Sparse", "ModelingToolkit")
3232
@time @safetestset "DAE jacobian2W sparse" include("dae_jacobian2w_sparse_tests.jl")
33+
@time @safetestset "prepare_user_sparsity mass matrix" include("prepare_user_sparsity_tests.jl")
3334
@time @safetestset "nzval helpers" include("nzval_helpers_tests.jl")
3435
@time @safetestset "prepare_sparse_jac!" include("prepare_sparse_jac_tests.jl")
3536
@time @safetestset "OOP J_t Tracking" include("oop_jt_tracking_test.jl")

0 commit comments

Comments
 (0)