From 1e3cb9c085f8d79efb361b66c7f6d1b545a7d179 Mon Sep 17 00:00:00 2001 From: araujoms Date: Fri, 17 Oct 2025 19:19:34 +0200 Subject: [PATCH 1/2] make eigvals(::Diagonal) accept sorteig --- src/diagonal.jl | 31 +++++++++++++++---------------- src/eigen.jl | 2 +- test/diagonal.jl | 6 +++--- 3 files changed, 19 insertions(+), 20 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index e60fb009..797d3f81 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -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) @@ -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] @@ -1051,12 +1051,11 @@ 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] @@ -1064,7 +1063,7 @@ function eigen(D::Diagonal{<:AbstractMatrix}; permute::Bool=true, scale::Bool=tr 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 @@ -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 diff --git a/src/eigen.jl b/src/eigen.jl index f565a547..57fefa96 100644 --- a/src/eigen.jl +++ b/src/eigen.jl @@ -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 diff --git a/test/diagonal.jl b/test/diagonal.jl index 712f426c..753a8cf9 100644 --- a/test/diagonal.jl +++ b/test/diagonal.jl @@ -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) @@ -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 @@ -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] From d1f41eb59e97c6103af2ea5b659570aad9067435 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mateus=20Ara=C3=BAjo?= Date: Fri, 17 Oct 2025 22:59:27 +0200 Subject: [PATCH 2/2] fix Furlong test --- test/unitful.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unitful.jl b/test/unitful.jl index 95feed97..87480cd7 100644 --- a/test/unitful.jl +++ b/test/unitful.jl @@ -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}}