Hi everyone,

I am trying to use the KLU method as linear solver for solving an ODE problem using Sundials.

I have tried to apply it to the CVODE_BDF solver for a simple problem, the rober equation. However, when I provide the sparsity pattern determined from the Jacobian function, that I also provide to the ODEFunction, it says:

“Sparsity Pattern in receiving SUNMatrix doesn’t match sending SparseMatrix”

The code is:

using Sundials

using IncompleteLU

using SparseArrays

using Plots

using LinearAlgebra################## Solving ODE ##################

ode repressentation

function rober!(du, u, p, t)

y₁, y₂, y₃ = u

k₁, k₂, k₃ = p

du[1] = -k₁ * y₁ + k₃ * y₂ * y₃

du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃

du[3] = k₂ * y₂^2

nothing

endu0 = [1.0; 0.0; 0.0]

tspan = (0.0, 1e4)

p = [0.04, 3e7, 1e4]

prob = ODEProblem(rober!, u0, tspan,p)

sol = solve(prob,CVODE_BDF())

p1 = plot();

plot!(p1,sol,label=“ODE”);################## Solving ODE With Jacobian ##################

#Jacobian

function rober_jac!(J, u, p, t)

y₁, y₂, y₃ = u

k₁, k₂, k₃ = p

J[1, 1] = k₁ * -1

J[2, 1] = k₁

J[3, 1] = 0

J[1, 2] = y₃ * k₃

J[2, 2] = y₂ * k₂ * -2 + y₃ * k₃ * -1

J[3, 2] = y₂ * 2 * k₂

J[1, 3] = k₃ * y₂

J[2, 3] = k₃ * y₂ * -1

J[3, 3] = 0

nothing

endu0 = [1.0, 0, 0]

tspan = (0.0, 1e4)

p = [0.04, 3e7, 1e4]

jac_proto = zeros(length(u0),length(u0))

rober_jac!(jac_proto, u0 .+ 1e-6, p, 0.0)

fun = ODEFunction(rober!; jac = rober_jac!, jac_prototype = sparse(jac_proto))

prob = ODEProblem(fun, u0, tspan,p)

sol = solve(prob,CVODE_BDF(; linear_solver = :KLU))

p1 = plot();

plot!(p1,sol,label=“ODE”);

As seen from the ODE definition and the Jacobian, the sparsity patters is indeed:

x x x

x x x

⋅ x ⋅

However, if I instead set the sparsity patters as a dense 3x3 matrix, it solves the problem without any warning.

fun = ODEFunction(rober!; jac = rober_jac!, jac_prototype = sparse(ones(3,3)))

prob = ODEProblem(fun, u0, tspan,p)

sol = solve(prob,CVODE_BDF(; linear_solver = :KLU))

p1 = plot();

plot!(p1,sol,label=“ODE”);

To me, that seems a bit odd as the elements [3,1] and [3,3] in the jacobian prototype are indeed zeros as also seen from the ODE definition (Rober!) and the Jacobian function (Rober_jac!).

Did I specify anything incorrectly?

I also looked at the examples from Sundials.jl/test/common_interface/jacobians.jl at master · SciML/Sundials.jl (github.com), however, the Jacobian for the Lotka problem is in this case a dense 2x2 matrix.

Did some of you guys try to use the KLU linear solver in a Sundials solver, either ODE or DAE, and specify a sparsity pattern that is not a dense matrix?