Add flag for type of shifts in general eigensolver and use Wilkinson

double shifts by default.

Add pertubation in case the trace of a 2x2 block is zero

Fixes #27

Author: andreasnoack

src/eigenGeneral.jl

@@ -70,75 +70,81 @@ 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
         end
 
         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

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