Skip to content

Commit 9d967cc

Browse files
andreasnoackandreasnoack
committed
Fix bug in symmetric eigensolvers for 1x1 matrices.
Convert many of the positional arguments in the symmetric solver into keyword arguments.
1 parent d735f8d commit 9d967cc

File tree

2 files changed

+84
-28
lines changed

2 files changed

+84
-28
lines changed

‎src/eigenSelfAdjoint.jl

Lines changed: 77 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -160,7 +160,7 @@ function eigQL2x2!(d::StridedVector, e::StridedVector, j::Integer, vectors::Matr
160160
return c, s
161161
end
162162

163-
function eigvalsPWK!(S::SymTridiagonal{T}, tol = eps(T), debug::Bool=false) where T<:Real
163+
function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), debug::Bool=false) where T<:Real
164164
d = S.dv
165165
e = S.ev
166166
n = length(d)
@@ -180,6 +180,7 @@ function eigvalsPWK!(S::SymTridiagonal{T}, tol = eps(T), debug::Bool=false) wher
180180
end
181181
blockend = n
182182
end
183+
183184
# Deflate?
184185
if blockstart == blockend
185186
# Yes
@@ -201,15 +202,15 @@ function eigvalsPWK!(S::SymTridiagonal{T}, tol = eps(T), debug::Bool=false) wher
201202

202203
debug && @printf("QL, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, μ: %e, rotations: %d\n", blockstart, blockend, e[blockstart], e[blockend-1], μ, iter += blockend - blockstart)
203204
end
204-
if blockstart == n
205+
if blockstart >= n
205206
break
206207
end
207208
end
208209
end
209210
sort!(d)
210211
end
211212

