%   [W,R] = house(A)
%
%   Description:
%       Compute implicit representation of a full QR factorization of the
%   matrix A = QR for an m x n matrix A using Householder reflections.
%
%   Inputs:
%       A       - m x n matrix
%   Outputs:
%       W       - m x n matrix of Householder reflectors
%       R       - upper triangular m x n matrix
%       Q (optional) - Q matrix
%          this version minimizes use of square roots
function [W,R,a,b] = house(A)
% function [W,R,Q] = house(A)
flops(0);tic
% save input
R = A;
% determine size
[m,n] = size(A);
% Initialize matrix to hold the Householder Reflectors
W = zeros(m,n);
% Perform upper triangularization while saving the reflectors
for k = 1:n,
   % Pick out 1st column of kth sub-matrix
   x = R(k:m,k);
   % Compute the reflector
   s=sign(x(1))*norm(x);
   x(1) = s + x(1);
% Normalize the reflector
   nn=norm(x);
%  nn=sqrt(W(k:m,k)'*W(k:m,k));
%  W(k:m,k) = x / nn;
   x = x / nn;
%% nn = x'*x; n2=2/nn;
   % Compute I-2P of kth sub-matrix
   if k~=n,
      R(k:m,k+1:n) = R(k:m,k+1:n) - 2*x*(x'*R(k:m,k+1:n));
%     R(k:m,k+1:n) = R(k:m,k+1:n) - 2*W(k:m,k)*(W(k:m,k)'*R(k:m,k+1:n));
%%    R(k:m,k+1:n) = R(k:m,k+1:n) - n2*x*(x'*R(k:m,k+1:n));
   end
%% n1=nn^(1/2);
%% W(k:m,k) = x/n1;
   W(k:m,k) = x;
   R(k:m,k)=0; R(k,k) = -s;
end;
%------------------- formQ logic
%Q = eye(m);
%flops(0);tic;
%for k = n:-1:1
%   Q(k:m,k:m) = Q(k:m,k:m) - 2*W(k:m,k) * (W(k:m,k)' * Q(k:m,k:m));
%end
%flops; toc;

a=flops; b=toc;
   return
