Redirect for ishermitian==true and match Base sorting (#56)

Author: ericphanson

src/eigenGeneral.jl

@@ -247,13 +247,41 @@ function doubleShiftQR!(HH::StridedMatrix, τ::Rotation, shiftTrace::Number, shi
 end
 
 _eigvals!(A::StridedMatrix; kwargs...)           = _eigvals!(_schur!(A; kwargs...))
-_eigvals!(H::HessenbergMatrix; kwargs...)        = _eigvals!(_schur!(H, kwargs...))
-_eigvals!(H::HessenbergFactorization; kwargs...) = _eigvals!(_schur!(H, kwargs...))
+_eigvals!(H::HessenbergMatrix; kwargs...)        = _eigvals!(_schur!(H; kwargs...))
+_eigvals!(H::HessenbergFactorization; kwargs...) = _eigvals!(_schur!(H; kwargs...))
 
 # Overload methods from LinearAlgebra to make them work generically
-LinearAlgebra.eigvals!(A::StridedMatrix; kwargs...)           = _eigvals!(A; kwargs...)
-LinearAlgebra.eigvals!(H::HessenbergMatrix; kwargs...)        = _eigvals!(H, kwargs...)
-LinearAlgebra.eigvals!(H::HessenbergFactorization; kwargs...) = _eigvals!(H, kwargs...)
+if VERSION > v"1.2.0-DEV.0"
+    function LinearAlgebra.eigvals!(A::StridedMatrix;
+                                    sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+                                    kwargs...)
+            if ishermitian(A)
+                return LinearAlgebra.sorteig!(eigvals!(Hermitian(A)), sortby)
+            end
+            LinearAlgebra.sorteig!(_eigvals!(A; kwargs...), sortby)
+    end
+
+    LinearAlgebra.eigvals!(H::HessenbergMatrix;
+                        sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+                        kwargs...) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby)
+
+    LinearAlgebra.eigvals!(H::HessenbergFactorization;
+                        sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
+                        kwargs...) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby)
+else
+
+    function LinearAlgebra.eigvals!(A::StridedMatrix; kwargs...)
+        if ishermitian(A)
+            return eigvals!(Hermitian(A))
+        end
+        _eigvals!(A; kwargs...)
+    end
+
+    LinearAlgebra.eigvals!(H::HessenbergMatrix; kwargs...) = _eigvals!(H; kwargs...)
+
+    LinearAlgebra.eigvals!(H::HessenbergFactorization; kwargs...) = _eigvals!(H; kwargs...)
+end
+
 
 function _eigvals!(S::Schur{T}; tol = eps(real(T))) where T
     HH = S.data

test/eigengeneral.jl

@@ -14,6 +14,9 @@ cplxord = t -> (real(t), imag(t))
     vBig    = eigvals(big.(A)) # not defined in LinearAlgebra so will dispatch to the version in GenericLinearAlgebra
     @test sort(vGLA, by = cplxord) ≈ sort(vLAPACK, by = cplxord)
     @test sort(vGLA, by = cplxord) ≈ sort(complex(eltype(A)).(vBig), by = cplxord)
+    if VERSION > v"1.2.0-DEV.0"
+        @test issorted(vBig, by = cplxord)
+    end
 end
 
 @testset "make sure that solver doesn't hang" begin
@@ -54,4 +57,4 @@ end
     @test sum(eigvals(schur(A).Schur)) ≈ sum(eigvals(Float64.(schur(big.(A)).Schur)))
 end
 
-end
+end

test/eigenselfadjoint.jl

@@ -9,6 +9,7 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
         # LinearAlgebra
         T = SymTridiagonal(big.(randn(n)), big.(randn(n - 1)))
         vals, vecs  = eigen(T)
+        @test issorted(vals)
         @testset "default" begin
             @test (vecs'*T)*vecs ≈ Diagonal(vals)
             @test eigvals(T) ≈ vals
@@ -17,13 +18,15 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
 
         @testset "eigen2" begin
             vals2, vecs2 = GenericLinearAlgebra.eigen2(T)
+            @test issorted(vals2)
             @test vals ≈ vals2
             @test vecs[[1,n],:] == vecs2
             @test vecs2*vecs2' ≈ Matrix(I, 2, 2)
         end
 
         @testset "QR version (QL is default)" begin
             vals, vecs = GenericLinearAlgebra.eigQR!(copy(T), vectors = Matrix{eltype(T)}(I, n, n))
+            @test issorted(vals)
             @test (vecs'*T)*vecs ≈ Diagonal(vals)
             @test eigvals(T) ≈ vals
             @test vecs'vecs ≈ Matrix(I, n, n)
@@ -37,13 +40,15 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
             @test vecs'*A*vecs ≈ diagm(0 => vals)
             @test eigvals(A) ≈ vals
             @test vecs'vecs ≈ Matrix(I, n, n)
+            @test issorted(vals)
         end
 
         @testset "eigen2" begin
             vals2, vecs2 = GenericLinearAlgebra.eigen2(A)
             @test vals ≈ vals2
             @test vecs[[1,n],:] ≈ vecs2
             @test vecs2*vecs2'  ≈ Matrix(I, 2, 2)
+            @test issorted(vals2)
         end
     end
 
@@ -52,6 +57,8 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
         λ = 10 .^ range(-8, stop=0, length=n)
         A = Hermitian(V*Diagonal(λ)*V' |> t -> (t + t')/2, uplo)
         vals, vecs = eigen(A)
+        @test issorted(vals)
+
         @testset "default" begin
             if uplo == :L # FixMe! Probably an conjugation is off somewhere. Don't have time to check now.
                 @test_broken vecs'*A*vecs ≈ diagm(0=> vals)
@@ -65,6 +72,7 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
 
         @testset "eigen2" begin
             vals2, vecs2 = GenericLinearAlgebra.eigen2(A)
+            @test issorted(vals2)
             @test vals ≈ vals2
             @test vecs[[1,n],:] ≈ vecs2
             @test vecs2*vecs2'  ≈ Matrix(I, 2, 2)
@@ -102,4 +110,10 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
         @test eigen(M1).values == GenericLinearAlgebra._eigen!(copy(M1)).values
         @test eigen(M2).values == GenericLinearAlgebra._eigen!(copy(M2)).values
     end
+
+    @testset "Sorting of `ishermitian(T)==true` matrices on pre-1.2" begin
+        T = big.(randn(5, 5))
+        T = T + T'
+        @test issorted(eigvals(T))
+    end
 end