Fix exceptional shift such that we can handle more matrices. (#69)

* Fix exceptional shift such that we can handle more matrices.
Also reduce the default number of iterations. Add complicated
test matrices.

* Drop Julia 1.0 support and delete azurepipelines.yml file

Author: andreasnoack

.travis.yml

@@ -4,7 +4,6 @@ os:
   - linux
 
 julia:
-  - 1.0
   - 1
   - nightly
 

Project.toml

@@ -15,4 +15,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 test = ["Quaternions", "Test"]
 
 [compat]
-julia = "1"
+julia = "1.3"

azure-pipelines.yml

@@ -85,7 +85,7 @@ function _schur!(
     tol = eps(real(T)),
     debug = false,
     shiftmethod = :Francis,
-    maxiter = 100*size(H, 1),
+    maxiter = 30*size(H, 1),
     kwargs...) where T
 
     n = size(H, 1)
@@ -99,6 +99,7 @@ function _schur!(
 
     @inbounds while true
         i += 1
+        debug && println("iteration: $i")
         if i > maxiter
             throw(ArgumentError("iteration limit $maxiter reached"))
         end
@@ -137,35 +138,20 @@ function _schur!(
         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(real(one(t))) : t # introduce a small pertubation for zero shifts
+            if iszero(i % 10)
+                # Use Eispack exceptional shift
+                β = abs(HH[iend, iend - 1]) + abs(HH[iend - 1, iend - 2])
+                d = (Hmm + β)^2 - Hmm*β/2
+                t = 2*Hmm + 3*β/2
+            else
+                d = Hm1m1*Hmm - HH[iend, iend - 1]*HH[iend - 1, iend]
+                t = Hm1m1 + Hmm
+            end
             debug && @printf("block start is: %6d, block end is: %6d, d: %10.3e, t: %10.3e\n", istart, iend, d, t)
 
             if shiftmethod == :Francis
-                # Run a bulge chase
-                if iszero(i % 10)
-                    # Vary the shift strategy to avoid dead locks
-                    # We use a Wilkinson-like shift as suggested in "Sandia technical report 96-0913J: How the QR algorithm fails to converge and how fix it".
-
-                    debug && @printf("Wilkinson-like shift! Subdiagonal is: %10.3e, last subdiagonal is: %10.3e\n", HH[iend, iend - 1], HH[iend - 1, iend - 2])
-                    _d = t*t - 4d
-                    if _d isa Real && _d >= 0
-                        # real eigenvalues
-                        a = t/2
-                        b = sqrt(_d)/2
-                        s = a > Hmm ? a - b : a + b
-                    else
-                        # complex case
-                        s = t/2
-                    end
-                    singleShiftQR!(HH, τ, s, istart, iend)
-                else
-                    # most of the time use Francis double shifts
-
-                    debug && @printf("Francis 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)
-                end
+                debug && @printf("Francis 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)
             elseif shiftmethod == :Rayleigh
                 debug && @printf("Single shift with Rayleigh shift! Subdiagonal is: %10.3e\n", HH[iend, iend - 1])
 

src/eigenGeneral.jl

@@ -93,4 +93,61 @@ end
     end
 end
 
+Demmel(η) = [0  1 0 0
+             1  0 η 0
+             0 -η 0 1
+             0  0 1 0]
+
+@testset "Demmel matrix" for t in (1e-10, 1e-9, 1e-8)
+    # See "Sandia technical report 96-0913J: How the QR algorithm fails to converge and how fix it"
+    A = Demmel(t)
+    vals = GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(A, maxiter=35))
+    @test abs.(vals) ≈ ones(4)
 end
+
+function Hevil2(θ, κ, α, γ)
+    # Eq (13) and (14)
+    β = ω = 0.0
+    ν = cos(θ)*cos(2γ) + cos(α + β + ω)*sin(2γ)*κ/2
+    σ = 1 + κ*sin(2γ)*cos(α + β + ω - θ) + κ^2*sin(γ)^2
+    μ = -sin(θ)*cos(2γ) - sin(α + β + ω)*sin(2γ)*κ/2
+    ρ = sqrt(σ - ν^2)
+
+    return [ν    (cos(2θ) - ν^2)/ρ              μ/ρ*(cos(2θ) - ν^2 + ρ^2)/sqrt(ρ^2 - μ^2)      (-2*μ*ν - sin(2*θ))/sqrt(ρ^2 - μ^2)
+            ρ   -ν - sin(2*θ)*μ/ρ^2            -μ/ρ^2*(μ*sin(2*θ) + 2*ν*ρ^2)/sqrt(ρ^2 - μ^2)    -μ/ρ*(cos(2*θ) - ν^2 + ρ^2)/sqrt(ρ^2 - μ^2)
+            0    sin(2*θ)*sqrt(ρ^2 - μ^2)/ρ^2   ν + sin(2*θ)*μ/ρ^2                              (cos(2*θ) - ν^2)/ρ
+            0    0                              ρ                                               -ν]
+end
+
+@testset "Complicate matrix from Sandia technical report" begin
+    H = Hevil2(0.111866322512629152, 1.08867072154101741, 0.338146383137297168, -0.313987810419091240)
+
+    @test Float64.(abs.(eigvals(big.(H)))) ≈ ones(4)
+end
+
+@testset "Issue 67" for (A, λs) in (
+    ([1 -2 1 -1 -1 0; 0 1 0 1 0 1; 1 -1 2 0 -1 0; 0 1 0 2 1 1; 1 0 1 0 0 0; 0 -1 1 -1 -2 0]    ,
+        [0.0, 0.0, (3 - √big(3)*im)/2, (3 + √big(3)*im)/2, (3 - √big(3)im)/2, (3 + √big(3)im)/2]),
+    ([1 0 -1 0 0 0; 0 1 1 1 -1 0; 1 0 -1 0 0 0; 1 -1 0 -1 1 1; 0 0 1 0 0 0; -1 0 1 0 0 0]      ,
+        zeros(6)),
+    ([1 -2 -1 1 -1 0; 0 1 0 -1 0 1; -1 1 2 0 1 0; 0 -1 -2 2 -1 -1; 1 0 -1 0 0 0; 0 -1 -1 1 0 0],
+        [(1 - √big(3)*im)/2, (1 + √big(3)*im)/2, (1 - √big(3)im)/2, (1 + √big(3)im)/2, 2, 2]),
+    ([2 0 -1 0 0 1; 0 2 0 -1 -1 0; -1 0 1 0 0 -1; 0 -1 0 1 1 0; 0 1 0 -1 0 0; -1 0 1 0 0 0]    ,
+        ones(6)),
+    ([1 0 1 0 -1 0; -2 1 2 -1 1 1; -1 0 -1 0 1 0; 2 1 2 -1 -2 1; 1 0 1 0 -1 0; 1 -1 -2 1 0 -1] ,
+        [-1, -1, 0, 0, 0, 0]))
+
+    @testset "shouldn't need excessive iterations (30*n) in the Float64 case" begin
+        GenericLinearAlgebra._schur!(float(A))
+    end
+
+    # For BigFloats, many iterations are required for convergence
+    # Improving this is probably a hard problem
+    vals = eigvals(big.(A), maxiter=1500)
+
+    # It's hard to compare complex conjugate pairs so we compare the real and imaginary parts separately
+    @test sort(real(vals)) ≈ sort(real(λs)) atol=1e-25
+    @test sort(imag(vals)) ≈ sort(imag(λs)) atol=1e-25
+end
+
+end