From f5412d7cfbafd06a2a3dccb82720d79ddf89c634 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 20 Feb 2025 16:00:12 +0530 Subject: [PATCH 1/2] Move `alg` to a keyword argument in symmetric eigen --- src/symmetriceigen.jl | 32 +++++++++++++++++++------------- test/symmetric.jl | 15 --------------- test/symmetriceigen.jl | 15 +++++++++++++++ 3 files changed, 34 insertions(+), 28 deletions(-) diff --git a/src/symmetriceigen.jl b/src/symmetriceigen.jl index 68a1b29f..1d92e4f8 100644 --- a/src/symmetriceigen.jl +++ b/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}; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing) if alg === DivideAndConquer() Eigen(sorteig!(LAPACK.syevd!('V', A.uplo, A.data)..., sortby)...) elseif alg === QRIteration() @@ -22,7 +28,7 @@ function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algo 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 @@ -45,19 +51,19 @@ 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) - _eigen(A, alg; sortby) +function eigen(A::RealHermSymComplexHerm; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing) + _eigen(A; alg, sortby) end # we dispatch on the eltype in an internal method to avoid ambiguities -function _eigen(A::RealHermSymComplexHerm, alg::Algorithm; sortby) +function _eigen(A::RealHermSymComplexHerm; alg::Algorithm, sortby) S = eigtype(eltype(A)) - eigen!(eigencopy_oftype(A, S), alg; sortby) + eigen!(eigencopy_oftype(A, S); alg, sortby) end -function _eigen(A::RealHermSymComplexHerm{Float16}, alg::Algorithm; sortby::Union{Function,Nothing}=nothing) +function _eigen(A::RealHermSymComplexHerm{Float16}; alg::Algorithm, sortby::Union{Function,Nothing}=nothing) S = eigtype(eltype(A)) - E = eigen!(eigencopy_oftype(A, S), alg, sortby=sortby) + E = eigen!(eigencopy_oftype(A, S); alg, sortby) values = convert(AbstractVector{Float16}, E.values) vectors = convert(AbstractMatrix{isreal(E.vectors) ? Float16 : Complex{Float16}}, E.vectors) return Eigen(values, vectors) @@ -114,7 +120,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}; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing) vals::Vector{real(eltype(A))} = if alg === DivideAndConquer() LAPACK.syevd!('N', A.uplo, A.data) elseif alg === QRIteration() @@ -129,7 +135,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::Algorithm = default_eigen_alg(A))) -> values Return the eigenvalues of `A`. @@ -143,9 +149,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; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing) S = eigtype(eltype(A)) - eigvals!(eigencopy_oftype(A, S), alg; sortby) + eigvals!(eigencopy_oftype(A, S); alg, sortby) end diff --git a/test/symmetric.jl b/test/symmetric.jl index edd3af48..68f0ff47 100644 --- a/test/symmetric.jl +++ b/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/test/symmetriceigen.jl b/test/symmetriceigen.jl index 71087ae4..0dc168fc 100644 --- a/test/symmetriceigen.jl +++ b/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 @@ -184,4 +185,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 From e25f2b90c77214762f9a1d139d1af832adbd085d Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 20 Feb 2025 18:26:30 +0530 Subject: [PATCH 2/2] Move alg to kwarg in test --- test/symmetriceigen.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/symmetriceigen.jl b/test/symmetriceigen.jl index 0dc168fc..7a2b88a5 100644 --- a/test/symmetriceigen.jl +++ b/test/symmetriceigen.jl @@ -174,7 +174,7 @@ end @test D.vectors ≈ D32.vectors # ensure that different algorithms dispatch correctly - λ, V = eigen(C, LinearAlgebra.QRIteration()) + λ, V = eigen(C; alg=LinearAlgebra.QRIteration()) @test λ isa Vector{Float16} @test C * V ≈ V * Diagonal(λ) end