From 7194038d84d98fabb8f78134e8a45b4a63471c2c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mateus=20Ara=C3=BAjo?= Date: Wed, 28 May 2025 17:17:40 +0200 Subject: [PATCH 1/2] clarify eigen docstring (#1364) Closes #1174 `eigen` is in fact guaranteed to return normalized eigenvectors in the general case, and orthonormal eigenvectors in the Hermitian case, and the docstring should state so. (cherry picked from commit c4a477fff61b1f72c7d867e3eebad8cc1eab4d8b) --- src/eigen.jl | 2 +- src/symmetriceigen.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/eigen.jl b/src/eigen.jl index e0124f2e..0837491a 100644 --- a/src/eigen.jl +++ b/src/eigen.jl @@ -182,7 +182,7 @@ end eigen(A; permute::Bool=true, scale::Bool=true, sortby) -> 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 +which contains the eigenvalues in `F.values` and the normalized eigenvectors in the columns of the matrix `F.vectors`. This corresponds to solving an eigenvalue problem of the form `Ax = λx`, where `A` is a matrix, `x` is an eigenvector, and `λ` is an eigenvalue. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.) diff --git a/src/symmetriceigen.jl b/src/symmetriceigen.jl index 1d92e4f8..08261f72 100644 --- a/src/symmetriceigen.jl +++ b/src/symmetriceigen.jl @@ -31,7 +31,7 @@ end 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 +which contains the eigenvalues in `F.values` and the orthonormal eigenvectors in the columns of the matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.) Iterating the decomposition produces the components `F.values` and `F.vectors`. @@ -76,7 +76,7 @@ eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, irange::UnitRange) eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> 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 +which contains the eigenvalues in `F.values` and the orthonormal eigenvectors in the columns of the matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.) Iterating the decomposition produces the components `F.values` and `F.vectors`. @@ -101,7 +101,7 @@ eigen!(A::RealHermSymComplexHerm{T,<:StridedMatrix}, vl::Real, vh::Real) where { eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> 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 +which contains the eigenvalues in `F.values` and the orthonormal eigenvectors in the columns of the matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.) Iterating the decomposition produces the components `F.values` and `F.vectors`. From 73ff52d47bade62b7dad8fc782b7f031cc9f2e2b Mon Sep 17 00:00:00 2001 From: Jae-Mo Lihm Date: Sun, 1 Jun 2025 13:48:02 +0200 Subject: [PATCH 2/2] Change default symmetriceigen algorithm back to RobustRepresentations (#1363) See #1313 A brief summary: RobustRepresentations (`LAPACK.syevr!`) was the only algorithm implemented until Julia v1.11. Then, https://github.com/JuliaLang/julia/pull/49355/ implemented interface to other algorithms, and changed the default to DivideAndConquer (`LAPACK.syevd!`) based on its better numerical accuracy and performance. But, it turned out that in some LAPACK implementation, the DivideAndConquer fails more frequently than RobustRepresentations (#1313). Based on the discussion in #1313 , this PR reverts the default algorithm to RobustRepresentations. Once the problem with the RobustRepresentations in the problematic LAPACK implementation has been sorted out, the default could be changed to DivideAndConquer. (cherry picked from commit b6ad8f9a05292266fe510d0d3be1373972cafb28) --- src/symmetriceigen.jl | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/symmetriceigen.jl b/src/symmetriceigen.jl index 08261f72..652698c0 100644 --- a/src/symmetriceigen.jl +++ b/src/symmetriceigen.jl @@ -10,9 +10,9 @@ 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!`. +Defaults to `LinearAlegbra.RobustRepresentations()`, which corresponds to the LAPACK function `LAPACK.syevr!`. """ -default_eigen_alg(@nospecialize(A)) = DivideAndConquer() +default_eigen_alg(@nospecialize(A)) = RobustRepresentations() # Eigensolvers for symmetric and Hermitian matrices function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing) @@ -37,9 +37,9 @@ matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vec Iterating the decomposition produces the components `F.values` and `F.vectors`. `alg` specifies which algorithm and LAPACK method to use for eigenvalue decomposition: -- `alg = DivideAndConquer()` (default): Calls `LAPACK.syevd!`. +- `alg = DivideAndConquer()`: Calls `LAPACK.syevd!`. - `alg = QRIteration()`: Calls `LAPACK.syev!`. -- `alg = RobustRepresentations()`: Multiple relatively robust representations method, Calls `LAPACK.syevr!`. +- `alg = RobustRepresentations()` (default): Multiple relatively robust representations method, Calls `LAPACK.syevr!`. See James W. Demmel et al, SIAM J. Sci. Comput. 30, 3, 1508 (2008) for a comparison of the accuracy and performance of different algorithms. @@ -140,14 +140,18 @@ end Return the eigenvalues of `A`. `alg` specifies which algorithm and LAPACK method to use for eigenvalue decomposition: -- `alg = DivideAndConquer()` (default): Calls `LAPACK.syevd!`. +- `alg = DivideAndConquer()`: Calls `LAPACK.syevd!`. - `alg = QRIteration()`: Calls `LAPACK.syev!`. -- `alg = RobustRepresentations()`: Multiple relatively robust representations method, Calls `LAPACK.syevr!`. +- `alg = RobustRepresentations()` (default): Multiple relatively robust representations method, Calls `LAPACK.syevr!`. See James W. Demmel et al, SIAM J. Sci. Comput. 30, 3, 1508 (2008) for a comparison of the accuracy and performance of different methods. The default `alg` used may change in the future. + +!!! compat "Julia 1.12" + The `alg` keyword argument requires Julia 1.12 or later. + """ function eigvals(A::RealHermSymComplexHerm; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing) S = eigtype(eltype(A))