Error with jac_prototype in Mass Matrix Form of ODEFunction (svd! on Sparse Matrix)

[SanjeevKhare]

Posted this on discourse but reposting here in the hope the someone might be able to help me.

I’m solving a DAE system using the mass matrix form. The issue arises when I include jac_prototype in the ODEFunction. If I define ODEFunction without jac_prototype, everything works fine. However, when I include jac_prototype, which is a sparse matrix, I get the following error during solve:

ERROR: MethodError: no method matching svd!(::SparseMatrixCSC{Float64, Int64}; full::Bool, alg::LinearAlgebra.DivideAndConquer)
The function `svd!` exists, but no method is defined for this combination of argument types.

This is my MWE

using LinearAlgebra
using SparseArrays
using ADTypes, SparseConnectivityTracer
# import DifferentialEquations as DE
using DifferentialEquations

function laplacian(n)
    main_diag = -2.0 * ones(n)
    off_diag = ones(n - 1)
    return spdiagm(0 => main_diag, -1 => off_diag, 1 => off_diag)
end

function landau_khalatnikov_debug!(du, u, p)
    poisson_residue_debug!(du, u, p)
    return nothing
end

function poisson_residue_debug!(du, u, p)
    l1, l2 = laplacian.(p)
    lap = kron(l2, I(p[1])) + kron(I(p[2]), l1)
    du, u = view(du, :), view(u, :)
    mul!(du, lap, u)
    return nothing
end

function system_residue_debug!(du, u, p, t)
    Ny, Nx, Nyfed = p
    num_nodes = Ny * Nx
    num_edges = Nyfed * Nx
    fill!(du, zero(eltype(u)))
    Py = view(u, 1:num_edges)
    dPy = view(du, 1:num_edges)
    ϕ = view(u, 1+num_edges:num_edges+num_nodes)
    dϕ = view(du, 1+num_edges:num_edges+num_nodes)

    Py = reshape(Py, Nyfed, Nx)
    dPy = reshape(dPy, Nyfed, Nx)
    ϕ = reshape(ϕ, Ny, Nx)
    dϕ = reshape(dϕ, Ny, Nx)

    landau_khalatnikov_debug!(dPy, Py, (Nyfed, Nx))
    poisson_residue_debug!(dϕ, ϕ, (Ny, Nx))
    return nothing
end

Nx, Ny, Nyfe = 32, 32, 13
params = Ny, Nx, Nyfe + 1

ϕ = zeros(Ny, Nx)
Py = rand(Nyfe + 1, Nx)
u0 = rand(length(Py) + length(ϕ))
du = similar(u0)

# Jacobian sparsity
system_residue_debug!(du, u0, params, 0.0)
jac_sparsity_system = ADTypes.jacobian_sparsity((du, u) -> system_residue_debug!(du, u, params, 0.0), similar(u0), copy(u0), TracerSparsityDetector())
jac_prototype = float.(jac_sparsity_system)

mass_matrix = Diagonal(vcat(ones(length(Py)), zeros(length(ϕ))))
ode_func = ODEFunction(system_residue_debug!, jac_prototype=jac_prototype, mass_matrix=mass_matrix)
ode_prob = ODEProblem(ode_func, u0, (0.0, 1e-19), params)
sol = solve(ode_prob, Rodas5(), reltol=1e-8, abstol=1e-8)

ode_func = ODEFunction(system_residue_debug!, mass_matrix=mass_matrix)
ode_prob = ODEProblem(ode_func, u0, (0.0, 1e-19), params)
sol = solve(ode_prob, Rodas5(), reltol=1e-8, abstol=1e-8)

I don't know how an ODE with mass matrix is solved, but from the stacktrace it seems that somewhere solver tries to invert Jacobian or tries to find pseudo inverse of Jacobian. And it required svd of Jacobian and julia doesn't have svd function for sparse matrix.

Providing sparse Jacobian/Jacobian sparsity to DAE should be a routine thing, so I feel like I'm calling the solve function improperly.

[ChrisRackauckas]

This was solved by SciML/NonlinearSolve.jl#577