Dimension mismatch with SymTridiagonal eigenvalue decomposition

[mforets]

julia> A = SymTridiagonal(fill(1.0, 3), fill(-1.0, 3))
3×3 SymTridiagonal{Float64, Vector{Float64}}:
  1.0  -1.0    ⋅ 
 -1.0   1.0  -1.0
   ⋅   -1.0   1.0

julia> eigen(Matrix(A))  # works
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
 -0.4142135623730927
  1.0000000000000013
  2.414213562373095
vectors:
3×3 Matrix{Float64}:
 0.5        0.707107      0.5
 0.707107  -7.85046e-16  -0.707107
 0.5       -0.707107      0.5

julia> eigen(A) # unexpected error.. (?)
ERROR: DimensionMismatch("ev has length 3 but needs one less than dv's length, 3)")
Stacktrace:
 [1] stegr!(jobz::Char, range::Char, dv::Vector{Float64}, ev::Vector{Float64}, vl::Float64, vu::Float64, il::Int64, iu::Int64)
   @ LinearAlgebra.LAPACK /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lapack.jl:3793
 [2] stegr!
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lapack.jl:3881 [inlined]
 [3] eigen!
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:280 [inlined]
 [4] eigen(A::SymTridiagonal{Float64, Vector{Float64}})
   @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:281
 [5] top-level scope
   @ REPL[28]:1

julia> versioninfo()
Julia Version 1.6.0-DEV.1699
Commit 8a5ac94e96* (2020-12-08 03:24 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-8705G CPU @ 3.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.0 (ORCJIT, skylake)

[mforets]

The question is whether eigen(A) should give a better error message? Or the constructor of SymTridiagonal should require that ev has dimension 2?

[dkarrasch]

I think this can be fixed "easily". In the LAPACK calling method, we have this:

https://github.com/JuliaLang/julia/blob/ce795bc73e3a7fcd4c338b4cf8edd5446748db40/stdlib/LinearAlgebra/src/lapack.jl#L3788-L3795

So we first check whether the length of ev is shorter by one than dv, and then append a zero in eev. I believe we should simply check whether ev already has the same length as dv, and make eev a copy of ev with zero last element in that case. I believe it is for this reason that we allow ev to have the same length as dv in the construction of SymTridiagonal, but then happen to be too strict in lapack.jl. Would you be able/willing to take a stab at this?

[dkarrasch]

As for the other symtriag methods near my quote, we could be more permissive there as well. It seems that those method do require length == n-1. If that's not satisfied, I think we should then simply pass ev[1:end-1] instead of ev.

[mforets]

Thanks for your detailed response! I'm willing to send a patch, though I've never contributed (and built Julia..).

[dkarrasch]

You may find useful information here: https://github.com/JuliaLang/julia/blob/master/CONTRIBUTING.md

[mforets]

It seems that those method do require length == n-1. If that's not satisfied, I think we should then simply pass ev[1:end-1] instead of ev.

Is it along these lines what you had in mind?

function eigvals!(A::SymTridiagonal{<:BlasReal})
    n = length(A.dv)
    ne = length(A.ev)
    if ne == n-1
        ev = A.ev
    elseif ne == n
        ev = view(A.ev, 1:n)
    else
        throw(...)
    end
    LAPACK.stev!('N', A.dv, ev)[1]
end