disp(' ');
disp(' QR_Householder_inv.m   ver 1.1   March 16, 2010');
disp(' By Tom Irvine  Email: tomirvine@aol.com');
disp('  ');
%
%  Matlab already has a function:  [Q,R] = qr(AA)
%
%  This script demonstates the underlying steps.
%
%  Reference: http://www.caam.rice.edu/~caam453/lecture4.pdf
%
%     
%      A = Q R
%
clear A;
clear Q;
clear R;
%
disp(' Enter the name of the matrix ');
A=input(' ')
%
[m,n]=(size(A));
Q=eye(m);
%
for(k=1:min(n-1,n))
%    
    ak=A(k:end,k);
%
    norm_ak=norm(ak);
    sign_ak=ak(1)/abs(ak(1));
    ek=zeros(length(ak),1);
    ek(1)=1;
%
    vk=ak + sign_ak*norm_ak*ek;
%
    ident=eye(m-k+1);
%    
    Hk=ident - 2*vk*vk'/(norm(vk)^2);  % construct reflector  
%
    Qk=[eye(k-1)  zeros(k-1,m-k+1);  zeros(m-k+1,k-1)  Hk];
%
    A=Qk*A;
    Q=Q*Qk;
end
Q
R=A
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
I=eye(n,n);
clear B;
B=Q'*I;
clear AA;
AA=R;
%
for(nv=1:n)
%
    B(n,nv)=B(n,nv)/AA(n,n);
%
    for(j=1:n-1)
        i=n-j;
%              
        for(k=i+1:n)
%
	        B(i,nv)=B(i,nv)-B(k,nv)*AA(i,k);
%            
        end
        B(i,nv)=B(i,nv)/AA(i,i);
%                 
    end
end
disp(' ');
disp(' The inverse matrix is ');
B