Skip to content

Commit

Permalink
Fix exceptional shift such that we can handle more matrices.
Browse files Browse the repository at this point in the history 
Also reduce the default number of iterations. Add complicated
test matrices.
  • Loading branch information
@andreasnoack
andreasnoack committed Aug 4, 2020
1 parent a3140d3 commit e66fda1
Show file tree
Hide file tree
Showing 2 changed files with 79 additions and 28 deletions.
50 changes: 22 additions & 28 deletions src/eigenGeneral.jl
View file
Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -137,40 +138,33 @@ 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
# t = iszero(t) ? eps(real(one(t))) : t # introduce a small pertubation for zero shifts
# if t*t >= 4*d
# # The real case
# # We make a double shift our of the Wilkinson shift
# λ̃ = wilkinson(Hmm, t, d)
# t = 2λ̃
# d = λ̃^2
# end
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])

# Run a bulge chase
singleShiftQR!(HH, τ, Hmm, istart, iend)
singleShiftQR!(HH, τ, σ, istart, iend)
else
throw(ArgumentError("only support supported shift methods are :Francis (default) and :Rayleigh. You supplied $shiftmethod"))
end
Expand Down
57 changes: 57 additions & 0 deletions test/eigengeneral.jl
View file
Expand Up @@ -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

0 comments on commit e66fda1

Please sign in to comment.