Generalized eigen for complex Hermitian BandedMatrix

src/lapack.jl

@@ -371,3 +371,76 @@ for (fname, elty) in ((:dsbgst_,:Float64),
         end
     end
 end
+
+# Convert a complex Hermitian Positive Definite generalized eigenvalue problem
+# to a symmetric eigenvalue problem assuming B has been processed by
+# a split-Cholesky factorization.
+for (fname, elty) in ((:zhbgst_, :ComplexF64),
+                      (:chbgst_, :ComplexF32))
+    @eval begin
+        #=
+        subroutine zhbgst   (
+        character   VECT,
+        character   UPLO,
+        integer     N,
+        integer     KA,
+        integer     KB,
+        complex*16, dimension( ldab, * )    AB,
+        integer     LDAB,
+        complex*16, dimension( ldbb, * )    BB,
+        integer     LDBB,
+        complex*16, dimension( ldx, * )     X,
+        integer     LDX,
+        complex*16, dimension( * )      WORK,
+        double precision, dimension( * )    RWORK,
+        integer     INFO
+        )
+        =#
+        local Relty = real($elty)
+        function hbgst!(vect::Char, uplo::Char, n::Int, ka::Int, kb::Int,
+                         AB::StridedMatrix{$elty}, BB::StridedMatrix{$elty},
+                         X::StridedMatrix{$elty}, work::StridedVector{$elty},
+                         rwork::StridedVector{Relty})
+            require_one_based_indexing(AB, BB, X, work)
+            chkstride1(AB, BB)
+            chkuplo(uplo)
+            chkvect(vect)
+            info  = Ref{BlasInt}()
+            ccall((@blasfunc($fname), liblapack), Nothing,
+                (
+                    Ref{UInt8},
+                    Ref{UInt8},
+                    Ref{BlasInt},
+                    Ref{BlasInt},
+                    Ref{BlasInt},
+                    Ptr{$elty},
+                    Ref{BlasInt},
+                    Ptr{$elty},
+                    Ref{BlasInt},
+                    Ptr{$elty},
+                    Ref{BlasInt},
+                    Ptr{$elty},
+                    Ptr{Relty},
+                    Ref{BlasInt},
+                ),
+                vect,
+                uplo,
+                n,
+                ka,
+                kb,
+                AB,
+                max(stride(AB,2),1),
+                BB,
+                max(stride(BB,2),1),
+                X,
+                max(stride(X,2),1),
+                work,
+                rwork,
+                info,
+                )
+
+            LAPACK.chklapackerror(info[])
+            AB
+        end
+    end
+end

src/symbanded/SplitCholesky.jl

@@ -33,18 +33,24 @@ SplitCholesky(factors::AbstractMatrix{T}, uplo::Char) where T = SplitCholesky{T,
 size(S::SplitCholesky) = size(S.factors)
 size(S::SplitCholesky, i::Integer) = size(S.factors, i)
 
-splitcholesky!(A::Symmetric{T,<:BandedMatrix{T}}) where T = splitcholesky!(A, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
-splitcholesky!(A::Symmetric{T,<:BandedMatrix{T}}, ::Type{Tr}) where {T,Tr} = splitcholesky!(MemoryLayout(typeof(A)), A, Tr)
+function splitcholesky!(A::HermOrSym{T,<:BandedMatrix{T}}) where T
+    splitcholesky!(A, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
+end
+function splitcholesky!(A::HermOrSym{T,<:BandedMatrix{T}}, Tr) where {T}
+    splitcholesky!(MemoryLayout(typeof(A)), A, Tr)
+end
 
 function splitcholesky!(::SymmetricLayout{<:BandedColumnMajor},
-                       A::AbstractMatrix{T}, ::Type{UpperTriangular}) where T<:BlasFloat
-    _, info = pbstf!('U', size(A, 1), bandwidth(A,2), symbandeddata(A))
+                       A::AbstractMatrix{T}, ::Type{LU}) where {T<:BlasFloat, LU}
+    uplo = LU == UpperTriangular ? 'U' : 'L'
+    pbstf!(uplo, size(A, 1), bandwidth(A,2), symbandeddata(A))
     SplitCholesky(A, A.uplo)
 end
 
-function splitcholesky!(::SymmetricLayout{<:BandedColumnMajor},
-                       A::AbstractMatrix{T}, ::Type{LowerTriangular}) where T<:BlasFloat
-    _, info = pbstf!('L', size(A, 1), bandwidth(A,1), symbandeddata(A))
+function splitcholesky!(::HermitianLayout{<:BandedColumnMajor},
+                       A::AbstractMatrix{T}, ::Type{LU}) where {T<:BlasFloat, LU}
+    uplo = LU == UpperTriangular ? 'U' : 'L'
+    pbstf!(uplo, size(A, 1), bandwidth(A,2), hermbandeddata(A))
     SplitCholesky(A, A.uplo)
 end
 
@@ -56,7 +62,7 @@ else
         throw(ArgumentError("transpose of SplitCholesky decomposition is not supported, consider using adjoint"))
 end
 
-function lmul!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
+function lmul!(S::SplitCholesky{T,<:HermOrSym{T,M}}, B::AbstractVecOrMat{T}) where {T<:Real,M<:BandedMatrix{T}}
     require_one_based_indexing(B)
     n, nrhs = size(B, 1), size(B, 2)
     if size(S, 1) != n
@@ -84,7 +90,7 @@ function lmul!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where
     B
 end
 
-function ldiv!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
+function ldiv!(S::SplitCholesky{T,<:HermOrSym{T,M}}, B::AbstractVecOrMat{T}) where {T<:Real,M<:BandedMatrix{T}}
     require_one_based_indexing(B)
     n, nrhs = size(B, 1), size(B, 2)
     if size(S, 1) != n
@@ -112,7 +118,7 @@ function ldiv!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where
     B
 end
 
-function lmul!(S::AdjointFact{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
+function lmul!(S::AdjointFact{T,<:SplitCholesky{T,<:HermOrSym{T,M}}}, B::AbstractVecOrMat{T}) where {T<:Real,M<:BandedMatrix{T}}
     require_one_based_indexing(B)
     n, nrhs = size(B, 1), size(B, 2)
     if size(S, 1) != n
@@ -147,7 +153,7 @@ function lmul!(S::AdjointFact{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVec
     B
 end
 
-function ldiv!(S::AdjointFact{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
+function ldiv!(S::AdjointFact{T,<:SplitCholesky{T,<:HermOrSym{T,M}}}, B::AbstractVecOrMat{T}) where {T<:Real,M<:BandedMatrix{T}}
     require_one_based_indexing(B)
     n, nrhs = size(B, 1), size(B, 2)
     if size(S, 1) != n

src/symbanded/bandedeigen.jl

@@ -1,17 +1,43 @@
 ## eigvals routine
 
 # the symmetric case uses lapack throughout
-eigvals(A::Symmetric{T,<:BandedMatrix{T}}) where T<:Union{Float32, Float64} = eigvals!(copy(A))
-eigvals(A::Symmetric{T,<:BandedMatrix{T}}, irange::UnitRange) where T<:Union{Float32, Float64} = eigvals!(copy(A), irange)
-eigvals(A::Symmetric{T,<:BandedMatrix{T}}, vl::Real, vu::Real) where T<:Union{Float32, Float64} = eigvals!(copy(A), vl, vu)
-
-eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}) = eigvals!(tridiagonalize(A))
-eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, irange::UnitRange) = eigvals!(tridiagonalize(A), irange)
-eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, vl::Real, vu::Real) = eigvals!(tridiagonalize(A), vl, vu)
+eigvals(A::Symmetric{T,<:BandedMatrix{T}}) where T<:BlasReal =
+    eigvals!(copy(A))
+eigvals(A::Symmetric{T,<:BandedMatrix{T}}, irange::UnitRange) where T<:BlasReal =
+    eigvals!(copy(A), irange)
+eigvals(A::Symmetric{T,<:BandedMatrix{T}}, vl::Real, vu::Real) where T<:BlasReal =
+    eigvals!(copy(A), vl, vu)
+
+eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}) =
+    eigvals!(tridiagonalize(A))
+eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, irange::UnitRange) =
+    eigvals!(tridiagonalize(A), irange)
+eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, vl::Real, vu::Real) =
+    eigvals!(tridiagonalize(A), vl, vu)
+
+# This isn't eigvals!(A, args...) to avoid incorrect dispatches
+# This is a cautious approach at the moment
+eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}) = _eigvals!(A)
+eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, irange::UnitRange) = _eigvals!(A, irange)
+eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, vl::Real, vu::Real) = _eigvals!(A, vl, vu)
+
+_eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, args...) =
+    eigvals!(tridiagonalize!(A), args...)
+
+function _copy_bandedsym(A, B)
+    if bandwidth(A) >= bandwidth(B)
+        copy(A)
+    else
+        copyto!(similar(B), A)
+    end
+end
 
-eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, args...) = eigvals!(tridiagonalize!(A), args...)
+function eigvals(A::HermOrSym{<:Any,<:BandedMatrix}, B::HermOrSym{<:Any,<:BandedMatrix})
+    AA = _copy_bandedsym(A, B)
+    eigvals!(AA, copy(B))
+end
 
-function eigvals!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T<:Real
+function eigvals!(A::HermOrSym{T,<:BandedMatrix{T}}, B::HermOrSym{T,<:BandedMatrix{T}}) where T<:BlasReal
     n = size(A, 1)
     @assert n == size(B, 1)
     @assert A.uplo == B.uplo
@@ -29,6 +55,25 @@ function eigvals!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatr
     eigvals!(A)
 end
 
+function eigvals!(A::Hermitian{T,<:BandedMatrix{T}}, B::Hermitian{T,<:BandedMatrix{T}}) where T<:BlasComplex
+    n = size(A, 1)
+    @assert n == size(B, 1)
+    @assert A.uplo == B.uplo
+    # compute split-Cholesky factorization of B.
+    kb = bandwidth(B)
+    B_data = hermbandeddata(B)
+    pbstf!(B.uplo, n, kb, B_data)
+    # convert to a regular symmetric eigenvalue problem.
+    ka = bandwidth(A)
+    A_data = hermbandeddata(A)
+    X = Array{T}(undef,0,0)
+    work = Vector{T}(undef,n)
+    rwork = Vector{real(T)}(undef,n)
+    hbgst!('N', A.uplo, n, ka, kb, A_data, B_data, X, work, rwork)
+    # compute eigenvalues of symmetric eigenvalue problem.
+    eigvals!(A)
+end
+
 abstract type AbstractBandedEigenvectors{T} <: AbstractMatrix{T} end
 
 getindex(B::AbstractBandedEigenvectors, i, j) = Matrix(B)[i,j]
@@ -57,7 +102,6 @@ struct BandedEigenvectors{T} <: AbstractBandedEigenvectors{T}
 end
 
 size(B::BandedEigenvectors) = size(B.Q)
-
 _get_scratch(B::BandedEigenvectors) = B.z1
 
 # V = S⁻¹ Q W
@@ -80,12 +124,16 @@ eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}) = eigen!(copy(A))
 eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, irange::UnitRange) = eigen!(copy(A), irange)
 eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, vl::Real, vu::Real) = eigen!(copy(A), vl, vu)
 
-function eigen(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
+function eigen(A::HermOrSym{T,<:BandedMatrix{T}}, B::HermOrSym{T,<:BandedMatrix{T}}) where T
     AA = _copy_bandedsym(A, B)
     eigen!(AA, copy(B))
 end
 
-function eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, args...) where T <: Real
+eigen!(A::HermOrSym{T,<:BandedMatrix{T}}) where T<:BlasFloat = _eigen!(A)
+eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, irange::UnitRange) where T<:BlasFloat = _eigen!(A, irange)
+eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, vl::Real, vu::Real) where T<:BlasFloat = _eigen!(A, vl, vu)
+
+function _eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, args...) where T<:BlasReal
     N = size(A, 1)
     KD = bandwidth(A)
     D = Vector{T}(undef, N)
@@ -98,7 +146,7 @@ function eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, args...) where T <: Real
     Eigen(Λ, BandedEigenvectors(G, Q, similar(Q, size(Q,1))))
 end
 
-function eigen!(A::Hermitian{T,<:BandedMatrix{T}}, args...) where T <: Complex
+function _eigen!(A::Hermitian{T,<:BandedMatrix{T}}, args...) where T <: BlasComplex
     N = size(A, 1)
     KD = bandwidth(A)
     D = Vector{real(T)}(undef, N)
@@ -111,7 +159,7 @@ function eigen!(A::Hermitian{T,<:BandedMatrix{T}}, args...) where T <: Complex
     Eigen(Λ, Q * W)
 end
 
-function eigen!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
+function eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, B::HermOrSym{T,<:BandedMatrix{T}}) where T <: BlasReal
     isdiag(A) || isdiag(B) || symmetricuplo(A) == symmetricuplo(B) || throw(ArgumentError("uplo of matrices do not match"))
     S = splitcholesky!(B)
     N = size(A, 1)
@@ -127,6 +175,23 @@ function eigen!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix
     GeneralizedEigen(Λ, BandedGeneralizedEigenvectors(S, Q, W))
 end
 
+function eigen!(A::Hermitian{T,<:BandedMatrix{T}}, B::Hermitian{T,<:BandedMatrix{T}}) where T <: BlasComplex
+    isdiag(A) || isdiag(B) || symmetricuplo(A) == symmetricuplo(B) || throw(ArgumentError("uplo of matrices do not match"))
+    splitcholesky!(B)
+    N = size(A, 1)
+    KA = bandwidth(A)
+    KB = bandwidth(B)
+    X = Matrix{T}(undef, N, N)
+    WORK = Vector{T}(undef, N)
+    RWORK = Vector{real(T)}(undef, N)
+    AB = hermbandeddata(A)
+    BB = hermbandeddata(B)
+    hbgst!('V', A.uplo, N, KA, KB, AB, BB, X, WORK, RWORK)
+    any(isnan, A) && throw(ArgumentError("NaN found in the standard form of A"))
+    Λ, W = eigen!(A)
+    GeneralizedEigen(Λ, X * W)
+end
+
 function Matrix(B::BandedEigenvectors)
     Q = copy(B.Q)
     G = B.G

