Skip to content

Commit 633e5bd

Browse files
jishnubjishnub
authored
Generalized eigen for complex Hermitian BandedMatrix (#360)
* in-place eigvals for complex hermitian BandedMatrix * implement hbgst * generalized complex hermitian eigen * generalized eigvals and tests * skip kwargs for now * test for symtridiag eigen on v1.9+ * fix eigen dispatch * limit BLAS eigvals to BlasReal/BlasComplex * Make method signatures more specific
1 parent 43764a2 commit 633e5bd

File tree

6 files changed

+188
-46
lines changed

6 files changed

+188
-46
lines changed

‎src/lapack.jl‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -371,3 +371,76 @@ for (fname, elty) in ((:dsbgst_,:Float64),
371371
end
372372
end
373373
end
374+
375+
# Convert a complex Hermitian Positive Definite generalized eigenvalue problem
376+
# to a symmetric eigenvalue problem assuming B has been processed by
377+
# a split-Cholesky factorization.
378+
for (fname, elty) in ((:zhbgst_, :ComplexF64),
379+
(:chbgst_, :ComplexF32))
380+
@eval begin
381+
#=
382+
subroutine zhbgst (
383+
character VECT,
384+
character UPLO,
385+
integer N,
386+
integer KA,
387+
integer KB,
388+
complex*16, dimension( ldab, * ) AB,
389+
integer LDAB,
390+
complex*16, dimension( ldbb, * ) BB,
391+
integer LDBB,
392+
complex*16, dimension( ldx, * ) X,
393+
integer LDX,
394+
complex*16, dimension( * ) WORK,
395+
double precision, dimension( * ) RWORK,
396+
integer INFO
397+
)
398+
=#
399+
local Relty = real($elty)
400+
function hbgst!(vect::Char, uplo::Char, n::Int, ka::Int, kb::Int,
401+
AB::StridedMatrix{$elty}, BB::StridedMatrix{$elty},
402+
X::StridedMatrix{$elty}, work::StridedVector{$elty},
403+
rwork::StridedVector{Relty})
404+
require_one_based_indexing(AB, BB, X, work)
405+
chkstride1(AB, BB)
406+
chkuplo(uplo)
407+
chkvect(vect)
408+
info = Ref{BlasInt}()
409+
ccall((@blasfunc($fname), liblapack), Nothing,
410+
(
411+
Ref{UInt8},
412+
Ref{UInt8},
413+
Ref{BlasInt},
414+
Ref{BlasInt},
415+
Ref{BlasInt},
416+
Ptr{$elty},
417+
Ref{BlasInt},
418+
Ptr{$elty},
419+
Ref{BlasInt},
420+
Ptr{$elty},
421+
Ref{BlasInt},
422+
Ptr{$elty},
423+
Ptr{Relty},
424+
Ref{BlasInt},
425+
),
426+
vect,
427+
uplo,
428+
n,
429+
ka,
430+
kb,
431+
AB,
432+
max(stride(AB,2),1),
433+
BB,
434+
max(stride(BB,2),1),
435+
X,
436+
max(stride(X,2),1),
437+
work,
438+
rwork,
439+
info,
440+
)
441+
442+
LAPACK.chklapackerror(info[])
443+
AB
444+
end
445+
end
446+
end

‎src/symbanded/SplitCholesky.jl‎

Lines changed: 17 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -33,18 +33,24 @@ SplitCholesky(factors::AbstractMatrix{T}, uplo::Char) where T = SplitCholesky{T,
3333
size(S::SplitCholesky) = size(S.factors)
3434
size(S::SplitCholesky, i::Integer) = size(S.factors, i)
3535

36-
splitcholesky!(A::Symmetric{T,<:BandedMatrix{T}}) where T = splitcholesky!(A, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
37-
splitcholesky!(A::Symmetric{T,<:BandedMatrix{T}}, ::Type{Tr}) where {T,Tr} = splitcholesky!(MemoryLayout(typeof(A)), A, Tr)
36+
function splitcholesky!(A::HermOrSym{T,<:BandedMatrix{T}}) where T
37+
splitcholesky!(A, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
38+
end
39+
function splitcholesky!(A::HermOrSym{T,<:BandedMatrix{T}}, Tr) where {T}
40+
splitcholesky!(MemoryLayout(typeof(A)), A, Tr)
41+
end
3842

3943
function splitcholesky!(::SymmetricLayout{<:BandedColumnMajor},
40-
A::AbstractMatrix{T}, ::Type{UpperTriangular}) where T<:BlasFloat
41-
_, info = pbstf!('U', size(A, 1), bandwidth(A,2), symbandeddata(A))
44+
A::AbstractMatrix{T}, ::Type{LU}) where {T<:BlasFloat, LU}
45+
uplo = LU == UpperTriangular ? 'U' : 'L'
46+
pbstf!(uplo, size(A, 1), bandwidth(A,2), symbandeddata(A))
4247
SplitCholesky(A, A.uplo)
4348
end
4449

45-
function splitcholesky!(::SymmetricLayout{<:BandedColumnMajor},
46-
A::AbstractMatrix{T}, ::Type{LowerTriangular}) where T<:BlasFloat
47-
_, info = pbstf!('L', size(A, 1), bandwidth(A,1), symbandeddata(A))
50+
function splitcholesky!(::HermitianLayout{<:BandedColumnMajor},
51+
A::AbstractMatrix{T}, ::Type{LU}) where {T<:BlasFloat, LU}
52+
uplo = LU == UpperTriangular ? 'U' : 'L'
53+
pbstf!(uplo, size(A, 1), bandwidth(A,2), hermbandeddata(A))
4854
SplitCholesky(A, A.uplo)
4955
end
5056

@@ -56,7 +62,7 @@ else
5662
throw(ArgumentError("transpose of SplitCholesky decomposition is not supported, consider using adjoint"))
5763
end
5864

59-
function lmul!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
65+
function lmul!(S::SplitCholesky{T,<:HermOrSym{T,M}}, B::AbstractVecOrMat{T}) where {T<:Real,M<:BandedMatrix{T}}
6066
require_one_based_indexing(B)
6167
n, nrhs = size(B, 1), size(B, 2)
6268
if size(S, 1) != n
@@ -84,7 +90,7 @@ function lmul!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where
8490
B
8591
end
8692

87-
function ldiv!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
93+
function ldiv!(S::SplitCholesky{T,<:HermOrSym{T,M}}, B::AbstractVecOrMat{T}) where {T<:Real,M<:BandedMatrix{T}}
8894
require_one_based_indexing(B)
8995
n, nrhs = size(B, 1), size(B, 2)
9096
if size(S, 1) != n
@@ -112,7 +118,7 @@ function ldiv!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where
112118
B
113119
end
114120

115-
function lmul!(S::AdjointFact{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
121+
function lmul!(S::AdjointFact{T,<:SplitCholesky{T,<:HermOrSym{T,M}}}, B::AbstractVecOrMat{T}) where {T<:Real,M<:BandedMatrix{T}}
116122
require_one_based_indexing(B)
117123
n, nrhs = size(B, 1), size(B, 2)
118124
if size(S, 1) != n
@@ -147,7 +153,7 @@ function lmul!(S::AdjointFact{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVec
147153
B
148154
end
149155

150-
function ldiv!(S::AdjointFact{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}}
156+
function ldiv!(S::AdjointFact{T,<:SplitCholesky{T,<:HermOrSym{T,M}}}, B::AbstractVecOrMat{T}) where {T<:Real,M<:BandedMatrix{T}}
151157
require_one_based_indexing(B)
152158
n, nrhs = size(B, 1), size(B, 2)
153159
if size(S, 1) != n

‎src/symbanded/bandedeigen.jl‎

Lines changed: 79 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,43 @@
11
## eigvals routine
22

33
# the symmetric case uses lapack throughout
4-
eigvals(A::Symmetric{T,<:BandedMatrix{T}}) where T<:Union{Float32, Float64} = eigvals!(copy(A))
5-
eigvals(A::Symmetric{T,<:BandedMatrix{T}}, irange::UnitRange) where T<:Union{Float32, Float64} = eigvals!(copy(A), irange)
6-
eigvals(A::Symmetric{T,<:BandedMatrix{T}}, vl::Real, vu::Real) where T<:Union{Float32, Float64} = eigvals!(copy(A), vl, vu)
7-
8-
eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}) = eigvals!(tridiagonalize(A))
9-
eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, irange::UnitRange) = eigvals!(tridiagonalize(A), irange)
10-
eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, vl::Real, vu::Real) = eigvals!(tridiagonalize(A), vl, vu)
4+
eigvals(A::Symmetric{T,<:BandedMatrix{T}}) where T<:BlasReal =
5+
eigvals!(copy(A))
6+
eigvals(A::Symmetric{T,<:BandedMatrix{T}}, irange::UnitRange) where T<:BlasReal =
7+
eigvals!(copy(A), irange)
8+
eigvals(A::Symmetric{T,<:BandedMatrix{T}}, vl::Real, vu::Real) where T<:BlasReal =
9+
eigvals!(copy(A), vl, vu)
10+
11+
eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}) =
12+
eigvals!(tridiagonalize(A))
13+
eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, irange::UnitRange) =
14+
eigvals!(tridiagonalize(A), irange)
15+
eigvals(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, vl::Real, vu::Real) =
16+
eigvals!(tridiagonalize(A), vl, vu)
17+
18+
# This isn't eigvals!(A, args...) to avoid incorrect dispatches
19+
# This is a cautious approach at the moment
20+
eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}) = _eigvals!(A)
21+
eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, irange::UnitRange) = _eigvals!(A, irange)
22+
eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, vl::Real, vu::Real) = _eigvals!(A, vl, vu)
23+
24+
_eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, args...) =
25+
eigvals!(tridiagonalize!(A), args...)
26+
27+
function _copy_bandedsym(A, B)
28+
if bandwidth(A) >= bandwidth(B)
29+
copy(A)
30+
else
31+
copyto!(similar(B), A)
32+
end
33+
end
1134

