Jutho / KrylovKit.jl

Krylov methods for linear problems, eigenvalues, singular values and matrix functions
Other
284 stars 37 forks source link

Access to undefined reference when using KrylovKit on a non-concrete vector #69

Open MasonProtter opened 1 year ago

MasonProtter commented 1 year ago

Someone on Slack ran into this issue when they accidentally ran

exponentiate(::SparseMatrixCSC{ComplexF32, Int64}, 1, ::Vector{Complex})

instead of using a concrete type.

Here's a MWE:

#+begin_src julia
using KrylovKit, SparseArrays
let H = sprand(1000, 1000, 0.001), v = zeros(Complex, 1000)
    v[1] = 1
    exponentiate(H, 1, v)
end
#+end_src

#+RESULTS:
UndefRefError: access to undefined reference

Stacktrace:
  [1] getindex
    @ ./essentials.jl:13 [inlined]
  [2] macro expansion
    @ ~/julia/usr/share/julia/stdlib/v1.9/LinearAlgebra/src/generic.jl:183 [inlined]
  [3] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [4] rmul!(X::Vector{Complex}, s::Bool)
    @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.9/LinearAlgebra/src/generic.jl:182
  [5] expintegrator
    @ ~/.julia/packages/KrylovKit/diNbc/src/matrixfun/expintegrator.jl:106 [inlined]
  [6] expintegrator(::SparseMatrixCSC{Float64, Int64}, ::Int64, ::Vector{Complex}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/matrixfun/expintegrator.jl:102
  [7] expintegrator
    @ ~/.julia/packages/KrylovKit/diNbc/src/matrixfun/expintegrator.jl:98 [inlined]
  [8] #exponentiate#82
    @ ~/.julia/packages/KrylovKit/diNbc/src/matrixfun/exponentiate.jl:79 [inlined]
  [9] exponentiate(A::SparseMatrixCSC{Float64, Int64}, t::Int64, v::Vector{Complex})
    @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/matrixfun/exponentiate.jl:79
 [10] top-level scope
    @ In[17]:4
Jutho commented 1 year ago

Do I infer correctly from the stack trace that this is actually a Julia Base (or actually LinearAlgebra) bug, namely with rmul!(X::Vector{Complex}, s::Bool)?

MasonProtter commented 1 year ago

Yeah I find that quite strange. I couldn't reproduce the bug by just calling rmul!(::Vector{Complex}, ::Bool), so I think that what’s happening is that KrylovKit is generating undefined junk somewhere and then trying to rmul! that junk

MasonProtter commented 1 year ago

Yeah, so if I make an undef vector of Complex then I get the error

julia> v = Vector{Complex}(undef, 2)
2-element Vector{Complex}:
 #undef
 #undef

julia> rmul!(v, false)
ERROR: UndefRefError: access to undefined reference
Stacktrace:
 [1] getindex
   @ ./essentials.jl:13 [inlined]
 [2] macro expansion
   @ ~/julia/usr/share/julia/stdlib/v1.9/LinearAlgebra/src/generic.jl:183 [inlined]
 [3] macro expansion
   @ ./simdloop.jl:77 [inlined]
 [4] rmul!(X::Vector{Complex}, s::Bool)
   @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.9/LinearAlgebra/src/generic.jl:182
 [5] top-level scope
   @ REPL[2]:1

So maybe the fix would just be to throw a more helpful error before actually hitting rmul!?

Jutho commented 1 year ago

The problem is that the implementation tries to be completely agnostic about the type of vector. So checking for undef or abstract type parameters is not really possible, unless by adding a special branch for Vector instances.