Comparing stuff

schur.jl

# Comparison of bulge chasing of real matrix with modified GenericLinearAlgebra vs https://github.com/haampie/IRAM.jl/pull/48
# Run with Julia Version 0.7.0-beta2.6.
# Modifications to GenericLinearAlgebra only include commenting out lines that apply / store rotations.

using LinearAlgebra, GenericLinearAlgebra, IRAM, BenchmarkTools

function compare_hessenberg_to_uppertriangular(n)
    # Random real Hessenberg matrix
    H = triu(rand(n, n), -1)

    H_iram = copy(H)
    H_gla = copy(H)

    GenericLinearAlgebra.schurfact!(H_gla)
    IRAM._local_schurfact!(H_iram, 1, n)

    # Make sure the number of nonzero off-diagonal values matches the number of conjugate eigvals
    f = x -> (real(x), imag(x))
    λs = sort!(eigvals(H), by = f)
    λs_iram = sort!(eigvals(H_iram), by = f)
    λs_gla = sort!(eigvals(H_gla), by = f)
    @assert 2 * count(!iszero, diag(H_iram, -1)) == count(!iszero ∘ imag, λs)
    
    # Make sure the eigenvalues are equal
    @show norm(λs_iram .- λs)
    @show norm(λs_gla .- λs)
    @show norm(λs_iram .- λs_gla)

    # Benchmark
    a = @benchmark GenericLinearAlgebra.schurfact!(HH) setup = (HH = GenericLinearAlgebra.hessfact!($H))
    b = @benchmark IRAM._local_schurfact!(HH, 1, $n) setup = (HH = copy($H))

    a, b
end

xrun.jl

julia> compare_hessenberg_to_uppertriangular(100)
norm(λs_iram .- λs) = 2.6629984796011266e-7
norm(λs_gla .- λs) = 2.662998761233004e-7
norm(λs_iram .- λs_gla) = 7.812790063035723e-14
(Trial(625.183 μs), Trial(395.125 ns))

julia> compare_hessenberg_to_uppertriangular(200)
norm(λs_iram .- λs) = 3.030021973997769e-5
norm(λs_gla .- λs) = 3.0300219748319796e-5
norm(λs_iram .- λs_gla) = 1.0270952878895238e-13
(Trial(28.105 ms), Trial(879.316 ns))