Skip to content

Commit 9e5cb8b

Browse files
andreasnoackandreasnoack
authored
Rely on LinearAlgebra.eigsort! for sorting eigenvalues of the (#151)
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.
1 parent bebae94 commit 9e5cb8b

File tree

2 files changed

+62
-33
lines changed

2 files changed

+62
-33
lines changed

‎src/eigenSelfAdjoint.jl‎

Lines changed: 61 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -162,7 +162,7 @@ function eig2x2!(d::StridedVector, e::StridedVector, j::Integer, vectors::Matrix
162162
return c, s
163163
end
164164

165-
function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) where {T<:Real}
165+
function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T)) where {T<:Real}
166166
d = S.dv
167167
e = S.ev
168168
n = length(d)
@@ -217,9 +217,7 @@ function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), sortby::Union{Function,
217217
end
218218
end
219219

220-
# LinearAlgebra.eigvals will pass sortby=nothing but LAPACK always sort the symmetric
221-
# eigenvalue problem so we'll will do the same here
222-
LinearAlgebra.sorteig!(d, sortby === nothing ? LinearAlgebra.eigsortby : sortby)
220+
return d
223221
end
224222

225223
function eigQL!(
@@ -287,8 +285,8 @@ function eigQL!(
287285
end
288286
end
289287
end
290-
p = sortperm(d)
291-
return d[p], vectors[:, p]
288+
289+
return d, vectors
292290
end
293291

294292
function eigQR!(
@@ -339,8 +337,8 @@ function eigQR!(
339337
end
340338
end
341339
end
342-
p = sortperm(d)
343-
return d[p], vectors[:, p]
340+
341+
return d, vectors
344342
end
345343

346344
# Own implementation based on Parlett's book
@@ -575,45 +573,71 @@ function symtriUpper!(
575573
)
576574
end
577575

576+
_eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
577+
LinearAlgebra.sorteig!(eigvalsPWK!(A; tol), sortby == nothing ? LinearAlgebra.eigsortby : sortby)
578578

579579
_eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
580-
eigvalsPWK!(A.diagonals; tol, sortby)
580+
_eigvals!(A.diagonals; tol, sortby)
581581

582-
_eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = eigvalsPWK!(A; tol, sortby)
582+
_eigvals!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
583+
_eigvals!(symtri!(A); tol, sortby)
583584

584-
_eigvals!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = eigvals!(symtri!(A); tol, sortby)
585-
586-
_eigvals!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = eigvals!(symtri!(A); tol, sortby)
585+
_eigvals!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
586+
_eigvals!(symtri!(A); tol, sortby)
587587

588588
LinearAlgebra.eigvals!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
589589
_eigvals!(A; tol, sortby)
590590

591591
LinearAlgebra.eigvals!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
592592
_eigvals!(A; tol, sortby)
593593

594-
LinearAlgebra.eigvals!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigvals!(A; tol, sortby)
594+
LinearAlgebra.eigvals!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
595+
_eigvals!(A; tol, sortby)
595596

596-
LinearAlgebra.eigvals!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigvals!(A; tol, sortby)
597+
LinearAlgebra.eigvals!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
598+
_eigvals!(A; tol, sortby)
597599

598600
_eigen!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
599-
LinearAlgebra.Eigen(LinearAlgebra.sorteig!(eigQL!(A.diagonals, vectors = Array(A.Q), tol = tol)..., sortby)...)
601+
LinearAlgebra.Eigen(
602+
LinearAlgebra.sorteig!(
603+
eigQL!(
604+
A.diagonals;
605+
vectors = Array(A.Q),
606+
tol
607+
)...,
608+
sortby == nothing ? LinearAlgebra.eigsortby : sortby
609+
)...
610+
)
600611

601-
_eigen!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = LinearAlgebra.Eigen(
602-
eigQL!(A, vectors = Matrix{eltype(A)}(I, size(A, 1), size(A, 1)), tol = tol)...,
603-
)
612+
_eigen!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
613+
LinearAlgebra.Eigen(
614+
LinearAlgebra.sorteig!(
615+
eigQL!(
616+
A;
617+
vectors = Matrix{eltype(A)}(I, size(A, 1), size(A, 1)),
618+
tol
619+
)...,
620+
sortby == nothing ? LinearAlgebra.eigsortby : sortby
621+
)...
622+
)
604623

605-
_eigen!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(symtri!(A), tol = tol)
624+
_eigen!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
625+
_eigen!(symtri!(A); tol, sortby)
606626

607-
_eigen!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(symtri!(A), tol = tol)
627+
_eigen!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
628+
_eigen!(symtri!(A); tol, sortby)
608629

609630
LinearAlgebra.eigen!(A::SymmetricTridiagonalFactorization; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
610631
_eigen!(A; tol, sortby)
611632

612-
LinearAlgebra.eigen!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(A; tol, sortby)
633+
LinearAlgebra.eigen!(A::SymTridiagonal; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
634+
_eigen!(A; tol, sortby)
613635

614-
LinearAlgebra.eigen!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(A; tol, sortby)
636+
LinearAlgebra.eigen!(A::Hermitian; tol = eps(real(eltype(A))), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
637+
_eigen!(A; tol, sortby)
615638

616-
LinearAlgebra.eigen!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) = _eigen!(A; tol, sortby)
639+
LinearAlgebra.eigen!(A::Symmetric{<:Real}; tol = eps(eltype(A)), sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
640+
_eigen!(A; tol, sortby)
617641

618642

619643
function eigen2!(
@@ -623,29 +647,35 @@ function eigen2!(
623647
V = zeros(eltype(A), 2, size(A, 1))
624648
V[1] = 1
625649
V[end] = 1
626-
eigQL!(A.diagonals, vectors = rmul!(V, A.Q), tol = tol)
650+
LinearAlgebra.sorteig!(
651+
eigQL!(A.diagonals; vectors = rmul!(V, A.Q), tol)...,
652+
LinearAlgebra.eigsortby
653+
)
627654
end
628655

629656
function eigen2!(A::SymTridiagonal; tol = eps(real(float(one(eltype(A))))))
630657
V = zeros(eltype(A), 2, size(A, 1))
631658
V[1] = 1
632659
V[end] = 1
633-
eigQL!(A, vectors = V, tol = tol)
660+
LinearAlgebra.sorteig!(
661+
eigQL!(A; vectors = V, tol)...,
662+
LinearAlgebra.eigsortby
663+
)
634664
end
635665

636666
eigen2!(A::Hermitian; tol = eps(float(real(one(eltype(A)))))) =
637-
eigen2!(symtri!(A), tol = tol)
667+
eigen2!(symtri!(A); tol)
638668

639669
eigen2!(A::Symmetric{<:Real}; tol = eps(float(one(eltype(A))))) =
640-
eigen2!(symtri!(A), tol = tol)
670+
eigen2!(symtri!(A); tol)
641671

642672

643673
eigen2(A::SymTridiagonal; tol = eps(float(real(one(eltype(A)))))) =
644-
eigen2!(copy(A), tol = tol)
674+
eigen2!(copy(A); tol)
645675

646-
eigen2(A::Hermitian, tol = eps(float(real(one(eltype(A)))))) = eigen2!(copy(A), tol = tol)
676+
eigen2(A::Hermitian; tol = eps(float(real(one(eltype(A)))))) = eigen2!(copy(A); tol)
647677

648-
eigen2(A::Symmetric{<:Real}, tol = eps(float(one(eltype(A))))) = eigen2!(copy(A), tol = tol)
678+
eigen2(A::Symmetric{<:Real}; tol = eps(float(one(eltype(A))))) = eigen2!(copy(A); tol)
649679

650680
# First method of each type here is identical to the method defined in
651681
# LinearAlgebra but is needed for disambiguation

‎test/eigenselfadjoint.jl‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -26,9 +26,8 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
2626
@testset "QR version (QL is default)" begin
2727
vals, vecs =
2828
GenericLinearAlgebra.eigQR!(copy(T), vectors = Matrix{eltype(T)}(I, n, n))
29-
@test issorted(vals)
3029
@test (vecs' * T) * vecs ≈ Diagonal(vals)
31-
@test eigvals(T) ≈ vals
30+
@test eigvals(T) ≈ sort(vals)
3231
@test vecs'vecs ≈ Matrix(I, n, n)
3332
end
3433
end

0 commit comments

Comments
 (0)