AdaptiveRadau produces error with KLUFactorization #2892
Description
Description of the bug 🐞
OrdinaryDiffEqFIRK docs say that:
"linsolve: Any LinearSolve.jl compatible linear solver. For example, to use KLU.jl, specify RadauIIA9(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver."
Which means that KLUFactorization should work without any problem, but produces an error.
Expected behavior
AdaptiveRadau(;linsolve=KLUFactorization()) works without any issue and produces correct result.
Minimal Reproducible Example 👇
KLUFactorization fails only for AdaptiveRadau and works for QNDF
using OrdinaryDiffEq
using OrdinaryDiffEqFIRK
using ADTypes
using SparseConnectivityTracer
using LinearSolve
using Sparspak
function test_sparse(ode_solver)
function f(du, u, p, t)
du .= [u[1], u[2]]
end
u0 = [1.0, 2.0]
p = ()
du0 = similar(u0)
jac_prototype = float.(ADTypes.jacobian_sparsity(
(du, u) -> f(du, u, p, 0.0),
du0,
u0, TracerSparsityDetector()))
ode_fun = ODEFunction(f, jac_prototype=jac_prototype)
prob = ODEProblem(ode_fun, u0, (0, 10))
sol = solve(prob, ode_solver)
end
test_sparse(AdaptiveRadau(;linsolve=LUFactorization())) # Success
test_sparse(AdaptiveRadau(;linsolve=SparspakFactorization())) # Success
test_sparse(QNDF(;linsolve=KLUFactorization())) # Success
test_sparse(AdaptiveRadau(;linsolve=KLUFactorization())) # ErrorError & Stacktrace
ERROR: type Nothing has no field colptr
Stacktrace:
[1] getproperty
@ .\Base.jl:49 [inlined]
[2] pattern_changed(fact::Nothing, A::SparseArrays.SparseMatrixCSC{ComplexF64, Int64})
@ LinearSolveSparseArraysExt C:\Users\user\.julia\packages\LinearSolve\4z4nI\ext\LinearSolveSparseArraysExt.jl:405
[3] solve!(cache::LinearSolve.LinearCache{…}, alg::KLUFactorization; kwargs::@Kwargs{…})
@ LinearSolveSparseArraysExt C:\Users\user\.julia\packages\LinearSolve\4z4nI\ext\LinearSolveSparseArraysExt.jl:288
[4] solve!(::LinearSolve.LinearCache{…}; kwargs::@Kwargs{…})
@ LinearSolve C:\Users\user\.julia\packages\LinearSolve\4z4nI\src\common.jl:424
[5] dolinsolve(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, linsolve::LinearSolve.LinearCache{…}; A::SparseArrays.SparseMatrixCSC{…}, linu::Vector{…}, b::Vector{…}, du::Nothing, u::Nothing, p::Nothing, t::Nothing, weight::Nothing, solverdata::Nothing, reltol::Float64)
@ OrdinaryDiffEqDifferentiation C:\Users\user\.julia\packages\OrdinaryDiffEqDifferentiation\TKWRC\src\linsolve_utils.jl:30
[6] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqFIRK.AdaptiveRadauCache{…}, repeat_step::Bool)
@ OrdinaryDiffEqFIRK C:\Users\user\.julia\packages\OrdinaryDiffEqFIRK\Dr9ga\src\firk_perform_step.jl:1877
[7] perform_step!
@ C:\Users\user\.julia\packages\OrdinaryDiffEqFIRK\Dr9ga\src\firk_perform_step.jl:1731 [inlined]
[8] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
@ OrdinaryDiffEqCore C:\Users\user\.julia\packages\OrdinaryDiffEqCore\GMkz9\src\solve.jl:611
[9] __solve(::ODEProblem{…}, ::AdaptiveRadau{…}; kwargs::@Kwargs{})
@ OrdinaryDiffEqCore C:\Users\user\.julia\packages\OrdinaryDiffEqCore\GMkz9\src\solve.jl:7
[10] __solve
@ C:\Users\user\.julia\packages\OrdinaryDiffEqCore\GMkz9\src\solve.jl:1 [inlined]
[11] #solve_call#33
@ C:\Users\user\.julia\packages\DiffEqBase\p82Yh\src\solve.jl:127 [inlined]
[12] solve_call(_prob::ODEProblem{…}, args::AdaptiveRadau{…})
@ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\p82Yh\src\solve.jl:84
[13] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::SciMLBase.NullParameters, args::AdaptiveRadau{…}; originator::SciMLBase.ChainRulesOriginator, kwargs::@Kwargs{})
@ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\p82Yh\src\solve.jl:563
[14] solve_up
@ C:\Users\user\.julia\packages\DiffEqBase\p82Yh\src\solve.jl:540 [inlined]
[15] solve(prob::ODEProblem{…}, args::AdaptiveRadau{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
@ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\p82Yh\src\solve.jl:530
[16] solve(prob::ODEProblem{…}, args::AdaptiveRadau{…})
@ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\p82Yh\src\solve.jl:520
[17] test_sparse(ode_solver::AdaptiveRadau{…})
@ Main c:\Code\AdaptiveRadauBug\src.jl:23
[18] top-level scope
@ c:\Code\AdaptiveRadauBug\src.jl:30
Some type information was truncated. Use `show(err)` to see complete types.Environment (please complete the following information):
- Output of
using Pkg; Pkg.status()
Status `C:\Code\AdaptiveRadauBug\Project.toml`
[47edcb42] ADTypes v1.18.0
[7ed4a6bd] LinearSolve v3.44.0
[1dea7af3] OrdinaryDiffEq v6.102.1
[5960d6e9] OrdinaryDiffEqFIRK v1.16.0
[9f842d2f] SparseConnectivityTracer v1.1.1
[e56a9233] Sparspak v0.3.14- Output of
versioninfo()
Julia Version 1.11.6
Commit 9615af0f26 (2025-07-09 12:58 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: 14 × Intel(R) Core(TM) Ultra 7 165U
WORD_SIZE: 64
LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
Threads: 10 default, 0 interactive, 5 GC (on 14 virtual cores)
Environment:
JULIA_EDITOR = code
JULIA_VSCODE_REPL = 1
JULIA_NUM_THREADS = 10