Bunchkaufman- and LU-decomposition based generalized eigenvalues and eigenvectors

stdlib/LinearAlgebra/src/bunchkaufman.jl

@@ -88,7 +88,7 @@ Base.iterate(S::BunchKaufman) = (S.D, Val(:UL))
 Base.iterate(S::BunchKaufman, ::Val{:UL}) = (S.uplo == 'L' ? S.L : S.U, Val(:p))
 Base.iterate(S::BunchKaufman, ::Val{:p}) = (S.p, Val(:done))
 Base.iterate(S::BunchKaufman, ::Val{:done}) = nothing
-
+copy(S::BunchKaufman) = BunchKaufman(copy(S.LD), copy(S.ipiv), S.uplo, S.symmetric, S.rook, S.info)
 
 """
     bunchkaufman!(A, rook::Bool=false; check = true) -> BunchKaufman

stdlib/LinearAlgebra/src/lapack.jl

@@ -6941,4 +6941,48 @@ Returns `X` (overwriting `C`) and `scale`.
 """
 trsyl!(transa::AbstractChar, transb::AbstractChar, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, isgn::Int=1)
 
+# Routines for column permutation
+for (fn, elty) in ((:slapmt_, :Float32),
+                (:dlapmt_, :Float64),
+                (:clapmt_, :ComplexF32),
+                (:zlapmt_, :ComplexF64))
+    @eval begin
+        function lapmt!(A::AbstractMatrix{$elty}, p::AbstractVector{<:BlasInt}, fb::Bool)
+            require_one_based_indexing(A)
+            chkstride1(A)
+            m, n = size(A)
+            lda = max(1, stride(A, 2))
+            lp = length(p)
+            if n != lp
+                throw(DimensionMismatch("The second dimension of A, ($m,$n), and P, ($lp), must match"))
+            end
+            ccall((@blasfunc($fn), libblastrampoline), Cvoid,
+            (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{BlasInt}),
+            fb, m, n, A, lda, p) ;
+        end
+    end
+end
+
+# Routines for row permutation
+for (fn, elty) in ((:slapmr_, :Float32),
+                (:dlapmr_, :Float64),
+                (:clapmr_, :ComplexF32),
+                (:zlapmr_, :ComplexF64))
+    @eval begin
+        function lapmr!(A::AbstractMatrix{$elty}, p::AbstractVector{<:BlasInt}, fb::Bool)
+            require_one_based_indexing(A)
+            chkstride1(A)
+            m, n = size(A)
+            lda = max(1, stride(A, 2))
+            lp = length(p)
+            if n != lp
+                throw(DimensionMismatch("The second dimension of A, ($m,$n), and P, ($lp), must match"))
+            end
+            ccall((@blasfunc($fn), libblastrampoline), Cvoid,
+            (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{BlasInt}),
+            fb, m, n, A, lda, p) ;
+        end
+    end
+end
+
 end # module

stdlib/LinearAlgebra/src/symmetriceigen.jl

