GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ordqz.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2020-2024 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 // Generalized eigenvalue reordering via LAPACK
27 
28 // Originally written by M. Koehler <koehlerm(AT)mpi-magdeburg.mpg.de>
29 
30 #undef DEBUG
31 
32 #if defined (HAVE_CONFIG_H)
33 # include "config.h"
34 #endif
35 
36 #include <cctype>
37 #include <cmath>
38 
39 #include "f77-fcn.h"
40 #include "lo-lapack-proto.h"
41 #include "qr.h"
42 #include "quit.h"
43 
44 #include "defun.h"
45 #include "error.h"
46 #include "errwarn.h"
47 #include "ovl.h"
48 
49 
50 #if defined (DEBUG)
51 # include "pager.h"
52 # include "pr-output.h"
53 #endif
54 
55 
56 OCTAVE_BEGIN_NAMESPACE(octave)
57 
58 DEFUN (ordqz, args, nargout,
59  doc: /* -*- texinfo -*-
60 @deftypefn {} {[@var{AR}, @var{BR}, @var{QR}, @var{ZR}] =} ordqz (@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{keyword})
61 @deftypefnx {} {[@var{AR}, @var{BR}, @var{QR}, @var{ZR}] =} ordqz (@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{select})
62 Reorder the QZ@tie{}decomposition of a generalized eigenvalue problem.
63 
64 The generalized eigenvalue problem is defined as
65 
66 @tex
67 $$A x = \lambda B x$$
68 @end tex
69 @ifnottex
70 
71 @math{A x = @var{lambda} B x}
72 
73 @end ifnottex
74 
75 Its generalized Schur decomposition is computed using the @code{qz} algorithm:
76 
77 @code{[@var{AA}, @var{BB}, @var{Q}, @var{Z}] = qz (@var{A}, @var{B})}
78 
79 where @var{AA}, @var{BB}, @var{Q}, and @var{Z} fulfill
80 @tex
81 $$ AA = Q \cdot A \cdot Z, BB = Q \cdot B \cdot Z $$
82 @end tex
83 @ifnottex
84 
85 @example
86 @group
87 
88 @var{AA} = @var{Q} * @var{A} * @var{Z}, @var{BB} = @var{Q} * @var{B} * @var{Z}
89 
90 @end group
91 @end example
92 
93 @end ifnottex
94 
95 The @code{ordqz} function computes a unitary transformation @var{QR} and
96 @var{ZR} such that the order of the eigenvalue on the diagonal of @var{AA} and
97 @var{BB} is changed. The resulting reordered matrices @var{AR} and @var{BR}
98 fulfill:
99 
100 @tex
101 $$ A_R = Q_R \cdot A \cdot Z_R, B_R = Q_R \cdot B \cdot Z_R $$
102 @end tex
103 @ifnottex
104 
105 @example
106 @group
107 
108 @var{AR} = @var{QR} * @var{A} * @var{ZR}, @var{BR} = @var{QR} * @var{B} * @var{ZR}
109 
110 @end group
111 @end example
112 
113 @end ifnottex
114 
115 The function can either be called with the @var{keyword} argument which
116 selects the eigenvalues in the top left block of @var{AR} and @var{BR} in the
117 following way:
118 
119 @table @asis
120 @item @qcode{"S"}, @nospell{@qcode{"udi"}}
121 small: leading block has all
122 @tex
123 $|\lambda| < 1$
124 @end tex
125 @ifnottex
126 |@var{lambda}| < 1
127 @end ifnottex
128 
129 @item @qcode{"B"}, @nospell{@qcode{"udo"}}
130 big: leading block has all
131 @tex
132 $|\lambda| \geq 1$
133 @end tex
134 @ifnottex
135 |@var{lambda}| @geq{} 1
136 @end ifnottex
137 
138 @item @qcode{"-"}, @nospell{@qcode{"lhp"}}
139 negative real part: leading block has all eigenvalues in the open left
140 half-plane
141 
142 @item @qcode{"+"}, @nospell{@qcode{"rhp"}}
143 non-negative real part: leading block has all eigenvalues in the closed right
144 half-plane
145 @end table
146 
147 If a logical vector @var{select} is given instead of a keyword the @code{ordqz}
148 function reorders all eigenvalues @code{k} to the left block for which
149 @code{select(k)} is true.
150 
151 Note: The keywords are compatible with the ones from @code{qr}.
152 
153 @seealso{eig, ordeig, qz, schur, ordschur}
154 @end deftypefn */)
155 {
156  enum { LHP, RHP, UDI, UDO, VEC, NONE } select_mode = NONE;
157 
158  if (args.length () != 5)
159  print_usage ();
160 
161  // Check select argument
162  if (args(4).is_string())
163  {
164  std::string opts = args(4).string_value ();
165  std::for_each (opts.begin (), opts.end (),
166  [] (char& c) { c = std::tolower (c); });
167  if (opts == "lhp" || opts == "-")
168  select_mode = LHP;
169  else if (opts == "rhp" || opts == "+")
170  select_mode = RHP;
171  else if (opts == "udi" || opts == "s")
172  select_mode = UDI;
173  else if (opts == "udo" || opts == "b")
174  select_mode = UDO;
175  else
176  error_with_id ("Octave:ordqz:unknown-keyword",
177  "ordqz: unknown KEYWORD, possible values: "
178  "lhp, rhp, udi, udo");
179  }
180  else if (args(4).isreal () || args(4).isinteger () || args(4).islogical ())
181  {
182  if (args(4).rows () > 1 && args(4).columns () > 1)
183  error_with_id ("Octave:ordqz:select-not-vector",
184  "ordqz: SELECT argument must be a vector");
185  select_mode = VEC;
186  }
187  else
188  error_with_id ("Octave:ordqz:unknown-arg",
189  "ordqz: OPT must be string or a logical vector");
190 
191  if (nargout > 4)
192  error_with_id ("Octave:ordqz:nargout",
193  "ordqz: at most four output arguments possible");
194 
195  // Matrix A: check dimensions.
196  F77_INT nn = to_f77_int (args(0).rows ());
197  F77_INT nc = to_f77_int (args(0).columns ());
198 
199  if (args(0).isempty ())
200  {
201  warn_empty_arg ("qz: A");
202  return octave_value_list (2, Matrix ());
203  }
204  else if (nc != nn)
205  err_square_matrix_required ("qz", "A");
206 
207  // Matrix A: get value.
208  Matrix aa;
209  ComplexMatrix caa;
210 
211  if (args(0).iscomplex ())
212  caa = args(0).complex_matrix_value ();
213  else
214  aa = args(0).matrix_value ();
215 
216  // Extract argument 2 (bb, or cbb if complex).
217  F77_INT b_nr = to_f77_int (args(1).rows ());
218  F77_INT b_nc = to_f77_int (args(1).columns ());
219 
220  if (nn != b_nc || nn != b_nr)
221  ::err_nonconformant ();
222 
223  Matrix bb;
224  ComplexMatrix cbb;
225 
226  if (args(1).iscomplex ())
227  cbb = args(1).complex_matrix_value ();
228  else
229  bb = args(1).matrix_value ();
230 
231  // Extract argument 3 (qq, or cqq if complex).
232  F77_INT q_nr = to_f77_int (args(2).rows ());
233  F77_INT q_nc = to_f77_int (args(2).columns ());
234 
235  if (nn != q_nc || nn != q_nr)
236  ::err_nonconformant ();
237 
238  Matrix qq;
239  ComplexMatrix cqq;
240 
241  if (args(2).iscomplex ())
242  cqq = args(2).complex_matrix_value ().hermitian ();
243  else
244  qq = args(2).matrix_value ().transpose ();
245 
246  // Extract argument 4 (zz, or czz if complex).
247  F77_INT z_nr = to_f77_int (args(3).rows ());
248  F77_INT z_nc = to_f77_int (args(3).columns ());
249 
250  if (nn != z_nc || nn != z_nr)
251  ::err_nonconformant ();
252 
253  Matrix zz;
254  ComplexMatrix czz;
255 
256  if (args(3).iscomplex ())
257  czz = args(3).complex_matrix_value ();
258  else
259  zz = args(3).matrix_value ();
260 
261  bool complex_case = (args(0).iscomplex () || args(1).iscomplex ()
262  || args(2).iscomplex () || args(3).iscomplex ());
263 
264  if (select_mode == VEC && args(4).rows () != nn && args(4).columns () != nn)
265  error_with_id ("Octave:ordqz:numel_select",
266  "ordqz: SELECT vector has the wrong number of elements");
267 
268  Array<double> select_array (dim_vector (nn, 1));
269  if (select_mode == VEC)
270  select_array = args(4).vector_value ();
271 
272  Array<F77_LOGICAL> select (dim_vector (nn, 1));
273 
274  if (complex_case)
275  {
276  // Complex
277  if (args(0).isreal ())
278  caa = ComplexMatrix (aa);
279  if (args(1).isreal ())
280  cbb = ComplexMatrix (bb);
281  if (args(2).isreal ())
282  cqq = ComplexMatrix (qq);
283  if (args(3).isreal ())
284  czz = ComplexMatrix (zz);
285 
286  ComplexRowVector alpha (dim_vector (nn, 1));
287  ComplexRowVector beta (dim_vector (nn, 1));
288  octave_idx_type k;
289 
290  for (k = 0; k < nn-1; k++)
291  {
292  if (caa(k+1, k) != 0.0)
293  error_with_id ("Octave:ordqz:unsupported_AA",
294  "ordqz: quasi upper triangular matrices are not "
295  "allowed with complex data");
296  }
297 
298  for (k = 0; k < nn; k++)
299  {
300  alpha(k) = caa(k, k);
301  beta(k) = cbb(k, k);
302  }
303 
304  for (k = 0; k < nn; k++)
305  {
306  switch (select_mode)
307  {
308  case LHP:
309  select(k) = real (alpha(k) * beta(k)) < 0;
310  break;
311  case RHP:
312  select(k) = real (alpha(k) * beta(k)) > 0;
313  break;
314  case UDI:
315  if (beta(k) != 0.0)
316  select(k) = abs (alpha(k)/beta(k)) < 1.0;
317  else
318  select(k) = false;
319  break;
320  case UDO:
321  if (beta(k) != 0.0)
322  select(k) = abs (alpha(k)/beta(k)) > 1.0;
323  else
324  select(k) = true;
325  break;
326  case VEC:
327  if (select_array(k) != 0.0)
328  select(k) = true;
329  else
330  select(k) = false;
331  break;
332  default:
333  // default: case just here to suppress compiler warning.
334  panic_impossible ();
335  }
336  }
337 
338  F77_LOGICAL wantq, wantz;
339  wantq = 1, wantz = 1;
340  F77_INT ijob, mm, lrwork3, liwork, info;
341  ijob = 0, lrwork3 = 1, liwork = 1;
342  F77_DBLE pl, pr;
343  ComplexRowVector work3 (lrwork3);
344  Array<F77_INT> iwork (dim_vector (liwork, 1));
345 
346  F77_XFCN (ztgsen, ZTGSEN,
347  (ijob, wantq, wantz,
348  select.fortran_vec (), nn,
349  F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
350  F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), nn,
351  F77_DBLE_CMPLX_ARG (alpha.fortran_vec ()),
352  F77_DBLE_CMPLX_ARG (beta.fortran_vec ()),
353  F77_DBLE_CMPLX_ARG (cqq.fortran_vec ()), nn,
354  F77_DBLE_CMPLX_ARG (czz.fortran_vec ()), nn,
355  mm,
356  pl, pr,
357  nullptr,
358  F77_DBLE_CMPLX_ARG (work3.fortran_vec ()), lrwork3,
359  iwork.fortran_vec (), liwork,
360  info));
361  if (info != 0)
362  error_with_id ("Octave:ordqz:ztgsen_failed",
363  "ordqz: failed to reorder eigenvalues");
364  }
365  else
366  {
367  // Extract eigenvalues
368  RowVector alphar (dim_vector (nn, 1));
369  RowVector alphai (dim_vector (nn, 1));
370  RowVector beta (dim_vector (nn, 1));
371 
372  octave_idx_type k;
373 
374  k = 0;
375  while (k < nn)
376  {
377 #ifdef DEBUG
378  octave_stdout << "ordqz: k = " << k << " nn = " << nn << " \n";
379 #endif
380  if ((k < nn-1 && aa(k+1, k) == 0.0) || k == nn-1)
381  {
382  alphar(k) = aa(k, k);
383  alphai(k) = 0.0;
384  beta(k) = bb(k, k);
385  k++;
386  }
387  else
388  {
389  double ar[2], ai[2], b[2], work[4];
390  char qz_job = 'E';
391  char comp_q = 'N';
392  char comp_z = 'N';
393  F77_INT nl = 2;
394  F77_INT ilo = 1;
395  F77_INT ihi = 2;
396  F77_INT lwork = 4;
397  F77_INT info = 0;
398  double *aa_vec = aa.fortran_vec ();
399  double *bb_vec = bb.fortran_vec ();
400 
401  F77_XFCN (dhgeqz, DHGEQZ,
402  (F77_CONST_CHAR_ARG2 (&qz_job, 1),
403  F77_CONST_CHAR_ARG2 (&comp_q, 1),
404  F77_CONST_CHAR_ARG2 (&comp_z, 1),
405  nl, ilo, ihi,
406  &aa_vec[k+k*nn], nn,
407  &bb_vec[k+k*nn], nn,
408  ar, ai, b,
409  nullptr, nn,
410  nullptr, nn, work, lwork, info
411  F77_CHAR_ARG_LEN (1)
412  F77_CHAR_ARG_LEN (1)
413  F77_CHAR_ARG_LEN (1)));
414  if (info != 0)
415  error("ordqz: failed to extract eigenvalues");
416 
417  alphar(k) = ar[0];
418  alphar(k+1) = ar[1];
419  alphai(k) = ai[0];
420  alphai(k+1) = ai[1];
421  beta(k) = b[0];
422  beta(k+1) = b[1];
423 
424  k += 2;
425  }
426 
427  }
428 
429  for (k = 0; k < nn; k++)
430  {
431  switch (select_mode)
432  {
433  case LHP:
434  select(k) = alphar(k) * beta(k) < 0;
435  break;
436  case RHP:
437  select(k) = alphar(k) * beta(k) > 0;
438  break;
439  case UDI:
440  select(k) = alphar(k)*alphar(k)
441  + alphai(k)*alphai(k) < beta(k)*beta(k);
442  break;
443  case UDO:
444  select(k) = alphar(k)*alphar(k)
445  + alphai(k)*alphai(k) > beta(k)*beta(k);
446  break;
447  case VEC:
448  if (select_array(k) != 0.0)
449  select(k) = true;
450  else
451  select(k) = false;
452  break;
453  default:
454  // default: case just here to suppress compiler warning.
455  panic_impossible();
456  }
457  }
458 
459  F77_LOGICAL wantq, wantz;
460  wantq = 1, wantz = 1;
461  F77_INT ijob, mm, lrwork3, liwork, info;
462  ijob = 0, lrwork3 = 4*nn+16, liwork = nn;
463  F77_DBLE pl, pr;
464  RowVector rwork3 (lrwork3);
465  Array<F77_INT> iwork (dim_vector (liwork, 1));
466 
467  F77_XFCN (dtgsen, DTGSEN,
468  (ijob, wantq, wantz,
469  select.fortran_vec (), nn,
470  aa.fortran_vec (), nn,
471  bb.fortran_vec (), nn,
472  alphar.fortran_vec (),
473  alphai.fortran_vec (),
474  beta.fortran_vec (),
475  qq.fortran_vec (), nn,
476  zz.fortran_vec (), nn,
477  mm,
478  pl, pr,
479  nullptr,
480  rwork3.fortran_vec (), lrwork3,
481  iwork.fortran_vec (), liwork,
482  info));
483  if (info != 0)
484  error("ordqz: failed to reorder eigenvalues");
485  }
486 
487  octave_value_list retval (nargout);
488  switch (nargout)
489  {
490  case 4:
491  if (complex_case)
492  retval(3) = czz;
493  else
494  retval(3) = zz;
495  OCTAVE_FALLTHROUGH;
496  case 3:
497  if (complex_case)
498  retval(2) = cqq.hermitian();
499  else
500  retval(2) = qq.transpose();
501  OCTAVE_FALLTHROUGH;
502  case 2:
503  if (complex_case)
504  retval(1) = cbb;
505  else
506  retval(1) = bb;
507  OCTAVE_FALLTHROUGH;
508  case 1:
509  if (complex_case)
510  retval(0) = caa;
511  else
512  retval(0) = aa;
513  break;
514  case 0:
515  if (complex_case)
516  retval(0) = caa;
517  else
518  retval(0) = aa;
519  break;
520  }
521 
522  return retval;
523 }
524 
525 
526 /*
527 %!shared A, B, AA, BB, QQ, ZZ, AC, BC, AAC, BBC, QQC, ZZC, select, selectc
528 %! A = [ -1.03428 0.24929 0.43205 -0.12860;
529 %! 1.16228 0.27870 2.12954 0.69250;
530 %! -0.51524 -0.34939 -0.77820 2.13721;
531 %! -1.32941 2.11870 0.72005 1.00835 ];
532 %! B = [ 1.407302 -0.632956 -0.360628 0.068534;
533 %! 0.149898 0.298248 0.991777 0.023652;
534 %! 0.169281 -0.405205 -1.775834 1.511730;
535 %! 0.717770 1.291390 -1.766607 -0.531352 ];
536 %! AC = [ 0.4577 + 0.7199i 0.1476 + 0.6946i 0.6202 + 0.2092i 0.7559 + 0.2759i;
537 %! 0.5868 + 0.7275i 0.9174 + 0.8781i 0.6741 + 0.1985i 0.4320 + 0.7023i;
538 %! 0.2408 + 0.6359i 0.2959 + 0.8501i 0.3904 + 0.5613i 0.5000 + 0.1428i;
539 %! 0.8177 + 0.8581i 0.2583 + 0.8970i 0.7706 + 0.5451i 0.1068 + 0.1650i];
540 %! BC = [ 0.089898 + 0.209257i 0.157769 + 0.311387i 0.018926 + 0.622517i 0.058825 + 0.374647i;
541 %! 0.009367 + 0.098211i 0.736087 + 0.095797i 0.973192 + 0.583765i 0.434018 + 0.461909i;
542 %! 0.880784 + 0.868215i 0.032839 + 0.569461i 0.873437 + 0.266081i 0.739426 + 0.362017i;
543 %! 0.121649 + 0.115111i 0.426695 + 0.492222i 0.247670 + 0.034414i 0.771629 + 0.078153i];
544 %! [AA, BB, QQ, ZZ] = qz (A, B);
545 %! [AAC, BBC, QQC, ZZC] = qz (AC, BC);
546 %! select = [0 0 1 1];
547 %! selectc = [0 0 0 1];
548 
549 %!test
550 %! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "rhp");
551 %! assert (all (real (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 0));
552 %! assert (all (real (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 0));
553 %! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
554 %! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
555 
556 %!test
557 %! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "+");
558 %! assert (all (real (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 0));
559 %! assert (all (real (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 0));
560 
561 %!test
562 %! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "lhp");
563 %! assert (all (real (eig (AAX(2:4,2:4), BBX(2:4,2:4))) >= 0));
564 %! assert (all (real (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 0));
565 %! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
566 %! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
567 
568 %!test
569 %! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "-");
570 %! assert (all (real (eig (AAX(2:4,2:4), BBX(2:4,2:4))) >= 0));
571 %! assert (all (real (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 0));
572 
573 %!test
574 %! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "udi");
575 %! assert (all (abs (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 1));
576 %! assert (all (abs (eig (AAX(2:4,2:4), BBX(2:4,2:4))) > 1));
577 %! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
578 %! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
579 
580 %!test
581 %! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "S");
582 %! assert (all (abs (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 1));
583 %! assert (all (abs (eig (AAX(2:4,2:4), BBX(2:4,2:4))) > 1));
584 
585 %!test
586 %! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "udo");
587 %! assert (all (abs (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 1));
588 %! assert (all (abs (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 1));
589 %! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
590 %! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
591 
592 %!test
593 %! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "B");
594 %! assert (all (abs (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 1));
595 %! assert (all (abs (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 1));
596 
597 %!test
598 %! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, select);
599 %! assert (all (iscomplex (eig (AAX(1:2,1:2), BBX(1:2,1:2)))));
600 %! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
601 %! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
602 
603 %!test
604 %! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "rhp");
605 %! assert (all (real (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) >= 0));
606 %! assert (all (real (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) < 0));
607 %! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
608 %! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
609 
610 %!test
611 %! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "lhp");
612 %! assert (all (real (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) < 0));
613 %! assert (all (real (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) >= 0));
614 %! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
615 %! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
616 
617 %!test
618 %! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "udi");
619 %! assert (all (abs (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) < 1));
620 %! assert (all (abs (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) >= 1));
621 %! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
622 %! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
623 
624 %!test
625 %! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "udo");
626 %! assert (all (abs (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) >= 1));
627 %! assert (all (abs (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) < 1));
628 %! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
629 %! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
630 
631 %!test
632 %! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, selectc);
633 %! ev = abs (eig (AACX(1:1,1:1), BBCX(1:1,1:1)));
634 %! assert(ev > 0.6 && ev < 0.7);
635 %! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
636 %! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
637 
638 %!test
639 %! A = toeplitz ([1,2,3,4]);
640 %! [B, A] = qr (A);
641 %! B = B';
642 %! [AA, BB, Q, Z] = qz (A, B);
643 %! [AAS, BBS, QS, ZS] = ordqz (AA, BB, Q, Z, "lhp");
644 %! E2 = ordeig (AAS, BBS);
645 %! ECOMP = [-3.414213562373092; -1.099019513592784;
646 %! -0.5857864376269046; 9.099019513592784];
647 %! assert (norm (ECOMP - E2, "Inf"), 0, 1e-8);
648 
649 ## Test input validation
650 %!error <Invalid call> ordqz ()
651 %!error <Invalid call> ordqz (eye (2))
652 %!error <Invalid call> ordqz (eye (2), eye (2))
653 %!error <Invalid call> ordqz (eye (2), eye (2), eye (2))
654 %!error <Invalid call> ordqz (eye (2), eye (2), eye (2), eye (2))
655 %!error id=Octave:ordqz:unknown-keyword
656 %! ordqz (eye (2), eye (2), eye (2), eye (2), "foobar");
657 %!error id=Octave:ordqz:select-not-vector
658 %! ordqz (eye (2), eye (2), eye (2), eye (2), eye (2));
659 %!error id=Octave:ordqz:unknown-arg
660 %! ordqz (eye (2), eye (2), eye (2), eye (2), {"foobar"});
661 %!error id=Octave:ordqz:nargout
662 %! [a,b,c,d,e] = ordqz (eye (2), eye (2), eye (2), eye (2), "udi");
663 %!warning <A: argument is empty matrix> ordqz ([], [], [], [], "udi");
664 %!error <A must be a square matrix> ordqz (ones (1,2), [], [], [], "udi");
665 %!error <nonconformant matrices>
666 %! ordqz (eye (3), eye (2), eye (2), eye (2), "udi");
667 %!error <nonconformant matrices>
668 %! ordqz (eye (2), eye (3), eye (2), eye (2), "udi");
669 %!error <nonconformant matrices>
670 %! ordqz (eye (2), eye (2), eye (3), eye (2), "udi");
671 %!error <nonconformant matrices>
672 %! ordqz (eye (2), eye (2), eye (2), eye (3), "udi");
673 %!error <SELECT vector .* wrong number of elements>
674 %! ordqz (eye (2), eye (2), eye (2), eye (2), ones (1,5));
675 %!error <quasi upper triangular matrices are not allowed with complex data>
676 %! AA = zeros (2, 2);
677 %! AA(2,1) = i;
678 %! ordqz (AA, eye (2), eye (2), eye (2), "udi");
679 
680 */
681 
682 OCTAVE_END_NAMESPACE(octave)
Array::fortran_vec
T * fortran_vec()
Size of the specified dimension.
Definition: Array-base.cc:1764
ComplexMatrix::hermitian
ComplexMatrix hermitian() const
Definition: CMatrix.h:170
Matrix
Definition: dMatrix.h:42
Matrix::transpose
Matrix transpose() const
Definition: dMatrix.h:140
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
octave_idx_type
octave_value_list
Definition: ovl.h:44
real
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
print_usage
void print_usage(void)
Definition: defun-int.h:72
DEFUN
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:56
error_with_id
void error_with_id(const char *id, const char *fmt,...)
Definition: error.cc:1033
error
void() error(const char *fmt,...)
Definition: error.cc:988
panic_impossible
#define panic_impossible()
Definition: error.h:503
err_square_matrix_required
void err_square_matrix_required(const char *fcn, const char *name)
Definition: errwarn.cc:122
warn_empty_arg
void warn_empty_arg(const char *name)
Definition: errwarn.cc:330
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:316
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:45
F77_DBLE
double F77_DBLE
Definition: f77-fcn.h:302
F77_LOGICAL
octave_f77_int_type F77_LOGICAL
Definition: f77-fcn.h:308
F77_INT
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
NONE
@ NONE
Definition: gl-render.cc:86
err_nonconformant
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
Definition: lo-array-errwarn.cc:71
isinteger
bool isinteger(double x)
Definition: lo-mappers.h:225
octave_stdout
#define octave_stdout
Definition: pager.h:309
abs
template int8_t abs(int8_t)