Exponential RK + Krylov + large, complex ODE system = segfault

[Lazersmoke]

When I try to solve my large (# of eqns = 128^2) ComplexF64-valued ODE with any krylov-enabled exponential RK method, using an in-place derivative function, I get a segfault. (and all of those are strict requirements, no segfault without them)
The segfault occurs for both hermitian (lanczos) and non-hermitian (arnoldi) linear parts in the semilinear problem.

It is somehow related to this line, which calls arnoldi! with a KrylovSubspace cache which was allocated here (A). This KrylovSubspace is slightly different from the one created by out-of-place arnoldi, see here (B).
In particular, (B) has KrylovSubspace.V as a sparse array, while (A) has it as a dense one instead.
Calling arnoldi! on its own with such an (A) constructed subspace object triggers the segfault.

The segfault would not trigger if the ConstantCache variants were used, but I don't know how to set this as an option to solve?

I don't know enough to say if this is really a DiffEq issue, or an ExponentialUtilities one, but it seems like it can be fixed by changing how the cache is set up somehow.

Gist with minimum reproducible code

I'm on the latest version:

(@v1.6) pkg> status OrdinaryDiffEq
      Status `C:\Users\Sam\.julia\environments\v1.6\Project.toml`
  [1dea7af3] OrdinaryDiffEq v5.60.0

Thanks!

[YingboMa]

That looks like JuliaLang/LinearAlgebra.jl#754 . Which Julia version are you using?

[Lazersmoke]

julia version 1.6.1

[Lazersmoke]

Can confirm

x = zeros(ComplexF64, 10001)
x' * x

reproduces

[YingboMa]

I think that should be fixed on 1.6.1. Are you using the Julia office binary from https://julialang.org/downloads/ ?

[Lazersmoke]

Yes, official binary
LinearAlgebra.BLAS.openblas_get_config() gives:
"OpenBLAS 0.3.10 USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=32"
Which is an earlier version than JuliaLang/julia#39216

[YingboMa]

I see. Do you mind to try Julia v1.7 beta out to see if the problem persists? I couldn't reproduce this segfault myself on 1.6.1.

[YingboMa]

BTW, I cannot run the MWE you posted. It errors with

julia> # krylov=true required
       # Also errors with ETDRK4
       sol = solve(prob,LawsonEuler(krylov=true);dt = 0.01)
ERROR: MethodError: no method matching eigen!(::SymTridiagonal{ComplexF64, Vector{ComplexF64}})
Closest candidates are:
  eigen!(::SymTridiagonal{var"#s832", V} where {var"#s832"<:Union{Float32, Float64}, V<:AbstractVector{var"#s832"}}) at /Users/scheme/.bin/julia-release/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:280
  eigen!(::SymTridiagonal{var"#s832", V} where {var"#s832"<:Union{Float32, Float64}, V<:AbstractVector{var"#s832"}}, ::UnitRange) at /Users/scheme/.bin/julia-release/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:283
  eigen!(::SymTridiagonal{var"#s832", V} where {var"#s832"<:Union{Float32, Float64}, V<:AbstractVector{var"#s832"}}, ::Real, ::Real) at /Users/scheme/.bin/julia-release/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:288
  ...
Stacktrace:
 [1] expv!(w::Vector{ComplexF64}, t::Float64, Ks::ExponentialUtilities; cache::ExponentialUtilities)
   @ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/vU95t/src/krylov_phiv.jl:85
 [2] perform_step!(integrator::OrdinaryDiffEq.ODEIn, cache::OrdinaryDiffEq.Lawso, repeat_step::Bool)
   @ OrdinaryDiffEq ~/src/julia/OrdinaryDiffEq/src/perform_step/exponential_rk_perform_step.jl:86
 [3] perform_step!
   @ ~/src/julia/OrdinaryDiffEq/src/perform_step/exponential_rk_perform_step.jl:71 [inlined]
 [4] solve!(integrator::OrdinaryDiffEq.ODEIn)
   @ OrdinaryDiffEq ~/src/julia/OrdinaryDiffEq/src/solve.jl:478
 [5] __solve(::ODEProblem{Vector{Co, ::LawsonEuler{DataType; kwargs::Base.Iterators.Pairs)
   @ OrdinaryDiffEq ~/src/julia/OrdinaryDiffEq/src/solve.jl:5
 [6] #solve_call#58
   @ ~/src/julia/DiffEqBase/src/solve.jl:61 [inlined]
 [7] #solve_up#60
   @ ~/src/julia/DiffEqBase/src/solve.jl:82 [inlined]
 [8] #solve#59
   @ ~/src/julia/DiffEqBase/src/solve.jl:70 [inlined]
 [9] top-level scope
   @ REPL[15]:3

[Lazersmoke]

Very interesting!
If I (in 1.6.1) add this line right after the using statements:
LinearAlgebra.BLAS.set_num_threads(1)
(this is the recommended work around, perhaps it is equivalent to just having a different BLAS version?)

I get (the same?) error:

ERROR: LoadError: MethodError: no method matching eigen!(::SymTridiagonal{ComplexF64, Vector{ComplexF64}})
Closest candidates are:
  eigen!(!Matched::SymTridiagonal{var"#s832", V} where {var"#s832"<:Union{Float32, Float64}, V<:AbstractVector{var"#s832"}}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\tridiag.jl:280
  eigen!(!Matched::SymTridiagonal{var"#s832", V} where {var"#s832"<:Union{Float32, Float64}, V<:AbstractVector{var"#s832"}}, !Matched::UnitRange) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\tridiag.jl:283
  eigen!(!Matched::SymTridiagonal{var"#s832", V} where {var"#s832"<:Union{Float32, Float64}, V<:AbstractVector{var"#s832"}}, !Matched::Real, !Matched::Real) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\tridiag.jl:288
  ...
Stacktrace:
 [1] expv!(w::Vector{ComplexF64}, t::Float64, Ks::ExponentialUtilities.KrylovSubspace{ComplexF64, ComplexF64, Float64, Matrix{ComplexF64}, Matrix{ComplexF64}}; cache::ExponentialUtilities.ExpvCache{ComplexF64})
   @ ExponentialUtilities C:\Users\Sam\.julia\packages\ExponentialUtilities\vU95t\src\krylov_phiv.jl:85
 [2] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{LawsonEuler{DataType}, true, Vector{ComplexF64}, Nothing, Float64, Int64, Float64, Float64, Float64, Float64, Vector{Vector{ComplexF64}}, ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Int64, SplitFunction{true, ODEFunction{true, DiffEqArrayOperator{ComplexF64, Diagonal{ComplexF64, Vector{ComplexF64}}, typeof(SciMLBase.DEFAULT_UPDATE_FUNC)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, typeof(nlfunc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SplitODEProblem{true}}, LawsonEuler{DataType}, OrdinaryDiffEq.InterpolationData{SplitFunction{true, ODEFunction{true, DiffEqArrayOperator{ComplexF64, Diagonal{ComplexF64, Vector{ComplexF64}}, typeof(SciMLBase.DEFAULT_UPDATE_FUNC)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, typeof(nlfunc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.LawsonEulerCache{Vector{ComplexF64}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Tuple{ExponentialUtilities.KrylovSubspace{ComplexF64, ComplexF64, Float64, Matrix{ComplexF64}, Matrix{ComplexF64}}, ExponentialUtilities.ExpvCache{ComplexF64}}}}, DiffEqBase.DEStats}, SplitFunction{true, ODEFunction{true, DiffEqArrayOperator{ComplexF64, Diagonal{ComplexF64, Vector{ComplexF64}}, typeof(SciMLBase.DEFAULT_UPDATE_FUNC)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, typeof(nlfunc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEq.LawsonEulerCache{Vector{ComplexF64}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Tuple{ExponentialUtilities.KrylovSubspace{ComplexF64, ComplexF64, Float64, Matrix{ComplexF64}, Matrix{ComplexF64}}, ExponentialUtilities.ExpvCache{ComplexF64}}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Int64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{ComplexF64}, ComplexF64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.LawsonEulerCache{Vector{ComplexF64}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Tuple{ExponentialUtilities.KrylovSubspace{ComplexF64, ComplexF64, Float64, Matrix{ComplexF64}, Matrix{ComplexF64}}, ExponentialUtilities.ExpvCache{ComplexF64}}}, repeat_step::Bool)
   @ OrdinaryDiffEq C:\Users\Sam\.julia\packages\OrdinaryDiffEq\a5qvy\src\perform_step\exponential_rk_perform_step.jl:86
 [3] perform_step!
   @ C:\Users\Sam\.julia\packages\OrdinaryDiffEq\a5qvy\src\perform_step\exponential_rk_perform_step.jl:71 [inlined]
 [4] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{LawsonEuler{DataType}, true, Vector{ComplexF64}, Nothing, Float64, Int64, Float64, Float64, Float64, Float64, Vector{Vector{ComplexF64}}, ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Int64, SplitFunction{true, ODEFunction{true, DiffEqArrayOperator{ComplexF64, Diagonal{ComplexF64, Vector{ComplexF64}}, typeof(SciMLBase.DEFAULT_UPDATE_FUNC)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, typeof(nlfunc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SplitODEProblem{true}}, LawsonEuler{DataType}, OrdinaryDiffEq.InterpolationData{SplitFunction{true, ODEFunction{true, DiffEqArrayOperator{ComplexF64, Diagonal{ComplexF64, Vector{ComplexF64}}, typeof(SciMLBase.DEFAULT_UPDATE_FUNC)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, typeof(nlfunc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.LawsonEulerCache{Vector{ComplexF64}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Tuple{ExponentialUtilities.KrylovSubspace{ComplexF64, ComplexF64, Float64, Matrix{ComplexF64}, Matrix{ComplexF64}}, ExponentialUtilities.ExpvCache{ComplexF64}}}}, DiffEqBase.DEStats}, SplitFunction{true, ODEFunction{true, DiffEqArrayOperator{ComplexF64, Diagonal{ComplexF64, Vector{ComplexF64}}, typeof(SciMLBase.DEFAULT_UPDATE_FUNC)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, typeof(nlfunc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEq.LawsonEulerCache{Vector{ComplexF64}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Tuple{ExponentialUtilities.KrylovSubspace{ComplexF64, ComplexF64, Float64, Matrix{ComplexF64}, Matrix{ComplexF64}}, ExponentialUtilities.ExpvCache{ComplexF64}}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Int64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{ComplexF64}, ComplexF64, Nothing, OrdinaryDiffEq.DefaultInit})
   @ OrdinaryDiffEq C:\Users\Sam\.julia\packages\OrdinaryDiffEq\a5qvy\src\solve.jl:478
 [5] __solve(::ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Int64, SplitFunction{true, ODEFunction{true, DiffEqArrayOperator{ComplexF64, Diagonal{ComplexF64, Vector{ComplexF64}}, typeof(SciMLBase.DEFAULT_UPDATE_FUNC)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, typeof(nlfunc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SplitODEProblem{true}}, ::LawsonEuler{DataType}; kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})
   @ OrdinaryDiffEq C:\Users\Sam\.julia\packages\OrdinaryDiffEq\a5qvy\src\solve.jl:5
 [6] #solve_call#58
   @ C:\Users\Sam\.julia\packages\DiffEqBase\oe7VF\src\solve.jl:61 [inlined]
 [7] #solve_up#60
   @ C:\Users\Sam\.julia\packages\DiffEqBase\oe7VF\src\solve.jl:82 [inlined]
 [8] #solve#59
   @ C:\Users\Sam\.julia\packages\DiffEqBase\oe7VF\src\solve.jl:70 [inlined]
 [9] top-level scope
   @ C:\Users\Sam\Documents\GitHub\julia\minrep.jl:21
in expression starting at C:\Users\Sam\Documents\GitHub\julia\minrep.jl:21

Trying 1.7 next

[Lazersmoke]

I get the same MethodError regardless of having the set_threads line or not in 1.7.0-beta3. So it looks like an ExponentialUtilities issue, in that they try to support complexes, but (in this version eigen! only works with real floats)