Skip to content
Merged
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
52 changes: 22 additions & 30 deletions src/eigen.jl
Original file line number Diff line number Diff line change
Expand Up @@ -112,8 +112,8 @@ end


@inline function _eig(s::Size, A::StaticMatrix, permute, scale)
# Only cover the hermitian branch, for now ast least
# This also solves some type-stability issues such as arise in Base
# Only cover the hermitian branch, for now at least
# This also solves some type-stability issues that arise in Base
if ishermitian(A)
return _eig(s, Hermitian(A), permute, scale)
else
Expand All @@ -136,19 +136,7 @@ end
TA = eltype(A)

@inbounds if A.uplo == 'U'
if a[3] == 0 # A is diagonal
A11 = a[1]
A22 = a[4]
if A11 < A22
vals = SVector(A11, A22)
vecs = @SMatrix [convert(TA, 1) convert(TA, 0);
convert(TA, 0) convert(TA, 1)]
else # A22 <= A11
vals = SVector(A22, A11)
vecs = @SMatrix [convert(TA, 0) convert(TA, 1);
convert(TA, 1) convert(TA, 0)]
end
else # A is not diagonal
if !iszero(a[3]) # A is not diagonal
t_half = real(a[1] + a[4]) / 2
d = real(a[1] * a[4] - a[3]' * a[3]) # Should be real

Expand All @@ -168,22 +156,11 @@ end

vecs = @SMatrix [ v11 v21 ;
v12 v22 ]

return (vals, vecs)
end
return (vals, vecs)
else # A.uplo == 'L'
if a[2] == 0 # A is diagonal
A11 = a[1]
A22 = a[4]
if A11 < A22
vals = SVector(A11, A22)
vecs = @SMatrix [convert(TA, 1) convert(TA, 0);
convert(TA, 0) convert(TA, 1)]
else # A22 <= A11
vals = SVector(A22, A11)
vecs = @SMatrix [convert(TA, 0) convert(TA, 1);
convert(TA, 1) convert(TA, 0)]
end
else # A is not diagonal
if !iszero(a[2]) # A is not diagonal
t_half = real(a[1] + a[4]) / 2
d = real(a[1] * a[4] - a[2]' * a[2]) # Should be real

Expand All @@ -203,9 +180,24 @@ end

vecs = @SMatrix [ v11 v21 ;
v12 v22 ]

return (vals,vecs)
end
return (vals,vecs)
end

# A must be diagonal if we reached this point; treatment of uplo 'L' and 'U' is then identical
A11 = real(a[1])
A22 = real(a[4])
if A11 < A22
vals = SVector(A11, A22)
vecs = @SMatrix [convert(TA, 1) convert(TA, 0);
convert(TA, 0) convert(TA, 1)]
else # A22 <= A11
vals = SVector(A22, A11)
vecs = @SMatrix [convert(TA, 0) convert(TA, 1);
convert(TA, 1) convert(TA, 0)]
end
return (vals,vecs)
end

# A small part of the code in the following method was inspired by works of David
Expand Down
8 changes: 8 additions & 0 deletions test/eigen.jl
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,14 @@ using StaticArrays, Test, LinearAlgebra
@test eigvals(m) ≈ sort(m_d)
@test eigvals(Hermitian(m)) ≈ sort(m_d)

m_c = complex(m) # issue 614 (diagonal complex Hermitian)
(vals, vecs) = eigen(Hermitian(m_c))
@test vals::SVector ≈ sort(m_d)
(vals, vecs) = eigen(Hermitian(m_c, :L))
@test vals::SVector ≈ sort(m_d)
@test eigvals(m_c) ≈ sort(m_d)
@test eigvals(Hermitian(m_c)) ≈ sort(m_d)

# issue #523
for (i, j) in ((1, 2), (2, 1)), uplo in (:U, :L)
A = SMatrix{2,2,Float64}((i, 0, 0, j))
Expand Down