Add interface to use symmetric eigendecomposition with different lapack method

stdlib/LinearAlgebra/src/lapack.jl

@@ -5595,11 +5595,6 @@ syevr!(jobz::AbstractChar, range::AbstractChar, uplo::AbstractChar, A::AbstractM
 Finds the eigenvalues (`jobz = N`) or eigenvalues and eigenvectors
 (`jobz = V`) of a symmetric matrix `A`. If `uplo = U`, the upper triangle
 of `A` is used. If `uplo = L`, the lower triangle of `A` is used.
-
-Use the divide-and-conquer method, instead of the QR iteration used by
-`syev!` or multiple relatively robust representations used by `syevr!`.
-See James W. Demmel et al, SIAM J. Sci. Comput. 30, 3, 1508 (2008) for
-a comparison of the accuracy and performatce of different methods.
 """
 syevd!(jobz::AbstractChar, uplo::AbstractChar, A::AbstractMatrix)
 

stdlib/LinearAlgebra/src/symmetriceigen.jl

@@ -5,14 +5,45 @@ eigencopy_oftype(A::Hermitian, S) = Hermitian(copy_similar(A, S), sym_uplo(A.upl
 eigencopy_oftype(A::Symmetric, S) = Symmetric(copy_similar(A, S), sym_uplo(A.uplo))
 
 # Eigensolvers for symmetric and Hermitian matrices
-eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) =
-    Eigen(sorteig!(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)..., sortby)...)
+function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, lapack_method::Symbol=:syevr; sortby::Union{Function,Nothing}=nothing)
+    if lapack_method === :syevr
+      Eigen(sorteig!(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)..., sortby)...)
+    elseif lapack_method === :syev
+      Eigen(sorteig!(LAPACK.syev!('V', A.uplo, A.data)..., sortby)...)
+    elseif lapack_method === :syevd
+      Eigen(sorteig!(LAPACK.syevd!('V', A.uplo, A.data)..., sortby)...)
+    else
+      throw(ArgumentError("Wrong lapack_method $lapack_method. Must be :syevr or :syev or :syevd."))
+    end
+end
+
+"""
+    eigen(A::Union{Hermitian, Symmetric}, lapack_method::Symbol=:syevr) -> 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
+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`.
+
+`lapack_method` specifies which LAPACK method to use for eigenvalue decomposition:
+- `lapack_method = :syevr` (default): Use multiple relatively robust representations
+- `lapack_method = :syev`: Use QR iteration
+- `lapack_method = :syevd`: Use divide-and-conquer method
+
+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 `lapack_method` used may change in the future.
 
-function eigen(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing)
+The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).
+"""
+function eigen(A::RealHermSymComplexHerm, lapack_method::Symbol=:syevr; sortby::Union{Function,Nothing}=nothing)
     S = eigtype(eltype(A))
-    eigen!(eigencopy_oftype(A, S), sortby=sortby)
+    eigen!(eigencopy_oftype(A, S), lapack_method; sortby)
 end
 
+
 eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, irange::UnitRange) =
     Eigen(LAPACK.syevr!('V', 'I', A.uplo, A.data, 0.0, 0.0, irange.start, irange.stop, -1.0)...)
 
@@ -63,17 +94,42 @@ function eigen(A::RealHermSymComplexHerm, vl::Real, vh::Real)
     eigen!(eigencopy_oftype(A, S), vl, vh)
 end
 
-function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing)
-    vals = LAPACK.syevr!('N', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)[1]
+
+function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, lapack_method::Symbol=:syevr; sortby::Union{Function,Nothing}=nothing)
+    if lapack_method === :syevr
+      vals = LAPACK.syevr!('N', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)[1]
+    elseif lapack_method === :syev
+      vals = LAPACK.syev!('N', A.uplo, A.data)
+    elseif lapack_method === :syevd
+      vals = LAPACK.syevd!('N', A.uplo, A.data)
+    else
+      throw(ArgumentError("Wrong lapack_method $lapack_method. Must be :syevr or :syev or :syevd."))
+    end
     !isnothing(sortby) && sort!(vals, by=sortby)
     return vals
 end
 
-function eigvals(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing)
+"""
+    eigvals(A::Union{Hermitian, Symmetric}, lapack_method::Symbol=:syevr) -> values
+
+Return the eigenvalues of `A`.
+
+`lapack_method` specifies which LAPACK method to use for eigenvalue decomposition.
+- `lapack_method = :syevr` (default): Use multiple relatively robust representations
+- `lapack_method = :syev`: Use QR iteration
+- `lapack_method = :syevd`: Use divide-and-conquer method
+
+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 `lapack_method` used may change in the future.
+"""
+function eigvals(A::RealHermSymComplexHerm, lapack_method::Symbol=:syevr; sortby::Union{Function,Nothing}=nothing)
     S = eigtype(eltype(A))
-    eigvals!(eigencopy_oftype(A, S), sortby=sortby)
+    eigvals!(eigencopy_oftype(A, S), lapack_method; sortby)
 end
 
+
 """
     eigvals!(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> values
 

stdlib/LinearAlgebra/test/symmetric.jl

@@ -571,6 +571,20 @@ end
     end
 end
 
+@testset "eigendecomposition with lapack_method" begin
+    for T in (Float64, ComplexF64)
+        n = 4
+        A = T <: Real ? Symmetric(randn(T, n,n)) : Hermitian(randn(T, n,n))
+        d, v = eigen(A)
+        for lapack_method in (:syev, :syevr, :syevd)
+            @test eigvals(A, lapack_method) ≈ d
+            d2, v2 = eigen(A, lapack_method)
+            @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