Skip to content

Commit 29d1990

Browse files
jaemolihmjaemolihm
authored
Add syevd LAPACK Hermitian eigensolver (#49262)
* Add LAPACK cheevd and zheevd wrapper * Add interface to real symmetric eigensolver syevd
1 parent 327da72 commit 29d1990

File tree

3 files changed

+119
-12
lines changed

3 files changed

+119
-12
lines changed

‎stdlib/LinearAlgebra/docs/src/index.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -838,6 +838,7 @@ LinearAlgebra.LAPACK.hetri!
838838
LinearAlgebra.LAPACK.hetrs!
839839
LinearAlgebra.LAPACK.syev!
840840
LinearAlgebra.LAPACK.syevr!
841+
LinearAlgebra.LAPACK.syevd!
841842
LinearAlgebra.LAPACK.sygvd!
842843
LinearAlgebra.LAPACK.bdsqr!
843844
LinearAlgebra.LAPACK.bdsdc!

‎stdlib/LinearAlgebra/src/lapack.jl

Lines changed: 111 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -5170,9 +5170,9 @@ solution `X`.
51705170
hetrs!(uplo::AbstractChar, A::AbstractMatrix, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat)
51715171

51725172
# Symmetric (real) eigensolvers
5173-
for (syev, syevr, sygvd, elty) in
5174-
((:dsyev_,:dsyevr_,:dsygvd_,:Float64),
5175-
(:ssyev_,:ssyevr_,:ssygvd_,:Float32))
5173+
for (syev, syevr, syevd, sygvd, elty) in
5174+
((:dsyev_,:dsyevr_,:dsyevd_,:dsygvd_,:Float64),
5175+
(:ssyev_,:ssyevr_,:ssyevd_,:ssygvd_,:Float32))
51765176
@eval begin
51775177
# SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
51785178
# * .. Scalar Arguments ..
@@ -5225,7 +5225,7 @@ for (syev, syevr, sygvd, elty) in
52255225
end
52265226
lda = stride(A,2)
52275227
m = Ref{BlasInt}()
5228-
w = similar(A, $elty, n)
5228+
W = similar(A, $elty, n)
52295229
ldz = n
52305230
if jobz == 'N'
52315231
Z = similar(A, $elty, ldz, 0)
@@ -5249,7 +5249,7 @@ for (syev, syevr, sygvd, elty) in
52495249
jobz, range, uplo, n,
52505250
A, max(1,lda), vl, vu,
52515251
il, iu, abstol, m,
5252-
w, Z, max(1,ldz), isuppz,
5252+
W, Z, max(1,ldz), isuppz,
52535253
work, lwork, iwork, liwork,
52545254
info, 1, 1, 1)
52555255
chklapackerror(info[])
@@ -5260,11 +5260,51 @@ for (syev, syevr, sygvd, elty) in
52605260
resize!(iwork, liwork)
52615261
end
52625262
end
5263-
w[1:m[]], Z[:,1:(jobz == 'V' ? m[] : 0)]
5263+
W[1:m[]], Z[:,1:(jobz == 'V' ? m[] : 0)]
52645264
end
52655265
syevr!(jobz::AbstractChar, A::AbstractMatrix{$elty}) =
52665266
syevr!(jobz, 'A', 'U', A, 0.0, 0.0, 0, 0, -1.0)
52675267

