hessfact2.jl

using BenchmarkTools
using Test
using LinearAlgebra
using Random
Random.seed!(1)

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

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

    @inbounds for i=1:n
        @inbounds for j=1:n
            Q[i,j] = i == j ? 1.0 : 0.0
            H[i,j] = A[i,j]
        end
    end

    vh = zeros(n, 1)
    hv = zeros(n, 1)
    qv = zeros(n, 1)

    @inbounds for j=1:n-2
        @inbounds for i=1:n
            u[i] = i < j+1 ? 0.0 : H[i, j]
        end
        u[j+1] = u[j+1] + sign(u[j+1]) * norm(u)
        normu = norm(u)
        @inbounds for i=1:n
            v[i] = u[i] / normu
        end

        mul!(vh, H', v)
        @inbounds for i=j+1:n
            @inbounds for k=1:n
                H[i,k] = H[i,k] - 2 * v[i] * vh[k]
            end
        end

        mul!(hv, H, v)
        mul!(qv, Q, v)
        @inbounds for i=1:n
            @inbounds for k=j+1:n
                H[i,k] = H[i,k] - 2 * hv[i] * v[k]
                Q[i,k] = Q[i,k] - 2 * qv[i] * v[k]
            end
        end
    end
end

const n = 100
const A = rand(n, n)
const H = zeros(n, n)
const Q = zeros(n, n)
const u = zeros(n)
const v = zeros(n)

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

# 2.246 ms (3 allocations: 2.63 KiB)