Materialize complex Symmetric matrices in eigen (#55348)

Currently, `eigen` for a complex Symmetric matrix fails, as there's no
specialized LAPACK function to handle such matrices. We may instead
materialize the matrix and use a generic solver. While a user may do it
by themselves, I think returning an answer is better than throwing an
error.

Author: jishnub

stdlib/LinearAlgebra/src/bunchkaufman.jl

@@ -127,6 +127,9 @@ function bunchkaufman!(A::StridedMatrix{<:BlasFloat}, rook::Bool = false; check:
     end
 end
 
+bkcopy_oftype(A, S) = eigencopy_oftype(A, S)
+bkcopy_oftype(A::Symmetric{<:Complex}, S) = Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo))
+
 """
     bunchkaufman(A, rook::Bool=false; check = true) -> S::BunchKaufman
 
@@ -206,7 +209,7 @@ julia> S.L*S.D*S.L' - A[S.p, S.p]
 ```
 """
 bunchkaufman(A::AbstractMatrix{T}, rook::Bool=false; check::Bool = true) where {T} =
-    bunchkaufman!(eigencopy_oftype(A, typeof(sqrt(oneunit(T)))), rook; check = check)
+    bunchkaufman!(bkcopy_oftype(A, typeof(sqrt(oneunit(T)))), rook; check = check)
 
 BunchKaufman{T}(B::BunchKaufman) where {T} =
     BunchKaufman(convert(Matrix{T}, B.LD), B.ipiv, B.uplo, B.symmetric, B.rook, B.info)
@@ -1540,7 +1543,7 @@ function bunchkaufman(A::AbstractMatrix{TS},
     rook::Bool = false;
     check::Bool = true
     ) where TS <: ClosedScalar{TR} where TR <: ClosedReal
-    return bunchkaufman!(eigencopy_oftype(A, TS), rook; check)
+    return bunchkaufman!(bkcopy_oftype(A, TS), rook; check)
 end
 
 function bunchkaufman(A::AbstractMatrix{TS},
@@ -1562,15 +1565,15 @@ function bunchkaufman(A::AbstractMatrix{TS},
     # We promote input to BigInt to avoid overflow problems
     if TA == Nothing
         if TS <: Integer
-            M = Rational{BigInt}.(eigencopy_oftype(A, TS))
+            M = Rational{BigInt}.(bkcopy_oftype(A, TS))
         else
-            M = Complex{Rational{BigInt}}.(eigencopy_oftype(A, TS))
+            M = Complex{Rational{BigInt}}.(bkcopy_oftype(A, TS))
         end
     else
         if TS <: Integer
-            M = TA(Rational{BigInt}.(eigencopy_oftype(A, TS)), Symbol(A.uplo))
+            M = TA(Rational{BigInt}.(bkcopy_oftype(A, TS)), Symbol(A.uplo))
         else
-            M = TA(Complex{Rational{BigInt}}.(eigencopy_oftype(A, TS)),
+            M = TA(Complex{Rational{BigInt}}.(bkcopy_oftype(A, TS)),
                 Symbol(A.uplo))
         end
     end

stdlib/LinearAlgebra/src/symmetriceigen.jl

@@ -4,6 +4,7 @@
 # Call `copytrito!` instead of `copy_similar` to only copy the matching triangular half
 eigencopy_oftype(A::Hermitian, S) = Hermitian(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo))
 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()
 

stdlib/LinearAlgebra/test/hessenberg.jl

@@ -272,4 +272,11 @@ end
     @test S[1,2] == S[Int8(1),UInt16(2)] == S[big(1), Int16(2)]
 end
 
+@testset "complex Symmetric" begin
+    D = diagm(0=>ComplexF64[1,2])
+    S = Symmetric(D)
+    H = hessenberg(S)
+    @test H.H == D
+end
+
 end # module TestHessenberg

stdlib/LinearAlgebra/test/symmetriceigen.jl

@@ -173,4 +173,10 @@ end
     @test D.vectors ≈ D32.vectors
 end
 
+@testset "complex Symmetric" begin
+    S = Symmetric(rand(ComplexF64,2,2))
+    λ, v = eigen(S)
+    @test S * v ≈ v * Diagonal(λ)
+end
+
 end # module TestSymmetricEigen