JuliaLinearAlgebra / BandedMatrices.jl

A Julia package for representing banded matrices
http://julialinearalgebra.github.io/BandedMatrices.jl/
Other
128 stars 38 forks source link

Unnecessary restiction of bandwidths in generalized eigenvalues #203

Closed KlausC closed 3 years ago

KlausC commented 4 years ago

In the generalized eigenvalue problem eigvals(A, B) with banded matrices, the bandwidth of B must not exceed that of A. There is no mathematical justification for that. It would be easy to extend the bandwidth in the used workspace to the maximum of both.

julia> A = Symmetric(BandedMatrix(0=>ones(3)))
3×3 Symmetric{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}:
 1.0   ⋅    ⋅ 
  ⋅   1.0   ⋅ 
  ⋅    ⋅   1.0

julia> B = Symmetric(BandedMatrix(0=>ones(3), 1=>ones(2)))
3×3 Symmetric{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}:
 1.0  1.0   ⋅ 
 1.0  1.0  1.0
  ⋅   1.0  1.0

julia> eigvals(A, B)
ERROR: LAPACKException(2)
Stacktrace:
 [1] chklapackerror at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lapack.jl:38 [inlined]
 [2] pbstf!(::Char, ::Int64, ::Int64, ::SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},Base.Slice{Base.OneTo{Int64}}},false}) at /home/crusius/.julia/dev/BandedMatrices/src/lapack.jl:282
 [3] eigvals!(::Symmetric{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}, ::Symmetric{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}) at /home/crusius/.julia/dev/BandedMatrices/src/symbanded/symbanded.jl:151
 [4] eigvals(::Symmetric{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}, ::Symmetric{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}) at /home/crusius/.julia/dev/BandedMatrices/src/symbanded/symbanded.jl:162
 [5] top-level scope at REPL[29]:1

julia> eigvals(B, A)
3-element Array{Float64,1}:
 -0.41421356237309503
  0.9999999999999999
  2.414213562373095
dlfivefifty commented 4 years ago

You are right, can you produce a PR? Probably in the eigvals overload