5268+
# SUBROUTINE DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK,
5269+
# $ IWORK, LIWORK, INFO )
5270+
# * .. Scalar Arguments ..
5271+
# CHARACTER JOBZ, UPLO
5272+
# INTEGER INFO, LDA, LIWORK, LWORK, N
5273+
# * ..
5274+
# * .. Array Arguments ..
5275+
# INTEGER IWORK( * )
5276+
# DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
5277+
function syevd!(jobz::AbstractChar, uplo::AbstractChar, A::AbstractMatrix{$elty})
5278+
chkstride1(A)
5279+
n = checksquare(A)
5280+
chkuplofinite(A, uplo)
5281+
lda = stride(A,2)
5282+
m = Ref{BlasInt}()
5283+
W = similar(A, $elty, n)
5284+
work = Vector{$elty}(undef, 1)
5285+
lwork = BlasInt(-1)
5286+
iwork = Vector{BlasInt}(undef, 1)
5287+
liwork = BlasInt(-1)
5288+
info = Ref{BlasInt}()
5289+
for i = 1:2 # first call returns lwork as work[1] and liwork as iwork[1]
5290+
ccall((@blasfunc($syevd), libblastrampoline), Cvoid,
5291+
(Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
5292+
Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
5293+
Ptr{BlasInt}, Clong, Clong),
5294+
jobz, uplo, n, A, max(1,lda),
5295+
W, work, lwork, iwork, liwork,
5296+
info, 1, 1)
5297+
chklapackerror(info[])
5298+
if i == 1
5299+
lwork = BlasInt(real(work[1]))
5300+
resize!(work, lwork)
5301+
liwork = iwork[1]
5302+
resize!(iwork, liwork)
5303+
end
5304+
end
5305+
jobz == 'V' ? (W, A) : W
5306+
end
5307+
52685308
# Generalized eigenproblem
52695309
# SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
52705310
# $ LWORK, IWORK, LIWORK, INFO )
@@ -5313,9 +5353,9 @@ for (syev, syevr, sygvd, elty) in
53135353
end
53145354
end
53155355
# Hermitian eigensolvers
5316-
for (syev, syevr, sygvd, elty, relty) in
5317-
((:zheev_,:zheevr_,:zhegvd_,:ComplexF64,:Float64),
5318-
(:cheev_,:cheevr_,:chegvd_,:ComplexF32,:Float32))
5356+
for (syev, syevr, syevd, sygvd, elty, relty) in
5357+
((:zheev_,:zheevr_,:zheevd_,:zhegvd_,:ComplexF64,:Float64),
5358+
(:cheev_,:cheevr_,:cheevd_,:chegvd_,:ComplexF32,:Float32))
53195359
@eval begin
53205360
# SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO )
53215361
# * .. Scalar Arguments ..
@@ -5376,7 +5416,7 @@ for (syev, syevr, sygvd, elty, relty) in
53765416
end
53775417
lda = max(1,stride(A,2))
53785418
m = Ref{BlasInt}()
5379-
w = similar(A, $relty, n)
5419+
W = similar(A, $relty, n)
53805420
if jobz == 'N'
53815421
ldz = 1
53825422
Z = similar(A, $elty, ldz, 0)
@@ -5404,7 +5444,7 @@ for (syev, syevr, sygvd, elty, relty) in
54045444
jobz, range, uplo, n,
54055445
A, lda, vl, vu,
54065446
il, iu, abstol, m,
5407-
w, Z, ldz, isuppz,
5447+
W, Z, ldz, isuppz,
54085448
work, lwork, rwork, lrwork,
54095449
iwork, liwork, info,
54105450
1, 1, 1)
@@ -5418,11 +5458,56 @@ for (syev, syevr, sygvd, elty, relty) in
54185458
resize!(iwork, liwork)
54195459
end
54205460
end
5421-
w[1:m[]], Z[:,1:(jobz == 'V' ? m[] : 0)]
5461+
W[1:m[]], Z[:,1:(jobz == 'V' ? m[] : 0)]
54225462
end
54235463
syevr!(jobz::AbstractChar, A::AbstractMatrix{$elty}) =
54245464
syevr!(jobz, 'A', 'U', A, 0.0, 0.0, 0, 0, -1.0)
54255465

