Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

RFC: sort eigenvalues in a canonical order #21598

Merged
merged 15 commits into from
Feb 5, 2019
Prev Previous commit
Next Next commit
more fixes
stevengj committed Feb 1, 2019
commit 8c6072e0ab97721d224d6e982b707cba43a029dc
4 changes: 2 additions & 2 deletions stdlib/LinearAlgebra/src/diagonal.jl
Original file line number Diff line number Diff line change
@@ -528,11 +528,11 @@ eigvals(D::Diagonal{<:Number}; permute::Bool=true, scale::Bool=true) = D.diag
eigvals(D::Diagonal; permute::Bool=true, scale::Bool=true) =
[eigvals(x) for x in D.diag] #For block matrices, etc.
eigvecs(D::Diagonal) = Matrix{eltype(D)}(I, size(D))
function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true)
function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=nothing)
if any(!isfinite, D.diag)
throw(ArgumentError("matrix contains Infs or NaNs"))
end
Eigen(eigvals(D), eigvecs(D))
Eigen(sorteig!(eigvals(D), eigvecs(D), sortby)...)
end

#Singular system
12 changes: 6 additions & 6 deletions stdlib/LinearAlgebra/src/eigen.jl
Original file line number Diff line number Diff line change
@@ -31,15 +31,15 @@ isposdef(A::Union{Eigen,GeneralizedEigen}) = isreal(A.values) && all(x -> x > 0,
# as is the LAPACK default (for complex λ — LAPACK sorts by λ for the Hermitian/Symmetric case)
eigsortby::Real) = λ
eigsortby::Complex) = (real(λ),imag(λ))
function sorteig!(λ, X, sortby=eigsortby)
function sorteig!::AbstractVector, X::AbstractMatrix, sortby::Union{Function,Nothing}=eigsortby)
if sortby !== nothing && !issorted(λ, by=sortby)
p = sortperm(λ; alg=QuickSort, by=sortby)
permute!(λ, p)
Base.permutecols!!(X, p)
end
return λ, X
end
sorteig!(λ, sortby=eigsortby) = sortby === nothing ? λ : sort!(λ, by=sortby)
sorteig!::AbstractVector, sortby::Union{Function,Nothing}=eigsortby) = sortby === nothing ? λ : sort!(λ, by=sortby)

"""
eigen!(A, [B]; permute, scale, sortby)
@@ -52,7 +52,7 @@ function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortb
n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
issymmetric(A) && return eigen!(Symmetric(A))
A, WR, WI, VL, VR, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)
iszero(WI) && return Eigen(WR, VR)
iszero(WI) && return Eigen(sorteig!(WR, VR, sortby)...)
evec = zeros(Complex{T}, n, n)
j = 1
while j <= n
@@ -133,9 +133,9 @@ julia> vals == F.values && vecs == F.vectors
true
```
"""
function eigen(A::StridedMatrix{T}; kws...) where T
function eigen(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T
AA = copy_oftype(A, eigtype(T))
isdiag(AA) && return eigen(Diagonal(AA); kws...)
isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
return eigen!(AA; kws...)
end
eigen(x::Number) = Eigen([x], fill(one(x), 1, 1))
@@ -330,7 +330,7 @@ function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function
issymmetric(A) && isposdef(B) && return eigen!(Symmetric(A), Symmetric(B))
n = size(A, 1)
alphar, alphai, beta, _, vr = LAPACK.ggev!('N', 'V', A, B)
iszero(alphai) && return GeneralizedEigen(alphar ./ beta, vr)
iszero(alphai) && return GeneralizedEigen(sorteig!(alphar ./ beta, vr, sortby)...)

vecs = zeros(Complex{T}, n, n)
j = 1
9 changes: 9 additions & 0 deletions stdlib/LinearAlgebra/test/diagonal.jl
Original file line number Diff line number Diff line change
@@ -493,4 +493,13 @@ end
end
end

@testset "eigenvalue sorting" begin
D = Diagonal([0.4, 0.2, -1.3])
@test eigvals(D) == eigen(D).values == [0.4, 0.2, -1.3] # not sorted by default
@test eigvals(Matrix(D)) == eigen(Matrix(D)).values == [-1.3, 0.2, 0.4] # sorted even if diagonal special case is detected
E = eigen(D, sortby=abs) # sortby keyword supported for eigen(::Diagonal)
@test E.values == [0.2, 0.4, -1.3]
@test E.vectors == [0 1 0; 1 0 0; 0 0 1]
end

end # module TestDiagonal
8 changes: 4 additions & 4 deletions stdlib/LinearAlgebra/test/lapack.jl
Original file line number Diff line number Diff line change
@@ -563,19 +563,19 @@ end
@testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
T = triu(rand(elty,10,10))
i = sortperm(diag(T), by=LinearAlgebra.eigsortby)[1]
S = copy(T)
v = eigvecs(Matrix(T))[:,1]
kmax = sortperm(v, by=abs)[end]
select = zeros(LinearAlgebra.BlasInt,10)
select[i] = 1
select,Vr = LAPACK.trevc!('R','S',select,copy(T))
@test Vr eigvecs(S)[:,1]
@test Vr v * Vr[kmax] / v[kmax]
select = zeros(LinearAlgebra.BlasInt,10)
select[i] = 1
select,Vl = LAPACK.trevc!('L','S',select,copy(T))
select = zeros(LinearAlgebra.BlasInt,10)
select[i] = 1
select,Vln,Vrn = LAPACK.trevc!('B','S',select,copy(T))
v = eigvecs(S)[:,1]
@test Vrn v * Vrn[i] / v[i]
@test Vrn v * Vrn[kmax] / v[kmax]
@test Vln Vl
@test_throws ArgumentError LAPACK.trevc!('V','S',select,copy(T))
@test_throws DimensionMismatch LAPACK.trrfs!('U','N','N',T,rand(elty,10,10),rand(elty,10,11))