12-
eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, args...) = eigvals!(tridiagonalize!(A), args...)
35+
function eigvals(A::HermOrSym{<:Any,<:BandedMatrix}, B::HermOrSym{<:Any,<:BandedMatrix})
36+
AA = _copy_bandedsym(A, B)
37+
eigvals!(AA, copy(B))
38+
end
1339

14-
function eigvals!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T<:Real
40+
function eigvals!(A::HermOrSym{T,<:BandedMatrix{T}}, B::HermOrSym{T,<:BandedMatrix{T}}) where T<:BlasReal
1541
n = size(A, 1)
1642
@assert n == size(B, 1)
1743
@assert A.uplo == B.uplo
@@ -29,6 +55,25 @@ function eigvals!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatr
2955
eigvals!(A)
3056
end
3157

58+
function eigvals!(A::Hermitian{T,<:BandedMatrix{T}}, B::Hermitian{T,<:BandedMatrix{T}}) where T<:BlasComplex
59+
n = size(A, 1)
60+
@assert n == size(B, 1)
61+
@assert A.uplo == B.uplo
62+
# compute split-Cholesky factorization of B.
63+
kb = bandwidth(B)
64+
B_data = hermbandeddata(B)
65+
pbstf!(B.uplo, n, kb, B_data)
66+
# convert to a regular symmetric eigenvalue problem.
67+
ka = bandwidth(A)
68+
A_data = hermbandeddata(A)
69+
X = Array{T}(undef,0,0)
70+
work = Vector{T}(undef,n)
71+
rwork = Vector{real(T)}(undef,n)
72+
hbgst!('N', A.uplo, n, ka, kb, A_data, B_data, X, work, rwork)
73+
# compute eigenvalues of symmetric eigenvalue problem.
74+
eigvals!(A)
75+
end
76+
3277
abstract type AbstractBandedEigenvectors{T} <: AbstractMatrix{T} end
3378

