disp(' ');
disp(' QR_Householder_eigen.m   ver 1.7   December 26, 2010');
disp(' By Tom Irvine  Email: tomirvine@aol.com');
disp('  ');
disp(' This script calculates the eigenvalues and eigenvectors ');
disp(' of a matrix via the QR Householder method. ');
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;
clear Qo;
clear Ro;
clear eig;
clear eigb;
clear Qt;
clear max;
%
disp(' Enter the name of the matrix ');
%
A=input(' ')
original_A=A;
[m,n]=(size(A));
%
Qt=eye(m);
III=eye(m);
%
for(ijk=1:100)
%
    clear Q;
    clear R;
%
    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];
%
       Qt=Qt*Qk;  
       A=Qk*A;
       Q=Q*Qk;
    end
%
    new_matrix=Q'*original_A*Q;
    A=new_matrix;
    R=A;
    original_A=new_matrix;
%
    if(ijk==1)
        Qo=Q;
        Ro=R;
    end
%
    for(i=1:n)
       eig(i)=A(i,i);
    end
%
    if(ijk>10)
       for(i=1:n) 
          if(abs(eigb(i))>1.0e-20)         
             error(i)=abs((eig(i)-eigb(i))/eigb(i));
          else
             error(i)=0.; 
          end
       end
       if( max(error) <0.0001)
          ijk
          break;
       end
    end    
%
    eigb=eig;
%
end   
%
disp(' ');
disp(' Original ');
Qo
Ro
%
disp(' Last Q matrix ');
Q
%
clear eig;
clear fn;
%
for(i=1:n)
    eig(i)=R(i,i);
end
eig=sort(eig);
disp(' ')
disp(' Eigenvalues ');
disp(' ');
for(i=1:n)
    out1=sprintf(' %10.5g ',eig(i));
    disp(out1);
end
fn=sqrt(eig)/(2*pi);
disp(' ')
disp(' Natural Frequencies (Hz) ');
disp(' ');
for(i=1:n)
    out1=sprintf(' %10.5g ',fn(i));
    disp(out1);
end
%
disp(' ');
disp(' Display eigenvectors?  1=yes 2=no ');
ie=input(' ');
disp(' ');
if(ie==1)
    disp(' Eigenvectors ');
    Qt
end