disp(' ');
disp('  Sturm_Sequence_alt.m   ver 1.2   September 8, 2012 ');
disp(' ');
disp('  by Tom Irvine  ');
disp(' ');
disp('  This script calculates the eigenvalues for a generalized ');
disp('  eigenvalue problem using the Sturm sequence method. ');
disp(' ');
disp('  The mass and stiffness matrices must be symmetric.  ');
disp('   ');
%
close all;
%
clear K;
clear M;
clear n;
clear fm;
clear ffm;
clear L;
clear D;
clear A;
clear Y;
clear Minv;
clear Kinv;
clear Xhat;
%
%
TP=2*pi;
%
% Sample Data Set 1
% K=[2 -1 0; -1 4 -1; 0 -1 2];
% M=[0.5 0 0; 0 1 0; 0 0 0.5];
% Sample Data Set 2
% K=[12 6 -12 6; 6 4 -6 2; -12 -6 12 -6; 6 2 -6 4];
% M=[156 22 54 -13; 22 4 13 -3; 54 13 156 -22; -13 -3 -22 4];
%
disp(' ');
M=input(' Enter the mass matrix name:  ');
disp(' ');
K=input(' Enter the stiffness matrix name:  ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=max(size(M));
%
disp(' ');
ntrials=input(' Enter the number of trials ');
%
disp(' ');
disp(' Enter number of modes of interest ');
nint=input(' ');
%
if(nint>n)
    n=nint;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:n
        for j=i+1:n
%
            if(abs(M(j,i))<1.0e-10)
                M(j,i)=0;
            end
%
            if(abs(K(j,i))<1.0e-10)
                K(j,i)=0;
            end            
%
            M(i,j)=M(j,i);
            K(i,j)=K(j,i);       
        end
    end
%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
tic
disp(' begin standard conversion ');
[Kstandard,lambdashift]=convert_general_standard(M,K);
toc
%
tic
disp(' begin tridiagonal conversion ');
[Ktridiagonal]=function_Householder_tridiagonal_alt(Kstandard);
toc
%
X=zeros(n,1);
Y=zeros(n,1);
%
num=10;
%
for i=1:n
    Y(i)=rand;
end
Y=Y-0.5;
%
Minv=pinv(M);
%
for i=1:num;
    Xhat=Minv*Y;
    Yhat=K*Xhat;
    rho=Xhat'*Yhat/(Xhat'*Y);
    Y=Yhat/sqrt(Xhat'*Y);
end
%
lambda_max=rho;
fnmax=sqrt(lambda_max)/TP;
%
% disp(' ');
% out1=sprintf(' The maximum eigenvalue = %8.4g ',rho);
% disp(out1);
% disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear Khat
%
iv=(1*2*pi)^2;
Khat=K+iv*M;
%
n=max(size(M));
%
X=zeros(n,1);
Y=zeros(n,1);
%
for i=1:n
    Y(i)=rand;
end
Y=Y-0.5;
%
Kinv=pinv(Khat);
%
for i=1:num;
    Xhat=Kinv*Y;
    Yhat=M*Xhat;
    rho=Xhat'*Yhat/(Xhat'*Y);
    Y=Yhat/sqrt(Xhat'*Y);
end
%
rho=1/rho;
%
rho=rho-iv;
%
lambda_min=rho;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
out1=sprintf(' %d eigenvalues below  lambda=%12.4g  or fn=%8.4g Hz',n,lambda_max,fnmax);
disp(out1);
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
M=eye(n,n);
%
K=Ktridiagonal;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
fm=zeros(n,2);
%
for i=1:n
    fm(i,2)=lambda_max;
end
%
mu=1;
%
tic
progressbar;
for ijk=1:ntrials
    progressbar(ijk/ntrials);
%
    ffm=max(fm(1:nint,2));
%
    lambda=(fm(mu,2)+fm(mu,1))/2;   
%
    if(ijk>4)
        mu=mu+1;
        if(mu>nint)
            mu=1;
        end
    end
%
%
    L=zeros(n,n);
    D=zeros(n,n);  
%    A=K-lambda*M;
    A=K;
    for i=1:n
        A(i,i)=A(i,i)-lambda;
    end
%
    for i=1:n-1
        A(i+1,i)=A(i,i+1);
    end
%
    [L, D] = ldlt_sytr(A);   
%
%
    num=0;
    nz=0;
%
    for i=1:n
        if(D(i,i)<0)
            num=num+1;
        end
        if(D(i,i)==0)
            nz=nz+1;
        end    
    end
%
%%    out1=sprintf(' %d eigenvalues below %12.4g  %8.4g',num,lambda,sqrt(lambda)/(2*pi));
%%    disp(out1);
%
    for i=1:num
        if(lambda<fm(i,2) && lambda>fm(i,1))
            fm(i,2)=lambda;
        end
    end
    for i=num+1:n
        if(lambda>fm(i,1) && lambda<fm(i,2))
            fm(i,1)=lambda;
        end
    end
%
end
%
progressbar(1);
toc
%
fm=fm-lambdashift;
%
disp(' ');
disp(' Natural Frequency Estimates ');
disp(' ');
disp(' n    lower    upper ');
%
for i=1:nint
    out1=sprintf(' %d  %8.4g %8.4g ',i,sqrt(fm(i,1))/TP,sqrt(fm(i,2))/TP);
    disp(out1);     
end