diff --git a/stdlib/LinearAlgebra/src/symmetriceigen.jl b/stdlib/LinearAlgebra/src/symmetriceigen.jl index fee524a702187..860a09a2c810e 100644 --- a/stdlib/LinearAlgebra/src/symmetriceigen.jl +++ b/stdlib/LinearAlgebra/src/symmetriceigen.jl @@ -6,10 +6,16 @@ eigencopy_oftype(A::Hermitian, S) = Hermitian(copytrito!(similar(parent(A), S, s eigencopy_oftype(A::Symmetric, S) = Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo)) eigencopy_oftype(A::Symmetric{<:Complex}, S) = copyto!(similar(parent(A), S), A) -default_eigen_alg(A) = DivideAndConquer() +""" + default_eigen_alg(A) + +Return the default algorithm used to solve the eigensystem `A v = λ v` for a symmetric matrix `A`. +Defaults to `LinearAlegbra.DivideAndConquer()`, which corresponds to the LAPACK function `LAPACK.syevd!`. +""" +default_eigen_alg(@nospecialize(A)) = DivideAndConquer() # Eigensolvers for symmetric and Hermitian matrices -function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing) +function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing, alg::Algorithm = default_eigen_alg(A)) if alg === DivideAndConquer() Eigen(sorteig!(LAPACK.syevd!('V', A.uplo, A.data)..., sortby)...) elseif alg === QRIteration() @@ -29,7 +35,7 @@ function eigen(A::RealHermSymComplexHerm{Float16}; sortby::Union{Function,Nothin end """ - eigen(A::Union{Hermitian, Symmetric}, alg::Algorithm = default_eigen_alg(A)) -> Eigen + eigen(A::Union{Hermitian, Symmetric}; alg::LinearAlgebra.Algorithm = LinearAlgebra.default_eigen_alg(A)) -> Eigen Compute the eigenvalue decomposition of `A`, returning an [`Eigen`](@ref) factorization object `F` which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the @@ -52,9 +58,9 @@ The default `alg` used may change in the future. The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref). """ -function eigen(A::RealHermSymComplexHerm, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing) +function eigen(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing, alg::Algorithm = default_eigen_alg(A)) S = eigtype(eltype(A)) - eigen!(eigencopy_oftype(A, S), alg; sortby) + eigen!(eigencopy_oftype(A, S); sortby, alg) end @@ -109,7 +115,7 @@ function eigen(A::RealHermSymComplexHerm, vl::Real, vh::Real) end -function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing) +function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing, alg::Algorithm = default_eigen_alg(A)) vals::Vector{real(eltype(A))} = if alg === DivideAndConquer() LAPACK.syevd!('N', A.uplo, A.data) elseif alg === QRIteration() @@ -124,7 +130,7 @@ function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Al end """ - eigvals(A::Union{Hermitian, Symmetric}, alg::Algorithm = default_eigen_alg(A))) -> values + eigvals(A::Union{Hermitian, Symmetric}; alg::LinearAlgebra.Algorithm = LinearAlgebra.default_eigen_alg(A))) -> values Return the eigenvalues of `A`. @@ -138,9 +144,9 @@ a comparison of the accuracy and performance of different methods. The default `alg` used may change in the future. """ -function eigvals(A::RealHermSymComplexHerm, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing) +function eigvals(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing, alg::Algorithm = default_eigen_alg(A)) S = eigtype(eltype(A)) - eigvals!(eigencopy_oftype(A, S), alg; sortby) + eigvals!(eigencopy_oftype(A, S); sortby, alg) end diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 89e9ca0d6a51d..f10540cf8e44c 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -720,21 +720,6 @@ end end end -@testset "eigendecomposition Algorithms" begin - using LinearAlgebra: DivideAndConquer, QRIteration, RobustRepresentations - for T in (Float64, ComplexF64, Float32, ComplexF32) - n = 4 - A = T <: Real ? Symmetric(randn(T, n, n)) : Hermitian(randn(T, n, n)) - d, v = eigen(A) - for alg in (DivideAndConquer(), QRIteration(), RobustRepresentations()) - @test (@inferred eigvals(A, alg)) ≈ d - d2, v2 = @inferred eigen(A, alg) - @test d2 ≈ d - @test A * v2 ≈ v2 * Diagonal(d2) - end - end -end - const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl")) using .Main.ImmutableArrays diff --git a/stdlib/LinearAlgebra/test/symmetriceigen.jl b/stdlib/LinearAlgebra/test/symmetriceigen.jl index d55d1deb6bf33..685a9789f1f99 100644 --- a/stdlib/LinearAlgebra/test/symmetriceigen.jl +++ b/stdlib/LinearAlgebra/test/symmetriceigen.jl @@ -3,6 +3,7 @@ module TestSymmetricEigen using Test, LinearAlgebra +using LinearAlgebra: DivideAndConquer, QRIteration, RobustRepresentations @testset "chol-eigen-eigvals" begin ## Cholesky decomposition based @@ -179,4 +180,18 @@ end @test S * v ≈ v * Diagonal(λ) end +@testset "eigendecomposition Algorithms" begin + for T in (Float64, ComplexF64, Float32, ComplexF32) + n = 4 + A = T <: Real ? Symmetric(randn(T, n, n)) : Hermitian(randn(T, n, n)) + d, v = eigen(A) + for alg in (DivideAndConquer(), QRIteration(), RobustRepresentations()) + @test (@inferred eigvals(A; alg)) ≈ d + d2, v2 = @inferred eigen(A; alg) + @test d2 ≈ d + @test A * v2 ≈ v2 * Diagonal(d2) + end + end +end + end # module TestSymmetricEigen