Fix bug in symtriLower/Upper! when diagonal contains non-real elements (#39)

Define a method to disambiguiate eigen(Hermitian)

Fixes #35

Author: andreasnoack

src/eigenSelfAdjoint.jl

@@ -438,7 +438,13 @@ function symtriLower!(AS::StridedMatrix{T},
                       u = Vector{T}(undef, size(AS, 1))) where T
     n = size(AS, 1)
 
+    # We ignore any non-real components of the diagonal
+    @inbounds for i in 1:n
+        AS[i,i] = real(AS[i,i])
+    end
+
     @inbounds for k = 1:(n - 2 + !(T<:Real))
+
         τk = LinearAlgebra.reflector!(view(AS, (k + 1):n, k))
         τ[k] = τk
 
@@ -492,6 +498,11 @@ function symtriUpper!(AS::StridedMatrix{T},
                       u = Vector{T}(undef, size(AS, 1))) where T
     n = LinearAlgebra.checksquare(AS)
 
+    # We ignore any non-real components of the diagonal
+    @inbounds for i in 1:n
+        AS[i,i] = real(AS[i,i])
+    end
+
     @inbounds for k = 1:(n - 2 + !(T<:Real))
         # This is a bit convoluted method to get the conjugated vector but conjugation is required for correctness of arrays of quaternions. Eventually, it should be sufficient to write vec(x') but it currently (July 10, 2018) hits a bug in LinearAlgebra
         τk = LinearAlgebra.reflector!(vec(transpose(view(AS, k, (k + 1):n)')))
@@ -618,6 +629,10 @@ function LinearAlgebra.eigvals(A::Hermitian{<:Real})
     T = typeof(sqrt(zero(eltype(A))))
     return eigvals!(LinearAlgebra.copy_oftype(A, T))
 end
+function LinearAlgebra.eigvals(A::Hermitian{<:Complex})
+    T = typeof(sqrt(zero(eltype(A))))
+    return eigvals!(LinearAlgebra.copy_oftype(A, T))
+end
 function LinearAlgebra.eigvals(A::Hermitian)
     T = typeof(sqrt(zero(eltype(A))))
     return eigvals!(LinearAlgebra.copy_oftype(A, T))
@@ -626,6 +641,10 @@ function LinearAlgebra.eigen(A::Hermitian{<:Real})
     T = typeof(sqrt(zero(eltype(A))))
     return eigen!(LinearAlgebra.copy_oftype(A, T))
 end
+function LinearAlgebra.eigen(A::Hermitian{<:Complex})
+    T = typeof(sqrt(zero(eltype(A))))
+    return eigen!(LinearAlgebra.copy_oftype(A, T))
+end
 function LinearAlgebra.eigen(A::Hermitian)
     T = typeof(sqrt(zero(eltype(A))))
     return eigen!(LinearAlgebra.copy_oftype(A, T))

test/eigenselfadjoint.jl

@@ -71,6 +71,13 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
         end
     end
 
+    @testset "big Hermitian{<:Complex}" begin
+        # This one used to cause an ambiguity error. See #35
+        A = complex.(randn(4,4), randn(4,4))
+        @test Float64.(eigen(Hermitian(big.(A))).values) ≈ eigen(Hermitian(copy(A))).values
+        @test Float64.(eigvals(Hermitian(big.(A))))      ≈ eigvals(Hermitian(copy(A)))
+    end
+
     @testset "generic Givens" begin
         x, y = randn(2)
         c, s, r = invoke(LinearAlgebra.givensAlgorithm, Tuple{Real,Real}, x, y)