Eigvecs for specific eigvals for Symmetric/Hermitian

Author: jishnub

src/symmetriceigen.jl

@@ -331,6 +331,12 @@ function eigvals!(A::Hermitian{T,S}, B::Hermitian{T,S}; sortby::Union{Function,N
 end
 eigvecs(A::HermOrSym) = eigvecs(eigen(A))
 
+function eigvecs(A::RealHermSymComplexHerm, eigvals::AbstractVector{<:Real})
+    F = hessenberg(A) # transform to SymTridiagonal form
+    X = eigvecs(F.H, eigvals)
+    return F.Q * X    # transform eigvecs of F.H back to eigvecs of A
+end
+
 function eigvals(A::AbstractMatrix, C::Cholesky; sortby::Union{Function,Nothing}=nothing)
     if ishermitian(A)
         eigvals!(eigencopy_oftype(Hermitian(A), eigtype(eltype(A))), C; sortby)

src/tridiag.jl

@@ -311,7 +311,7 @@ eigmin(A::SymTridiagonal) = eigvals(A, 1:1)[1]
 
 #Compute selected eigenvectors only corresponding to particular eigenvalues
 """
-    eigvecs(A::SymTridiagonal[, eigvals]) -> Matrix
+    eigvecs(A::Union{Symmetric, Hermitian, SymTridiagonal}[, eigvals])::Matrix
 
 Return a matrix `M` whose columns are the eigenvectors of `A`. (The `k`th eigenvector can
 be obtained from the slice `M[:, k]`.)
@@ -346,7 +346,7 @@ julia> eigvecs(A, [1.])
  -0.5547001962252291
 ```
 """
-eigvecs(A::SymTridiagonal{<:BlasFloat,<:StridedVector}, eigvals::Vector{<:Real}) = LAPACK.stein!(A.dv, A.ev, eigvals)
+eigvecs(A::SymTridiagonal{<:BlasFloat,<:StridedVector}, eigvals::StridedVector{<:Real}) = LAPACK.stein!(A.dv, A.ev, eigvals)
 
 function svdvals!(A::SymTridiagonal)
     vals = eigvals!(A)

test/symmetriceigen.jl

@@ -199,4 +199,23 @@ end
     end
 end
 
+@testset "eigvecs for specific eigvals" begin
+    function testeigvecs(S, vals)
+        V = eigvecs(S, vals)
+        @test S * V ≈ V * Diagonal(vals)
+    end
+    for T in (Symmetric, Hermitian)
+        S = T(rand(3,3))
+        vals = eigvals(S)
+        testeigvecs(S, vals)
+        testeigvecs(S, vals[1:2])
+        testeigvecs(S, @view vals[2:2])
+    end
+    H = Hermitian(rand(ComplexF64,3,3))
+    vals = eigvals(H)
+    testeigvecs(H, vals)
+    testeigvecs(H, vals[1:2])
+    testeigvecs(H, @view vals[2:2])
+end
+
 end # module TestSymmetricEigen