function [Qt R]=myQR(A)

% This program finds the full QR factorization of a mxn matrix A,
%    A = QR
% where Q is mxm orthogonal, and R is mxn upper-triangular, 
% using Householder reflections the big expensive way, actually 
% computing the Householder reflection matrices

[m n]=size(A);


R=A;
for i=1:(m-1)
  [Qhatx Qhat]=ZeroOut(R(i:end,i),1);
  if i==1
      Q=Qhat
      Qt=Q;
      
  else
      Q=[eye(i-1) zeros(1,n-i+1);  zeros(n-i+1,1)  Qhat]
      Qt=Qt*Q;
  end
  
  R=Q*R   % Building the upper triangular R
  
  
end


fprintf('\n ||A - QR|| = %e\n',norm(A-Qt*R))