error with SymTridiagonal Jacobian

[isaacsas]

On OrdinaryDiffEq v6.58.0

using OrdinaryDiffEq, LinearAlgebra

M = 10
d = -2*ones(M)
u = ones(M-1)
L = SymTridiagonal(d,u)

function f(du,u,p,t)
       mul!(du, p.L, u)
       nothing
end

function jac(J,u,p,t)
       J .= p.L
       nothing
end

p = (; L)
x = (1:M) * (1 / (M+1))
u0 = sin.(pi*x)

odefun = ODEFunction(f; jac, jac_prototype = L)
oprob = ODEProblem(odefun, u0, (0.0, 1.0), p)
sol = solve(oprob, Trapezoid())

gives the following error

ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.
Stacktrace:
  [1] ldlt!(S::SymTridiagonal{Float64, Vector{Float64}})
    @ LinearAlgebra /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/LinearAlgebra/src/ldlt.jl:122
  [2] ldlt(M::SymTridiagonal{Float64, Vector{Float64}}; shift::Bool)
    @ LinearAlgebra /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/LinearAlgebra/src/ldlt.jl:169
  [3] ldlt
    @ LinearSolve /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/LinearAlgebra/src/ldlt.jl:163 [inlined]
  [4] ldlt_instance
    @ LinearSolve ~/.julia/packages/ArrayInterface/a4Wwt/src/ArrayInterface.jl:569 [inlined]
  [5] 
    @ LinearSolve ~/.julia/packages/LinearSolve/tPdfJ/src/factorization.jl:340 [inlined]
  [6] macro expansion
    @ LinearSolve ~/.julia/packages/LinearSolve/tPdfJ/src/default.jl:299 [inlined]
  [7] init_cacheval(alg::LinearSolve.DefaultLinearSolver, A::SymTridiagonal{…}, b::Vector{…}, u::Vector{…}, Pl::LinearSolve.InvPreconditioner{…}, Pr::Diagonal{…}, maxiters::Int64, abstol::Float64, reltol::Float64, verbose::Bool, assump::OperatorAssumptions{…})
    @ LinearSolve ~/.julia/packages/LinearSolve/tPdfJ/src/default.jl:282
  [8] #init#3
    @ OrdinaryDiffEq ~/.julia/packages/LinearSolve/tPdfJ/src/common.jl:153 [inlined]
  [9] init
    @ OrdinaryDiffEq ~/.julia/packages/LinearSolve/tPdfJ/src/common.jl:114 [inlined]
 [10] build_nlsolver(alg::Trapezoid{…}, nlalg::NLNewton{…}, u::Vector{…}, uprev::Vector{…}, p::@NamedTuple{…}, t::Float64, dt::Float64, f::ODEFunction{…}, rate_prototype::Vector{…}, ::Type{…}, ::Type{…}, ::Type{…}, γ::Rational{…}, c::Int64, α::Int64, ::Val{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yppG9/src/nlsolve/utils.jl:198
 [11] build_nlsolver(alg::Trapezoid{…}, nlalg::NLNewton{…}, u::Vector{…}, uprev::Vector{…}, p::@NamedTuple{…}, t::Float64, dt::Float64, f::ODEFunction{…}, rate_prototype::Vector{…}, ::Type{…}, ::Type{…}, ::Type{…}, γ::Rational{…}, c::Int64, α::Int64, ::Val{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yppG9/src/nlsolve/utils.jl:146 [inlined]
 [12] build_nlsolver(alg::Trapezoid{…}, nlalg::NLNewton{…}, u::Vector{…}, uprev::Vector{…}, p::@NamedTuple{…}, t::Float64, dt::Float64, f::ODEFunction{…}, rate_prototype::Vector{…}, ::Type{…}, ::Type{…}, ::Type{…}, γ::Rational{…}, c::Int64, α::Int64, ::Val{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yppG9/src/nlsolve/utils.jl:136 [inlined]
 [13] alg_cache(alg::Trapezoid{…}, u::Vector{…}, rate_prototype::Vector{…}, ::Type{…}, ::Type{…}, ::Type{…}, uprev::Vector{…}, uprev2::Vector{…}, f::ODEFunction{…}, t::Float64, dt::Float64, reltol::Float64, p::@NamedTuple{…}, calck::Bool, ::Val{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yppG9/src/caches/sdirk_caches.jl:111
 [14] __init(prob::ODEProblem{…}, alg::Trapezoid{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Rational{…}, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yppG9/src/solve.jl:324
 [15] __init (repeats 5 times)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yppG9/src/solve.jl:10 [inlined]
 [16] __solve(::ODEProblem{…}, ::Trapezoid{…}; kwargs::@Kwargs{})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yppG9/src/solve.jl:5
 [17] __solve
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yppG9/src/solve.jl:1 [inlined]
 [18] #solve_call#34
    @ DiffEqBase ~/.julia/packages/DiffEqBase/xSmHR/src/solve.jl:571 [inlined]
 [19] solve_call(_prob::ODEProblem{…}, args::Trapezoid{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/xSmHR/src/solve.jl:537
 [20] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::@NamedTuple{…}, args::Trapezoid{…}; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/xSmHR/src/solve.jl:1033
 [21] solve_up
    @ DiffEqBase ~/.julia/packages/DiffEqBase/xSmHR/src/solve.jl:1006 [inlined]
 [22] solve(prob::ODEProblem{…}, args::Trapezoid{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/xSmHR/src/solve.jl:943
 [23] solve(prob::ODEProblem{…}, args::Trapezoid{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/xSmHR/src/solve.jl:933
 [24] top-level scope
    @ REPL[43]:1
Some type information was truncated. Use `show(err)` to see complete types.

Making the matrix Tridiagonal works fine however.

[ChrisRackauckas]

I think it's a Base Julia issue on v1.9 since its fine in v1.10, and just something inside of ldlt! that changed.

[isaacsas]

This errors for me on 1.10-beta3. Also, the error is non-deterministic. I just started multiple fresh REPLS (~10?) and it only errored in like 1/3 of them.

[isaacsas]

Should I open an issue on JuliaLang?

[isaacsas]

Have confirmed I see the same (non-deterministic error on 1.9 too).

[ChrisRackauckas]

Interesting. Yeah open on Julialang with just the linear solve.

[isaacsas]

I don't see this error though when directly applying ldlt! to the discrete Laplacian here. It is only when going through OrdinaryDiffEq (which I assume is applying ldlt! to some kind of M matrix).

[ChrisRackauckas]

is this with check=false?

[isaacsas]

I'm not setting such an argument to any function and/or when starting Julia from the command line (but I'm not sure what you are referring to with check=false).

[ChrisRackauckas]

JuliaArrays/ArrayInterface.jl#423 will make this always work.

[isaacsas]

Awesome, thanks!