src/symbanded/symbanded.jl

@@ -114,6 +114,8 @@ function materialize!(M::BlasMatMulVecAdd{<:HermitianLayout{<:BandedColumnMajor}
     _banded_hbmv!(symmetricuplo(A), α, A, x, β, y)
 end
 
+## eigvals routine
+
 function copyto!(A::Symmetric{<:Number,<:BandedMatrix}, B::Symmetric{<:Number,<:BandedMatrix})
     size(A) == size(B) || throw(ArgumentError("sizes of A and B must match"))
     bandwidth(A) >= bandwidth(B) || throw(ArgumentError("bandwidth of A must exceed that of B"))
@@ -140,16 +142,3 @@ function copyto!(A::Symmetric{<:Number,<:BandedMatrix}, B::Symmetric{<:Number,<:
     end
     return A
 end
-
-function _copy_bandedsym(A, B)
-    if bandwidth(A) >= bandwidth(B)
-        copy(A)
-    else
-        copyto!(similar(B), A)
-    end
-end
-
-function eigvals(A::Symmetric{<:Any,<:BandedMatrix}, B::Symmetric{<:Any,<:BandedMatrix})
-    AA = _copy_bandedsym(A, B)
-    eigvals!(AA, copy(B))
-end

src/symbanded/tridiagonalize.jl

@@ -84,7 +84,7 @@ function copybands(A::AbstractMatrix{T}, d::Integer) where T
     B
 end
 
-function _tridiagonalize!(A::AbstractMatrix{T}, ::SymmetricLayout{<:BandedColumnMajor}) where T<:Union{Float32,Float64}
+function _tridiagonalize!(A::AbstractMatrix{T}, ::SymmetricLayout{<:BandedColumnMajor}) where T<:BlasReal
     n=size(A, 1)
     d = Vector{T}(undef,n)
     e = Vector{T}(undef,n-1)
@@ -94,7 +94,7 @@ function _tridiagonalize!(A::AbstractMatrix{T}, ::SymmetricLayout{<:BandedColumn
     SymTridiagonal(d,e)
 end
 
-function _tridiagonalize!(A::AbstractMatrix{T}, ::HermitianLayout{<:BandedColumnMajor}) where T<:Union{ComplexF32,ComplexF64}
+function _tridiagonalize!(A::AbstractMatrix{T}, ::HermitianLayout{<:BandedColumnMajor}) where T<:BlasComplex
     n=size(A, 1)
     d = Vector{real(T)}(undef,n)
     e = Vector{real(T)}(undef,n-1)
@@ -104,4 +104,4 @@ function _tridiagonalize!(A::AbstractMatrix{T}, ::HermitianLayout{<:BandedColumn
     SymTridiagonal(d,e)
 end
 
-tridiagonalize!(A::AbstractMatrix) = _tridiagonalize!(A, MemoryLayout(typeof(A)))
+tridiagonalize!(A::AbstractMatrix{<:BlasFloat}) = _tridiagonalize!(A, MemoryLayout(typeof(A)))

test/test_symbanded.jl

@@ -313,10 +313,10 @@ end
 @testset "Generalized eigenvalues $W{$T}($Ua,$Ub)($n,$wa-$wb)" for (T,W) in (
                                     (Float32, Symmetric),
                                     (Float64, Symmetric),
-                                    #(Float32, Hermitian),
-                                    #(Float64, Hermitian),
-                                    #(ComplexF32, Hermitian),
-                                    #(ComplexF64, Hermitian),
+                                    (Float32, Hermitian),
+                                    (Float64, Hermitian),
+                                    (ComplexF32, Hermitian),
+                                    (ComplexF64, Hermitian),
                                    ),
     (Ua, Ub) in  ((:L,:L), (:U,:U)),
     (wa, wb) in ((2, 3), (3, 2)), n in (4,)
@@ -328,7 +328,16 @@ end
     end
     A = sbmatrix(W, T, Ua, wa, n)
     B = sbmatrix(W, T, Ub, wb, n)
-    @test eigvals(A, B) ≈ eigvals(Matrix(A), Matrix(B))
+    AM, BM = Matrix.((A,B))
+    @test eigvals(A, B) ≈ eigvals(AM, BM)
+    if VERSION >= v"1.9"
+        Λ, V = eigen(A, B)
+        VM = Matrix(V)
+        Λ2, V2 = eigen(AM, BM)
+        @test Λ ≈ Λ2
+        @test VM' * AM * VM ≈ V2' * AM * V2
+        @test VM' * AM * VM ≈ VM' * BM * VM * Diagonal(Λ)
+    end
 end
 
 @testset "Cholesky" begin