Skip to content

Commit 2ce3f73

Browse files
jishnubViralBShah
authored andcommitted
Move alg to a keyword argument in symmetric eigen (#1214)
This reopens JuliaLang/julia#55481 One argument there was that packages may wish to define custom eigenvalue algorithms. I'm unsure if this is common, but in principle, we may provide an internal method `eigen(A, alg)` that packages extend, and users may call the method with `alg` as a keyword argument. Fixes #1208 --------- Co-authored-by: Viral B. Shah <viral@juliacomputing.com>
1 parent 4f203c3 commit 2ce3f73

File tree

3 files changed

+35
-29
lines changed

3 files changed

+35
-29
lines changed

src/symmetriceigen.jl

Lines changed: 19 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -6,10 +6,16 @@ eigencopy_oftype(A::Hermitian, S) = Hermitian(copytrito!(similar(parent(A), S, s
66
eigencopy_oftype(A::Symmetric, S) = Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo))
77
eigencopy_oftype(A::Symmetric{<:Complex}, S) = copyto!(similar(parent(A), S), A)
88

9-
default_eigen_alg(A) = DivideAndConquer()
9+
"""
10+
default_eigen_alg(A)
11+
12+
Return the default algorithm used to solve the eigensystem `A v = λ v` for a symmetric matrix `A`.
13+
Defaults to `LinearAlegbra.DivideAndConquer()`, which corresponds to the LAPACK function `LAPACK.syevd!`.
14+
"""
15+
default_eigen_alg(@nospecialize(A)) = DivideAndConquer()
1016

1117
# Eigensolvers for symmetric and Hermitian matrices
12-
function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
18+
function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
1319
if alg === DivideAndConquer()
1420
Eigen(sorteig!(LAPACK.syevd!('V', A.uplo, A.data)..., sortby)...)
1521
elseif alg === QRIteration()
@@ -22,7 +28,7 @@ function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algo
2228
end
2329

2430
"""
25-
eigen(A::Union{Hermitian, Symmetric}, alg::Algorithm = default_eigen_alg(A)) -> Eigen
31+
eigen(A::Union{Hermitian, Symmetric}; alg::LinearAlgebra.Algorithm = LinearAlgebra.default_eigen_alg(A)) -> Eigen
2632
2733
Compute the eigenvalue decomposition of `A`, returning an [`Eigen`](@ref) factorization object `F`
2834
which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the
@@ -45,19 +51,19 @@ The default `alg` used may change in the future.
4551
4652
The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).
4753
"""
48-
function eigen(A::RealHermSymComplexHerm, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
49-
_eigen(A, alg; sortby)
54+
function eigen(A::RealHermSymComplexHerm; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
55+
_eigen(A; alg, sortby)
5056
end
5157

5258
# we dispatch on the eltype in an internal method to avoid ambiguities
53-
function _eigen(A::RealHermSymComplexHerm, alg::Algorithm; sortby)
59+
function _eigen(A::RealHermSymComplexHerm; alg::Algorithm, sortby)
5460
S = eigtype(eltype(A))
55-
eigen!(eigencopy_oftype(A, S), alg; sortby)
61+
eigen!(eigencopy_oftype(A, S); alg, sortby)
5662
end
5763

58-
function _eigen(A::RealHermSymComplexHerm{Float16}, alg::Algorithm; sortby::Union{Function,Nothing}=nothing)
64+
function _eigen(A::RealHermSymComplexHerm{Float16}; alg::Algorithm, sortby::Union{Function,Nothing}=nothing)
5965
S = eigtype(eltype(A))
60-
E = eigen!(eigencopy_oftype(A, S), alg, sortby=sortby)
66+
E = eigen!(eigencopy_oftype(A, S); alg, sortby)
6167
values = convert(AbstractVector{Float16}, E.values)
6268
vectors = convert(AbstractMatrix{isreal(E.vectors) ? Float16 : Complex{Float16}}, E.vectors)
6369
return Eigen(values, vectors)
@@ -114,7 +120,7 @@ function eigen(A::RealHermSymComplexHerm, vl::Real, vh::Real)
114120
end
115121

116122

117-
function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
123+
function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
118124
vals::Vector{real(eltype(A))} = if alg === DivideAndConquer()
119125
LAPACK.syevd!('N', A.uplo, A.data)
120126
elseif alg === QRIteration()
@@ -129,7 +135,7 @@ function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Al
129135
end
130136

131137
"""
132-
eigvals(A::Union{Hermitian, Symmetric}, alg::Algorithm = default_eigen_alg(A))) -> values
138+
eigvals(A::Union{Hermitian, Symmetric}; alg::Algorithm = default_eigen_alg(A))) -> values
133139
134140
Return the eigenvalues of `A`.
135141
@@ -143,9 +149,9 @@ a comparison of the accuracy and performance of different methods.
143149
144150
The default `alg` used may change in the future.
145151
"""
146-
function eigvals(A::RealHermSymComplexHerm, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
152+
function eigvals(A::RealHermSymComplexHerm; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
147153
S = eigtype(eltype(A))
148-
eigvals!(eigencopy_oftype(A, S), alg; sortby)
154+
eigvals!(eigencopy_oftype(A, S); alg, sortby)
149155
end
150156

151157

test/symmetric.jl

Lines changed: 0 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -720,21 +720,6 @@ end
720720
end
721721
end
722722

723-
@testset "eigendecomposition Algorithms" begin
724-
using LinearAlgebra: DivideAndConquer, QRIteration, RobustRepresentations
725-
for T in (Float64, ComplexF64, Float32, ComplexF32)
726-
n = 4
727-
A = T <: Real ? Symmetric(randn(T, n, n)) : Hermitian(randn(T, n, n))
728-
d, v = eigen(A)
729-
for alg in (DivideAndConquer(), QRIteration(), RobustRepresentations())
730-
@test (@inferred eigvals(A, alg)) d
731-
d2, v2 = @inferred eigen(A, alg)
732-
@test d2 d
733-
@test A * v2 v2 * Diagonal(d2)
734-
end
735-
end
736-
end
737-
738723
const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
739724
isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl"))
740725
using .Main.ImmutableArrays

test/symmetriceigen.jl

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
module TestSymmetricEigen
44

55
using Test, LinearAlgebra
6+
using LinearAlgebra: DivideAndConquer, QRIteration, RobustRepresentations
67

78
@testset "chol-eigen-eigvals" begin
89
## Cholesky decomposition based
@@ -173,7 +174,7 @@ end
173174
@test D.vectors D32.vectors
174175

175176
# ensure that different algorithms dispatch correctly
176-
λ, V = eigen(C, LinearAlgebra.QRIteration())
177+
λ, V = eigen(C; alg=LinearAlgebra.QRIteration())
177178
@test λ isa Vector{Float16}
178179
@test C * V V * Diagonal(λ)
179180
end
@@ -184,4 +185,18 @@ end
184185
@test S * v v * Diagonal(λ)
185186
end
186187

188+
@testset "eigendecomposition Algorithms" begin
189+
for T in (Float64, ComplexF64, Float32, ComplexF32)
190+
n = 4
191+
A = T <: Real ? Symmetric(randn(T, n, n)) : Hermitian(randn(T, n, n))
192+
d, v = eigen(A)
193+
for alg in (DivideAndConquer(), QRIteration(), RobustRepresentations())
194+
@test (@inferred eigvals(A; alg)) d
195+
d2, v2 = @inferred eigen(A; alg)
196+
@test d2 d
197+
@test A * v2 v2 * Diagonal(d2)
198+
end
199+
end
200+
end
201+
187202
end # module TestSymmetricEigen

0 commit comments

Comments
 (0)