3479
getindex(B::AbstractBandedEigenvectors, i, j) = Matrix(B)[i,j]
@@ -57,7 +102,6 @@ struct BandedEigenvectors{T} <: AbstractBandedEigenvectors{T}
57102
end
58103

59104
size(B::BandedEigenvectors) = size(B.Q)
60-
61105
_get_scratch(B::BandedEigenvectors) = B.z1
62106

63107
# V = S⁻¹ Q W
@@ -80,12 +124,16 @@ eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}) = eigen!(copy(A))
80124
eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, irange::UnitRange) = eigen!(copy(A), irange)
81125
eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, vl::Real, vu::Real) = eigen!(copy(A), vl, vu)
82126

83-
function eigen(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
127+
function eigen(A::HermOrSym{T,<:BandedMatrix{T}}, B::HermOrSym{T,<:BandedMatrix{T}}) where T
84128
AA = _copy_bandedsym(A, B)
85129
eigen!(AA, copy(B))
86130
end
87131

88-
function eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, args...) where T <: Real
132+
eigen!(A::HermOrSym{T,<:BandedMatrix{T}}) where T<:BlasFloat = _eigen!(A)
133+
eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, irange::UnitRange) where T<:BlasFloat = _eigen!(A, irange)
134+
eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, vl::Real, vu::Real) where T<:BlasFloat = _eigen!(A, vl, vu)
135+
136+
function _eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, args...) where T<:BlasReal
89137
N = size(A, 1)
90138
KD = bandwidth(A)
91139
D = Vector{T}(undef, N)
@@ -98,7 +146,7 @@ function eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, args...) where T <: Real
98146
Eigen(Λ, BandedEigenvectors(G, Q, similar(Q, size(Q,1))))
99147
end
100148

