Skip to content

prepare_user_sparsity: unwrap operator mass matrices before findall (#2929)#3944

Open
singhharsh1708 wants to merge 1 commit into
SciML:masterfrom
singhharsh1708:fix-2929-matrixoperator-massmatrix-sparsity
Open

prepare_user_sparsity: unwrap operator mass matrices before findall (#2929)#3944
singhharsh1708 wants to merge 1 commit into
SciML:masterfrom
singhharsh1708:fix-2929-matrixoperator-massmatrix-sparsity

Conversation

@singhharsh1708

Copy link
Copy Markdown
Contributor

Partially addresses #2929.

Bug

prepare_user_sparsity (lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl) seeds the sparsity pattern from the mass matrix with:

idxs = findall(!iszero, prob.f.mass_matrix)
for idx in idxs
    sparsity[idx] = prob.f.mass_matrix[idx]
    ...

when the mass matrix isn't a UniformScaling. findall/getindex need a concrete matrix; a SciMLOperator mass matrix (e.g. MatrixOperator) supports neither directly:

using OrdinaryDiffEq, SparseArrays, SciMLOperators
M = sparse([1.0 2.0; 0.0 1.0])
mass_matrix = MatrixOperator(M; update_func = (A, u, p, t) -> M)
jac_prototype = sparse([1.0 0.0; 0.0 1.0])
prob = ODEProblem(ODEFunction((du,u,p,t)->du.=u; mass_matrix, jac_prototype), ones(2), (0.0, 1.0))
solve(prob, ImplicitEuler())
# MethodError: no method matching keys(::MatrixOperator{...})

matching #2929's original report exactly.

Fix

Unwrap the operator to its concrete matrix first: convert(AbstractMatrix, mm). This only reads the operator's currently-cached array (SciMLOperators.matrix.jl's convert(::Type{AbstractMatrix}, L::MatrixOperator) = convert(AbstractMatrix, L.A)), not a live re-evaluation, so no update_coefficients! call is involved.

What this doesn't fix

While tracing #2929's repro I found the baseline case (no jac_prototype at all) has separately regressed since the issue was filed — it now fails with a different error, setfield!: immutable struct of type MatrixOperator cannot be changed. That's entirely inside SciMLOperators.jl: update_coefficients! on a MatrixOperator built with only an out-of-place update_func (no update_func!) does a raw L.A = ... field assignment on what's declared an immutable struct — update_coefficients (no bang) handles the same case correctly via @reset. Root-caused to src/matrix.jl:227 in SciMLOperators v1.24.3, confirmed with a minimal SciMLOperators-only reproducer (no OrdinaryDiffEq involved). That's a separate upstream bug I can't fix here; will report it in SciMLOperators.jl. #2929's repro won't fully pass end-to-end until that's also fixed, but this PR fixes the part that's actually ours.

Tests

New prepare_user_sparsity_tests.jl: MatrixOperator mass matrix (fails on master with the keys MethodError, passes here with the correct merged sparsity pattern), plus regression checks that the UniformScaling and plain-sparse-matrix paths are unchanged.

OrdinaryDiffEqDifferentiation Core suite passes (39/39, including 6 new).

AI Disclosure

Claude assisted with this work.

@singhharsh1708
singhharsh1708 force-pushed the fix-2929-matrixoperator-massmatrix-sparsity branch from 39aa6be to dd6b78a Compare July 18, 2026 07:33
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 SciML#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 SciML#2929's repro is a separate SciMLOperators.jl
bug, not fixable here).
@singhharsh1708
singhharsh1708 force-pushed the fix-2929-matrixoperator-massmatrix-sparsity branch from dd6b78a to 9410f9b Compare July 21, 2026 14:31
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant