disp(' ');
disp(' Householder_tridiagonal_alt.m   ver 1.6   February 21, 2011');
disp(' By Tom Irvine  Email: tomirvine@aol.com');
disp('  ');
%
%  References: 
%  Meirovitch p156
%  http://www.mpi-hd.mpg.de/astrophysik/HEA/internal/Numerical_Recipes/f11-
%  2.pdf
%
clear M;
clear K;
clear Kstandard;
clear v;
clear P;
clear MMA;
clear KKA;
clear n;
%
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));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%
%   LDLT decomposition to form standardized matrix
%
    clear A;
    clear L;
    clear U;
    clear D;
    A=MMA;
%
    L=zeros(n,n);
    U=zeros(n,n);
    D=zeros(n,n);
%
%  Triangular Decomposition
%
    L=eye(n,n);
%
    i=1;
%
    U(1,1)=A(1,1);
    for i=2:n
        L(i,1)=A(i,1)/U(1,1);
        U(1,i)=A(1,i);
    end
%
    for i=2:n
%
        sum=0.;
        for k=1:i-1
            sum=sum+U(k,i)*L(i,k);
        end
        U(i,i)=A(i,i)-sum;
%
        for j=i+1:n
%        
            sum=0.;
            for k=1:j-1
                sum=sum+U(k,i)*L(j,k);
            end
            L(j,i)=(A(j,i)-sum)/U(i,i);
%
            sum=0.;
            for k=1:j-1
                sum=sum+U(k,j)*L(i,k);
            end
            U(i,j)=(A(i,j)-sum);
        end
    end
%
% Solve for Diagonal Matrix D
%
    for i=1:n
        D(i,i)=U(i,i);
    end
    LT=L';
%
    disp(' ');
%
    n=max(size(MMA));
    QQQ=eye(n,n);
%
    Lm=L*sqrt(D);
    Lminv=Lm\QQQ;
    Kstandard=(Lminv*KKA*Lminv');
    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
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Symmetry check
%
   disp(' ');
   disp(' Symmetry check ');
   jflag=0;
   for(i=1:n)
       for(j=i+1:n)
           if(abs(K(i,j)-K(j,i))>abs((1.0e-06)*(K(i,j)+K(j,i))))
               out1=sprintf('\n K(%d,%d)=%8.4g   K(%d,%d)=%8.4g \n',i,j,K(i,j),j,i,K(j,i));
               disp(out1);
               diff=abs(K(i,j)-K(j,i));
               out2=sprintf(' diff=%8.4g \n',diff);
               disp(out2);
               jflag=1;
               break;
           end
       end
       if(jflag==1)
           break;
       end
   end
%
   if(jflag==1)
       disp(' Warning:  symmetry test fails ');
       disp(' ');
       qqqq=input(' press Enter to continue ');
   else
       disp(' Symmetry test passes ');
   end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Determine whether the matrix is already tridiagonal
%
disp(' 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 
%
    Kstandard=K
%      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Numerical Recipes method
%
%% tic
%
disp(' ');
disp(' Calculate eigenvalues?  1=yes  2=no ');
ic=input(' ');
%
disp(' ');
disp(' Calculate eigenvectors?  1=yes  2=no ');
ims=input(' ');
if(ims==1)
    ic=1;
end
%
IM=eye(n);
Q=IM;
%
A=K;
%
if(iflag==1)
    for(ijk=1:n-2)
%
        out1=sprintf(' step %d \n',ijk);
        disp(out1);
        clear P;
        clear r;
        sum=0;
%
        for(j=ijk+1:n)
            sum=sum+A(ijk,j)^2;
        end
        alpha=sqrt(sum);
%%        out1=sprintf('ijk=%d alpha= %8.4g ',ijk,alpha);
%%        disp(out1);   
%
        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
        for(j=ijk+2:n)
            v(j)=A(ijk,j)/(2*alpha*v(ijk+1));
        end
        v;
%
%%        P=IM-2*(v'*v)
%%        K=P*K*P
%
        if(ims==1)
            P=IM-2*(v'*v);
            Q=Q*P;
        end
%   
        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';
%
    end
    tol=max(max(abs(A)))*1.0e-10;
    for(i=1:n)
        for(j=1:n)
            if(abs(A(i,j))<tol && i~=j)
                A(i,j)=0;
            end
        end
    end
    A
    clear K;
    K=A;
else
    K
end
%% toc
%
if(ic==1)
    [fn,omega,ModeShapes,MST]=Standard_Eigen(K,3);
%
    if(im==1)
        if(ims==1)
            ModeShapes=Lminv'*(Q*ModeShapes);
        end
    end
end
%
clear Ktridiagonal;
Ktridiagonal=K;
fn=fn';
disp(' ');
disp(' Output arrays: ');
if(ic==1)
    disp('   fn ');
end
if(ims==1)
    disp('   ModeShapes ');
end
disp('   Kstandard ');
disp('   Ktridiagonal ');