Fix and speed up sparsity pattern check

[hersle]

See SciML/OrdinaryDiffEq.jl#2924. Now this works:

using ModelingToolkit, OrdinaryDiffEqRosenbrock, LinearSolve, BenchmarkTools
D, t = ModelingToolkit.D_nounits, ModelingToolkit.t_nounits
@variables x(t) y(t)
@named M = System([D(x) ~ x + y, D(y) ~ x], t)
M = mtkcompile(M)
prob = ODEProblem(M, [x => 1.0, y => 1.0], (0.0, 1.0); jac = true, sparse = true)

@btime sol = solve(prob, Rodas5P(linsolve = KLUFactorization(check_pattern = true)); adaptive = false, dt = 1e-6, save_everystep = false)
@btime sol = solve(prob, Rodas5P(linsolve = KLUFactorization(check_pattern = false)); adaptive = false, dt = 1e-6, save_everystep = false)

  761.541 ms (13000303 allocations: 305.19 MiB)
  675.513 ms (9000307 allocations: 137.34 MiB)

[hersle]

This early length test makes sure a new factorization is produced when the ODE solver does the first factorization, even when check_pattern = false.

There could very well be cleaner ways to accomplish the same, just not obvious to me. Give me a hint if you see a better solution :)

[hersle]

Now also for UMFPACK.

And with a non-allocating pattern check, so the default check_pattern = true also performs better (same benchmark as in OP):

  685.764 ms (9000307 allocations: 137.34 MiB)
  668.947 ms (9000307 allocations: 137.34 MiB)

[hersle]

Finally adding this specialization on the common field in KLUFactorization. Without it, all the Ref(klu.common) calls in KLU allocate (also seen in yellow in the flame graph in SciML/OrdinaryDiffEq.jl#2924). Then most allocations in the example are gone:

  644.805 ms (306 allocations: 14.77 KiB)
  647.733 ms (306 allocations: 14.77 KiB)

[ChrisRackauckas]

Yeah that's a nice shortcut that should make this relatively cheap.

[hersle]

Looks like one real test failure, will look into it. Looks like it should be trivial to fix.

https://github.com/SciML/LinearSolve.jl/actions/runs/20082406087/job/57612355234?pr=845#step:6:1244

[ChrisRackauckas]

yes that looks real

[hersle]

Sorry @ChrisRackauckas I was mistaken, that test fails also on master, see e.g. here. That is the same failure with fact being nothing.

I am not sure what has to be done to fix the test, but I think this PR looks unrelated.

[ChrisRackauckas]

Yeah I just saw that. I think @oscardssmith that might be related to the WOperator move? Though it's not entirely clear why.

[oscardssmith]

I don't think so. I don't believe we've merged the PR to actually use the moved WOperator.

[hersle]

Looks to me like that OrdinaryDiffEq test starts to fail with 5abf86b:

shell> git checkout 5abf86b4~
Previous HEAD position was 5abf86b4 Make verbose parameter default more robust for Julia LTS
HEAD is now at c4ce1de1 Update JET tests that now pass with improved type stability

julia> sol = solve(prob, Rosenbrock23(linsolve = KLUFactorization()))
┌ Warning: At t=0.0, dt was forced below floating point epsilon 5.0e-324, and step error estimate = NaN. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase C:\Users\herma\.julia\dev\SciMLBase\src\integrator_interface.jl:671
retcode: Unstable
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

shell> git checkout 5abf86b4                                 
Previous HEAD position was c4ce1de1 Update JET tests that now pass with improved type stability
HEAD is now at 5abf86b4 Make verbose parameter default more robust for Julia LTS

julia> sol = solve(prob, Rosenbrock23(linsolve = KLUFactorization()))
ERROR: FieldError: type Nothing has no field `colptr`; Nothing has no fields at all.
Stacktrace:
  [1] getproperty
    @ .\Base_compiler.jl:54 [inlined]
  [2] pattern_changed(fact::Nothing, A::SparseMatrixCSC{Float64, Int64})
    @ LinearSolveSparseArraysExt C:\Users\herma\.julia\dev\LinearSolve\ext\LinearSolveSparseArraysExt.jl:397
  [3] solve!(cache::LinearSolve.LinearCache{…}, alg::KLUFactorization; kwargs::@Kwargs{…})
    @ LinearSolveSparseArraysExt C:\Users\herma\.julia\dev\LinearSolve\ext\LinearSolveSparseArraysExt.jl:283
  [4] solve!(::LinearSolve.LinearCache{…}; kwargs::@Kwargs{…})
    @ LinearSolve C:\Users\herma\.julia\dev\LinearSolve\src\common.jl:441
  [5] dolinsolve(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, linsolve::LinearSolve.LinearCache{…}; A::OrdinaryDiffEqDifferentiation.WOperator{…}, linu::Nothing, b::Vector{…}, du::Vector{…}, u::Vector{…}, p::SciMLBase.NullParameters, t::Float64, weight::Vector{…}, solverdata::@NamedTuple{…}, reltol::Float64)
    @ OrdinaryDiffEqDifferentiation c:\Users\herma\.julia\dev\OrdinaryDiffEq\lib\OrdinaryDiffEqDifferentiation\src\linsolve_utils.jl:30
  [6] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqRosenbrock.Rosenbrock23Cache{…},
 repeat_step::Bool)
    @ OrdinaryDiffEqRosenbrock c:\Users\herma\.julia\dev\OrdinaryDiffEq\lib\OrdinaryDiffEqRosenbrock\src\rosenbrock_perform_step.jl:58
  [7] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqRosenbrock.Rosenbrock23Cache{…})
    @ OrdinaryDiffEqRosenbrock c:\Users\herma\.julia\dev\OrdinaryDiffEq\lib\OrdinaryDiffEqRosenbrock\src\rosenbrock_perform_step.jl:27
  [8] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ OrdinaryDiffEqCore c:\Users\herma\.julia\dev\OrdinaryDiffEq\lib\OrdinaryDiffEqCore\src\solve.jl:611
  [9] __solve(::ODEProblem{…}, ::Rosenbrock23{…}; kwargs::@Kwargs{})
    @ OrdinaryDiffEqCore c:\Users\herma\.julia\dev\OrdinaryDiffEq\lib\OrdinaryDiffEqCore\src\solve.jl:7
 [10] __solve
    @ c:\Users\herma\.julia\dev\OrdinaryDiffEq\lib\OrdinaryDiffEqCore\src\solve.jl:1 [inlined]
 [11] #solve_call#23
    @ C:\Users\herma\.julia\dev\DiffEqBase\src\solve.jl:127 [inlined]
 [12] solve_call
    @ C:\Users\herma\.julia\dev\DiffEqBase\src\solve.jl:84 [inlined]
 [13] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::SciMLBase.NullParameters, args::Rosenbrock23{…}; originator::SciMLBase.ChainRulesOriginator, kwargs::@Kwargs{})
    @ DiffEqBase C:\Users\herma\.julia\dev\DiffEqBase\src\solve.jl:563
 [14] solve_up
    @ C:\Users\herma\.julia\dev\DiffEqBase\src\solve.jl:540 [inlined]
 [15] #solve#29
    @ C:\Users\herma\.julia\dev\DiffEqBase\src\solve.jl:530 [inlined]
 [16] solve(prob::ODEProblem{…}, args::Rosenbrock23{…})      
    @ DiffEqBase C:\Users\herma\.julia\dev\DiffEqBase\src\solve.jl:520
 [17] top-level scope
    @ REPL[34]:1
Some type information was truncated. Use `show(err)` to see complete types.

But like you see, it looks like it was broken/unstable before it errored/failed, too.