Fix for JuliaLang#49533 (2): Remove Cholesky decomposition based eige…

…n! functions for generalised Normal+Symmetric matrices and forward to the generic functions in eigen.jl
aravindh-krishnamoorthy · May 6, 2023 · 4701fbc · 4701fbc
1 parent 3b564f2
commit 4701fbc
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 12 deletions.
diff --git a/stdlib/LinearAlgebra/src/symmetriceigen.jl b/stdlib/LinearAlgebra/src/symmetriceigen.jl
@@ -164,17 +164,11 @@ function eigen!(A::Hermitian{T,S}, B::Hermitian{T,S}; sortby::Union{Function,Not
     vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))
     GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
 end
-function eigen!(A::RealHermSymComplexHerm{T,<:StridedMatrix}, B::AbstractMatrix{T}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
-    return _choleigen!(A, B, sortby)
+function eigen!(A::RealHermSymComplexHerm{T,S}, B::AbstractMatrix{T}; sortby::Union{Function,Nothing}=nothing) where {T<:Number,S<:StridedMatrix}
+    return eigen!(Matrix{T}(A), B; sortby) ;
 end
-function eigen!(A::StridedMatrix{T}, B::Union{RealHermSymComplexHerm{T},Diagonal{T}}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
-    return _choleigen!(A, B, sortby)
-end
-function _choleigen!(A, B, sortby)
-    U = cholesky(B).U
-    vals, w = eigen!(UtiAUi!(A, U))
-    vecs = U \ w
-    GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
+function eigen!(A::AbstractMatrix{T}, B::RealHermSymComplexHerm{T,S}; sortby::Union{Function,Nothing}=nothing) where {T<:Number,S<:StridedMatrix}
+    return eigen!(A, Matrix{T}(B); sortby) ;
 end
 
 # Perform U' \ A / U in-place, where U::Union{UpperTriangular,Diagonal}
@@ -197,10 +191,10 @@ end
 eigvecs(A::HermOrSym) = eigvecs(eigen(A))
 
 # Note: No specilized LAPACK routines for Matrix+Symmetric and Symmetric+Matrix exist. Hence, the calls are forwarded to conventional Matrix eigvals!
-function eigvals!(A::AbstractMatrix{T}, B::HermOrSym{T,S}; sortby::Union{Function,Nothing}=nothing) where {T<:Union{BlasReal,BlasComplex},S<:StridedMatrix}
+function eigvals!(A::AbstractMatrix{T}, B::RealHermSymComplexHerm{T,S}; sortby::Union{Function,Nothing}=nothing) where {T<:Number,S<:StridedMatrix}
     return eigvals!(A, Matrix{T}(B); sortby) ;
 end
 
-function eigvals!(A::HermOrSym{T,S}, B::AbstractMatrix{T}; sortby::Union{Function,Nothing}=nothing) where {T<:Union{BlasReal,BlasComplex},S<:StridedMatrix}
+function eigvals!(A::RealHermSymComplexHerm{T,S}, B::AbstractMatrix{T}; sortby::Union{Function,Nothing}=nothing) where {T<:Number,S<:StridedMatrix}
     return eigvals!(Matrix{T}(A), B; sortby) ;
 end
diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl
@@ -827,6 +827,13 @@ end
 @testset "issue #49533" begin
     A = Float64[1 1 0 0; 1 2 1 0; 0 1 3 1; 0 0 1 4] ;
     B = Matrix(Diagonal(Float64[1:4;])) ;
+    # eigen
+    e0,v0 = eigen(A, B)
+    e1,v1 = eigen(A, Symmetric(B))
+    e2,v2 = eigen(Symmetric(A), B)
+    @test e0 ≈ e1 && v0 ≈ v1
+    @test e0 ≈ e2 && v0 ≈ v2
+    # eigvals
     @test eigvals(A, B) ≈ eigvals(A, Symmetric(B))
     @test eigvals(A, B) ≈ eigvals(Symmetric(A), B)
 end