212-
function eigQL!(S::SymTridiagonal{T},
213+
function eigQL!(S::SymTridiagonal{T};
213214
vectors::Matrix = zeros(T, 0, size(S, 1)),
214215
tol = eps(T),
215216
debug::Bool=false) where T<:Real
@@ -254,7 +255,7 @@ function eigQL!(S::SymTridiagonal{T},
254255
singleShiftQL!(S, μ, blockstart, blockend, vectors)
255256
debug && @printf("QL, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, μ: %f\n", blockstart, blockend, e[blockstart], e[blockend-1], μ)
256257
end
257-
if blockstart == n
258+
if blockstart >= n
258259
break
259260
end
260261
end
@@ -263,7 +264,7 @@ function eigQL!(S::SymTridiagonal{T},
263264
return d[p], vectors[:,p]
264265
end
265266

266-
function eigQR!(S::SymTridiagonal{T},
267+
function eigQR!(S::SymTridiagonal{T};
267268
vectors::Matrix = zeros(T, 0, size(S, 1)),
268269
tol = eps(T),
269270
debug::Bool=false) where T<:Real
@@ -301,7 +302,7 @@ function eigQR!(S::SymTridiagonal{T},
301302

302303
debug && @printf("QR, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, d[blockend]: %f, μ: %f\n", blockstart, blockend, e[blockstart], e[blockend-1], d[blockend], μ)
303304
end
304-
if blockstart == n
305+
if blockstart >= n
305306
break
306307
end
307308
end
@@ -527,40 +528,89 @@ function symtriUpper!(AS::StridedMatrix{T},
527528
SymmetricTridiagonalFactorization('U', AS, τ, SymTridiagonal(real(diag(AS)), real(diag(AS, 1))))
528529
end
529530

530-
_eigvals!(A::SymmetricTridiagonalFactorization) = eigvalsPWK!(A.diagonals, eps(eltype(A.diagonals)), false)
531-
_eigvals!(A::SymTridiagonal ) = eigvalsPWK!(A, eps(eltype(A)) , false)
532-
_eigvals!(A::Hermitian ) = eigvals!(symtri!(A))
533531

534-
LinearAlgebra.eigvals!(A::SymmetricTridiagonalFactorization) = _eigvals!(A)
535-
LinearAlgebra.eigvals!(A::SymTridiagonal ) = _eigvals!(A)
536-
LinearAlgebra.eigvals!(A::Hermitian ) = _eigvals!(A)
532+
_eigvals!(A::SymmetricTridiagonalFactorization;
533+
tol = eps(real(eltype(A))),
534+
debug = false) = eigvalsPWK!(A.diagonals, tol = tol, debug = debug)
535+
536+
_eigvals!(A::SymTridiagonal;
537+
tol = eps(real(eltype(A))),
538+
debug = false) = eigvalsPWK!(A, tol = tol, debug = debug)
539+
540+
_eigvals!(A::Hermitian;
541+
tol = eps(real(eltype(A))),
542+
debug = false) = eigvals!(symtri!(A), tol = tol, debug = debug)
543+
544+
545+
LinearAlgebra.eigvals!(A::SymmetricTridiagonalFactorization;
546+
tol = eps(real(eltype(A))),
547+
debug = false) = _eigvals!(A, tol = tol, debug = debug)
548+
549+
LinearAlgebra.eigvals!(A::SymTridiagonal;
550+
tol = eps(real(eltype(A))),
551+
debug = false) = _eigvals!(A, tol = tol, debug = debug)
552+
553+
LinearAlgebra.eigvals!(A::Hermitian;
554+
tol = eps(real(eltype(A))),
555+
debug = false) = _eigvals!(A, tol = tol, debug = debug)
556+
557+
558+
_eigen!(A::SymmetricTridiagonalFactorization;
559+
tol = eps(real(eltype(A))),
560+
debug = false) =
561+
LinearAlgebra.Eigen(eigQL!(A.diagonals, vectors = Array(A.Q), tol = tol, debug = debug)...)
537562

538-
_eigen!(A::SymmetricTridiagonalFactorization) =
539-
LinearAlgebra.Eigen(eigQL!(A.diagonals, Array(A.Q), eps(eltype(A.diagonals)), false)...)
540-
_eigen!(A::SymTridiagonal) =
541-
LinearAlgebra.Eigen(eigQL!(A, Matrix{eltype(A)}(I, size(A, 1), size(A, 1)), eps(eltype(A)), false)...)
542-
_eigen!(A::Hermitian) = _eigen!(symtri!(A))
563+
_eigen!(A::SymTridiagonal;
564+
tol = eps(real(eltype(A))),
565+
debug = false) =
566+
LinearAlgebra.Eigen(eigQL!(A, vectors = Matrix{eltype(A)}(I, size(A, 1), size(A, 1)), tol = tol, debug = debug)...)
543567

544-
LinearAlgebra.eigen!(A::SymmetricTridiagonalFactorization) = _eigen!(A)
545-
LinearAlgebra.eigen!(A::SymTridiagonal ) = _eigen!(A)
546-
LinearAlgebra.eigen!(A::Hermitian ) = _eigen!(A)
568+
_eigen!(A::Hermitian;
569+
tol = eps(real(eltype(A))),
570+
debug = false) = _eigen!(symtri!(A), tol = tol, debug = debug)
547571

548-
function eigen2!(A::SymmetricTridiagonalFactorization, tol = eps(real(float(one(eltype(A))))), debug = false)
572+
573+
LinearAlgebra.eigen!(A::SymmetricTridiagonalFactorization;
574+
tol = eps(real(eltype(A))),
575+
debug = false) = _eigen!(A, tol = tol, debug = debug)
576+
577+
LinearAlgebra.eigen!(A::SymTridiagonal;
578+
tol = eps(real(eltype(A))),
579+
debug = false) = _eigen!(A, tol = tol, debug = debug)
580+
581+
LinearAlgebra.eigen!(A::Hermitian;
582+
tol = eps(real(eltype(A))),
583+
debug = false) = _eigen!(A, tol = tol, debug = debug)
584+
585+
586+
function eigen2!(A::SymmetricTridiagonalFactorization;
587+
tol = eps(real(float(one(eltype(A))))),
588+
debug = false)
549589
V = zeros(eltype(A), 2, size(A, 1))
550590
V[1] = 1
551591
V[end] = 1
552-
eigQL!(A.diagonals, rmul!(V, A.Q), tol, debug)
592+
eigQL!(A.diagonals, vectors = rmul!(V, A.Q), tol = tol, debug = debug)
553593
end
554-
function eigen2!(A::SymTridiagonal, tol = eps(real(float(one(eltype(A))))), debug = false)
594+
595+
function eigen2!(A::SymTridiagonal;
596+
tol = eps(real(float(one(eltype(A))))),
597+
debug = false)
555598
V = zeros(eltype(A), 2, size(A, 1))
556599
V[1] = 1
557600
V[end] = 1
558-
eigQL!(A, V, tol, debug)
601+
eigQL!(A, vectors = V, tol = tol, debug = debug)
559602
end
560-
eigen2!(A::Hermitian, tol = eps(float(real(one(eltype(A))))), debug = false) = eigen2!(symtri!(A), tol, debug)
561603

562-
eigen2(A::SymTridiagonal, tol = eps(float(real(one(eltype(A))))), debug = false) = eigen2!(copy(A), tol, debug)
563-
eigen2(A::Hermitian , tol = eps(float(real(one(eltype(A))))), debug = false) = eigen2!(copy(A), tol, debug)
604+
eigen2!(A::Hermitian;
605+
tol = eps(float(real(one(eltype(A))))),
606+
debug = false) = eigen2!(symtri!(A), tol = tol, debug = debug)
607+
608+
609+
eigen2(A::SymTridiagonal;
610+
tol = eps(float(real(one(eltype(A))))),
611+
debug = false) = eigen2!(copy(A), tol = tol, debug = debug)
612+
613+
eigen2(A::Hermitian , tol = eps(float(real(one(eltype(A))))), debug = false) = eigen2!(copy(A), tol = tol, debug = debug)
564614

565615
# First method of each type here is identical to the method defined in
566616
# LinearAlgebra but is needed for disambiguation

‎test/eigenselfadjoint.jl

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
2323
end
2424

2525
@testset "QR version (QL is default)" begin
26-
vals, vecs = GenericLinearAlgebra.eigQR!(copy(T), Matrix{eltype(T)}(I, n, n))
26+
vals, vecs = GenericLinearAlgebra.eigQR!(copy(T), vectors = Matrix{eltype(T)}(I, n, n))
2727
@test (vecs'*T)*vecs ≈ Diagonal(vals)
2828
@test eigvals(T) ≈ vals
2929
@test vecs'vecs ≈ Matrix(I, n, n)
@@ -77,4 +77,10 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
7777
@test c*x + s*y ≈ r
7878
@test c*y ≈ s*x
7979
end
80+
81+
@testset "out-of-bounds issue in 1x1 case" begin
82+
@test GenericLinearAlgebra._eigvals!(SymTridiagonal([1.0], Float64[])) == [1.0]
83+
@test GenericLinearAlgebra._eigen!(SymTridiagonal([1.0], Float64[])).values == [1.0]
84+
@test GenericLinearAlgebra._eigen!(SymTridiagonal([1.0], Float64[])).vectors == fill(1.0, 1, 1)
85+
end
8086
end

0 commit comments

Comments
 (0)