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.

Sign up for GitHub

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
Merged

RFC: sort eigenvalues in a canonical order #21598

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
726b309
sort eigenvalues in a canonical order
stevengj Jan 10, 2019
c493100
no sortby for special matrix types
stevengj Jan 10, 2019
75cc29d
simplification
stevengj Jan 10, 2019
e639171
fixes
stevengj Jan 24, 2019
c0576ea
pass through keyword args from eigen to eigen! and similar
stevengj Jan 31, 2019
8c6072e
more fixes
stevengj Feb 1, 2019
bc1ce98
whoops
stevengj Feb 1, 2019
167424e
fix bidiag test
stevengj Feb 1, 2019
c3fc98d
simplify lapack test fixes
stevengj Feb 1, 2019
79fb517
don't discard sortby if the matrix happens to be hermitian
stevengj Feb 1, 2019
cf08bda
missing check in permutecols
stevengj Feb 1, 2019
e37586e
Merge branch 'master' into eigsort
stevengj Feb 1, 2019
d6acbb0
test fixes
stevengj Feb 4, 2019
118bcf5
Merge branch 'eigsort' of git://github.com/stevengj/julia into eigsort
stevengj Feb 4, 2019
1b913d5
NEWS [ci skip]
stevengj Feb 5, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
don't discard sortby if the matrix happens to be hermitian
  • Loading branch information
@stevengj
stevengj committed Feb 1, 2019
commit 79fb517ebcda4998a878664c3407aca11ed6bcd6
16 changes: 8 additions & 8 deletions stdlib/LinearAlgebra/src/eigen.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ Same as [`eigen`](@ref), but saves space by overwriting the input `A` (and
function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T<:BlasReal
n = size(A, 2)
n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
issymmetric(A) && return eigen!(Symmetric(A))
issymmetric(A) && return eigen!(Symmetric(A), sortby=sortby)
A, WR, WI, VL, VR, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)
iszero(WI) && return Eigen(sorteig!(WR, VR, sortby)...)
evec = zeros(Complex{T}, n, n)
Expand All @@ -73,7 +73,7 @@ end
function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex
n = size(A, 2)
n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
ishermitian(A) && return eigen!(Hermitian(A))
ishermitian(A) && return eigen!(Hermitian(A), sortby=sortby)
eval, evec = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)[[2,4]]
return Eigen(sorteig!(eval, evec, sortby)...)
end
Expand Down Expand Up @@ -191,12 +191,12 @@ julia> A
```
"""
function eigvals!(A::StridedMatrix{<:BlasReal}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby)
issymmetric(A) && return eigvals!(Symmetric(A))
issymmetric(A) && return sorteig!(eigvals!(Symmetric(A)), sortby)
_, valsre, valsim, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'N', 'N', A)
return sorteig!(iszero(valsim) ? valsre : complex.(valsre, valsim), sortby)
end
function eigvals!(A::StridedMatrix{<:BlasComplex}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby)
ishermitian(A) && return eigvals(Hermitian(A))
ishermitian(A) && return sorteig!(eigvals(Hermitian(A)), sortby)
return sorteig!(LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'N', 'N', A)[2], sortby)
end

Expand Down Expand Up @@ -327,7 +327,7 @@ det(A::Eigen) = prod(A.values)

# Generalized eigenproblem
function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasReal
issymmetric(A) && isposdef(B) && return eigen!(Symmetric(A), Symmetric(B))
issymmetric(A) && isposdef(B) && return eigen!(Symmetric(A), Symmetric(B), sortby=sortby)
n = size(A, 1)
alphar, alphai, beta, _, vr = LAPACK.ggev!('N', 'V', A, B)
iszero(alphai) && return GeneralizedEigen(sorteig!(alphar ./ beta, vr, sortby)...)
Expand All @@ -350,7 +350,7 @@ function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function
end

function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex
ishermitian(A) && isposdef(B) && return eigen!(Hermitian(A), Hermitian(B))
ishermitian(A) && isposdef(B) && return eigen!(Hermitian(A), Hermitian(B), sortby=sortby)
alpha, beta, _, vr = LAPACK.ggev!('N', 'V', A, B)
return GeneralizedEigen(sorteig!(alpha./beta, vr, sortby)...)
end
Expand Down Expand Up @@ -444,12 +444,12 @@ julia> B
```
"""
function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasReal
issymmetric(A) && isposdef(B) && return eigvals!(Symmetric(A), Symmetric(B))
issymmetric(A) && isposdef(B) && return sorteig!(eigvals!(Symmetric(A), Symmetric(B)), sortby)
alphar, alphai, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B)
return sorteig!((iszero(alphai) ? alphar : complex.(alphar, alphai))./beta, sortby)
end
function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex
ishermitian(A) && isposdef(B) && return eigvals!(Hermitian(A), Hermitian(B))
ishermitian(A) && isposdef(B) && return sorteig!(eigvals!(Hermitian(A), Hermitian(B)), sortby)
alpha, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B)
return sorteig!(alpha./beta, sortby)
end
Expand Down
14 changes: 7 additions & 7 deletions stdlib/LinearAlgebra/src/symmetric.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -506,12 +506,12 @@ end
inv(A::Hermitian{<:Any,<:StridedMatrix}) = Hermitian(_inv(A), sym_uplo(A.uplo))
inv(A::Symmetric{<:Any,<:StridedMatrix}) = Symmetric(_inv(A), sym_uplo(A.uplo))

eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)...)
eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) = Eigen(sorteig!(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)..., sortby)...)

function eigen(A::RealHermSymComplexHerm)
function eigen(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing)
T = eltype(A)
S = eigtype(T)
eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A))
eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), sortby=sortby)
end

eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, irange::UnitRange) = Eigen(LAPACK.syevr!('V', 'I', A.uplo, A.data, 0.0, 0.0, irange.start, irange.stop, -1.0)...)
Expand Down Expand Up @@ -656,13 +656,13 @@ end
eigmax(A::RealHermSymComplexHerm{<:Real,<:StridedMatrix}) = eigvals(A, size(A, 1):size(A, 1))[1]
eigmin(A::RealHermSymComplexHerm{<:Real,<:StridedMatrix}) = eigvals(A, 1:1)[1]

function eigen!(A::HermOrSym{T,S}, B::HermOrSym{T,S}) where {T<:BlasReal,S<:StridedMatrix}
function eigen!(A::HermOrSym{T,S}, B::HermOrSym{T,S}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasReal,S<:StridedMatrix}
vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))
GeneralizedEigen(vals, vecs)
GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
end
function eigen!(A::Hermitian{T,S}, B::Hermitian{T,S}) where {T<:BlasComplex,S<:StridedMatrix}
function eigen!(A::Hermitian{T,S}, B::Hermitian{T,S}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasComplex,S<:StridedMatrix}
vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))
GeneralizedEigen(vals, vecs)
GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
end

eigvals!(A::HermOrSym{T,S}, B::HermOrSym{T,S}) where {T<:BlasReal,S<:StridedMatrix} =
Expand Down