From d9e8fbaad28c422a5907bfe3e36ce33b9911d41a Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 13 Aug 2024 16:58:50 +0530 Subject: [PATCH 1/4] LinearAlgebra: move alg to a keyword argument in symmetric eigen --- stdlib/LinearAlgebra/src/symmetriceigen.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetriceigen.jl b/stdlib/LinearAlgebra/src/symmetriceigen.jl index fee524a702187..8cc693ec4180d 100644 --- a/stdlib/LinearAlgebra/src/symmetriceigen.jl +++ b/stdlib/LinearAlgebra/src/symmetriceigen.jl @@ -6,10 +6,10 @@ 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(@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 +29,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 +52,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 +109,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 +124,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 +138,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 From 9f5b852ad773db5fa0f55216d8fc832885361d2d Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 13 Aug 2024 17:26:48 +0530 Subject: [PATCH 2/4] Tests for alg --- stdlib/LinearAlgebra/test/symmetriceigen.jl | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/stdlib/LinearAlgebra/test/symmetriceigen.jl b/stdlib/LinearAlgebra/test/symmetriceigen.jl index d55d1deb6bf33..1c540d3768ec4 100644 --- a/stdlib/LinearAlgebra/test/symmetriceigen.jl +++ b/stdlib/LinearAlgebra/test/symmetriceigen.jl @@ -179,4 +179,21 @@ end @test S * v ≈ v * Diagonal(λ) end +@testset "eigvals/eigen alg" begin + n = 4 + S = Symmetric(rand(n,n)) + λ1 = eigvals(S, alg=LinearAlgebra.DivideAndConquer()) + λ2 = eigvals(S, alg=LinearAlgebra.RobustRepresentations()) + λ3 = eigvals(S, alg=LinearAlgebra.QRIteration()) + @test λ1 ≈ λ2 ≈ λ3 + + λ1, v1 = eigen(S, alg=LinearAlgebra.DivideAndConquer()) + λ2, v2 = eigen(S, alg=LinearAlgebra.RobustRepresentations()) + λ3, v3 = eigen(S, alg=LinearAlgebra.QRIteration()) + @test λ1 ≈ λ2 ≈ λ3 + @test S * v1 ≈ v1 * Diagonal(λ1) + @test S * v2 ≈ v2 * Diagonal(λ1) + @test S * v3 ≈ v3 * Diagonal(λ1) +end + end # module TestSymmetricEigen From c679ec14147cd0a3fffe65b9caaca815a76a71e1 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 13 Aug 2024 21:51:24 +0530 Subject: [PATCH 3/4] Update tests and move to test/symmetriceigen.jl --- stdlib/LinearAlgebra/test/symmetric.jl | 15 ----------- stdlib/LinearAlgebra/test/symmetriceigen.jl | 28 ++++++++++----------- 2 files changed, 13 insertions(+), 30 deletions(-) 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 1c540d3768ec4..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,21 +180,18 @@ end @test S * v ≈ v * Diagonal(λ) end -@testset "eigvals/eigen alg" begin - n = 4 - S = Symmetric(rand(n,n)) - λ1 = eigvals(S, alg=LinearAlgebra.DivideAndConquer()) - λ2 = eigvals(S, alg=LinearAlgebra.RobustRepresentations()) - λ3 = eigvals(S, alg=LinearAlgebra.QRIteration()) - @test λ1 ≈ λ2 ≈ λ3 - - λ1, v1 = eigen(S, alg=LinearAlgebra.DivideAndConquer()) - λ2, v2 = eigen(S, alg=LinearAlgebra.RobustRepresentations()) - λ3, v3 = eigen(S, alg=LinearAlgebra.QRIteration()) - @test λ1 ≈ λ2 ≈ λ3 - @test S * v1 ≈ v1 * Diagonal(λ1) - @test S * v2 ≈ v2 * Diagonal(λ1) - @test S * v3 ≈ v3 * Diagonal(λ1) +@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 69daf2a545fcc0b19d871b99f0c1d6a9035d33f2 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Fri, 25 Oct 2024 16:57:55 +0530 Subject: [PATCH 4/4] Add a docstring to default_eigen_alg --- stdlib/LinearAlgebra/src/symmetriceigen.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/stdlib/LinearAlgebra/src/symmetriceigen.jl b/stdlib/LinearAlgebra/src/symmetriceigen.jl index 8cc693ec4180d..860a09a2c810e 100644 --- a/stdlib/LinearAlgebra/src/symmetriceigen.jl +++ b/stdlib/LinearAlgebra/src/symmetriceigen.jl @@ -6,6 +6,12 @@ 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) + +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