function [R, U] = hqr(A)
%HQR Householder reflections for QR decomposition.
% [R, U] = hqr(A) returns
% R, the upper triangular factor, and
% U, the reflector generators.    
    H = @(u,x) x - 2*u*(u'*x);
    [n,~] = size(A);
    U = zeros(n,n);
    R = A;
    for j = 1:n
        u = hreflect(R(j:n, j));
        U(j:n, j) = u;
        R(j:n, j:n) = H(u, R(j:n,j:n));
        R(j+1:n, j) = 0;
    end
end