% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-12
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

Y1 = [ 1 1 1;
       0 1 1 ;
       2 2 -2; 
       2 0 1; 
       1 -2 -1];
M1 = size(Y1,1); 

Y2 = [ 2 -5 -2;
       1 -4 1;
       -1 -5 0;
       2 -1 2;];
M2 = size(Y2,1); 

% The combined patterns:
Y = [ Y1; Y2 ]; 
M = M1 + M2;

% The original dimension of the pattern space:
m = size(Y1,2); 

% Compute matrix B (there must be a better way): 
B = zeros(m,m); 
for jj=1:m,
  for kk=1:m,
    
    tmp=0;
    for q=1:M1,
      for p=1:M2,
	tmp=tmp + (Y1(q,jj)-Y2(p,jj))*(Y1(q,kk)-Y2(p,kk));
      end
    end
    
    B(jj,kk)=tmp/(M1*M2); 
    
  end
end

% Compute the matrix C (there ...): 
C = zeros(m,m); 
for jj=1:m,
  for kk=1:m,
    
    tmp=0;
    for q=1:M,
      for p=1:M,
	tmp=tmp + (Y(q,jj)-Y(p,jj))*(Y(q,kk)-Y(p,kk));
      end
    end
    
    C(jj,kk)=2*tmp/(M*(M-1)); 
    
  end
end

% Compute the matrix B * C^{-1}:
A = B * inv(C); 

% Find the eigenstructure: 
[V,D] = eig(A)