Strange error with stiff DiffEq solvers

[arlk]

I'm not sure if this issue belongs here or somewhere downstream on the DiffEq side. But this is the MWE:

f!(D, z, p, t) = nothing
prob = ODEProblem(f!, ComponentArray(x=zeros(4)), (0.0, 1.0), 0.0) 
solve(prob, Rodas4())

and I get this error:

ERROR: LoadError: LinearAlgebra.SingularException(1)                                                                                                                                                                                                                              
Stacktrace:                                                                                                                                                                                                                                                                       
 [1] checknonsingular at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/factorization.jl:19 [inlined]                                                                                                                                    
 [2] generic_lufact!(::ComponentMatrix{Float64}, ::Val{true}; check::Bool) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:178                                                                                                   
 [3] #lu!#2 at /home/arun/.julia/packages/RecursiveFactorization/wAEz3/src/lu.jl:12 [inlined]                                                                                                                                                                                     
 [4] lu! at /home/arun/.julia/packages/RecursiveFactorization/wAEz3/src/lu.jl:9 [inlined] (repeats 2 times)                                                                                                                                                                       
 [5] (::DefaultLinSolve)(::ComponentVector{Float64}, ::ComponentMatrix{Float64}, ::ComponentVector{Float64}, ::Bool; tol::Nothing, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/arun/.julia/packages/ComponentArrays/zGwDt/src/if_requir
ed/diffeqbase.jl:44                                                                                                                                                                                                                                                               
 [6] DefaultLinSolve at /home/arun/.julia/packages/ComponentArrays/zGwDt/src/if_required/diffeqbase.jl:35 [inlined]                                                                                                                                                               
 [7] perform_step!(::OrdinaryDiffEq.ODEIntegrator{Rodas4{0,true,DefaultLinSolve,DataType},true,ComponentVector{Float64},Nothing,Float64,Float64,Float64,Float64,Float64,Array{ComponentVector{Float64},1},ODESolution{Float64,2,Array{ComponentVector{Float64},1},Nothing,Nothing,
Array{Float64,1},Array{Array{ComponentVector{Float64},1},1},ODEProblem{ComponentVector{Float64},Tuple{Float64,Float64},true,Float64,ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing
,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Rodas4{0,true,DefaultLinSolve,DataType},OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Noth
ing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{ComponentVector{Float64},1},Array{Float64,1},Array{Array{ComponentVector{Float64},1},1},OrdinaryDiffEq.Rodas4Cache{ComponentVector{Float64},ComponentVector{Float64},ComponentVector{Float64},C
omponentMatrix{Float64},ComponentMatrix{Float64},OrdinaryDiffEq.RodasTableau{Float64,Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothin
g,Nothing},ComponentVector{Float64},Float64},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64},DefaultLinSolve,SparseDi
ffTools.ForwardColorJacCache{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,
Float64},Float64},Float64,4}},ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64
,Float64},Float64},Float64,4}},ComponentVector{Float64},Array{Array{NTuple{4,Bool},1},1},UnitRange{Int64},Nothing},ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,
Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ComponentVector{Float64},Float64},Float64},Float64,1}}}},DiffEqBase.DEStats},ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothi
ng,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Rodas4Cache{ComponentVector{Float64},ComponentVector{Float64},ComponentVector{Float64},ComponentMatrix{Float64},ComponentMatrix{Float64},OrdinaryDiffEq.RodasTableau{Float64,Float64},DiffEqBase.TimeGradientWrapper{OD
EFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ComponentVector{Float64},Float64},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScali
ng{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64},DefaultLinSolve,SparseDiffTools.ForwardColorJacCache{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,type
of(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64},Float64},Float64,4}},ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typ
eof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64},Float64},Float64,4}},ComponentVector{Float64},Array{Array{NTuple{4,Bool},1},1},UnitRange{Int64},Nothing},ComponentVec
tor{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ComponentVector{Float64},Float64},Float64},Flo
at64,1}}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DE
FAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,Base.Order.ForwardOrdering},DataStructures.BinaryHeap{Float64,Base.Order.ForwardOrdering},Nothing,Nothing,Int64,Tuple{},Tuple{},Tuple{}},ComponentVector{Float64},Float64,Nothing,OrdinaryDiffEq.DefaultInit}, ::OrdinaryD
iffEq.Rodas4Cache{ComponentVector{Float64},ComponentVector{Float64},ComponentVector{Float64},ComponentMatrix{Float64},ComponentMatrix{Float64},OrdinaryDiffEq.RodasTableau{Float64,Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScalin
g{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ComponentVector{Float64},Float64},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Noth
ing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64},DefaultLinSolve,SparseDiffTools.ForwardColorJacCache{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Not
hing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64},Float64},Float64,4}},ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,No
thing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64},Float64},Float64,4}},ComponentVector{Float64},Array{Array{NTuple{4,Bool},1},1},UnitRange{Int64},Nothing},ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.TimeG
radientWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ComponentVector{Float64},Float64},Float64},Float64,1}}}, ::Bool) at /home/arun/.julia/packages/Ordi
naryDiffEq/pH3tP/src/perform_step/rosenbrock_perform_step.jl:779                                                                                                                                                                                                                  
 [8] perform_step! at /home/arun/.julia/packages/OrdinaryDiffEq/pH3tP/src/perform_step/rosenbrock_perform_step.jl:745 [inlined] 
