function d = GenEig(A,B)
% A and B are n-by-n nonsingular matrices.
% d is a column n-vector consisting of the eigenvalues of the matrix
% M = (A*inv(B) + B*inv(A))/2

[TA,TB] = qz(A,B,'real');
[n,n] = size(A);
k = 1;
d = [];
while k<n
if TA(k+1,k)==0
mu = TA(k,k)/TB(k,k);
d = [d ; (mu + 1/mu)/2];
k = k+1;
else
mu = eig(TA(k:k+1,k:k+1),TB(k:k+1,k:k+1));
d = [d ; (mu + 1./mu)/2];
k = k+2;
end
end
if k==n
d = [d; (TA(k,k)/TB(k,k) + TB(k,k)/TA(k,k))/2];
end