function [X,E] = polyeig(varargin)
%POLYEIG Polynomial eigenvalue problem.
%   [X,E] = POLYEIG(A0,A1,..,Ap) solves the polynomial eigenvalue problem
%   of degree p:
%       (A0 + lambda*A1 + ... + lambda^p*Ap)*x = 0.
%   The input is p+1 square matrices, A0, A1, ..., Ap, all of the same
%   order, n.  The output is an n-by-n*p matrix, X, whose columns   
%   are the eigenvectors, and a vector of length n*p, E, whose
%   elements are the eigenvalues.
%       for j = 1:n*p
%          lambda = E(j)
%          x = X(:,j)
%          (A0 + lambda*A1 + ... + lambda^p*Ap)*x is approximately 0.
%       end
% 
%   E = POLYEIG(A0,A1,..,Ap) is a vector of length n*p whose
%   elements are the eigenvalues of the polynomial eigenvalue problem.
%
%   Special cases:
%       p = 0, polyeig(A), the standard eigenvalue problem, eig(A).
%       p = 1, polyeig(A,B), the generalized eigenvalue problem, eig(A,-B).
%       n = 1, polyeig(a0,a1,..,ap), for scalars a0, ..., ap, 
%       is the standard polynomial problem, roots([ap .. a1 a0])
%
%   If both A0 and Ap are singular the problem is potentially ill-posed.
%   Theoretically, the solutions might not exist or might not be unique.
%   Computationally, the computed solutions may be inaccurate.  An attempt
%   is made to detect this situation, and a warning message may result.
%   If one, but not both, of A0 and Ap is singular, the problem is well
%   posed, but some of the eigenvalues may be zero or "infinite".

%   Mofified: Nicholas J. Higham and Francoise Tisseur, 3-3-00.
%   C. Moler, 5-5-93, 12-29-93.
%   Copyright 1984-2000 The MathWorks, Inc.
%   $Revision: 5.9 $  $Date: 2000/06/01 02:04:19 $
 
% Build two n*p-by-n*p matrices:
%    A = [A0   0   0   0]   B = [-A1 -A2 -A3 -A4]
%        [ 0   I   0   0]       [  I   0   0   0]
%        [ 0   0   I   0]       [  0   I   0   0]
%        [ 0   0   0   I]       [  0   0   I   0]

n = length(varargin{1});
p = nargin-1;
A = eye(n*p);
A(1:n,1:n) = varargin{1};
if p == 0
   B = eye(n);
   p = 1;
else
   B = diag(ones(n*(p-1),1),-n);
   j = 1:n;
   for k = 1:p
      B(1:n,j) = - varargin{k+1};
      j = j+n;
   end
end

% Use the QZ algorithm on the big matrix pair (A,B).
if nargout > 1
   [alpha,beta,Q,Z,X] = qz(A,B,'complex');
else
   [alpha,beta] = qz(A,B,'complex');
end

% Extract the eigenvalues.
alpha = diag(alpha);
beta = diag(beta);
atol = 100*n*max(abs(alpha))*eps;
btol = 100*n*max(abs(beta))*eps;
if any(abs(alpha) < atol & abs(beta) < btol)
   wrnstr = sprintf(['Rank deficient generalized eigenvalue problem.\n' ...
            '         Eigenvalues are not well determined.  Results may be inaccurate.']);
   warning(wrnstr);
end
E = alpha./beta;

if nargout <= 1, X = E; return, end
if p == 1, return; end

% For each eigenvalue, extract the eigenvector from whichever portion
% of the big eigenvector matrix X gives the smallest normalized residual.
V = zeros(n,p);
% Division by zero possible for zero eigenvalues, so turn off warnings.
warns = warning; warning('off');
for j = 1:p*n
   V(:) = X(:,j);
   R = varargin{p+1};
   for k = p:-1:1
       R = varargin{k} + E(j)*R;
   end
   R = R*V;
   res = sum(abs(R))./ sum(abs(V));  % Normalized residuals.
   [s,ind] = min(res);
   X(1:n,j) = V(:,ind)/norm(V(:,ind));  % Eigenvector with unit 2-norm.
end
X = X(1:n,:);
warning(warns);