101-
function eigen!(A::Hermitian{T,<:BandedMatrix{T}}, args...) where T <: Complex
149+
function _eigen!(A::Hermitian{T,<:BandedMatrix{T}}, args...) where T <: BlasComplex
102150
N = size(A, 1)
103151
KD = bandwidth(A)
104152
D = Vector{real(T)}(undef, N)
@@ -111,7 +159,7 @@ function eigen!(A::Hermitian{T,<:BandedMatrix{T}}, args...) where T <: Complex
111159
Eigen(Λ, Q * W)
112160
end
113161

114-
function eigen!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
162+
function eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, B::HermOrSym{T,<:BandedMatrix{T}}) where T <: BlasReal
115163
isdiag(A) || isdiag(B) || symmetricuplo(A) == symmetricuplo(B) || throw(ArgumentError("uplo of matrices do not match"))
116164
S = splitcholesky!(B)
117165
N = size(A, 1)
@@ -127,6 +175,23 @@ function eigen!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix
127175
GeneralizedEigen(Λ, BandedGeneralizedEigenvectors(S, Q, W))
128176
end
129177

178+
function eigen!(A::Hermitian{T,<:BandedMatrix{T}}, B::Hermitian{T,<:BandedMatrix{T}}) where T <: BlasComplex
179+
isdiag(A) || isdiag(B) || symmetricuplo(A) == symmetricuplo(B) || throw(ArgumentError("uplo of matrices do not match"))
180+
splitcholesky!(B)
181+
N = size(A, 1)
182+
KA = bandwidth(A)
183+
KB = bandwidth(B)
184+
X = Matrix{T}(undef, N, N)
185+
WORK = Vector{T}(undef, N)
186+
RWORK = Vector{real(T)}(undef, N)
187+
AB = hermbandeddata(A)
188+
BB = hermbandeddata(B)
189+
hbgst!('V', A.uplo, N, KA, KB, AB, BB, X, WORK, RWORK)
190+
any(isnan, A) && throw(ArgumentError("NaN found in the standard form of A"))
191+
Λ, W = eigen!(A)
192+
GeneralizedEigen(Λ, X * W)
193+
end
194+
130195
function Matrix(B::BandedEigenvectors)
131196
Q = copy(B.Q)
132197
G = B.G

‎src/symbanded/symbanded.jl‎

