Hessenberg factorizations using Householder transformations

hessfact1.jl

using BenchmarkTools
using Test
using LinearAlgebra

"""
    hessfact!(H, Q, A)

Calculate Hessenberg factorization H = Q' * A * Q using Householder
transformations.
"""
function hessfact!(H, Q, A)
    n, m = size(H)
    n == m || error("Not a square matrix")

 @. H = A
    fill!(Q, 0.0)
    for i in 1:n
        Q[i,i] = 1.0
    end

    u = zeros(n)
    for j in 1:n-2
        u = H[j+1:end, j]
        u[1] = u[1] + sign(u[1]) * norm(u)
        v = u / norm(u)
        H[j+1:end, :] = H[j+1:end, :] - 2*v*(v'*H[j+1:end,:])
        H[:, j+1:end] = H[:, j+1:end] - (H[:,j+1:end] * (2*v)) * v'
        Q[:, j+1:end] = Q[:, j+1:end] - (Q[:,j+1:end] * (2*v)) * v'
    end
end

const A = rand(100,100)
const H = zeros(size(A))
const Q = zeros(size(A))

hessfact!(H, Q, A)
@test isapprox(H, Q'A*Q)
@btime hessfact!($H, $Q, $A)

# 11.572 ms (3203 allocations: 45.92 MiB)