Bunchkaufman- and LU-decomposition based generalized eigenvalues and eigenvectors

base/combinatorics.jl

@@ -136,6 +136,44 @@ function permutecols!!(a::AbstractMatrix, p::AbstractVector{<:Integer})
     a
 end
 
+# Row and column permutations for AbstractMatrix
+permutecols!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
+    _permute!(a, p, Base.swapcols!)
+permuterows!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
+    _permute!(a, p, Base.swaprows!)
+@inline function _permute!(a::AbstractMatrix, p::AbstractVector{<:Integer}, swapfun!::F) where {F}
+    require_one_based_indexing(a, p)
+    p .= .-p
+    for i in 1:length(p)
+        p[i] > 0 && continue
+        j = i
+        in = p[j] = -p[j]
+        while p[in] < 0
+            swapfun!(a, in, j)
+            j = in
+            in = p[in] = -p[in]
+        end
+    end
+    a
+end
+invpermutecols!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
+    _invpermute!(a, p, Base.swapcols!)
+invpermuterows!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
+    _invpermute!(a, p, Base.swaprows!)
+@inline function _invpermute!(a::AbstractMatrix, p::AbstractVector{<:Integer}, swapfun!::F) where {F}
+    require_one_based_indexing(a, p)
+    p .= .-p
+    for i in 1:length(p)
+        p[i] > 0 && continue
+        j = p[i] = -p[i]
+        while j != i
+           swapfun!(a, j, i)
+           j = p[j] = -p[j]
+        end
+     end
+    a
+end
+
 """
     permute!(v, p)
 

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -11,11 +11,11 @@ import Base: \, /, *, ^, +, -, ==
 import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, asec, asech,
     asin, asinh, atan, atanh, axes, big, broadcast, cbrt, ceil, cis, collect, conj, convert,
     copy, copyto!, copymutable, cos, cosh, cot, coth, csc, csch, eltype, exp, fill!, floor,
-    getindex, hcat, getproperty, imag, inv, isapprox, isequal, isone, iszero, IndexStyle,
-    kron, kron!, length, log, map, ndims, one, oneunit, parent, permutedims,
-    power_by_squaring, promote_rule, real, sec, sech, setindex!, show, similar, sin,
-    sincos, sinh, size, sqrt, strides, stride, tan, tanh, transpose, trunc, typed_hcat,
-    vec, view, zero
+    getindex, hcat, getproperty, imag, inv, invpermuterows!, isapprox, isequal, isone, iszero,
+    IndexStyle, kron, kron!, length, log, map, ndims, one, oneunit, parent, permutecols!,
+    permutedims, permuterows!, power_by_squaring, promote_rule, real, sec, sech, setindex!,
+    show, similar, sin, sincos, sinh, size, sqrt, strides, stride, tan, tanh, transpose, trunc,
+    typed_hcat, vec, view, zero
 using Base: IndexLinear, promote_eltype, promote_op, promote_typeof, print_matrix,
     @propagate_inbounds, reduce, typed_hvcat, typed_vcat, require_one_based_indexing,
     splat

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
@@ -278,6 +278,27 @@ end
 Base.propertynames(B::BunchKaufman, private::Bool=false) =
     (:p, :P, :L, :U, :D, (private ? fieldnames(typeof(B)) : ())...)
 
+function getproperties!(B::BunchKaufman{T,<:StridedMatrix}) where {T<:BlasFloat}
+    # NOTE: Unlike in the 'getproperty' function, in this function L/U and D are computed in place.
+    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]
+        # Avoid aliasing dl and du.
+        dl = B.symmetric ? du : conj.(du)
+    else
+        M = UnitLowerTriangular(LUD)
+        dl = od[1:end-1]
+        # Avoid aliasing dl and du.
+        du = B.symmetric ? dl : conj.(dl)
+    end
+    return (M, Tridiagonal(dl, diag(LUD), du), B.p)
+end
+
 issuccess(B::BunchKaufman) = B.info == 0
 
 function adjoint(B::BunchKaufman)

stdlib/LinearAlgebra/src/lapack.jl

@@ -7045,4 +7045,48 @@ julia> LAPACK.lacpy!(B, A, 'U')
 """
 lacpy!(B::AbstractMatrix, A::AbstractMatrix, uplo::AbstractChar)
 
