disp(' ');
disp(' QR_Householder.m   ver 1.3   April 8, 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(' ')
original_A=A;
%
[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
new_matrix=Q'*original_A*Q
%
disp(' ');
disp(' Repeat with new_matrix for eigenvalue analysis ');