[9] solve!(::OrdinaryDiffEq.ODEIntegrator{Rodas4{0,true,DefaultLinSolve,DataType},true,ComponentVector{Float64},Nothing,Float64,Float64,Float64,Float64,Float64,Array{ComponentVector{Float64},1},ODESolution{Float64,2,Array{ComponentVector{Float64},1},Nothing,Nothing,Array{F
loat64,1},Array{Array{ComponentVector{Float64},1},1},ODEProblem{ComponentVector{Float64},Tuple{Float64,Float64},true,Float64,ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothin
g,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Rodas4{0,true,DefaultLinSolve,DataType},OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Not
hing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{ComponentVector{Float64},1},Array{Float64,1},Array{Array{ComponentVector{Float64},1},1},OrdinaryDiffEq.Rodas4Cache{ComponentVector{Float64},ComponentVector{Float64},ComponentVector{Float64},Componen
tMatrix{Float64},ComponentMatrix{Float64},OrdinaryDiffEq.RodasTableau{Float64,Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothi
ng},ComponentVector{Float64},Float64},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64},DefaultLinSolve,SparseDiffTools
.ForwardColorJacCache{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64
},Float64},Float64,4}},ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float6
4},Float64},Float64,4}},ComponentVector{Float64},Array{Array{NTuple{4,Bool},1},1},UnitRange{Int64},Nothing},ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing
,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ComponentVector{Float64},Float64},Float64},Float64,1}}}},DiffEqBase.DEStats},ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Noth
ing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Rodas4Cache{ComponentVector{Float64},ComponentVector{Float64},ComponentVector{Float64},ComponentMatrix{Float64},ComponentMatrix{Float64},OrdinaryDiffEq.RodasTableau{Float64,Float64},DiffEqBase.TimeGradientWrapper{ODEFuncti
on{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ComponentVector{Float64},Float64},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64},DefaultLinSolve,SparseDiffTools.ForwardColorJacCache{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64},Float64},Float64,4}},ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64},Float64},Float64,4}},ComponentVector{Float64},Array{Array{NTuple{4,Bool},1},1},UnitRange{Int64},Nothing},ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ComponentVector{Float64},Float64},Float64},Float64,1}}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,Base.Order.ForwardOrdering},DataStructures.BinaryHeap{Float64,Base.Order.ForwardOrdering},Nothing,Nothing,Int64,Tuple{},Tuple{},Tuple{}},ComponentVector{Float64},Float64,Nothing,OrdinaryDiffEq.DefaultInit}) at /home/arun/.julia/packages/OrdinaryDiffEq/pH3tP/src/solve.jl:450                   
 [10] __solve(::ODEProblem{ComponentVector{Float64},Tuple{Float64,Float64},true,Float64,ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Rodas4{0,true,DefaultLinSolve,DataType}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/arun/.julia/packages/OrdinaryDiffEq/pH3tP/src/solve.jl:5
 [11] __solve at /home/arun/.julia/packages/OrdinaryDiffEq/pH3tP/src/solve.jl:4 [inlined]                                                
 [12] #solve_call#456 at /home/arun/.julia/packages/DiffEqBase/3iigH/src/solve.jl:65 [inlined]                                           
 [13] solve_call at /home/arun/.julia/packages/DiffEqBase/3iigH/src/solve.jl:52 [inlined]                                                
 [14] #solve_up#458 at /home/arun/.julia/packages/DiffEqBase/3iigH/src/solve.jl:86 [inlined]                                             
 [15] solve_up at /home/arun/.julia/packages/DiffEqBase/3iigH/src/solve.jl:79 [inlined]                                                  
 [16] #solve#457 at /home/arun/.julia/packages/DiffEqBase/3iigH/src/solve.jl:74 [inlined]                                                
 [17] solve(::ODEProblem{ComponentVector{Float64},Tuple{Float64,Float64},true,Float64,ODEFunction{true,typeof(f!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{}
,Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Rodas4{0,true,DefaultLinSolve,DataType}) at /home/arun/.julia/packages/DiffEqBase/3iigH/src/solve.jl:72
 [18] top-level scope at /tmp/test.jl:9
 [19] include(::String) at ./client.jl:457
 [20] top-level scope at REPL[2]:1                                  
in expression starting at /tmp/test.jl:9

Instead, if I downgrade to v.0.7.4 (or not use ComponentArrays altogether) this error disappears and everything starts working as intended.

[jonniedie]

Huh. This is a weird one. lu! is pretty much doing what I would expect it to, so that's not the issue.

julia> using ComponentArrays, RecursiveFactorization

julia> A = ComponentArray(rand(4,4), Axis(x=1:4), Axis(x=1:4))
4×4 ComponentMatrix{Float64} with axes Axis(x = 1:4,) × Axis(x = 1:4,)
 0.783502  0.842274  0.430332  0.194907
 0.211523  0.137423  0.975525  0.342043
 0.702106  0.233538  0.791773  0.254408
 0.401241  0.773116  0.414625  0.549047

julia> B = ComponentArray(zeros(4,4), Axis(x=1:4), Axis(x=1:4))
4×4 ComponentMatrix{Float64} with axes Axis(x = 1:4,) × Axis(x = 1:4,)
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> RecursiveFactorization.lu!(A)
LU{Float64,ComponentMatrix{Float64}}
L factor:
4×4 Array{Float64,2}:
 1.0        0.0       0.0       0.0
 0.896112   1.0       0.0       0.0
 0.269971   0.172604  1.0       0.0
 0.512112  -0.655709  0.583546  1.0
U factor:
4×4 Array{Float64,2}:
 0.783502   0.842274  0.430332  0.194907
 0.0       -0.521234  0.406147  0.0797491
 0.0        0.0       0.789246  0.275659
 0.0        0.0       0.0       0.340665

julia> RecursiveFactorization.lu!(B)
ERROR: SingularException(1)
Stacktrace:
 [1] checknonsingular at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\factorization.jl:19 [inlined]
 [2] generic_lufact!(::ComponentMatrix{Float64}, ::Val{true}; check::Bool) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\lu.jl:178
 [3] #lu!#2 at C:\Users\jdiegelm\.julia\packages\RecursiveFactorization\wAEz3\src\lu.jl:12 [inlined]
 [4] lu! at C:\Users\jdiegelm\.julia\packages\RecursiveFactorization\wAEz3\src\lu.jl:9 [inlined] (repeats 2 times)
 [5] top-level scope at REPL[45]:1

[jonniedie]

I'll have to dig a little deeper tonight to see what's going on. Because it's throwing a SingularException even for cases where nothing should be singular

f!(D, z, p, t) = (D .= z)
prob = ODEProblem(f!, ComponentArray(x=ones(4)), (0.0, 1.0), 0.0) 
solve(prob, Rodas4())

[jonniedie]

Oooh, that was a good catch. Thanks!

ComponentArrays used to return views for all indexing, even if it was plain indexing like x[3:4]. I recently decided that this was a bad thing as most people would expect x[3:4] to return a new array, so I changed it to do so. However, the old method of view would just forward to the getindex since getindex used to always just work like view. When I made the change to getindex, the forwarded view no longer viewed, but returned a new array. This made broadcasted setting like x[3:4] .= 0 not work because dotview eventually calls view here. When DifferentialEquations was making the Jacobian, it was setting by row like J[idx, :] .= ..., which was using the broken view method.