From ec8c7a3ae060dd935ca46de8631d48a5c12ed346 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 20 Apr 2025 22:36:31 +0200 Subject: [PATCH] Rely on LinearAlgebra.eigsort! for sorting eigenvalues of the symmetric solver Some profiling revealed that the sorting used until this PR had significant overhead for smaller and medium sized problems. In addition, the values would be sorted twice if the user supplied a non-standard sorting function. --- src/eigenSelfAdjoint.jl | 92 ++++++++++++++++++++++++++-------------- test/eigenselfadjoint.jl | 3 +- 2 files changed, 62 insertions(+), 33 deletions(-) diff --git a/src/eigenSelfAdjoint.jl b/src/eigenSelfAdjoint.jl index 402da7d..c98a90e 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), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) where {T<:Real} +function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T)) where {T<:Real} d = S.dv e = S.ev n = length(d) @@ -217,9 +217,7 @@ function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), sortby::Union{Function, end end - # LinearAlgebra.eigvals will pass sortby=nothing but LAPACK always sort the symmetric - # eigenvalue problem so we'll will do the same here - LinearAlgebra.sorteig!(d, sortby === nothing ? LinearAlgebra.eigsortby : sortby) + return d end function eigQL!( @@ -287,8 +285,8 @@ function eigQL!( end end end - p = sortperm(d) - return d[p], vectors[:, p] + + return d, vectors end function eigQR!( @@ -339,8 +337,8 @@ function eigQR!( end end end - p = sortperm(d) - return d[p], vectors[:, p] + + return d, vectors end # Own implementation based on Parlett's book @@ -575,15 +573,17 @@ function symtriUpper!( ) end +_eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + LinearAlgebra.sorteig!(eigvalsPWK!(A; tol), sortby == nothing ? LinearAlgebra.eigsortby : sortby) _eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = - eigvalsPWK!(A.diagonals; tol, sortby) + _eigvals!(A.diagonals; tol, sortby) -_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))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + _eigvals!(symtri!(A); tol, sortby) -_eigvals!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = eigvals!(symtri!(A); tol, sortby) - -_eigvals!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = eigvals!(symtri!(A); tol, sortby) +_eigvals!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + _eigvals!(symtri!(A); tol, sortby) LinearAlgebra.eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigvals!(A; tol, sortby) @@ -591,29 +591,53 @@ LinearAlgebra.eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(elty 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))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigvals!(A; tol, sortby) +LinearAlgebra.eigvals!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + _eigvals!(A; tol, sortby) -LinearAlgebra.eigvals!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigvals!(A; tol, sortby) +LinearAlgebra.eigvals!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + _eigvals!(A; tol, sortby) _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)...) + LinearAlgebra.Eigen( + LinearAlgebra.sorteig!( + eigQL!( + A.diagonals; + vectors = Array(A.Q), + tol + )..., + sortby == nothing ? LinearAlgebra.eigsortby : sortby + )... + ) -_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::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + LinearAlgebra.Eigen( + LinearAlgebra.sorteig!( + eigQL!( + A; + vectors = Matrix{eltype(A)}(I, size(A, 1), size(A, 1)), + tol + )..., + sortby == nothing ? LinearAlgebra.eigsortby : sortby + )... + ) -_eigen!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(symtri!(A), tol = tol) +_eigen!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + _eigen!(symtri!(A); tol, sortby) -_eigen!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(symtri!(A), tol = tol) +_eigen!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + _eigen!(symtri!(A); tol, sortby) 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))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(A; tol, sortby) +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))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(A; tol, sortby) +LinearAlgebra.eigen!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + _eigen!(A; tol, sortby) -LinearAlgebra.eigen!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(A; tol, sortby) +LinearAlgebra.eigen!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = + _eigen!(A; tol, sortby) function eigen2!( @@ -623,29 +647,35 @@ function eigen2!( V = zeros(eltype(A), 2, size(A, 1)) V[1] = 1 V[end] = 1 - eigQL!(A.diagonals, vectors = rmul!(V, A.Q), tol = tol) + LinearAlgebra.sorteig!( + eigQL!(A.diagonals; vectors = rmul!(V, A.Q), tol)..., + LinearAlgebra.eigsortby + ) end function eigen2!(A::SymTridiagonal; tol = eps(real(float(one(eltype(A)))))) V = zeros(eltype(A), 2, size(A, 1)) V[1] = 1 V[end] = 1 - eigQL!(A, vectors = V, tol = tol) + LinearAlgebra.sorteig!( + eigQL!(A; vectors = V, tol)..., + LinearAlgebra.eigsortby + ) end eigen2!(A::Hermitian; tol = eps(float(real(one(eltype(A)))))) = - eigen2!(symtri!(A), tol = tol) + eigen2!(symtri!(A); tol) eigen2!(A::Symmetric{<:Real}; tol = eps(float(one(eltype(A))))) = - eigen2!(symtri!(A), tol = tol) + eigen2!(symtri!(A); tol) eigen2(A::SymTridiagonal; tol = eps(float(real(one(eltype(A)))))) = - eigen2!(copy(A), tol = tol) + eigen2!(copy(A); tol) -eigen2(A::Hermitian, tol = eps(float(real(one(eltype(A)))))) = eigen2!(copy(A), tol = tol) +eigen2(A::Hermitian; tol = eps(float(real(one(eltype(A)))))) = eigen2!(copy(A); tol) -eigen2(A::Symmetric{<:Real}, tol = eps(float(one(eltype(A))))) = eigen2!(copy(A), tol = tol) +eigen2(A::Symmetric{<:Real}; tol = eps(float(one(eltype(A))))) = eigen2!(copy(A); tol) # First method of each type here is identical to the method defined in # LinearAlgebra but is needed for disambiguation diff --git a/test/eigenselfadjoint.jl b/test/eigenselfadjoint.jl index 3fe8eb7..4c515fb 100644 --- a/test/eigenselfadjoint.jl +++ b/test/eigenselfadjoint.jl @@ -26,9 +26,8 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions @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 eigvals(T) ≈ sort(vals) @test vecs'vecs ≈ Matrix(I, n, n) end end