function [Q, R] = hqr(A)
%HQR Compute the full QR factorization of A using
%Householder reflections.
%
%   [Q R] = hqr(A)computes the QR factorization of the
%   m x n matrix A. There are no restrictions on the relative
%   sizes of m and n. Q is an orthogonal m x m matrix, and R
%   is an m x n upper triangular matrix R such that A = QR

[m, n] = size(A);

R = A;
Q = eye(m);
k = min(m-1, n);

for i=1:k
   [R(i:m,i:n), u] = hzero1(R(i:m,i:n));
   % post multiply by Hu
   if norm(u) ~= 0
      Q(1:m,i:m) = Q(1:m,i:m)-(2/norm(u)^2)*Q(1:m,i:m)*(u*u');
   end
end
