function [QQ,RR] = HouseholderQR(A)
% Householder QR factorization to
% transform an m by n matrix A with m >= n
% to upper triangular form.
% OUTPUT: an m by m orthogonal matrix QQ
% & an m by n upper triangular matrix RR
% such that QQ*RR = A.
% Algorithm 3.1 on p.124 Heath
% Example 3.12 & 13 on p.135 & 136 Heath
% K. Ming Leung, 04/6/03

[m,n] = size(A); QQ=eye(m);
for k = 1:n
    v = zeros(m,1);
    alpha = -sign(A(k,k))*norm(A(k:m,k));
    v(k:m) = A(k:m,k); v(k) = v(k) - alpha;
    beta = (v'*v)/2;
    H=eye(m)-v*v.'/beta;
    QQ = QQ*H;
    if (beta~=0)
        % for j=k
        gamma = (v'*A(:,k))/beta;
        A(k,k) = A(k,k) - gamma*v(k);
        A(k+1:m,k) = zeros(m-k,1);
        % for the rest of the columns
        for j = k+1:n
            gamma = (v'*A(:,j))/beta;
            A(:,j) = A(:,j) - gamma*v;
        end
    end    % end if
end        % end for

RR=A;