MKLLUFactorization broken with v3.25.0?

[hersle]

When updating only LinearSolve from v3.24.0 to v3.25.0, this example breaks:

using MKL, ModelingToolkit, OrdinaryDiffEq, LinearSolve
using ModelingToolkit: t_nounits as t, D_nounits as D

@variables x(t) y(t)
eqs = [D(x) ~  x + y, D(y) ~ x - y]
@mtkbuild sys = System(eqs, t)
prob = ODEProblem(sys, [x => 1.0, y => 1.0], (0.0, 1.0))
sol = solve(prob, KenCarp4(linsolve = MKLLUFactorization()))

┌ Warning: At t=0.0, dt was forced below floating point epsilon 5.0e-324, and step error estimate = 1.0. 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 ~/.julia/packages/SciMLBase/dld4W/src/integrator_interface.jl:623
retcode: Unstable
Interpolation: 3rd order Hermite
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [1.0, 1.0]

The ODE solves successfully when downgrading to v3.24.0 without any errors and warnings.
It also solves successfully with linsolve = LUFactorization() on v3.25.0.

(@v1.11) pkg> status
Status `~/.julia/environments/v1.11/Project.toml`
  [6e4b80f9] BenchmarkTools v1.6.0
  [0e736298] Bessels v0.2.8
  [d360d2e6] ChainRulesCore v1.25.2
⌅ [f68482b8] Cthulhu v2.16.5
  [82cc6244] DataInterpolations v8.3.1
  [163ba53b] DiffResults v1.1.0
  [e30172f5] Documenter v1.14.1
  [6a86dc24] FiniteDiff v2.27.0
⌅ [f6369f11] ForwardDiff v0.10.38
  [c9556dd2] ForwardDiffChainRules v0.3.0
  [7073ff75] IJulia v1.29.2
⌅ [a98d9a8b] Interpolations v0.15.1
  [98e50ef6] JuliaFormatter v2.1.6
  [23fbe1c1] Latexify v0.16.8
  [7ed4a6bd] LinearSolve v3.25.0
  [2fda8390] LsqFit v0.15.1
  [33e6dc65] MKL v0.9.0
  [961ee093] ModelingToolkit v10.14.0
  [77ba4419] NaNMath v1.1.3
  [8913a72c] NonlinearSolve v4.10.0
  [e7bfaba1] NumericalIntegration v0.3.4
  [429524aa] Optim v1.13.2
  [1dea7af3] OrdinaryDiffEq v6.101.0
  [5ad8b20f] PhysicalConstants v0.2.4
  [14b8a8f1] PkgTemplates v0.7.56
  [91a5bcdd] Plots v1.40.17
  [c3e4b0f8] Pluto v0.20.13
  [d236fae5] PreallocationTools v0.4.29
  [efd6af41] ProfileCanvas v0.1.6
  [c46f51b8] ProfileView v1.10.1
  [1fd47b50] QuadGK v2.11.2
⌃ [295af30f] Revise v3.8.0
  [0bca4576] SciMLBase v2.107.0
  [53ae85a6] SciMLStructures v1.7.0
  [efcf1570] Setfield v1.1.2
  [276daf66] SpecialFunctions v2.5.1
  [f3b207a7] StatsPlots v0.15.7
  [64c5a815] SymBoltz v0.6.0 `../../dev/SymBoltz`
  [2efcf032] SymbolicIndexingInterface v0.3.42
  [d1185830] SymbolicUtils v3.29.0
  [0c5d862f] Symbolics v6.48.1
  [0c614874] TerminalPager v0.6.4
  [a759f4b9] TimerOutputs v0.5.29
  [fce5fe82] Turing v0.39.10
  [1986cc42] Unitful v1.24.0
  [6112ee07] UnitfulAstro v1.2.2
  [51b3953f] VimBindings v0.3.15
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

[ChrisRackauckas]

The only change I can think of is https://github.com/SciML/LinearSolve.jl/pull/667/files#diff-c304873533e7a8d67eb08f6fec74dded4f147924d0f0b7d3962d365600341359L230-L247. Try reverting that and see? It is tested here though so I'm surprised CI didn't catch it.

Maybe try this branch? #680

[hersle]

Thanks! Can confirm that #680 fixes it:

(@v1.11) pkg> add https://github.com/SciML/LinearSolve.jl#ChrisRackauckas-patch-2

...

(@v1.11) pkg> status
Status `~/.julia/environments/v1.11/Project.toml`
  ...
  [7ed4a6bd] LinearSolve v3.25.0 `https://github.com/SciML/LinearSolve.jl#ChrisRackauckas-patch-2`
  ...

julia> sol = solve(prob, KenCarp4(linsolve = MKLLUFactorization()))
retcode: Success
Interpolation: 3rd order Hermite
t: 5-element Vector{Float64}:
 0.0
 0.047298897789691714
 0.3659213952903471
 0.7746067188716641
 1.0
u: 5-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [1.0022380199497678, 1.0969063739859737]
 [1.1369051297182677, 1.9018270568034834]
 [1.6623791798651655, 3.5404109662598438]
 [2.17804408986812, 4.914418755592536]