m=8;n=6;
H=hilb(m); A=H(:,1:n); 
[AA,d]=HouseholderQR(A)
Q=[];
for i=eye(m),                                  % compute Q explicit
  z=HouseholderQy(AA,i); Q=[Q z];
end
R=triu(AA);
for i=1:n, R(i,i)=d(i); end                    % add diagonal to R
[norm(Q'*Q-eye(m)) norm(A-Q*R)] 
[q,r]=qr(A);                                   % compare with Matlab qr
[norm(q'*q-eye(m)) norm(A-q*r)]
