function [X,E,s] = polyeig(varargin)
%MATLAB Code Generation Library Function

%   Limitations:
%   * The basis of eigenvectors may be different, and the eigenvalues for
%     real inputs will not generally be sorted so that complex conjugate
%     pairs are adjacent in E.
%   * Differences in eigenvectors and ordering of eigenvalues may also lead
%     to differences in condition number outputs, s.

%   Copyright 1984-2020 The MathWorks, Inc.
%#codegen

coder.columnMajor;
coder.inline('never');
narginchk(1,inf);
coder.internal.assert(nargout <= 2 || nargin >= 2, ...
    'MATLAB:polyeig:tooFewInputs');
ONE = coder.internal.indexInt(1);
TWO = coder.internal.indexInt(2);
n = coder.internal.indexInt(size(varargin{1},1));
coder.internal.assert(size(varargin{1},2) == n,'MATLAB:square');
p = coder.internal.indexInt(nargin - 1);
np = n*p;
xZERO = coder.internal.scalarEg(varargin{:});
xONE = cast(1,'like',xZERO);
rZERO = real(xZERO);
%   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]
if nargin == 1
    A = varargin{1};
    B = eye(n,'like',xZERO);
    p = ONE;
else
    A = zeros(np,'like',xZERO);
    for j = 1:n
        for i = 1:n
            A(i,j) = varargin{1}(i,j);
        end
    end
    for j = n+1:np
        A(j,j) = 1;
    end
    B = zeros(np,'like',xZERO);
    for k = coder.unroll(TWO:nargin)
        coder.internal.assert(size(varargin{k},1) == n && ...
            size(varargin{k},2) == n,'MATLAB:dimagree');
        offset = (k - 2)*n;
        for j = 1:n
            for i = 1:n
                B(i,j + offset) = -varargin{k}(i,j);
            end
        end
        for j = 1:np-n
            B(j+n,j) = 1;
        end
    end
end
% Use the QZ algorithm on the big matrix pair (A,B).
if nargout == 1
    X = eig(A,B);
    return
end
if p >= 2
    [XTall,E] = eig(A,B,'vector');
    % For each eigenvalue, extract the eigenvector from whichever portion
    % of the big eigenvector matrix X gives the smallest normalized
    % residual.
    X = coder.nullcopy(zeros(n,np,'like',XTall));
    for j = 1:np
        % Evaluate the polynomial at lambda = E(j).
        R = evaluatePolynomial(E(j),varargin{:});
        % Find the offset to the eigenvector with the minimum residual.
        offset = findMinResidOffset(XTall,j,R,p);
        % Normalize the eigenvector.
        nrmV = coder.internal.blas.xnrm2(n,XTall, ...
            (j - 1)*np + offset + ONE,ONE);
        % Eigenvector with unit 2-norm.
        for i = 1:n
            X(i,j) = XTall(offset + i,j)/nrmV;
        end
    end
else
    [X,E] = eig(A,B,'vector');
