function [W,R]=house_fact(X)
%
%  Householder factorization of a matrix X
%   
%  Input
%   X- matrix to factored
%
%  Output
%   W - matrix of Householder vectors
%   R - upper triangular factor
%
   [m,n]=size(X);
   p=min([m n]);
%
%  Although I did not require it, this method allows for m < n
%
    W=zeros(m,p); R=zeros(p,n);
    for k=1:p
%
%  Compute Householder vector and diagonal element
%
      [W(k:m,k),X(k,k)]=house(X(k:m,k));
%
%  Alter X
%
      f=2*X(k:m,k+1:n)'*W(k:m,k);
       X(k:m,k+1:n)=X(k:m,k+1:n)-W(k:m,k)*f';
%
%  Update finished part of R
%
       R(k,k:n)=X(k,k:n);       
    end;