function x = householder(A,b)

% function x = householder(A,b)
% solution of the matrix equation Ax=b 
% via QR-factorisation

[Q,R,I]=qr_factor(A); 
x=r_solve(R,Q'*b);
x(:,I)=x;