end
if nargout == 3
    % Construct matrix Y whose rows are conjugates of left eigenvectors.
    rcond_p = rcond(varargin{p + 1});
    rcond_0 = rcond(varargin{1});
    coder.internal.assert(eml_max(rcond_p,rcond_0) > eps(class(xZERO)), ...
        'MATLAB:polyeig:nonSingularCoeffMatrix');
    if rcond_p >= rcond_0
        V = cast(varargin{p + 1},'like',X);
        E1 = E;
    else
        V = cast(varargin{1},'like',X);
        E1 = 1./E;
    end
    % Y is a stack of p copies of X with the entries in column j scaled by
    % successive powers of E1(j), starting with E1(j)^0 and ending with
    % E1(j)^(p-1).
    Y = coder.nullcopy(zeros(np,'like',X));
    for j = 1:np
        for i = 1:n
            Y(i,j) = X(i,j);
        end
        for k = 1:p - 1
            offset = n*k;
            for i = 1:n
                Y(offset + i,j) = Y(offset - n + i,j)*E1(j);
            end
        end
    end
    % W is an n*p-by-n matrix of all zeros except that the lowermost
    % n-by-n submatrix is the identity.
    W = zeros(np,n,'like',rZERO);
    offset = np - n;
    for j = 1:n
        W(offset + j,j) = 1;
    end
    % Compute the left eigenvectors.
    W = Y\W/V;
    % Normalize the left eigenvectors.
    for i = 1:np
        nrm = coder.internal.blas.xnrm2(n,W,i,np);
        for j = 1:n
            W(i,j) = W(i,j)/nrm;
        end
    end
    % Reconstruct alpha-beta representation of eigenvalues: E(i) =
    % a(i)/b(i).
    a = E;
    b = ones(size(a),'like',E);
    for k = 1:n
        if isinf(a(k))
            a(k) = 1;
            b(k) = 0;
        end
    end
    % Calculate the Frobenius Norm of the input matrices.
    nu = zeros(p + 1,1,'like',rZERO);
    for i = 1:p + 1
        nu(i) = norm(varargin{i},'fro');
    end
    s = zeros(np,1,'like',rZERO);
    % Compute condition numbers.
    ab = coder.nullcopy(zeros(1,p + 1,'like',E));
    dp = cast(p,'like',rZERO);
    nn = n*n;
    for j = 1:np
        % ab = (a(j).^(0:p-1)).*(b(j).^(p-1:-1:0));
        for k = 1:p
            dk = cast(k,'like',a);
            ab(k) = (a(j).^(dk - 1)).*(b(j).^(dp - dk));
        end
        Da = ab(1)*varargin{2};
        Db = dp*ab(1)*varargin{1};
        for k = TWO:p
            dk = cast(k,'like',rZERO);
            Da = Da + dk*ab(k)*varargin{k + 1};
            Db = Db + (dp - dk + 1)*ab(k)*varargin{k};
        end
        % nab = norm((a(j).^(0:p)).*(b(j).^(p:-1:0)).* nu');
        for k = 1:p + 1
            dk = cast(k,'like',rZERO);
            ab(k) = (a(j)^(dk - 1))*(b(j)^(dp + 1 - dk))*nu(k);
        end
        % Now calculate W(j,:)*(conj(b(j))*Da - conj(a(j))*Db)*X(:,j) using
        % BLAS calls.
        caj = conj(a(j));
        cbj = conj(b(j));
        % Compute conj(b(j))*Da - conj(a(j))*Db and store it in Da.
        for k = 1:nn
            Da(k) = cbj*Da(k) - caj*Db(k);
        end
        % Compute W(j,:)*Da and store the result in Db(1:n).
        % Here we formulate this as Db(1:n) = Da.'*W(j,:).'
        Db = coder.internal.blas.xgemv('T',n,n,xONE,Da,ONE,n, ...
            W,j,np,xZERO,Db,ONE,ONE);
        % Compute d = Db(1:n).'*X(:,j);
        offset = 1 + (j - 1)*n;
        d = coder.internal.blas.xdot(n,Db,ONE,ONE,X,offset,ONE);
        s(j) = norm(ab)/abs(d);
    end
end

%--------------------------------------------------------------------------

function R = evaluatePolynomial(lambda,varargin)
p = coder.internal.indexInt(length(varargin) - 1);
R = cast(varargin{p + 1},'like',lambda);
if ~isinf(lambda)
    for k = coder.unroll(p:-1:1)
        R = varargin{k} + lambda*R;
    end
end

%--------------------------------------------------------------------------

function offset = findMinResidOffset(X,col,R,p)
% X(idx + (1:n),col) is the segment of X(:,col) with the smallest residual.
% Assumptions: X is an n*p-by-n*p matrix, R is n-by-n, and p is a
% non-negative integer.
coder.internal.prefer_const(p);
ZERO = coder.internal.indexInt(0);
ONE = coder.internal.indexInt(1);
xZERO = zeros('like',X);
rZERO = real(xZERO);
xONE = ones('like',X);
n = coder.internal.indexInt(size(R,1));
np = coder.internal.indexInt(size(X,1));
% Compute RX = R*reshape(X(:,col),n,p).
RX = zeros(n,p,'like',X);
RX = coder.internal.blas.xgemm('N','N',n,p,n,xONE,R,ONE,n, ...
    X,(col - 1)*np + ONE,n,xZERO,RX,ONE,n);
% Compute the residuals for each length n segment of X(:,col) and note the
% location if it is the smallest yet.
firstres = true;
minres = rZERO;
offset = ZERO;
for j = ONE:p
    num = rZERO;
    for i = ONE:n
        num = num + abs(RX(i,j));
    end
    den = rZERO;
    rowOffset = (j - 1)*n;
    for i = ONE:n
        den = den + abs(X(rowOffset + i,col));
    end
    res = num/den;
    if (firstres || res < minres) && ~isnan(res)
        firstres = false;
        minres = res;
        offset = rowOffset;
    end
end

%--------------------------------------------------------------------------