5466+
# SUBROUTINE ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
5467+
# $ LRWORK, IWORK, LIWORK, INFO )
5468+
# * .. Scalar Arguments ..
5469+
# CHARACTER JOBZ, UPLO
5470+
# INTEGER INFO, LDA, LIWORK, LRWORK, LWORK, N
5471+
# * ..
5472+
# * .. Array Arguments ..
5473+
# INTEGER IWORK( * )
5474+
# DOUBLE PRECISION RWORK( * )
5475+
# COMPLEX*16 A( LDA, * ), WORK( * )
5476+
function syevd!(jobz::AbstractChar, uplo::AbstractChar, A::AbstractMatrix{$elty})
5477+
chkstride1(A)
5478+
chkuplofinite(A, uplo)
5479+
n = checksquare(A)
5480+
lda = max(1, stride(A,2))
5481+
m = Ref{BlasInt}()
5482+
W = similar(A, $relty, n)
5483+
work = Vector{$elty}(undef, 1)
5484+
lwork = BlasInt(-1)
5485+
rwork = Vector{$relty}(undef, 1)
5486+
lrwork = BlasInt(-1)
5487+
iwork = Vector{BlasInt}(undef, 1)
5488+
liwork = BlasInt(-1)
5489+
info = Ref{BlasInt}()
5490+
for i = 1:2 # first call returns lwork as work[1], lrwork as rwork[1] and liwork as iwork[1]
5491+
ccall((@blasfunc($syevd), liblapack), Cvoid,
5492+
(Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
5493+
Ptr{$relty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$relty}, Ref{BlasInt},
5494+
Ptr{BlasInt}, Ref{BlasInt}, Ptr{BlasInt}, Clong, Clong),
5495+
jobz, uplo, n, A, stride(A,2),
5496+
W, work, lwork, rwork, lrwork,
5497+
iwork, liwork, info, 1, 1)
5498+
chklapackerror(info[])
5499+
if i == 1
5500+
lwork = BlasInt(real(work[1]))
5501+
resize!(work, lwork)
5502+
lrwork = BlasInt(rwork[1])
5503+
resize!(rwork, lrwork)
5504+
liwork = iwork[1]
5505+
resize!(iwork, liwork)
5506+
end
5507+
end
5508+
jobz == 'V' ? (W, A) : W
5509+
end
5510+
54265511
# SUBROUTINE ZHEGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
54275512
# $ LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
54285513
# * .. Scalar Arguments ..
@@ -5504,6 +5589,20 @@ The eigenvalues are returned in `W` and the eigenvectors in `Z`.
55045589
syevr!(jobz::AbstractChar, range::AbstractChar, uplo::AbstractChar, A::AbstractMatrix,
55055590
vl::AbstractFloat, vu::AbstractFloat, il::Integer, iu::Integer, abstol::AbstractFloat)
55065591

5592+
"""
5593+
syevd!(jobz, uplo, A)
5594+
5595+
Finds the eigenvalues (`jobz = N`) or eigenvalues and eigenvectors
5596+
(`jobz = V`) of a symmetric matrix `A`. If `uplo = U`, the upper triangle
5597+
of `A` is used. If `uplo = L`, the lower triangle of `A` is used.
5598+
5599+
Use the divide-and-conquer method, instead of the QR iteration used by
5600+
`syev!` or multiple relatively robust representations used by `syevr!`.
5601+
See James W. Demmel et al, SIAM J. Sci. Comput. 30, 3, 1508 (2008) for
5602+
a comparison of the accuracy and performatce of different methods.
5603+
"""
5604+
syevd!(jobz::AbstractChar, uplo::AbstractChar, A::AbstractMatrix)
5605+
55075606
"""
55085607
sygvd!(itype, jobz, uplo, A, B) -> (w, A, B)
55095608

‎stdlib/LinearAlgebra/test/lapack.jl

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,13 @@ using LinearAlgebra: BlasInt
2727
@test LAPACK.syevr!('N', 'V', 'U', copy(Asym), 0.0, 1.0, 4, 5, -1.0)[1] ≈ vals[vals .< 1.0]
2828
@test LAPACK.syevr!('N', 'I', 'U', copy(Asym), 0.0, 1.0, 4, 5, -1.0)[1] ≈ vals[4:5]
2929
@test vals ≈ LAPACK.syev!('N', 'U', copy(Asym))
30+
@test vals ≈ LAPACK.syevd!('N', 'U', copy(Asym))
31+
vals_test, Z_test = LAPACK.syev!('V', 'U', copy(Asym))
32+
@test vals_test ≈ vals
33+
@test Z_test*(Diagonal(vals)*Z_test') ≈ Asym
34+
vals_test, Z_test = LAPACK.syevd!('V', 'U', copy(Asym))
35+
@test vals_test ≈ vals
36+
@test Z_test*(Diagonal(vals)*Z_test') ≈ Asym
3037
@test_throws DimensionMismatch LAPACK.sygvd!(1, 'V', 'U', copy(Asym), zeros(elty, 6, 6))
3138
end
3239
end

0 commit comments

Comments
 (0)