From 81eee0d5e9c00f50a0fe79c761d99ef7373dff2c Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 10 Aug 2024 09:34:47 +0100 Subject: [PATCH] Fix #133 (Support sortby) (#134) --- .gitignore | 2 +- Project.toml | 2 +- src/eigenSelfAdjoint.jl | 38 +++++++++++++++++++------------------- test/eigenselfadjoint.jl | 10 +++++++++- 4 files changed, 30 insertions(+), 22 deletions(-) diff --git a/.gitignore b/.gitignore index c10b7b2..1699e29 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,3 @@ Manifest.toml docs/build/ -.github/.DS_Store +.DS_Store diff --git a/Project.toml b/Project.toml index 6a42cf0..1efa285 100644 --- a/Project.toml +++ b/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" diff --git a/src/eigenSelfAdjoint.jl b/src/eigenSelfAdjoint.jl index c0b6a2a..d8008d3 100644 --- a/src/eigenSelfAdjoint.jl +++ b/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!( diff --git a/test/eigenselfadjoint.jl b/test/eigenselfadjoint.jl index a5a7508..8f097d9 100644 --- a/test/eigenselfadjoint.jl +++ b/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