From a11835ae4a36aa94238271fcdf8d4cac1f12723d Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 10 Jan 2026 13:34:01 +0100 Subject: [PATCH 1/6] diagonal `eig` outputs sorted values --- src/implementations/eig.jl | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/implementations/eig.jl b/src/implementations/eig.jl index 5a0dd679..eca5d89a 100644 --- a/src/implementations/eig.jl +++ b/src/implementations/eig.jl @@ -129,10 +129,20 @@ end # Diagonal logic # -------------- -function eig_full!(A::Diagonal, (D, V)::Tuple{Diagonal, Diagonal}, alg::DiagonalAlgorithm) - check_input(eig_full!, A, (D, V), alg) - D === A || copy!(D, A) - one!(V) +function eig_full!(A::Diagonal, DV, alg::DiagonalAlgorithm) + check_input(eig_full!, A, DV, alg) + D, V = DV + diagA = diagview(A) + I = sortperm(diagA; by = real) + if D === A + permute!(diagA, I) + else + diagview(D) .= view(diagA, I) + end + zero!(V) + n = size(A, 1) + I .+= (0:(n - 1)) .* n + V[I] .= Ref(one(eltype(V))) return D, V end @@ -140,6 +150,7 @@ function eig_vals!(A::Diagonal, D::AbstractVector, alg::DiagonalAlgorithm) check_input(eig_vals!, A, D, alg) Ad = diagview(A) D === Ad || copy!(D, Ad) + sort!(D; by = real) return D end From b9add605269ff9c8409b2cfce2749472ea19f15a Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 12 Jan 2026 14:04:20 +0100 Subject: [PATCH 2/6] correct output types --- src/implementations/eig.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/implementations/eig.jl b/src/implementations/eig.jl index eca5d89a..a038e174 100644 --- a/src/implementations/eig.jl +++ b/src/implementations/eig.jl @@ -32,12 +32,12 @@ function check_input(::typeof(eig_full!), A::AbstractMatrix, DV, ::DiagonalAlgor m, n = size(A) @assert m == n && isdiag(A) D, V = DV - @assert D isa Diagonal && V isa Diagonal + @assert D isa Diagonal @check_size(D, (m, m)) @check_size(V, (m, m)) # Diagonal doesn't need to promote to complex scalartype since we know it is diagonalizable @check_scalar(D, A) - @check_scalar(V, A) + @check_scalar(V, A, real) return nothing end function check_input(::typeof(eig_vals!), A::AbstractMatrix, D, ::DiagonalAlgorithm) @@ -70,7 +70,7 @@ function initialize_output(::Union{typeof(eig_trunc!), typeof(eig_trunc_no_error end function initialize_output(::typeof(eig_full!), A::Diagonal, ::DiagonalAlgorithm) - return A, similar(A) + return A, similar(A, real(eltype(A)), size(A)) end function initialize_output(::typeof(eig_vals!), A::Diagonal, ::DiagonalAlgorithm) return diagview(A) From ac64d3e7ec232eeb5f91e142a5029064a2266cdc Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 12 Jan 2026 14:04:27 +0100 Subject: [PATCH 3/6] correct sorting order --- src/implementations/eig.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/implementations/eig.jl b/src/implementations/eig.jl index a038e174..d685f740 100644 --- a/src/implementations/eig.jl +++ b/src/implementations/eig.jl @@ -129,11 +129,12 @@ end # Diagonal logic # -------------- +eig_sortby(x::T) where {T <: Number} = T <: Complex ? (real(x), imag(X)) : x function eig_full!(A::Diagonal, DV, alg::DiagonalAlgorithm) check_input(eig_full!, A, DV, alg) D, V = DV diagA = diagview(A) - I = sortperm(diagA; by = real) + I = sortperm(diagA; by = eig_sortby) if D === A permute!(diagA, I) else @@ -150,7 +151,7 @@ function eig_vals!(A::Diagonal, D::AbstractVector, alg::DiagonalAlgorithm) check_input(eig_vals!, A, D, alg) Ad = diagview(A) D === Ad || copy!(D, Ad) - sort!(D; by = real) + sort!(D; by = eig_sortby) return D end From b8327b11ab059b64173c20649351fd9456c1b4ae Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 12 Jan 2026 15:09:57 +0100 Subject: [PATCH 4/6] Fix typo --- src/implementations/eig.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/implementations/eig.jl b/src/implementations/eig.jl index d685f740..450bbf7a 100644 --- a/src/implementations/eig.jl +++ b/src/implementations/eig.jl @@ -129,7 +129,7 @@ end # Diagonal logic # -------------- -eig_sortby(x::T) where {T <: Number} = T <: Complex ? (real(x), imag(X)) : x +eig_sortby(x::T) where {T <: Number} = T <: Complex ? (real(x), imag(x)) : x function eig_full!(A::Diagonal, DV, alg::DiagonalAlgorithm) check_input(eig_full!, A, DV, alg) D, V = DV From a29f73a3fe22791b9d6386a1942aa27e3a373cb3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 15 Jan 2026 10:05:24 +0100 Subject: [PATCH 5/6] standardize `eig` output types --- src/implementations/eig.jl | 22 ++++++++++++---------- test/testsuite/eig.jl | 4 ++-- 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/src/implementations/eig.jl b/src/implementations/eig.jl index 450bbf7a..9b785898 100644 --- a/src/implementations/eig.jl +++ b/src/implementations/eig.jl @@ -30,23 +30,21 @@ end function check_input(::typeof(eig_full!), A::AbstractMatrix, DV, ::DiagonalAlgorithm) m, n = size(A) - @assert m == n && isdiag(A) + ((m == n) && isdiag(A)) || throw(DimensionMismatch("diagonal input matrix expected")) D, V = DV - @assert D isa Diagonal + @assert D isa Diagonal && V isa AbstractMatrix @check_size(D, (m, m)) + @check_scalar(D, A, complex) @check_size(V, (m, m)) - # Diagonal doesn't need to promote to complex scalartype since we know it is diagonalizable - @check_scalar(D, A) - @check_scalar(V, A, real) + @check_scalar(V, A, complex) return nothing end function check_input(::typeof(eig_vals!), A::AbstractMatrix, D, ::DiagonalAlgorithm) m, n = size(A) - @assert m == n && isdiag(A) + ((m == n) && isdiag(A)) || throw(DimensionMismatch("diagonal input matrix expected")) @assert D isa AbstractVector @check_size(D, (n,)) - # Diagonal doesn't need to promote to complex scalartype since we know it is diagonalizable - @check_scalar(D, A) + @check_scalar(D, A, complex) return nothing end @@ -70,10 +68,14 @@ function initialize_output(::Union{typeof(eig_trunc!), typeof(eig_trunc_no_error end function initialize_output(::typeof(eig_full!), A::Diagonal, ::DiagonalAlgorithm) - return A, similar(A, real(eltype(A)), size(A)) + T = eltype(A) + Tc = complex(T) + D = T <: Complex ? A : Diagonal(similar(A, Tc, size(A, 1))) + return D, similar(A, Tc, size(A)) end function initialize_output(::typeof(eig_vals!), A::Diagonal, ::DiagonalAlgorithm) - return diagview(A) + T = eltype(A) + return T <: Complex ? diagview(A) : similar(A, complex(T), size(A, 1)) end # Implementation diff --git a/test/testsuite/eig.jl b/test/testsuite/eig.jl index 61ed1fc8..6ddfef5d 100644 --- a/test/testsuite/eig.jl +++ b/test/testsuite/eig.jl @@ -28,7 +28,7 @@ function test_eig_full( return @testset "eig_full! $summary_str" begin A = instantiate_matrix(T, sz) Ac = deepcopy(A) - Tc = isa(A, Diagonal) ? eltype(T) : complex(eltype(T)) + Tc = complex(eltype(T)) D, V = @testinferred eig_full(A) @test eltype(D) == eltype(V) == Tc @test A * V ≈ V * D @@ -51,7 +51,7 @@ function test_eig_full_algs( return @testset "eig_full! algorithm $alg $summary_str" for alg in algs A = instantiate_matrix(T, sz) Ac = deepcopy(A) - Tc = isa(A, Diagonal) ? eltype(T) : complex(eltype(T)) + Tc = complex(eltype(T)) D, V = @testinferred eig_full(A; alg) @test eltype(D) == eltype(V) == Tc @test A * V ≈ V * D From a98c6888df1bc365abaf26834cd2225b35da443b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 15 Jan 2026 11:40:45 +0100 Subject: [PATCH 6/6] update changelog --- docs/src/changelog.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 331b691f..2e295231 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -30,6 +30,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Fixed +- Eigenvalue decompositions of diagonal inputs are sorted and have the same type as non-diagonal inputs ([#151](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/151) + ## [0.6.2](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/compare/v0.6.1...v0.6.2) - 2026-01-08 ### Added