function [Q,R] = qrfact(A)
% QRFACT  QR factorization by Householder reflections.
% (demo only--not efficient)
% Input:
%   A      m-by-n matrix
% Output:
%   Q,R    A=QR, Q m-by-m orthogonal, R m-by-n upper triangular

[m,n] = size(A);
Q = eye(m);
R = A;
for k = 1:n
  z = R(k:m,k);
  v = [ -sign(z(1))*norm(z) - z(1); -z(2:end) ];
  P = eye(m-k+1) - 2*(v*v')/(v'*v);   % HH reflection
  R(k:m,:) = P*R(k:m,:);
  Q(k:m,:) = P*Q(k:m,:);
end

Q = Q';
R = triu(R);                          % enforce exact triangularity