@@ -184,6 +184,50 @@ function eigen!(A::AbstractMatrix, C::Cholesky; sortby::Union{Function,Nothing}=
     GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
 end
 
+# Bunch-Kaufmann (LDLT) based solution for generalized eigenvalues and eigenvectors
+function eigen(A::AbstractMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    eigen!(copy(A), copy(B); sortby)
+end
+function eigen!(A::AbstractMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+   # NOTE: getproperty(B, :D) and getproperty(B, :L/U) are not in-place. Hence, in the following, sysconvf/_rook is
+   # directly utilized to obtain M, du, d = diag(LUD), and dl. See bunchkaufman.jl/getproperty for details.
+   if B.rook
+        LUD, od = LAPACK.syconvf_rook!(B.uplo, 'C', B.LD, B.ipiv)
+    else
+        LUD, od = LAPACK.syconv!(B.uplo, B.LD, B.ipiv)
+    end
+    if B.uplo == 'U'
+        M = UnitUpperTriangular(LUD)
+        du = od[2:end]
+        # Aliasing of dl and du is not allowed for gtsv! below.
+        dl = B.symmetric ? copy(du) : conj.(copy(du))
+    else
+        M = UnitLowerTriangular(LUD)
+        dl = od[1:end-1]
+        # Aliasing of dl and du is not allowed for gtsv! below.
+        du = B.symmetric ? copy(dl) : conj.(copy(dl))
+    end
+    # In-place permutation of A, A .= A[B.p, B.p]
+    LAPACK.lapmt!(A, B.p, true)
+    LAPACK.lapmr!(A, B.p, true)
+    # In-place A .= inv(M)*A
+    LAPACK.trtrs!(B.uplo, 'N', 'U', M.data, A)
+    # In-place A .= A*inv(M^H)
+    BLAS.trsm!('R', B.uplo , 'C', 'U', one(T), M.data, A)
+    # In-place A .= inv(Tridiagonal(dl, d = diag(LUD), du))*A
+    LAPACK.gtsv!(dl, diag(LUD), du, A)
+    # Compute the eigenvalues and eigenvectors of A, which are the generalized eigenvalues
+    # and eigenvectors.
+    # NOTE: That vecs overwrites A for matrices with real-valued generalized eigenvalues
+    # and complex-valued eigenvalues. In other cases, a copy is created, see eigen.jl for details.
+    vals, vecs = eigen!(A; sortby)
+    # In-place vecs .= inv(M^H)*vecs
+    LAPACK.trtrs!(B.uplo, 'C', 'U', convert(typeof(vecs), M.data), vecs)
+    # In-place vecs .= vecs[invperm(B.p), :]
+    LAPACK.lapmr!(vecs, invperm(B.p), true)
+    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 +262,39 @@ 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}
+    eigvals!(copy(A), copy(B); sortby)
+end
+function eigvals!(A::AbstractMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    # NOTE: getproperty(B, :D) and getproperty(B, :L/U) are not in-place. Hence, in the following, sysconvf/_rook is
+    # directly utilized to obtain M, du, d = diag(LUD), and dl. See bunchkaufman.jl/getproperty for details.
+    if B.rook
+        LUD, od = LAPACK.syconvf_rook!(B.uplo, 'C', B.LD, B.ipiv)
+    else
+        LUD, od = LAPACK.syconv!(B.uplo, B.LD, B.ipiv)
+    end
+    if B.uplo == 'U'
+        M = UnitUpperTriangular(LUD)
+        du = od[2:end]
+        # Aliasing of dl and du is not allowed for gtsv! below.
+        dl = B.symmetric ? copy(du) : conj.(copy(du))
+    else
+        M = UnitLowerTriangular(LUD)
+        dl = od[1:end-1]
+        # Aliasing of dl and du is not allowed for gtsv! below.
+        du = B.symmetric ? copy(dl) : conj.(copy(dl))
+    end
+    # In-place permutation of A, A .= A[B.p, B.p]
+    LAPACK.lapmt!(A, B.p, true)
+    LAPACK.lapmr!(A, B.p, true)
+    # In-place A .= inv(M)*A
+    LAPACK.trtrs!(B.uplo, 'N', 'U', M.data, A)
+    # In-place A .= A*inv(M^H)
+    BLAS.trsm!('R', B.uplo , 'C', 'U', one(T), M.data, A)
+    # In-place A .= inv(Tridiagonal(dl, d = diag(LUD), du))*A
+    LAPACK.gtsv!(dl, diag(LUD), du, A)
+    # Compute the eigenvalues of A, which are the generalized eigenvalues.
+    return eigvals!(A; sortby)
+end

stdlib/LinearAlgebra/test/symmetriceigen.jl

@@ -8,7 +8,7 @@ using Test, LinearAlgebra
     ## Cholesky decomposition based
 
     # eigenvalue sorting
-    sf = x->(real(x),imag(x))
+    sf = x->(imag(x),real(x))
 
     ## Real valued
     A = Float64[1 1 0 0; 1 2 1 0; 0 1 3 1; 0 0 1 4]
@@ -40,6 +40,9 @@ using Test, LinearAlgebra
 end
 
 @testset "issue #49533" begin
+    # eigenvalue sorting
+    sf = x->(imag(x),real(x))
+
     ## Real valued
     A = Float64[1 1 0 0; 1 2 1 0; 0 1 3 1; 0 0 1 4]
     B = Matrix(Diagonal(Float64[1:4;]))
@@ -62,7 +65,6 @@ end
     B =  [2.0+2.0im 1.0+1.0im 4.0+4.0im 3.0+3.0im; 0 3.0+2.0im 1.0+1.0im 3.0+4.0im; 3.0+3.0im 1.0+4.0im 0 0; 0 1.0+2.0im 3.0+1.0im 1.0+1.0im]
     BH = B'B
     # eigen
-    sf = x->(real(x),imag(x))
     e1,v1 = eigen(A, Hermitian(BH))
     @test A*v1 ≈ Hermitian(BH)*v1*Diagonal(e1)
     e2,v2 = eigen(Hermitian(AH), B)
@@ -75,4 +77,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->(imag(x),real(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