1########################################################################
2##
3## Copyright (C) 2012-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25
26## -*- texinfo -*-
27## @deftypefn  {} {@var{z} =} polyeig (@var{C0}, @var{C1}, @dots{}, @var{Cl})
28## @deftypefnx {} {[@var{v}, @var{z}] =} polyeig (@var{C0}, @var{C1}, @dots{}, @var{Cl})
29##
30## Solve the polynomial eigenvalue problem of degree @var{l}.
31##
32## Given an @var{n}x@var{n} matrix polynomial
33##
34## @code{@var{C}(@var{s}) = @var{C0} + @var{C1} @var{s} + @dots{} + @var{Cl}
35## @var{s}^@var{l}}
36##
37## @code{polyeig} solves the eigenvalue problem
38##
39## @code{(@var{C0} + @var{C1} @var{z} + @dots{} + @var{Cl} @var{z}^@var{l})
40## @var{v} = 0}.
41##
42## Note that the eigenvalues @var{z} are the zeros of the matrix polynomial.
43## @var{z} is a row vector with @code{@var{n}*@var{l}} elements.  @var{v} is a
44## matrix (@var{n} x @var{n}*@var{l}) with columns that correspond to the
45## eigenvectors.
46##
47## @seealso{eig, eigs, compan}
48## @end deftypefn
49
50function [z, v] = polyeig (varargin)
51
52  if (nargin < 1 || nargout > 2)
53    print_usage ();
54  endif
55
56  n = rows (varargin{1});
57
58  for i = 1 : nargin
59    if (! issquare (varargin{i}))
60      error ("polyeig: coefficients must be square matrices");
61    endif
62    if (rows (varargin{i}) != n)
63      error ("polyeig: coefficients must have the same dimensions");
64    endif
65  endfor
66
67  ## matrix polynomial degree
68  l = nargin - 1;
69
70  ## form needed matrices
71  C = [ zeros(n * (l - 1), n), eye(n * (l - 1));
72       -cell2mat(varargin(1:end-1)) ];
73
74  D = [ eye(n * (l - 1)), zeros(n * (l - 1), n);
75        zeros(n, n * (l - 1)), varargin{end} ];
76
77  ## solve generalized eigenvalue problem
78  if (nargout < 2)
79    z = eig (C, D);
80  else
81    [z, v] = eig (C, D);
82    v = diag (v);
83    ## return n-element eigenvectors normalized so that the infinity-norm = 1
84    z = z(1:n,:);
85    t = max (z);    # max() takes the abs if complex.
86    z ./= t;
87  endif
88
89endfunction
90
91
92%!shared C0, C1
93%! C0 = [8, 0; 0, 4];
94%! C1 = [1, 0; 0, 1];
95
96%!test
97%! z = polyeig (C0, C1);
98%! assert (z, [-8; -4]);
99
100%!test
101%! [v,z] = polyeig (C0, C1);
102%! assert (z, [-8; -4]);
103%! z = diag (z);
104%! d = C0*v + C1*v*z;
105%! assert (norm (d), 0.0);
106
107## Test input validation
108%!error polyeig ()
109%!error [a,b,c] = polyeig (1)
110%!error <coefficients must be square matrices> polyeig (ones (3,2))
111%!error <coefficients must have the same dimensions>
112%! polyeig (ones (3,3), ones (2,2))
113