Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
31 changes: 15 additions & 16 deletions src/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1007,16 +1007,16 @@ end
_ortho_eltype(T) = Base.promote_op(/, T, T)
_ortho_eltype(T::Type{<:Number}) = typeof(one(T)/one(T))

# TODO Docstrings for eigvals, eigvecs, eigen all mention permute, scale, sortby as keyword args
# but not all of them below provide them. Do we need to fix that?
#Eigensystem
eigvals(D::Diagonal{<:Number}; permute::Bool=true, scale::Bool=true) = copy(D.diag)
eigvals(D::Diagonal; permute::Bool=true, scale::Bool=true) =
reduce(vcat, eigvals(x) for x in D.diag) #For block matrices, etc.
function eigvecs(D::Diagonal{T}) where {T<:AbstractMatrix}
diag_vecs = [ eigvecs(x) for x in D.diag ]
eigvals(D::Diagonal{<:Number}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) = sorteig!(copy(D.diag), sortby)
eigvals(D::Diagonal; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) =
sorteig!(reduce(vcat, eigvals(x; sortby=nothing) for x in D.diag), sortby) #For block matrices, etc.
function _eigen(D::Diagonal{T}) where {T<:AbstractMatrix}
facts = [eigen(x; sortby=nothing) for x in D.diag]
λ = reduce(vcat, f.values for f in facts)
diag_vecs = [f.vectors for f in facts]
matT = promote_type(map(typeof, diag_vecs)...)
ncols_diag = [ size(x, 2) for x in D.diag ]
ncols_diag = [size(x, 2) for x in D.diag]
nrows = size(D, 1)
vecs = Matrix{Vector{eltype(matT)}}(undef, nrows, sum(ncols_diag))
for j in axes(D, 2), i in axes(D, 1)
Expand All @@ -1031,14 +1031,14 @@ function eigvecs(D::Diagonal{T}) where {T<:AbstractMatrix}
end
end
end
return vecs
return λ, vecs
end
function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=nothing)
function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby)
if any(!isfinite, D.diag)
throw(ArgumentError("matrix contains Infs or NaNs"))
end
Td = _ortho_eltype(eltype(D))
λ = eigvals(D)
λ = eigvals(D; sortby=nothing)
if !isnothing(sortby)
p = sortperm(λ; alg=QuickSort, by=sortby)
λ = λ[p]
Expand All @@ -1051,20 +1051,19 @@ function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true, sortby::Union{
end
Eigen(λ, evecs)
end
function eigen(D::Diagonal{<:AbstractMatrix}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=nothing)
function eigen(D::Diagonal{<:AbstractMatrix}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby)
if any(any(!isfinite, x) for x in D.diag)
throw(ArgumentError("matrix contains Infs or NaNs"))
end
λ = eigvals(D)
evecs = eigvecs(D)
λ, evecs = _eigen(D)
if !isnothing(sortby)
p = sortperm(λ; alg=QuickSort, by=sortby)
λ = λ[p]
evecs = evecs[:,p]
end
Eigen(λ, evecs)
end
function eigen(Da::Diagonal, Db::Diagonal; sortby::Union{Function,Nothing}=nothing)
function eigen(Da::Diagonal, Db::Diagonal; sortby::Union{Function,Nothing}=eigsortby)
if any(!isfinite, Da.diag) || any(!isfinite, Db.diag)
throw(ArgumentError("matrices contain Infs or NaNs"))
end
Expand All @@ -1073,7 +1072,7 @@ function eigen(Da::Diagonal, Db::Diagonal; sortby::Union{Function,Nothing}=nothi
end
return GeneralizedEigen(eigen(Db \ Da; sortby)...)
end
function eigen(A::AbstractMatrix, D::Diagonal; sortby::Union{Function,Nothing}=nothing)
function eigen(A::AbstractMatrix, D::Diagonal; sortby::Union{Function,Nothing}=eigsortby)
if any(iszero, D.diag)
throw(ArgumentError("right-hand side diagonal matrix is singular"))
end
Expand Down
2 changes: 1 addition & 1 deletion src/eigen.jl
Original file line number Diff line number Diff line change
Expand Up @@ -199,7 +199,7 @@ 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`](@ref) or [`SymTridiagonal`](@ref)) may implement their own sorting convention and not
(e.g. [`SymTridiagonal`](@ref)) may implement their own sorting convention and not
accept a `sortby` keyword.

# Examples
Expand Down
6 changes: 3 additions & 3 deletions test/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -413,7 +413,7 @@ LinearAlgebra.istril(N::NotDiagonal) = istril(N.a)
@test factorize(D) == D

@testset "Eigensystem" begin
eigD = eigen(D)
eigD = eigen(D, sortby=nothing)
@test Diagonal(eigD.values) == D
@test eigD.vectors == Matrix(I, size(D))
eigsortD = eigen(D, sortby=LinearAlgebra.eigsortby)
Expand Down Expand Up @@ -564,7 +564,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 sort([eigvals(D)...;], by=LinearAlgebra.eigsortby) ≈ eigvals([D.diag[1] zeros(3,2); zeros(2,3) D.diag[2]])
@test eigvals(D, sortby=LinearAlgebra.eigsortby) ≈ eigvals([D.diag[1] zeros(3,2); zeros(2,3) D.diag[2]])
end

@testset "eigvals should return a copy of the diagonal" begin
Expand Down Expand Up @@ -1065,7 +1065,7 @@ 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(D) == eigen(D).values == [-1.3, 0.2, 0.4] # 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]
Expand Down
2 changes: 1 addition & 1 deletion test/unitful.jl
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ end
@test U isa AbstractMatrix{<:Union{Real,Complex}}
@test V isa AbstractMatrix{<:Union{Real,Complex}}
@test s isa AbstractVector{<:Furlong{1}}
E = eigen(Du)
E = eigen(Du; sortby=nothing)
vals, vecs = E
@test Matrix(E) == Du
@test vals isa AbstractVector{<:Furlong{1}}
Expand Down