diff --git a/src/eigenGeneral.jl b/src/eigenGeneral.jl index 34ef581..49b71b6 100644 --- a/src/eigenGeneral.jl +++ b/src/eigenGeneral.jl @@ -70,67 +70,73 @@ module EigenGeneral R::Rotation end - function schurfact!{T<:Real}(H::HessenbergFactorization{T}; tol = eps(T), debug = false) + function wilkinson(Hmm, t, d) + λ1 = (t + sqrt(t*t - 4d))/2 + λ2 = (t - sqrt(t*t - 4d))/2 + return ifelse(abs(Hmm - λ1) < abs(Hmm - λ2), λ1, λ2) + end + + + function schurfact!{T<:Real}(H::HessenbergFactorization{T}; tol = eps(T), debug = false, shiftmethod = :Wilkinson, maxiter = 100*size(H, 1)) n = size(H, 1) istart = 1 iend = n HH = H.data τ = Rotation(Givens{T}[]) + # iteration count + i = 0 + @inbounds while true + i += 1 + if i > maxiter + throw(ArgumentError("iteration limit $maxiter reached")) + end + # Determine if the matrix splits. Find lowest positioned subdiagonal "zero" for istart = iend - 1:-1:1 - # debug && @printf("istart: %6d, iend %6d\n", istart, iend) - # istart == minstart && break if abs(HH[istart + 1, istart]) < tol*(abs(HH[istart, istart]) + abs(HH[istart + 1, istart + 1])) istart += 1 - debug && @printf("Top deflation! Subdiagonal element is: %10.3e and istart now %6d\n", HH[istart, istart - 1], istart) + debug && @printf("Split! Subdiagonal element is: %10.3e and istart now %6d\n", HH[istart, istart - 1], istart) break elseif istart > 1 && abs(HH[istart, istart - 1]) < tol*(abs(HH[istart - 1, istart - 1]) + abs(HH[istart, istart])) - debug && @printf("Top deflation! Next subdiagonal element is: %10.3e and istart now %6d\n", HH[istart, istart - 1], istart) + debug && @printf("Split! Next subdiagonal element is: %10.3e and istart now %6d\n", HH[istart, istart - 1], istart) break end end # if block size is one we deflate if istart >= iend + debug && @printf("Bottom deflation! Block size is one. New iend is %6d\n", iend - 1) iend -= 1 # and the same for a 2x2 block elseif istart + 1 == iend + debug && @printf("Bottom deflation! Block size is two. New iend is %6d\n", iend - 2) iend -= 2 - # if we don't deflate we'll run either a single or double shift bulge chase + # run a QR iteration + # shift method is specified with shiftmethod kw argument else Hmm = HH[iend, iend] Hm1m1 = HH[iend - 1, iend - 1] d = Hm1m1*Hmm - HH[iend, iend - 1]*HH[iend - 1, iend] t = Hm1m1 + Hmm + t = iszero(t) ? eps(one(t)) : t # introduce a small pertubation for zero shifts debug && @printf("block start is: %6d, block end is: %6d, d: %10.3e, t: %10.3e\n", istart, iend, d, t) - # For small (sub) problems use Raleigh quotion shift and single shift - if iend <= istart + 2 - σ = HH[iend, iend] + if shiftmethod == :Wilkinson + debug && @printf("Double shift with Wilkinson shift! Subdiagonal is: %10.3e, last subdiagonal is: %10.3e\n", HH[iend, iend - 1], HH[iend - 1, iend - 2]) # Run a bulge chase - singleShiftQR!(HH, τ, σ, istart, iend) - - # If the eigenvales of the 2x2 block are real use single shift - elseif t*t > 4d - debug && @printf("Single shift! subdiagonal is: %10.3e\n", HH[iend, iend - 1]) - - # Calculate the Wilkinson shift - λ1 = (t + sqrt(t*t - 4d))/2 - λ2 = (t - sqrt(t*t - 4d))/2 - σ = abs(Hmm - λ1) < abs(Hmm - λ2) ? λ1 : λ2 + doubleShiftQR!(HH, τ, t, d, istart, iend) + elseif shiftmethod == :Rayleigh + debug && @printf("Single shift with Rayleigh shift! Subdiagonal is: %10.3e\n", HH[iend, iend - 1]) # Run a bulge chase - singleShiftQR!(HH, τ, σ, istart, iend) - - # else use double shift + singleShiftQR!(HH, τ, Hmm, istart, iend) else - debug && @printf("Double shift! subdiagonal is: %10.3e, last subdiagonal is: %10.3e\n", HH[iend, iend - 1], HH[iend - 1, iend - 2]) - doubleShiftQR!(HH, τ, t, d, istart, iend) + throw(ArgumentError("only support supported shift methods are :Wilkinson (default) and :Rayleigh. You supplied $shiftmethod")) end end if iend <= 2 break end @@ -138,7 +144,7 @@ module EigenGeneral return Schur{T,typeof(HH)}(HH, τ) end - schurfact!(A::StridedMatrix; tol = eps(float(real(eltype(A)))), debug = false) = schurfact!(hessfact!(A), tol = tol, debug = debug) + schurfact!(A::StridedMatrix; kwargs...) = schurfact!(hessfact!(A); kwargs...) function singleShiftQR!(HH::StridedMatrix, τ::Rotation, shift::Number, istart::Integer, iend::Integer) m = size(HH, 1) @@ -172,16 +178,21 @@ module EigenGeneral H21 = HH[istart + 1, istart] Htmp11 = HH[istart + 2, istart] HH[istart + 2, istart] = 0 - Htmp21 = HH[istart + 3, istart] - HH[istart + 3, istart] = 0 - Htmp22 = HH[istart + 3, istart + 1] - HH[istart + 3, istart + 1] = 0 + if istart + 3 <= m + Htmp21 = HH[istart + 3, istart] + HH[istart + 3, istart] = 0 + Htmp22 = HH[istart + 3, istart + 1] + HH[istart + 3, istart + 1] = 0 + else + # values doen't matter in this case but variables should be initialized + Htmp21 = Htmp22 = Htmp11 + end G1, r = givens(H11*H11 + HH[istart, istart + 1]*H21 - shiftTrace*H11 + shiftDeterminant, H21*(H11 + HH[istart + 1, istart + 1] - shiftTrace), istart, istart + 1) G2, _ = givens(r, H21*HH[istart + 2, istart + 1], istart, istart + 2) vHH = view(HH, :, istart:m) A_mul_B!(G1, vHH) A_mul_B!(G2, vHH) - vHH = view(HH, 1:istart + 3, :) + vHH = view(HH, 1:min(istart + 3, m), :) A_mul_Bc!(vHH, G1) A_mul_Bc!(vHH, G2) A_mul_B!(G1, τ) @@ -209,9 +220,9 @@ module EigenGeneral return HH end - eigvals!(A::StridedMatrix; tol = eps(one(A[1])), debug = false) = eigvals!(schurfact!(A, tol = tol, debug = debug)) - eigvals!(H::HessenbergMatrix; tol = eps(one(A[1])), debug = false) = eigvals!(schurfact!(H, tol = tol, debug = debug)) - eigvals!(H::HessenbergFactorization; tol = eps(one(A[1])), debug = false) = eigvals!(schurfact!(H, tol = tol, debug = debug)) + eigvals!(A::StridedMatrix; kwargs...) = eigvals!(schurfact!(A; kwargs...)) + eigvals!(H::HessenbergMatrix; kwargs...) = eigvals!(schurfact!(H, kwargs...)) + eigvals!(H::HessenbergFactorization; kwargs...) = eigvals!(schurfact!(H, kwargs...)) function eigvals!{T}(S::Schur{T}; tol = eps(T)) HH = S.data diff --git a/test/eigengeneral.jl b/test/eigengeneral.jl index ef2d6ce..11ef538 100644 --- a/test/eigengeneral.jl +++ b/test/eigengeneral.jl @@ -10,4 +10,11 @@ using LinearAlgebra @test sort(imag(v1)) ≈ sort(imag(v2)) @test sort(real(v1)) ≈ sort(real(map(Complex{Float64}, vBig))) @test sort(imag(v1)) ≈ sort(imag(map(Complex{Float64}, vBig))) +end + +@testset "make sure that solver doesn't hang" begin + for i in 1:1000 + A = randn(8, 8) + sort(abs.(LinearAlgebra.EigenGeneral.eigvals!(copy(A)))) ≈ sort(abs.(eigvals(A))) + end end \ No newline at end of file