disp(' ');
disp(' QR_Householder_eigen_mo.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=sparse(A);
[m,n]=(size(A));
%
Qt=eye(m);
%
R_rec=zeros(n,n);
A_rec=zeros(n,n);
Q_rec=zeros(n,n);
%
error_rec=1.0e+20;
%
for(ijk=1:100)
tic
%
    disp(' ref 1');
    clear Q;
    clear R;
%
    A=full(A);
    Am=max(abs(diag(A)));
        Am=Am*1.0e-12;
        nq=max(size(A));
        for(i=1:nq)
            for(j=i+1:nq)
                if(abs(A(i,j))<Am)
                    A(i,j)=0;
                    A(j,i)=0;
                end
            end
        end
    A=sparse(A);
%
    disp(' ref 2');
    [Q,R] = qr(A);
    disp(' ref 3');
%
    Q=full(Q);
    Am=max(abs(diag(Q)));
        Am=Am*1.0e-12;
        nq=max(size(Q));
        for(i=1:nq)
            for(j=i+1:nq)
                if(abs(Q(i,j))<Am)
                    Q(i,j)=0;
                    Q(j,i)=0;
                end
            end
        end    
%
    Q=sparse(Q);
       disp(' ref 4');
%
    new_matrix=Q'*original_A*Q;
%
%
    new_matrix=full(new_matrix);
    Am=max(abs(diag(new_matrix)));
        Am=Am*1.0e-12;
        nq=max(size(new_matrix));
        for(i=1:nq)
            for(j=1:nq)
                if(abs(new_matrix(i,j))<Am)
                    new_matrix(i,j)=0;
                end
            end
        end 
        new_matrix=sparse(new_matrix);
%
    A=new_matrix;
    R=A;
    original_A=new_matrix;
%
   disp(' ref 5');
    if(ijk==1)
        Qo=Q;
        Ro=R;
    end
%
    eig=full(diag(full(A)));
%
    if(ijk>=2)
           disp(' ref 5');
       for(i=1:n) 
          if(abs(eigb(i))>1.0e-20)         
             error(i)=full(abs((eig(i)-eigb(i))/eigb(i)));
          else
             error(i)=0.; 
          end
       end
%
       if(max(full(error))<error_rec)
          error_rec=max(full(error));
          eig_rec=eig;
          R_rec=R;
          A_rec=A;
          Q_rec=Q;
       end
%
       out1=sprintf('%ld  max error=%8.4g',ijk,max(full(error)));
       disp(out1);
%
       if( max(error) <0.00001)
          break;
       end 
    end
    eigb=eig;
%
toc
disp(' ');
end 
%
%%disp(' ');
%%disp(' Original ');
%%Qo;
%%Ro'
%
%%disp(' Last Q matrix ');
%%Q;
%
clear eig;
clear fn;
%
Qt=full(Qt);
R=full(R_rec);
eig=abs(diag(R_rec));
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear EEE;
clear max;
szz=max(size(eig));
for(i=1:szz)
    EEE(i,1)=eig(i);
end
EEE(:,2:szz+1)=Qt';
EEE=sortrows(EEE);
Eigenvalues=EEE(:,1);
Qt=EEE(:,2:szz+1)';
%
omega = full(sqrt(Eigenvalues));
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fn=omega/(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