+# 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 m != lp
+                throw(DimensionMismatch("The first 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,49 @@ 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::StridedMatrix{T}, B::BunchKaufman{T,<:AbstractMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasFloat}
+    eigen!(copy(A), copy(B); sortby)
+end
+function eigen!(A::StridedMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasFloat}
+    M, TD, p = getproperties!(B)
+    # Compute generalized eigenvalues of equivalent matrix:
+    #    A' = inv(Tridiagonal(dl,d,du))*inv(M)*P*A*P'*inv(M')
+    # See: https://github.com/JuliaLang/julia/pull/50471#issuecomment-1627836781
+    permutecols!(A, p)
+    permuterows!(A, p)
+    ldiv!(M, A)
+    rdiv!(A, M')
+    ldiv!(TD, A)
+    vals, vecs = eigen!(A; sortby)
+    # Compute generalized eigenvectors from 'vecs':
+    #   vecs = P'*inv(M')*vecs
+    # See: https://github.com/JuliaLang/julia/pull/50471#issuecomment-1627836781
+    M = B.uplo == 'U' ? UnitUpperTriangular{eltype(vecs)}(M) : UnitLowerTriangular{eltype(vecs)}(M) ;
+    ldiv!(M', vecs)
+    invpermuterows!(vecs, p)
+    GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
+end
+
+# LU based solution for generalized eigenvalues and eigenvectors
+function eigen(A::StridedMatrix{T}, F::LU{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T}
+    return eigen!(copy(A), copy(F); sortby)
+end
+function eigen!(A::StridedMatrix{T}, F::LU{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T}
+    L = UnitLowerTriangular(F.L)
+    U = UpperTriangular(F.U)
+    permuterows!(A, F.p)
+    ldiv!(L, A)
+    rdiv!(A, U)
+    vals, vecs = eigen!(A; sortby)
+    # Compute generalized eigenvectors from 'vecs':
+    #   vecs = P'*inv(M')*vecs
+    # See: https://github.com/JuliaLang/julia/pull/50471#issuecomment-1627836781
+    U = UpperTriangular{eltype(vecs)}(U)
+    ldiv!(U, vecs)
+    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 +261,36 @@ 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::StridedMatrix{T}, B::BunchKaufman{T,<:AbstractMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasFloat}
+    eigvals!(copy(A), copy(B); sortby)
+end
+function eigvals!(A::StridedMatrix{T}, B::BunchKaufman{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasFloat}
+    M, TD, p = getproperties!(B)
+    # Compute generalized eigenvalues of equivalent matrix:
+    #    A' = inv(Tridiagonal(dl,d,du))*inv(M)*P*A*P'*inv(M')
+    # See: https://github.com/JuliaLang/julia/pull/50471#issuecomment-1627836781
+    permutecols!(A, p)
+    permuterows!(A, p)
+    ldiv!(M, A)
+    rdiv!(A, M')
+    ldiv!(TD, A)
+    return eigvals!(A; sortby)
+end
+
+# LU based solution for generalized eigenvalues
+function eigvals(A::StridedMatrix{T}, F::LU{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T}
+    return eigvals!(copy(A), copy(F); sortby)
+end
+function eigvals!(A::StridedMatrix{T}, F::LU{T,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) where {T}
+    L = UnitLowerTriangular(F.L)
+    U = UpperTriangular(F.U)
+    # Compute generalized eigenvalues of equivalent matrix:
+    #    A' = inv(L)*(P*A)*inv(U)
+    # See: https://github.com/JuliaLang/julia/pull/50471#issuecomment-1627836781
+    permuterows!(A, F.p)
+    ldiv!(L, A)
+    rdiv!(A, U)
+    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,72 @@ end
     @test eigvals(AH, BH; sortby=sf) ≈ eigvals(Hermitian(AH), Hermitian(BH); sortby=sf)
 end
 
+@testset "bk-lu-eigen-eigvals" 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)
+    e,v = eigen(A,lu(Hermitian(BH,:L)); sortby=sf)
+    @test e0 ≈ e
+    @test A*v ≈ BH*v*Diagonal(e)
+    e,v = eigen(A,lu(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
+    el = eigvals(A,lu(Hermitian(BH,:L)); sortby=sf)
+    eu = eigvals(A,lu(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)
+    e,v = eigen(A,lu(Hermitian(BH,:L)); sortby=sf)
+    @test e0 ≈ e
+    @test A*v ≈ BH*v*Diagonal(e)
+    e,v = eigen(A,lu(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
+    el = eigvals(A,lu(Hermitian(BH,:L)); sortby=sf)
+    eu = eigvals(A,lu(Hermitian(BH,:U)); sortby=sf)
+    @test e0 ≈ el
+    @test e0 ≈ eu
+end
+
 end # module TestSymmetricEigen