GNU Octave  8.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
ordqz.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2020-2023 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)
OCTAVE_END_NAMESPACE
OCTAVE_END_NAMESPACE(octave)
nn
static F77_INT nn
Definition: DASPK.cc:66
Array::fortran_vec
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array-base.cc:1766
ComplexMatrix::hermitian
ComplexMatrix hermitian(void) const
Definition: CMatrix.h:170
Matrix
Definition: dMatrix.h:42
Matrix::transpose
Matrix transpose(void) 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
OCTINTERP_API 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:1024
error
void error(const char *fmt,...)
Definition: error.cc:979
panic_impossible
#define panic_impossible()
Definition: error.h:508
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:85
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
ComplexMatrix
class OCTAVE_API ComplexMatrix
Definition: mx-fwd.h:32
octave_stdout
#define octave_stdout
Definition: pager.h:314
abs
static T abs(T x)
Definition: pr-output.cc:1678