diff --git a/NEWS.md b/NEWS.md index 4912b6f3358c7..6f8fb82d0cd5c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -37,6 +37,7 @@ Standard library changes * Added keyword arguments `rtol`, `atol` to `pinv` and `nullspace` ([#29998]). * `UniformScaling` instances are now callable such that e.g. `I(3)` will produce a `Diagonal` matrix ([#30298]). +* Eigenvalues λ of general matrices are now sorted lexicographically by (Re λ, Im λ) ([#21598]). #### SparseArrays diff --git a/base/combinatorics.jl b/base/combinatorics.jl index f19a2938b9b5b..af2240dc93cc1 100644 --- a/base/combinatorics.jl +++ b/base/combinatorics.jl @@ -63,6 +63,37 @@ isperm(p::Tuple{}) = true isperm(p::Tuple{Int}) = p[1] == 1 isperm(p::Tuple{Int,Int}) = ((p[1] == 1) & (p[2] == 2)) | ((p[1] == 2) & (p[2] == 1)) +# swap columns i and j of a, in-place +function swapcols!(a::AbstractMatrix, i, j) + i == j && return + cols = axes(a,2) + @boundscheck i in cols || throw(BoundsError(a, (:,i))) + @boundscheck j in cols || throw(BoundsError(a, (:,j))) + for k in axes(a,1) + @inbounds a[k,i],a[k,j] = a[k,j],a[k,i] + end +end +# like permute!! applied to each row of a, in-place in a (overwriting p). +function permutecols!!(a::AbstractMatrix, p::AbstractVector{<:Integer}) + require_one_based_indexing(a, p) + count = 0 + start = 0 + while count < length(p) + ptr = start = findnext(!iszero, p, start+1)::Int + next = p[start] + count += 1 + while next != start + swapcols!(a, ptr, next) + p[ptr] = 0 + ptr = next + next = p[next] + count += 1 + end + p[ptr] = 0 + end + a +end + function permute!!(a, p::AbstractVector{<:Integer}) require_one_based_indexing(a, p) count = 0 diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 0b245a5143716..01969e1104d51 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -38,13 +38,13 @@ julia> A = [-4. -17.; 2. 2.] julia> eigvals(A) 2-element Array{Complex{Float64},1}: - -1.0 + 5.0im -1.0 - 5.0im + -1.0 + 5.0im julia> eigvecs(A) 2×2 Array{Complex{Float64},2}: - 0.945905+0.0im 0.945905-0.0im - -0.166924-0.278207im -0.166924+0.278207im + 0.945905-0.0im 0.945905+0.0im + -0.166924+0.278207im -0.166924-0.278207im ``` In addition, Julia provides many [factorizations](@ref man-linalg-factorizations) which can be used to diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 6ea49ad50479e..58126a45d7d9f 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -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 diff --git a/stdlib/LinearAlgebra/src/eigen.jl b/stdlib/LinearAlgebra/src/eigen.jl index fc6cfaa858aea..841abd22a80c8 100644 --- a/stdlib/LinearAlgebra/src/eigen.jl +++ b/stdlib/LinearAlgebra/src/eigen.jl @@ -27,18 +27,32 @@ Base.iterate(S::Union{Eigen,GeneralizedEigen}, ::Val{:done}) = nothing isposdef(A::Union{Eigen,GeneralizedEigen}) = isreal(A.values) && all(x -> x > 0, A.values) +# pick a canonical ordering to avoid returning eigenvalues in "random" order +# as is the LAPACK default (for complex λ — LAPACK sorts by λ for the Hermitian/Symmetric case) +eigsortby(λ::Real) = λ +eigsortby(λ::Complex) = (real(λ),imag(λ)) +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!(λ::AbstractVector, sortby::Union{Function,Nothing}=eigsortby) = sortby === nothing ? λ : sort!(λ, by=sortby) + """ - eigen!(A, [B]) + eigen!(A, [B]; permute, scale, sortby) Same as [`eigen`](@ref), but saves space by overwriting the input `A` (and `B`), instead of creating a copy. """ -function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T<:BlasReal +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(WR, VR) + iszero(WI) && return Eigen(sorteig!(WR, VR, sortby)...) evec = zeros(Complex{T}, n, n) j = 1 while j <= n @@ -53,18 +67,19 @@ function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where end j += 1 end - return Eigen(complex.(WR, WI), evec) + return Eigen(sorteig!(complex.(WR, WI), evec, sortby)...) end -function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T<:BlasComplex +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)) - return Eigen(LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)[[2,4]]...) + 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 """ - eigen(A; permute::Bool=true, scale::Bool=true) -> Eigen + eigen(A; permute::Bool=true, scale::Bool=true, sortby) -> Eigen Computes the eigenvalue decomposition of `A`, returning an `Eigen` factorization object `F` which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the @@ -79,6 +94,12 @@ before the eigenvector calculation. The option `permute=true` permutes the matri closer to upper triangular, and `scale=true` scales the matrix by its diagonal elements to make rows and columns more equal in norm. The default is `true` for both options. +By default, the eigenvalues and vectors are sorted lexicographically by `(real(λ),imag(λ))`. +A different comparison function `by(λ)` can be passed to `sortby`, or you can pass +`sortby=nothing` to leave the eigenvalues in an arbitrary order. Some special matrix types +(e.g. `Diagonal` or `SymTridiagonal`) may implement their own sorting convention and not +accept a `sortby` keyword. + # Examples ```jldoctest julia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0]) @@ -112,18 +133,18 @@ julia> vals == F.values && vecs == F.vectors true ``` """ -function eigen(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) 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), permute = permute, scale = scale) - return eigen!(AA, permute = permute, scale = scale) + isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby) + return eigen!(AA; permute=permute, scale=scale, sortby=sortby) end eigen(x::Number) = Eigen([x], fill(one(x), 1, 1)) """ - eigvecs(A; permute::Bool=true, scale::Bool=true) -> Matrix + eigvecs(A; permute::Bool=true, scale::Bool=true, `sortby`) -> Matrix Return a matrix `M` whose columns are the eigenvectors of `A`. (The `k`th eigenvector can -be obtained from the slice `M[:, k]`.) The `permute` and `scale` keywords are the same as +be obtained from the slice `M[:, k]`.) The `permute`, `scale`, and `sortby` keywords are the same as for [`eigen`](@ref). # Examples @@ -135,19 +156,17 @@ julia> eigvecs([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0]) 0.0 0.0 1.0 ``` """ -eigvecs(A::Union{Number, AbstractMatrix}; permute::Bool=true, scale::Bool=true) = - eigvecs(eigen(A, permute=permute, scale=scale)) +eigvecs(A::Union{Number, AbstractMatrix}; kws...) = + eigvecs(eigen(A; kws...)) eigvecs(F::Union{Eigen, GeneralizedEigen}) = F.vectors eigvals(F::Union{Eigen, GeneralizedEigen}) = F.values """ - eigvals!(A; permute::Bool=true, scale::Bool=true) -> values + eigvals!(A; permute::Bool=true, scale::Bool=true, sortby) -> values Same as [`eigvals`](@ref), but saves space by overwriting the input `A`, instead of creating a copy. -The option `permute=true` permutes the matrix to become -closer to upper triangular, and `scale=true` scales the matrix by its diagonal elements to -make rows and columns more equal in norm. +The `permute`, `scale`, and `sortby` keywords are the same as for [`eigen`](@ref). !!! note The input matrix `A` will not contain its eigenvalues after `eigvals!` is @@ -171,29 +190,27 @@ julia> A 0.0 5.37228 ``` """ -function eigvals!(A::StridedMatrix{<:BlasReal}; permute::Bool=true, scale::Bool=true) - issymmetric(A) && return eigvals!(Symmetric(A)) +function eigvals!(A::StridedMatrix{<:BlasReal}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) + issymmetric(A) && return sorteig!(eigvals!(Symmetric(A)), sortby) _, valsre, valsim, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'N', 'N', A) - return iszero(valsim) ? valsre : complex.(valsre, valsim) + return sorteig!(iszero(valsim) ? valsre : complex.(valsre, valsim), sortby) end -function eigvals!(A::StridedMatrix{<:BlasComplex}; permute::Bool=true, scale::Bool=true) - ishermitian(A) && return eigvals(Hermitian(A)) - return LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'N', 'N', A)[2] +function eigvals!(A::StridedMatrix{<:BlasComplex}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) + 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 # promotion type to use for eigenvalues of a Matrix{T} eigtype(T) = promote_type(Float32, typeof(zero(T)/sqrt(abs2(one(T))))) """ - eigvals(A; permute::Bool=true, scale::Bool=true) -> values + eigvals(A; permute::Bool=true, scale::Bool=true, sortby) -> values Return the eigenvalues of `A`. For general non-symmetric matrices it is possible to specify how the matrix is balanced -before the eigenvalue calculation. The option `permute=true` permutes the matrix to -become closer to upper triangular, and `scale=true` scales the matrix by its diagonal -elements to make rows and columns more equal in norm. The default is `true` for both -options. +before the eigenvalue calculation. The `permute`, `scale`, and `sortby` keywords are +the same as for [`eigen!`](@ref). # Examples ```jldoctest @@ -208,8 +225,8 @@ julia> eigvals(diag_matrix) 4.0 ``` """ -eigvals(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T = - eigvals!(copy_oftype(A, eigtype(T)), permute = permute, scale = scale) +eigvals(A::StridedMatrix{T}; kws...) where T = + eigvals!(copy_oftype(A, eigtype(T)); kws...) """ For a scalar input, `eigvals` will return a scalar. @@ -309,11 +326,11 @@ inv(A::Eigen) = A.vectors * inv(Diagonal(A.values)) / A.vectors det(A::Eigen) = prod(A.values) # Generalized eigenproblem -function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasReal - issymmetric(A) && isposdef(B) && return eigen!(Symmetric(A), Symmetric(B)) +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), sortby=sortby) 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 @@ -329,13 +346,13 @@ function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasReal end j += 1 end - return GeneralizedEigen(complex.(alphar, alphai)./beta, vecs) + return GeneralizedEigen(sorteig!(complex.(alphar, alphai)./beta, vecs, sortby)...) end -function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasComplex - ishermitian(A) && isposdef(B) && return eigen!(Hermitian(A), Hermitian(B)) +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), sortby=sortby) alpha, beta, _, vr = LAPACK.ggev!('N', 'V', A, B) - return GeneralizedEigen(alpha./beta, vr) + return GeneralizedEigen(sorteig!(alpha./beta, vr, sortby)...) end """ @@ -348,6 +365,9 @@ Computes the generalized eigenvalue decomposition of `A` and `B`, returning a Iterating the decomposition produces the components `F.values` and `F.vectors`. +Any keyword arguments passed to `eigen` are passed through to the lower-level +[`eigen!`](@ref) function. + # Examples ```jldoctest julia> A = [1 0; 0 -1] @@ -364,13 +384,13 @@ julia> F = eigen(A, B); julia> F.values 2-element Array{Complex{Float64},1}: - 0.0 + 1.0im 0.0 - 1.0im + 0.0 + 1.0im julia> F.vectors 2×2 Array{Complex{Float64},2}: - 0.0-1.0im 0.0+1.0im - -1.0-0.0im -1.0+0.0im + 0.0+1.0im 0.0-1.0im + -1.0+0.0im -1.0-0.0im julia> vals, vecs = F; # destructuring via iteration @@ -378,15 +398,15 @@ julia> vals == F.values && vecs == F.vectors true ``` """ -function eigen(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) where {TA,TB} +function eigen(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB} S = promote_type(eigtype(TA),TB) - return eigen!(copy_oftype(A, S), copy_oftype(B, S)) + eigen!(copy_oftype(A, S), copy_oftype(B, S); kws...) end eigen(A::Number, B::Number) = eigen(fill(A,1,1), fill(B,1,1)) """ - eigvals!(A, B) -> values + eigvals!(A, B; sortby) -> values Same as [`eigvals`](@ref), but saves space by overwriting the input `A` (and `B`), instead of creating copies. @@ -409,8 +429,8 @@ julia> B = [0. 1.; 1. 0.] julia> eigvals!(A, B) 2-element Array{Complex{Float64},1}: - 0.0 + 1.0im 0.0 - 1.0im + 0.0 + 1.0im julia> A 2×2 Array{Float64,2}: @@ -423,15 +443,15 @@ julia> B 0.0 1.0 ``` """ -function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasReal - issymmetric(A) && isposdef(B) && return eigvals!(Symmetric(A), Symmetric(B)) +function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasReal + issymmetric(A) && isposdef(B) && return sorteig!(eigvals!(Symmetric(A), Symmetric(B)), sortby) alphar, alphai, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B) - return (iszero(alphai) ? alphar : complex.(alphar, alphai))./beta + return sorteig!((iszero(alphai) ? alphar : complex.(alphar, alphai))./beta, sortby) end -function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasComplex - ishermitian(A) && isposdef(B) && return eigvals!(Hermitian(A), Hermitian(B)) +function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex + ishermitian(A) && isposdef(B) && return sorteig!(eigvals!(Hermitian(A), Hermitian(B)), sortby) alpha, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B) - alpha./beta + return sorteig!(alpha./beta, sortby) end """ @@ -453,13 +473,13 @@ julia> B = [0 1; 1 0] julia> eigvals(A,B) 2-element Array{Complex{Float64},1}: - 0.0 + 1.0im 0.0 - 1.0im + 0.0 + 1.0im ``` """ -function eigvals(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) where {TA,TB} +function eigvals(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB} S = promote_type(eigtype(TA),TB) - return eigvals!(copy_oftype(A, S), copy_oftype(B, S)) + return eigvals!(copy_oftype(A, S), copy_oftype(B, S); kws...) end """ @@ -482,11 +502,11 @@ julia> B = [0 1; 1 0] julia> eigvecs(A, B) 2×2 Array{Complex{Float64},2}: - 0.0-1.0im 0.0+1.0im - -1.0-0.0im -1.0+0.0im + 0.0+1.0im 0.0-1.0im + -1.0+0.0im -1.0-0.0im ``` """ -eigvecs(A::AbstractMatrix, B::AbstractMatrix) = eigvecs(eigen(A, B)) +eigvecs(A::AbstractMatrix, B::AbstractMatrix; kws...) = eigvecs(eigen(A, B; kws...)) function show(io::IO, mime::MIME{Symbol("text/plain")}, F::Union{Eigen,GeneralizedEigen}) summary(io, F); println(io) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 3a08824f66fb6..977768f4478f8 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -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)...) @@ -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} = diff --git a/stdlib/LinearAlgebra/test/bidiag.jl b/stdlib/LinearAlgebra/test/bidiag.jl index 0c127ff8dbd48..1a1dcd4fa0ba7 100644 --- a/stdlib/LinearAlgebra/test/bidiag.jl +++ b/stdlib/LinearAlgebra/test/bidiag.jl @@ -252,7 +252,7 @@ Random.seed!(1) @testset "Eigensystems" begin if relty <: AbstractFloat d1, v1 = eigen(T) - d2, v2 = eigen(map(elty<:Complex ? ComplexF64 : Float64,Tfull)) + d2, v2 = eigen(map(elty<:Complex ? ComplexF64 : Float64,Tfull), sortby=nothing) @test (uplo == :U ? d1 : reverse(d1)) ≈ d2 if elty <: Real test_approx_eq_modphase(v1, uplo == :U ? v2 : v2[:,n:-1:1]) diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index d7cf098e42fe7..0d2123db8931e 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -304,8 +304,7 @@ end @testset "svdvals and eigvals (#11120/#11247)" begin D = Diagonal(Matrix{Float64}[randn(3,3), randn(2,2)]) @test sort([svdvals(D)...;], rev = true) ≈ svdvals([D.diag[1] zeros(3,2); zeros(2,3) D.diag[2]]) - @test [eigvals(D)...;] ≈ eigvals([D.diag[1] zeros(3,2); zeros(2,3) D.diag[2]]) - + @test sort([eigvals(D)...;], by=LinearAlgebra.eigsortby) ≈ eigvals([D.diag[1] zeros(3,2); zeros(2,3) D.diag[2]]) end @testset "eigmin (#27847)" begin @@ -501,5 +500,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 diff --git a/stdlib/LinearAlgebra/test/lapack.jl b/stdlib/LinearAlgebra/test/lapack.jl index 610d7f0e5a96d..47037ced3a665 100644 --- a/stdlib/LinearAlgebra/test/lapack.jl +++ b/stdlib/LinearAlgebra/test/lapack.jl @@ -231,7 +231,7 @@ end @testset for elty in (ComplexF32, ComplexF64) A = rand(elty,10,10) Aw, Avl, Avr = LAPACK.geev!('N','V',copy(A)) - fA = eigen(A) + fA = eigen(A, sortby=nothing) @test fA.values ≈ Aw @test fA.vectors ≈ Avr end @@ -561,18 +561,18 @@ end @testset "trrfs & trevc" begin @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) T = triu(rand(elty,10,10)) - S = copy(T) + v = eigvecs(T, sortby=nothing)[:,1] select = zeros(LinearAlgebra.BlasInt,10) select[1] = 1 select,Vr = LAPACK.trevc!('R','S',select,copy(T)) - @test Vr ≈ eigvecs(S)[:,1] + @test Vr ≈ v select = zeros(LinearAlgebra.BlasInt,10) select[1] = 1 select,Vl = LAPACK.trevc!('L','S',select,copy(T)) select = zeros(LinearAlgebra.BlasInt,10) select[1] = 1 select,Vln,Vrn = LAPACK.trevc!('B','S',select,copy(T)) - @test Vrn ≈ eigvecs(S)[:,1] + @test Vrn ≈ v @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))