AdaptiveRadau produces error with KLUFactorization

[mrarat]

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()))      # Error

Error & 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

[ChrisRackauckas]

The bigger thing here is just the Radau needs to support sparse matrices, given it's a FIRK it needs to do a few things to the sparsity pattern but it should just do what is required. It may be simple enough to actually Claude it though.

[ChrisRackauckas-Claude]

I've identified the root cause of this issue and have a partial fix in progress.

Root Cause: When (and other Radau methods) initialize complex sparse matrices for the linear solver, they use similar(J, Complex{eltype(W1)}) followed by zero-filling. This creates a sparse matrix with the correct sparsity pattern but all-zero values. When this is passed to LinearSolve.init with KLUFactorization, KLU cannot properly initialize its factorization cache, resulting in the type Nothing has no field colptr error.

Why QNDF works: QNDF doesn't use complex sparse matrices in the same way - it works through the nonlinear solver infrastructure differently.

Fix in Progress: I'm working on a fix to properly preserve the jacobian prototype's sparsity pattern when creating complex W matrices. The challenge is ensuring KLU can initialize with these matrices.

I'll open a PR with the current progress and continue iterating on the solution.

[ChrisRackauckas-Claude]

I've identified the root cause of this issue and have a partial fix in progress.

Root Cause: When AdaptiveRadau (and other Radau methods) initialize complex sparse matrices for the linear solver, they use similar(J, Complex) followed by zero-filling. This creates a sparse matrix with the correct sparsity pattern but all-zero values. When passed to LinearSolve.init with KLUFactorization, KLU cannot properly initialize its factorization cache, resulting in the 'type Nothing has no field colptr' error.

Why QNDF works: QNDF doesn't use complex sparse matrices in the same way.

Fix in Progress: Working on properly preserving the jacobian prototype's sparsity pattern when creating complex W matrices. Opening a draft PR to get feedback.