%
%   subspace_iteration.m    ver 1.o  July 14, 2004
%   by Tom Irvine   Email: tomirvine@aol.com
%  
%   Solution of the Generalized Eigenvalue Problem
%
clear K;
clear M;
clear nnn;
%
%  Sample stiffness matrix K and mass matrix M
%
K=[ 2 -1 0 ; -1 3 -2 ; 0 -2 2 ]
M=[ 1 -0 0 ; 0 1 0   ; 0 0 2 ]
%
%  Sample trial vectors
%
U1(:,1)=[1 2 3]';
U1(:,2)=[1 2 -1]';
%
nnn=size(U1);  % number of columns
n1=nnn(1,1);
n2=nnn(1,2);
n=n2;
%
% Normalize trial eigenvectors with respect to mass matrix
%
    nn=U1'*M*U1;
%
for i=1:2
    U1(:,i)= U1(:,i)/sqrt(nn(i,i));
end
%
U2=K\(M*U1);
%
for i=1:10
%
    K2=U2'*K*U2;
    M2=U2'*M*U2;
%
    [P2,A2]=eig(K2,M2);
%   
    nn=P2'*M2*P2;
    for i=1:2
        P2(:,i)= P2(:,i)/sqrt(nn(i,i));
    end    
%
    Uhat2=U2*P2;
%
    U3=K\M*Uhat2;
    U2=U3;
end
disp(' Eigenvectors ')
Uhat2
disp(' Eigenvalues ')
A2