function [W,A] = householder(A)
% HOUSEHOLDER computes a QR factorization upper triangular matrix R 

[m,n] = size(A);

ident = eye(m);
W = zeros(m,n);

for k=1:n
    e1=ident(k:m,k);
    x = A(k:m,k);
    W(k:m,k)=( (x(1)>=0) - (x(1)<0) )*norm(x)*e1+x;
    W(k:m,k)=W(k:m,k)/norm(W(k:m,k));
    A(k:m,k:n)=A(k:m,k:n)-2*W(k:m,k)*W(k:m,k)'*A(k:m,k:n);
end