disp(' ');
disp('   subspace_iteration.m    ver 1.2  Sep 5, 2012 ');
disp('   by Tom Irvine   Email: tomirvine@aol.com ');
disp('  ');
disp('   Solution of the Generalized Eigenvalue Problem. ');
%
clear f;
clear eee;
clear A2;
clear K;
clear M;
clear U1;
clear Uhat2;
clear P2;
clear nnn;
clear Ev;
clear fn;
%
disp(' ');
disp(' Enter the mass matrix name ');
M=input(' ');
disp(' ');
disp(' Enter the stiffness matrix name ');
K=input(' ');
%
disp(' ');
disp(' Enter number of vectors ');
nv=input(' ');
%
disp(' ');
disp(' Enter number of iterations ');
niter=input(' ');
%
%  Sample trial vectors
%
%
ndof=max(size(M));
%
% form trial vectors
%
U1=4*(0.5-rand(ndof,nv)); 
%
% Normalize trial eigenvectors with respect to mass matrix
%
    nn=U1'*M*U1;
%
for i=1:nv
    U1(:,i)= U1(:,i)/sqrt(nn(i,i));
end
%
Kinv=pinv(K);
KinvM=Kinv*M;
%
U2=KinvM*U1;
%
for ijk=1:niter
%
    K2=U2'*K*U2;
    M2=U2'*M*U2;
%
    [P2,A2]=eig(K2,M2);
%   
%   normalize with respect to mass matrix
%
    nn=P2'*M2*P2;
    for i=1:nv
        P2(:,i)= P2(:,i)/sqrt(nn(i,i));
    end    
%
    Uhat2=U2*P2;
%
    U3=KinvM*Uhat2;
    U2=U3;
end
%
%% disp(' Eigenvalues ')
%% A2
%
Ev=zeros(nv,1);
fn=zeros(nv,1);
%
[Ev,Uhat2]=sort_eigenvalues_eigenvectors(A2,Uhat2);
%
for i=1:nv
    fn(i)=sqrt(Ev(i))/(2*pi);
end
%
disp(' ');
disp(' n  eigenvalue  fn(Hz) ');
for i=1:nv
    out1=sprintf('%d  %8.4g  %8.4g ',i,Ev(i),fn(i));
    disp(out1);
end
%
disp(' ');
disp(' Display eigenvectors?  1=yes 2=no ');
ide=input(' ');
if(ide==1)
    disp(' ');
    disp(' Eigenvectors ')
    Uhat2
end