disp(' ');
disp(' Householder_tridiagonal_mo.m   ver 1.6   December 25, 2010');
disp(' By Tom Irvine  Email: tomirvine@aol.com');
disp('  ');
%
%  Reference: Meirovitch p156
%
clear M;
clear K;
clear v;
clear P;
clear MMA;
clear KKA;
clear n;
clear Kstandard;
clear Q;
%
disp(' Select input method ');
disp('  1=mass & stiffness matrices ');
disp('  2=standardized stiffness matrix ');
im=input(' ');
%
if(im==1)
    disp(' The input matrices must be real, symmetric. ');
    disp(' ');
    disp(' Enter the name of the mass matrix ');
    MMA=input(' ');
    disp(' Enter the name of the stiffness matrix ');
    KKA=input(' ');  
    n=max(size(KKA));
    MMS=sparse(MMA);
    KKS=sparse(KKA);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%
%   Cholesky decomposition to form standardized matrix
%
tic
   disp(' ');
   disp(' Cholesky ');
%
%
    clear R;
    clear Rinv;
    R=sparse(chol(MMS));
    Rinv=sparse(inv(R'));
    Kstandard=(Rinv*KKS*Rinv');    
    K=Kstandard;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
    disp(' The input matrix must be real, symmetric. ');
    disp(' ');
    disp(' Enter the name of the standardized matrix ');
    K=input(' ');
    n=max(size(K));
end
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Symmetry check
%
   disp(' ');
   disp(' Symmetry check ');
   jflag=0;
   Kmax=max(max(abs(K)));
   KK=abs(K-K');
   KKmax=max(max(KK));
%
   if(KKmax>1.0e-12*Kmax)
       disp(' Warning:  symmetry test fails ');
       disp(' ');
       qqqq=input(' press Enter to continue ');
   else
       disp(' Symmetry test passes ');
   end
   toc
   tic
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Determine whether the matrix is already tridiagonal
%
		   iflag=0;
%        	
           if(iflag==0)
%           
		        for(i=1:n)
%		        
                    for(j=i+2:n)
%                                                  
                         if(abs(K(i,j))>1.0e-18)
 %                        
                              iflag=1;
                              break;                       
                         end      
                    end
                    if(iflag==1)                  
                         break;
                    end                       
                end      	
           end  
%      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
toc
tic
Kstandard=K;
K=sparse(K);
IM=eye(n);
Q=sparse(IM);
%
A=sparse(K);
%
        A=full(A);
        Am=max(abs(diag(A)));
        Am=Am*1.0e-12;
        progressbar(1);
        for(i=1:n)
            progressbar(i/n)
            for(j=i:n)
                if(abs(A(i,j))<Am && i~=j)
                    A(i,j)=0;
                end
            end
        end
        A=sparse(A);
%
disp(' begin tridiagonal ');
progressbar(1);
if(iflag==1)
    for(ijk=1:n-2)
        progressbar(ijk/(n-2)) % Update figure  
        clear P;
        clear r;
%
        alpha=norm(A(ijk,ijk+1:n));
%
        v(1:ijk)=0;
        ccc=1+ (A(ijk,ijk+1)/alpha);
        v(ijk+1)=sqrt(0.5*ccc);
        if(abs(ccc)<1.0e-20)
            out1=sprintf('\n Warning:  ccc=%8.5g \n',ccc);
            disp(out1);
        end
%
        v(ijk+2:n)=A(ijk,ijk+2:n)/(2*alpha*v(ijk+1));
        v;
%
%%        P=IM-2*(v'*v)
%%        K=P*K*P
%   
        u=sqrt(2)*v';
        H=0.5*norm(u)^2;
%        
        p=A*u/H;
        k=u'*p/(2*H);
        q=p-k*u;
%       
        A=A-(q*u'+u*q');
%
        Am=max(abs(diag(A)));
        Am=Am*1.0e-12;
        nq=max(size(A));
        for(i=ijk:ijk)
            for(j=i+1:nq)
                if(abs(A(i,j))<Am)
                    A(i,j)=0;
                    A(j,i)=0;
                end
            end
        end
%
end
    K=A;
else
    K;
end
Ktridiagonal=full(K);
toc
%
disp(' ');
disp(' Calculate eigenvalues?  1=yes  2=no ');
ic=input(' ');
%
if(ic==1)
    K=full(K);
    [fn,omega,ModeShapes,MST]=Standard_Eigen(K,3);
    %
    if(im==1)
        disp(' ');
        disp(' Display Mode Shapes?  1=yes  2=no ');
        ims=input(' ');
        if(ims==1)
            ModeShapes=Lminv'*(Q*ModeShapes)
        end
    end
end