From ece7669bc8f3d3620a08741463480013cf9bc9b0 Mon Sep 17 00:00:00 2001
From: Aravindh Krishnamoorthy <aravindh.krishnamoorthy@fau.de>
Date: Wed, 5 Jul 2023 13:48:20 +0200
Subject: [PATCH] Carry forward implementation from PR JuliaLang#49673

---
 stdlib/LinearAlgebra/src/symmetriceigen.jl  | 21 +++++++++
 stdlib/LinearAlgebra/test/symmetriceigen.jl | 48 +++++++++++++++++++++
 2 files changed, 69 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/symmetriceigen.jl b/stdlib/LinearAlgebra/src/symmetriceigen.jl
index 279577c31d664..45d71830471a1 100644
--- a/stdlib/LinearAlgebra/src/symmetriceigen.jl
+++ b/stdlib/LinearAlgebra/src/symmetriceigen.jl
@@ -184,6 +184,18 @@ function eigen!(A::AbstractMatrix, C::Cholesky; sortby::Union{Function,Nothing}=
     GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
 end
 
+function eigen(A::AbstractMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    # Bunchkaufman decomposition based eigenvalues and eigenvectors
+    if B.uplo == 'U'
+        vals, w = eigen!(ldiv(lu!(copy(B.D)), UtiAUi!(A[B.p,B.p], B.U')); sortby)
+        vecs = (B.U' \ w)[invperm(B.p),:]
+    else # B.uplo == 'L'
+        vals, w = eigen!(ldiv(lu!(copy(B.D)), UtiAUi!(A[B.p,B.p], B.L')); sortby)
+        vecs = (B.L' \ w)[invperm(B.p),:]
+    end
+    GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
+end
+
 # Perform U' \ A / U in-place, where U::Union{UpperTriangular,Diagonal}
 UtiAUi!(A, U) = _UtiAUi!(A, U)
 UtiAUi!(A::Symmetric, U) = Symmetric(_UtiAUi!(copytri!(parent(A), A.uplo), U), sym_uplo(A.uplo))
@@ -218,3 +230,12 @@ function eigvals!(A::AbstractMatrix{T}, C::Cholesky{T, <:AbstractMatrix}; sortby
     # Cholesky decomposition based eigenvalues
     return eigvals!(UtiAUi!(A, C.U); sortby)
 end
+
+# Bunch-Kaufmann (LDLT) based solution for generalized eigenvalues
+function eigvals(A::AbstractMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    if B.uplo == 'U'
+        return eigvals!(ldiv(lu!(copy(B.D)), UtiAUi!(A[B.p,B.p], B.U')); sortby)
+    else # B.uplo == 'L'
+        return eigvals!(ldiv(lu!(copy(B.D)), UtiAUi!(A[B.p,B.p], B.L')); sortby)
+    end
+end
diff --git a/stdlib/LinearAlgebra/test/symmetriceigen.jl b/stdlib/LinearAlgebra/test/symmetriceigen.jl
index c28c17255c222..0d10ef301e565 100644
--- a/stdlib/LinearAlgebra/test/symmetriceigen.jl
+++ b/stdlib/LinearAlgebra/test/symmetriceigen.jl
@@ -75,4 +75,52 @@ end
     @test eigvals(AH, BH; sortby=sf) ≈ eigvals(Hermitian(AH), Hermitian(BH); sortby=sf)
 end
 
+@testset "bk-eigen-eigals" begin
+    # Bunchkaufman decomposition based
+
+    # eigenvalue sorting
+    sf = x->(real(x),imag(x))
+
+    # Real-valued random matrix
+    N = 10
+    A = randn(N,N)
+    B = randn(N,N)
+    BH = (B+B')/2
+    # eigen
+    e0 = eigvals(A,BH; sortby=sf)
+    e,v = eigen(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
+    @test e0 ≈ e
+    @test A*v ≈ BH*v*Diagonal(e)
+    e,v = eigen(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
+    @test e0 ≈ e
+    @test A*v ≈ BH*v*Diagonal(e)
+    # eigvals
+    e0 = eigvals(A,BH; sortby=sf)
+    el = eigvals(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
+    eu = eigvals(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
+    @test e0 ≈ el
+    @test e0 ≈ eu
+
+    # Complex-valued random matrix
+    N = 10
+    A = complex.(randn(N,N),randn(N,N))
+    B = complex.(randn(N,N),randn(N,N))
+    BH = (B+B')/2
+    sf = x->(real(x),imag(x))
+    # eigen
+    e0 = eigvals(A,BH; sortby=sf)
+    e,v = eigen(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
+    @test e0 ≈ e
+    @test A*v ≈ BH*v*Diagonal(e)
+    e,v = eigen(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
+    @test e0 ≈ e
+    @test A*v ≈ BH*v*Diagonal(e)
+    # eigvals
+    e0 = eigvals(A,BH; sortby=sf)
+    el = eigvals(A,bunchkaufman(Hermitian(BH,:L)); sortby=sf)
+    eu = eigvals(A,bunchkaufman(Hermitian(BH,:U)); sortby=sf)
+    @test e0 ≈ el
+    @test e0 ≈ eu
+end
+
 end # module TestSymmetricEigen