RFC: Methods for eigen on complex hermitian tridiagonal matrices (J…

…uliaLang#49546)

Co-authored-by: Daniel Karrasch <[email protected]>
Co-authored-by: Steven G. Johnson <[email protected]>

NEWS.md

@@ -115,6 +115,7 @@ Standard library changes
 * `eigvals/eigen(A, bunchkaufman(B))` and `eigvals/eigen(A, lu(B))`, which utilize the Bunchkaufman (LDL) and LU decomposition of `B`,
  respectively, now efficiently compute the generalized eigenvalues (`eigen`: and eigenvectors) of `A` and `B`. Note: The second
  argument is the output of `bunchkaufman` or `lu` ([#50471]).
+* There is now a specialized dispatch for `eigvals/eigen(::Hermitian{<:Tridiagonal})` which performs a similarity transformation to create a real symmetrix triagonal matrix, and solve that using the LAPACK routines ([#49546]).
 * Structured matrices now retain either the axes of the parent (for `Symmetric`/`Hermitian`/`AbstractTriangular`/`UpperHessenberg`), or that of the principal diagonal (for banded matrices) ([#52480]).
 * `bunchkaufman` and `bunchkaufman!` now work for any `AbstractFloat`, `Rational` and their complex variants. `bunchkaufman` now supports `Integer` types, by making an internal conversion to `Rational{BigInt}`. Added new function `inertia` that computes the inertia of the diagonal factor given by the `BunchKaufman` factorization object of a real symmetric or Hermitian matrix. For complex symmetric matrices, `inertia` only computes the number of zero eigenvalues of the diagonal factor ([#51487]).
 

stdlib/LinearAlgebra/src/symmetriceigen.jl

@@ -294,3 +294,34 @@ function eigvals!(A::StridedMatrix{T}, F::LU{T,<:StridedMatrix}; sortby::Union{F
  rdiv!(A, U)
  return eigvals!(A; sortby)
 end
+
+
+function eigen(A::Hermitian{Complex{T}, <:Tridiagonal}; kwargs...) where {T}
+ (; dl, d, du) = parent(A)
+ N = length(d)
+ if N <= 1
+ eigen(parent(A); kwargs...)
+ else
+ if A.uplo == 'U'
+ E = du'
+ Er = abs.(du)
+ else
+ E = dl
+ Er = abs.(E)
+ end
+ S = Vector{Complex{T}}(undef, N)
+ S[1] = 1
+ for i ∈ 1:N-1
+ S[i+1] = iszero(Er[i]) ? one(Complex{T}) : S[i] * sign(E[i])
+ end
+ B = SymTridiagonal(real(d), Er)
+ Λ, Φ = eigen(B; kwargs...)
+ return Eigen(Λ, Diagonal(S) * Φ)
+ end
+end
+
+function eigvals(A::Hermitian{Complex{T}, <:Tridiagonal}; kwargs...) where {T}
+ (; dl, d, du) = parent(A)
+ E = A.uplo == 'U' ? abs.(du) : abs.(dl)
+ eigvals(SymTridiagonal(real.(d), E); kwargs...)
+end

stdlib/LinearAlgebra/test/eigen.jl

@@ -45,12 +45,13 @@ aimg = randn(n,n)/2
  @test eigvecs(f) === f.vectors
  @test Array(f) ≈ a
 
- for T in (Tridiagonal(a), Hermitian(Tridiagonal(a)))
+ for T in (Tridiagonal(a), Hermitian(Tridiagonal(a), :U), Hermitian(Tridiagonal(a), :L))
  f = eigen(T)
  d, v = f
  for i in 1:size(a,2)
  @test T*v[:,i] ≈ d[i]*v[:,i]
  end
+ @test eigvals(T) ≈ d
  @test det(T) ≈ det(f)
  @test inv(T) ≈ inv(f)
  end