function [Q,R]=house(A);
% HOUSE computes QR factors by Householder triangularization.
%    Essentially Trefethen & Bau Algorithm 10.1 along with 10.3 to form Q.
%    See CLGS, MGS for Gram-Schmidt versions and QR for Matlab built in version.
%
%    [Q,R] = house(A)  produces Q with orthonormal columns and R upper
%    triangular so that A = QR.  For m by n A with m>=n with full rank.
%    Returned Q is m by n, R is n by n, that is, "economy size".
%ELB 10/03

myeps=2*eps; % use machine eps times 2
[m n]=size(A);
scal=max(max(abs(A)));
if scal==0, Q=eye(m,n); R=zeros(n,n); return, end % case of zero matrix

%alg 10.1
for k=1:n
    clear v
    v=A(k:m,k); %size(v)=[m-k+1 1]
    if v(1,1)<0
        v(1,1)=-norm(v,2)+v(1,1);
    else
        v(1,1)=norm(v,2)+v(1,1);
    end
    r=norm(v,2);
    if abs(r)/scal<myeps, error('Householder reflector nearly zero'), end
    v=v/r;
    V(:,k)=[zeros(k-1,1); v]; 
    A(k:m,k:n)=A(k:m,k:n)-2*v*(v'*A(k:m,k:n));
end    
R=A(1:n,:);

%alg 10.3 to get Q which is m by n
for i=1:n
    e=zeros(m,1); e(i)=1;
    for k=n:-1:1
        v=V(k:m,k);
        e(k:m)=e(k:m)-2*v*(v'*e(k:m));
    end
    Q(:,i)=e;
end