# Error with jac\_prototype in Mass Matrix Form of ODEFunction (svd! on Sparse Matrix)
**URL:** https://discourse.julialang.org/t/error-with-jac-prototype-in-mass-matrix-form-of-odefunction-svd-on-sparse-matrix/127513
**Category:** Numerics
**Tags:** differentialequation
**Created:** [March 30, 2025, 1:20pm UTC](https://discourse.julialang.org/t/error-with-jac-prototype-in-mass-matrix-form-of-odefunction-svd-on-sparse-matrix/127513 "2025-03-30T13:20:13Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@sanjeev](https://discourse.julialang.org/u/sanjeev)
#### Post date: [March 30, 2025, 1:20pm UTC](https://discourse.julialang.org/t/error-with-jac-prototype-in-mass-matrix-form-of-odefunction-svd-on-sparse-matrix/127513/1 "2025-03-30T13:20:13Z")
Hi all,
I’m solving a DAE system using DifferentialEquations.jl in 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:
```julia-repl
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
```julia
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)
```
The documentation suggests that this should be possible with Rodas5() algorithm.
Any help is appreciated.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 14, 2025, 1:33am UTC](https://discourse.julialang.org/t/error-with-jac-prototype-in-mass-matrix-form-of-odefunction-svd-on-sparse-matrix/127513/2 "2025-04-14T01:33:02Z")
Can you open an issue for this one? It ended up being weirder than expected. The solution for now is to force a different nonlinear solver during initialization, but that shouldn’t be required. I believe @Oscar_Smith is looking into this but because it’s not in the bug tracker it’s not being tracked and will get lost.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [April 14, 2025, 5:07pm UTC](https://discourse.julialang.org/t/error-with-jac-prototype-in-mass-matrix-form-of-odefunction-svd-on-sparse-matrix/127513/3 "2025-04-14T17:07:17Z")
So the specific problem is an issue with solving a nonlinear problem with `Broyden(; init_jacobian = Val(:true_jacobian), autodiff=AutoFiniteDiff())`.
Specifically, the MWE is
```julia
using LinearAlgebra, SparseArrays, NonlinearSolveQuasiNewton, ForwardDiff
function f(du, u,p)
@. du = u^2 - 2
du[1] += u[4]
end
prob = NonlinearProblem(NonlinearFunction(f; jac_prototype=sparse(ones(4, 4))), rand(4))
solve(prob, Broyden(; init_jacobian = Val(:true_jacobian)))
```
---
### Author:  [@sanjeev](https://discourse.julialang.org/u/sanjeev)
#### Post date: [April 15, 2025, 9:45pm UTC](https://discourse.julialang.org/t/error-with-jac-prototype-in-mass-matrix-form-of-odefunction-svd-on-sparse-matrix/127513/4 "2025-04-15T21:45:24Z")
I have opened the issue [here](https://github.com/SciML/DifferentialEquations.jl/issues/1083).
Thanks, I will try to set the nonlinear solvers as NewtonRaphson and see if I can make it work.
I also read up on this, and it seems that solving mass matrix ODE is same as normal ODE, it is just that effectively during the stage computation of an implicit method, coefficent of y\_next is modified from (I - γhJ) to (M - γhJ). So just for educational purposes I am going to implement it.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 17, 2025, 8:20pm UTC](https://discourse.julialang.org/t/error-with-jac-prototype-in-mass-matrix-form-of-odefunction-svd-on-sparse-matrix/127513/5 "2025-04-17T20:20:14Z")
This should be solved with latest releases now.
---
### Author:  [@sanjeev](https://discourse.julialang.org/u/sanjeev)
#### Post date: [April 18, 2025, 3:31am UTC](https://discourse.julialang.org/t/error-with-jac-prototype-in-mass-matrix-form-of-odefunction-svd-on-sparse-matrix/127513/6 "2025-04-18T03:31:04Z")
Thank you!