disp(' ');
disp(' Hessenberg_tridiagonal.m   ver 1.0   December 27, 2010');
disp(' By Tom Irvine  Email: tomirvine@aol.com');
disp('  ');
%
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
%
   tic
   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);
[Q,Ktridiagonal] = hess(full(K));
K=Ktridiagonal;
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