Skip to content

Commit 88a4a22

Browse files
jishnubjishnub
authored
In-place eigvals for complex Hermitian BandedMatrix (#359)
* in-place eigvals for complex hermitian BandedMatrix * bump version to 0.17.24 * fix pointer types * complex hermitian inplace eigen * accept ranges in eigen * reorder tests
1 parent f922182 commit 88a4a22

File tree

6 files changed

+186
-55
lines changed

6 files changed

+186
-55
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "BandedMatrices"
22
uuid = "aae01518-5342-5314-be14-df237901396f"
3-
version = "0.17.23"
3+
version = "0.17.24"
44

55
[deps]
66
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

‎src/lapack.jl‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -115,6 +115,57 @@ for (fname, elty) in ((:dsbtrd_,:Float64),
115115
end
116116
end
117117

118+
for (fname, elty) in ((:zhbtrd_,:ComplexF64),
119+
(:chbtrd_,:ComplexF32))
120+
@eval begin
121+
local Relty = real($elty)
122+
function hbtrd!(vect::Char, uplo::Char,
123+
m::Int, k::Int, A::AbstractMatrix{$elty},
124+
d::AbstractVector{Relty}, e::AbstractVector{Relty}, Q::AbstractMatrix{$elty},
125+
work::AbstractVector{$elty})
126+
require_one_based_indexing(A)
127+
chkstride1(A)
128+
chkuplo(uplo)
129+
chkvect(vect)
130+
info = Ref{BlasInt}()
131+
n = size(A,2)
132+
n ≠ m && throw(ArgumentError("Matrix must be square"))
133+
size(A,1) < k+1 && throw(ArgumentError("Not enough bands"))
134+
info = Ref{BlasInt}()
135+
ccall((@blasfunc($fname), liblapack), Nothing,
136+
(
137+
Ref{UInt8},
138+
Ref{UInt8},
139+
Ref{BlasInt},
140+
Ref{BlasInt},
141+
Ptr{$elty},
142+
Ref{BlasInt},
143+
Ptr{Relty},
144+
Ptr{Relty},
145+
Ptr{$elty},
146+
Ref{BlasInt},
147+
Ptr{$elty},
148+
Ptr{BlasInt},
149+
),
150+
vect,
151+
uplo,
152+
n,
153+
k,
154+
A,
155+
max(1,stride(A,2)),
156+
d,
157+
e,
158+
Q,
159+
max(1,stride(Q,2)),
160+
work,
161+
info
162+
)
163+
LAPACK.chklapackerror(info[])
164+
d, e, Q
165+
end
166+
end
167+
end
168+
118169
## Bidiagonalization
119170
for (fname, elty) in ((:dgbbrd_,:Float64),
120171
(:sgbbrd_,:Float32))

‎src/symbanded/bandedeigen.jl‎

Lines changed: 50 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,34 @@
1+
## eigvals routine
2+
3+
# 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)
11+
12+
eigvals!(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, args...) = eigvals!(tridiagonalize!(A), args...)
13+
14+
function eigvals!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T<:Real
15+
n = size(A, 1)
16+
@assert n == size(B, 1)
17+
@assert A.uplo == B.uplo
18+
# compute split-Cholesky factorization of B.
19+
kb = bandwidth(B)
20+
B_data = symbandeddata(B)
21+
pbstf!(B.uplo, n, kb, B_data)
22+
# convert to a regular symmetric eigenvalue problem.
23+
ka = bandwidth(A)
24+
A_data = symbandeddata(A)
25+
X = Array{T}(undef,0,0)
26+
work = Vector{T}(undef,2n)
27+
sbgst!('N', A.uplo, n, ka, kb, A_data, B_data, X, work)
28+
# compute eigenvalues of symmetric eigenvalue problem.
29+
eigvals!(A)
30+
end
31+
132
# V = G Q
233
struct BandedEigenvectors{T} <: AbstractMatrix{T}
334
G::Vector{Givens{T}}
@@ -44,13 +75,16 @@ convert(::Type{GeneralizedEigen{T, T, Matrix{T}, Vector{T}}}, F::GeneralizedEige
4475
compress(F::Eigen{T, T, BandedEigenvectors{T}, Vector{T}}) where T = convert(Eigen{T, T, Matrix{T}, Vector{T}}, F)
4576
compress(F::GeneralizedEigen{T, T, BandedGeneralizedEigenvectors{T}, Vector{T}}) where T = convert(GeneralizedEigen{T, T, Matrix{T}, Vector{T}}, F)
4677

47-
eigen(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real = eigen!(copy(A))
78+
eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}) = eigen!(copy(A))
79+
eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, irange::UnitRange) = eigen!(copy(A), irange)
80+
eigen(A::RealHermSymComplexHerm{<:Real,<:BandedMatrix}, vl::Real, vu::Real) = eigen!(copy(A), vl, vu)
81+
4882
function eigen(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
4983
AA = _copy_bandedsym(A, B)
5084
eigen!(AA, copy(B))
5185
end
5286

53-
function eigen!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
87+
function eigen!(A::HermOrSym{T,<:BandedMatrix{T}}, args...) where T <: Real
5488
N = size(A, 1)
5589
KD = bandwidth(A)
5690
D = Vector{T}(undef, N)
@@ -59,10 +93,23 @@ function eigen!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
5993
WORK = Vector{T}(undef, N)
6094
AB = symbandeddata(A)
6195
sbtrd!('V', A.uplo, N, KD, AB, D, E, G, WORK)
62-
Λ, Q = eigen(SymTridiagonal(D, E))
96+
Λ, Q = eigen!(SymTridiagonal(D, E), args...)
6397
Eigen(Λ, BandedEigenvectors(G, Q, similar(Q, size(Q,1))))
6498
end
6599

100+
function eigen!(A::Hermitian{T,<:BandedMatrix{T}}, args...) where T <: Complex
101+
N = size(A, 1)
102+
KD = bandwidth(A)
103+
D = Vector{real(T)}(undef, N)
104+
E = Vector{real(T)}(undef, N-1)
105+
Q = Matrix{T}(undef, N, N)
106+
WORK = Vector{T}(undef, N)
107+
AB = hermbandeddata(A)
108+
hbtrd!('V', A.uplo, N, KD, AB, D, E, Q, WORK)
109+
Λ, W = eigen!(SymTridiagonal(D, E), args...)
110+
Eigen(Λ, Q * W)
111+
end
112+
66113
function eigen!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T <: Real
67114
isdiag(A) || isdiag(B) || symmetricuplo(A) == symmetricuplo(B) || throw(ArgumentError("uplo of matrices do not match"))
68115
S = splitcholesky!(B)

‎src/symbanded/symbanded.jl‎

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

117-
118-
## eigvals routine
119-
120-
121-
function _tridiagonalize!(A::AbstractMatrix{T}, ::SymmetricLayout{<:BandedColumnMajor}) where T
122-
n=size(A, 1)
123-
d = Vector{T}(undef,n)
124-
e = Vector{T}(undef,n-1)
125-
Q = Matrix{T}(undef,0,0)
126-
work = Vector{T}(undef,n)
127-
sbtrd!('N', symmetricuplo(A), size(A,1), bandwidth(A), symbandeddata(A), d, e, Q, work)
128-
SymTridiagonal(d,e)
129-
end
130-
131-
tridiagonalize!(A::AbstractMatrix) = _tridiagonalize!(A, MemoryLayout(typeof(A)))
132-
133-
eigvals(A::Symmetric{T,<:BandedMatrix{T}}) where T<:Base.IEEEFloat = eigvals!(copy(A))
134-
eigvals(A::Symmetric{T,<:BandedMatrix{T}}) where T<:Real = eigvals!(tridiagonalize(A))
135-
eigvals(A::Hermitian{T,<:BandedMatrix{T}}) where T<:Real = eigvals!(tridiagonalize(A))
136-
eigvals(A::Hermitian{T,<:BandedMatrix{T}}) where T<:Complex = eigvals!(tridiagonalize(A))
137-
138-
eigvals!(A::Symmetric{T,<:BandedMatrix{T}}) where T <: Real = eigvals!(tridiagonalize!(A))
139-
function eigvals!(A::Symmetric{T,<:BandedMatrix{T}}, B::Symmetric{T,<:BandedMatrix{T}}) where T<:Real
140-
n = size(A, 1)
141-
@assert n == size(B, 1)
142-
@assert A.uplo == B.uplo
143-
# compute split-Cholesky factorization of B.
144-
kb = bandwidth(B)
145-
B_data = symbandeddata(B)
146-
pbstf!(B.uplo, n, kb, B_data)
147-
# convert to a regular symmetric eigenvalue problem.
148-
ka = bandwidth(A)
149-
A_data = symbandeddata(A)
150-
X = Array{T}(undef,0,0)
151-
work = Vector{T}(undef,2n)
152-
sbgst!('N', A.uplo, n, ka, kb, A_data, B_data, X, work)
153-
# compute eigenvalues of symmetric eigenvalue problem.
154-
eigvals!(A)
155-
end
156-
157117
function copyto!(A::Symmetric{<:Number,<:BandedMatrix}, B::Symmetric{<:Number,<:BandedMatrix})
158118
size(A) == size(B) || throw(ArgumentError("sizes of A and B must match"))
159119
bandwidth(A) >= bandwidth(B) || throw(ArgumentError("bandwidth of A must exceed that of B"))

‎src/symbanded/tridiagonalize.jl‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -84,3 +84,24 @@ 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}
88+
n=size(A, 1)
89+
d = Vector{T}(undef,n)
90+
e = Vector{T}(undef,n-1)
91+
Q = Matrix{T}(undef,0,0)
92+
work = Vector{T}(undef,n)
93+
sbtrd!('N', symmetricuplo(A), size(A,1), bandwidth(A), symbandeddata(A), d, e, Q, work)
94+
SymTridiagonal(d,e)
95+
end
96+
97+
function _tridiagonalize!(A::AbstractMatrix{T}, ::HermitianLayout{<:BandedColumnMajor}) where T<:Union{ComplexF32,ComplexF64}
98+
n=size(A, 1)
99+
d = Vector{real(T)}(undef,n)
100+
e = Vector{real(T)}(undef,n-1)
101+
Q = Matrix{T}(undef,0,0)
102+
work = Vector{T}(undef,n)
103+
hbtrd!('N', symmetricuplo(A), size(A,1), bandwidth(A), hermbandeddata(A), d, e, Q, work)
104+
SymTridiagonal(d,e)
105+
end
106+
107+
tridiagonalize!(A::AbstractMatrix) = _tridiagonalize!(A, MemoryLayout(typeof(A)))

‎test/test_symbanded.jl‎

Lines changed: 63 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -47,13 +47,20 @@ import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedCol
4747
# (generalized) eigen & eigvals
4848
Random.seed!(0)
4949

50-
A = brand(Float64, 100, 100, 2, 4)
51-
std = eigvals(Symmetric(Matrix(A)))
52-
@test eigvals(Symmetric(A)) ≈ std
53-
@test eigvals(Hermitian(A)) ≈ std
54-
@test eigvals(Hermitian(big.(A))) ≈ std
50+
@testset for T in (Float32, Float64)
51+
A = brand(T, 10, 10, 2, 4)
52+
std = eigvals(Symmetric(Matrix(A)))
53+
@test eigvals(Symmetric(A)) ≈ std
54+
@test eigvals(Symmetric(A), 2:4) ≈ std[2:4]
55+
@test eigvals(Symmetric(A), 1, 2) ≈ std[1 .<= std .<= 2]
56+
@test eigvals!(copy(Symmetric(A))) ≈ std
57+
@test eigvals!(copy(Symmetric(A)), 2:4) ≈ std[2:4]
58+
@test eigvals!(copy(Symmetric(A)), 1, 2) ≈ std[1 .<= std .<= 2]
59+
60+
@test eigvals(Symmetric(A, :L)) ≈ eigvals(Symmetric(Matrix(A), :L))
61+
end
5562

56-
A = brand(ComplexF64, 100, 100, 4, 0)
63+
A = brand(ComplexF64, 20, 20, 4, 0)
5764
@test Symmetric(A)[2:10,1:9] isa BandedMatrix
5865
@test Hermitian(A)[2:10,1:9] isa BandedMatrix
5966
@test isempty(Hermitian(A)[1:0,1:9])
@@ -62,11 +69,7 @@ import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedCol
6269
@test [Symmetric(A)[k,j] for k=2:10, j=1:9] == Symmetric(A)[2:10,1:9]
6370
@test [Hermitian(A)[k,j] for k=2:10, j=1:9] == Hermitian(A)[2:10,1:9]
6471

65-
std = eigvals(Hermitian(Matrix(A), :L))
66-
@test eigvals(Hermitian(A, :L)) ≈ std
67-
@test eigvals(Hermitian(big.(A), :L)) ≈ std
68-
69-
A = Symmetric(brand(Float64, 100, 100, 2, 4))
72+
A = Symmetric(brand(Float64, 10, 10, 2, 4))
7073
F = eigen(A)
7174
Λ, Q = F
7275
@test Q'Matrix(A)*Q ≈ Diagonal(Λ)
@@ -79,6 +82,16 @@ import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedCol
7982
@test Q[:,3] ≈ MQ[:,3]
8083
@test Q[1,2] ≈ MQ[1,2]
8184

85+
F = eigen(A, 2:4)
86+
Λ, Q = F
87+
QM = Matrix(Q)
88+
@test QM' * (Matrix(A)*QM) ≈ Diagonal(Λ)
89+
90+
F = eigen(A, 1, 2)
91+
Λ, Q = F
92+
QM = Matrix(Q)
93+
@test QM' * (Matrix(A)*QM) ≈ Diagonal(Λ)
94+
8295
function An(::Type{T}, N::Int) where {T}
8396
A = Symmetric(BandedMatrix(Zeros{T}(N,N), (0, 2)))
8497
for n = 0:N-1
@@ -177,6 +190,45 @@ end
177190
@test all(A*x .=== muladd!(one(T),A,x,zero(T),copy(x)) .===
178191
materialize!(MulAdd(one(T),A,x,zero(T),copy(x))) .===
179192
BLAS.hbmv!('L', 1, one(T), view(parent(A).data,3:4,:), x, zero(T), similar(x)))
193+
194+
@testset "eigen" begin
195+
@testset for T in (Float32, Float64, ComplexF64, ComplexF32), uplo in (:L, :U)
196+
A = brand(T, 20, 20, 4, 2)
197+
std = eigvals(Hermitian(Matrix(A), uplo))
198+
@test eigvals(Hermitian(A, uplo)) ≈ std
199+
@test eigvals(Hermitian(A, uplo), 2:4) ≈ std[2:4]
200+
@test eigvals(Hermitian(A, uplo), 1, 2) ≈ std[1 .<= std .<= 2]
201+
@test eigvals(Hermitian(big.(A), uplo)) ≈ std
202+
@test eigvals!(copy(Hermitian(A, uplo))) ≈ std
203+
@test eigvals!(copy(Hermitian(A, uplo)), 2:4) ≈ std[2:4]
204+
@test eigvals!(copy(Hermitian(A, uplo)), 1, 2) ≈ std[1 .<= std .<= 2]
205+
end
206+
207+
@testset for T in (Float32, Float64, ComplexF32, ComplexF64), uplo in (:L, :U)
208+
A = Hermitian(brand(T, 20, 20, 2, 4), uplo)
209+
MA = Hermitian(Matrix(A), uplo)
210+
λ1, v1 = eigen!(copy(A))
211+
λ2, v2 = eigen!(copy(MA))
212+
@test λ1 ≈ λ2
213+
@test v1' * MA * v1 ≈ Diagonal(λ1)
214+
215+
λ3, v3 = eigen!(copy(A), 2:4)
216+
@test λ3 ≈ λ1[2:4]
217+
@test v3' * (MA * v3) ≈ Diagonal(λ3)
218+
219+
λ4, v4 = eigen(A, 2:4)
220+
@test λ4 ≈ λ3
221+
@test v4' * (MA * v4) ≈ Diagonal(λ4)
222+
223+
λ3, v3 = eigen!(copy(A), 1, 2)
224+
@test λ3 ≈ λ1[1 .<= λ1 .<= 2]
225+
@test v3' * (MA * v3) ≈ Diagonal(λ3)
226+
227+
λ4, v4 = eigen!(copy(A), 1, 2)
228+
@test λ4 ≈ λ1[1 .<= λ1 .<= 2]
229+
@test v4' * (MA * v4) ≈ Diagonal(λ4)
230+
end
231+
end
180232
end
181233

182234
@testset "LDLᵀ" begin

0 commit comments

Comments
 (0)