relation

stuff.jl

using LinearAlgebra
using Test

function generate_real_H_with_imaginary_eigs(n, T::Type = Float64)
    while true
        H = triu(rand(T, n + 1, n), -1)
        λs = sort!(eigvals(view(H, 1 : n, 1 : n)), by = abs)

        for i = 1 : n
            μ = λs[i]
            if imag(μ) != 0
                # Remove conjugate pair
                deleteat!(λs, (i, i + 1))
                return H, λs, μ
            end
        end
    end
end

function helloworld(n = 5)
    H, λs, μ = generate_real_H_with_imaginary_eigs(n)

    # Use complex shifts matrix for simplicity
    Hc = complex(H)

    Q = Matrix{ComplexF64}(I, n + 1, n + 1)

    # First QR step

    for i = 1 : n
        Hc[i, i] -= μ
    end

    Q₁′, R₁ = qr(Hc)
    Q₁ = Matrix(Q₁′)
    
    Q *= Q₁
    Hc = R₁ * Q₁[1 : n, 1 : n - 1]
    for i = 1 : n - 1
        Hc[i,i] += μ
    end

    # Second QR step

    for i = 1 : n - 1
        Hc[i, i] -= μ'
    end

    Q₂′, R₂ = qr(Hc)
    Q₂ = Matrix(Q₂′)

    Q *= Q₂
    Hc = R₂ * Q₂[1 : n - 1, 1 : n - 2]
    for i = 1 : n - 2
        Hc[i,i] += μ'
    end

    # Imaginary part should be 0
    Hc = real(Hc)




    ### Relations that should hold:
    @test norm(Q' * H * Q[1 : n, 1 : n - 2] - Hc) < 1e-14

    return H, Hc, Q
end