Lines changed: 2 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -114,6 +114,8 @@ function materialize!(M::BlasMatMulVecAdd{<:HermitianLayout{<:BandedColumnMajor}
114114
_banded_hbmv!(symmetricuplo(A), α, A, x, β, y)
115115
end
116116

117+
## eigvals routine
118+
117119
function copyto!(A::Symmetric{<:Number,<:BandedMatrix}, B::Symmetric{<:Number,<:BandedMatrix})
118120
size(A) == size(B) || throw(ArgumentError("sizes of A and B must match"))
119121
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,<:
140142
end
141143
return A
142144
end
143-
144-
function _copy_bandedsym(A, B)
145-
if bandwidth(A) >= bandwidth(B)
146-
copy(A)
147-
else
148-
copyto!(similar(B), A)
149-
end
150-
end
151-
152-
function eigvals(A::Symmetric{<:Any,<:BandedMatrix}, B::Symmetric{<:Any,<:BandedMatrix})
153-
AA = _copy_bandedsym(A, B)
154-
eigvals!(AA, copy(B))
155-
end

‎src/symbanded/tridiagonalize.jl‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -84,7 +84,7 @@ function copybands(A::AbstractMatrix{T}, d::Integer) where T
8484
B
8585
end
8686

87-
function _tridiagonalize!(A::AbstractMatrix{T}, ::SymmetricLayout{<:BandedColumnMajor}) where T<:Union{Float32,Float64}
87+
function _tridiagonalize!(A::AbstractMatrix{T}, ::SymmetricLayout{<:BandedColumnMajor}) where T<:BlasReal
8888
n=size(A, 1)
8989
d = Vector{T}(undef,n)
9090
e = Vector{T}(undef,n-1)
@@ -94,7 +94,7 @@ function _tridiagonalize!(A::AbstractMatrix{T}, ::SymmetricLayout{<:BandedColumn
9494
SymTridiagonal(d,e)
9595
end
9696

97-
function _tridiagonalize!(A::AbstractMatrix{T}, ::HermitianLayout{<:BandedColumnMajor}) where T<:Union{ComplexF32,ComplexF64}
97+
function _tridiagonalize!(A::AbstractMatrix{T}, ::HermitianLayout{<:BandedColumnMajor}) where T<:BlasComplex
9898
n=size(A, 1)
9999
d = Vector{real(T)}(undef,n)
100100
e = Vector{real(T)}(undef,n-1)
@@ -104,4 +104,4 @@ function _tridiagonalize!(A::AbstractMatrix{T}, ::HermitianLayout{<:BandedColumn
104104
SymTridiagonal(d,e)
105105
end
106106

107-
tridiagonalize!(A::AbstractMatrix) = _tridiagonalize!(A, MemoryLayout(typeof(A)))
107+
tridiagonalize!(A::AbstractMatrix{<:BlasFloat}) = _tridiagonalize!(A, MemoryLayout(typeof(A)))

‎test/test_symbanded.jl‎

Lines changed: 14 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -313,10 +313,10 @@ end
313313
@testset "Generalized eigenvalues $W{$T}($Ua,$Ub)($n,$wa-$wb)" for (T,W) in (
314314
(Float32, Symmetric),
315315
(Float64, Symmetric),
316-
#(Float32, Hermitian),
317-
#(Float64, Hermitian),
318-
#(ComplexF32, Hermitian),
319-
#(ComplexF64, Hermitian),
316+
(Float32, Hermitian),
317+
(Float64, Hermitian),
318+
(ComplexF32, Hermitian),
319+
(ComplexF64, Hermitian),
320320
),
321321
(Ua, Ub) in ((:L,:L), (:U,:U)),
322322
(wa, wb) in ((2, 3), (3, 2)), n in (4,)
@@ -328,7 +328,16 @@ end
328328
end
329329
A = sbmatrix(W, T, Ua, wa, n)
330330
B = sbmatrix(W, T, Ub, wb, n)
331-
@test eigvals(A, B) ≈ eigvals(Matrix(A), Matrix(B))
331+
AM, BM = Matrix.((A,B))
332+
@test eigvals(A, B) ≈ eigvals(AM, BM)
333+
if VERSION >= v"1.9"
334+
Λ, V = eigen(A, B)
335+
VM = Matrix(V)
336+
Λ2, V2 = eigen(AM, BM)
337+
@test Λ ≈ Λ2
338+
@test VM' * AM * VM ≈ V2' * AM * V2
339+
@test VM' * AM * VM ≈ VM' * BM * VM * Diagonal(Λ)
340+
end
332341
end
333342

334343
@testset "Cholesky" begin

0 commit comments

Comments
 (0)