Fix #133 (Support sortby) (#134)

Author: dlfivefifty

.gitignore

@@ -1,3 +1,3 @@
 Manifest.toml
 docs/build/
-.github/.DS_Store
+.DS_Store

Project.toml

@@ -1,6 +1,6 @@
 name = "GenericLinearAlgebra"
 uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
-version = "0.3.11"
+version = "0.3.12"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/eigenSelfAdjoint.jl

@@ -162,7 +162,7 @@ function eig2x2!(d::StridedVector, e::StridedVector, j::Integer, vectors::Matrix
     return c, s
 end
 
-function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T)) where {T<:Real}
+function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) where {T<:Real}
     d = S.dv
     e = S.ev
     n = length(d)
@@ -216,7 +216,7 @@ function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T)) where {T<:Real}
             end
         end
     end
-    sort!(d)
+    LinearAlgebra.sorteig!(d, sortby)
 end
 
 function eigQL!(
@@ -570,39 +570,39 @@ function symtriUpper!(
 end
 
 
-_eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A)))) =
-    eigvalsPWK!(A.diagonals, tol = tol)
+_eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+    eigvalsPWK!(A.diagonals; tol, sortby)
 
-_eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A)))) = eigvalsPWK!(A, tol = tol)
+_eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = eigvalsPWK!(A; tol, sortby)
 
-_eigvals!(A::Hermitian; tol = eps(real(eltype(A)))) = eigvals!(symtri!(A), tol = tol)
+_eigvals!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = eigvals!(symtri!(A); tol, sortby)
 
 
-LinearAlgebra.eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A)))) =
-    _eigvals!(A, tol = tol)
+LinearAlgebra.eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+    _eigvals!(A; tol, sortby)
 
-LinearAlgebra.eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A)))) =
-    _eigvals!(A, tol = tol)
+LinearAlgebra.eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+    _eigvals!(A; tol, sortby)
 
-LinearAlgebra.eigvals!(A::Hermitian; tol = eps(real(eltype(A)))) = _eigvals!(A, tol = tol)
+LinearAlgebra.eigvals!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigvals!(A; tol, sortby)
 
 
-_eigen!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A)))) =
-    LinearAlgebra.Eigen(eigQL!(A.diagonals, vectors = Array(A.Q), tol = tol)...)
+_eigen!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+    LinearAlgebra.Eigen(LinearAlgebra.sorteig!(eigQL!(A.diagonals, vectors = Array(A.Q), tol = tol)..., sortby)...)
 
-_eigen!(A::SymTridiagonal; tol = eps(real(eltype(A)))) = LinearAlgebra.Eigen(
+_eigen!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = LinearAlgebra.Eigen(
     eigQL!(A, vectors = Matrix{eltype(A)}(I, size(A, 1), size(A, 1)), tol = tol)...,
 )
 
-_eigen!(A::Hermitian; tol = eps(real(eltype(A)))) = _eigen!(symtri!(A), tol = tol)
+_eigen!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(symtri!(A), tol = tol)
 
 
-LinearAlgebra.eigen!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A)))) =
-    _eigen!(A, tol = tol)
+LinearAlgebra.eigen!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+    _eigen!(A; tol, sortby)
 
-LinearAlgebra.eigen!(A::SymTridiagonal; tol = eps(real(eltype(A)))) = _eigen!(A, tol = tol)
+LinearAlgebra.eigen!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(A; tol, sortby)
 
-LinearAlgebra.eigen!(A::Hermitian; tol = eps(real(eltype(A)))) = _eigen!(A, tol = tol)
+LinearAlgebra.eigen!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(A; tol, sortby)
 
 
 function eigen2!(

test/eigenselfadjoint.jl

@@ -1,5 +1,4 @@
 using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
-Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
 
 @testset "The selfadjoint eigen problem" begin
     n = 50
@@ -157,4 +156,13 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
             @test hypot(c, s) ≈ 1
         end
     end
+
+    @testset "#133" begin
+        A = SymTridiagonal{BigFloat}(randn(5), randn(4))
+        T = Tridiagonal(A)
+        @test eigvals(A) == eigvals(T) == eigvals(A; sortby=LinearAlgebra.eigsortby) == eigvals(T; sortby=LinearAlgebra.eigsortby) ==  eigvals!(deepcopy(A); sortby=LinearAlgebra.eigsortby)
+        @test eigen(A).values == eigen(T).values == eigen(A; sortby=LinearAlgebra.eigsortby).values == eigen(T; sortby=LinearAlgebra.eigsortby).values
+        # compare abs to avoid sign issues
+        @test abs.(eigen(A).vectors) == abs.(eigen(T).vectors) == abs.(eigen(A; sortby=LinearAlgebra.eigsortby).vectors) == abs.(eigen(T; sortby=LinearAlgebra.eigsortby).vectors)
+    end
 end