GNU Octave 7.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-2022 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
56OCTAVE_NAMESPACE_BEGIN
57
58DEFUN (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})
62Reorder the QZ@tie{}decomposition of a generalized eigenvalue problem.
63
64The 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
75Its 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
79where @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
95The @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}
98fulfill:
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
115The function can either be called with the @var{keyword} argument which
116selects the eigenvalues in the top left block of @var{AR} and @var{BR} in the
117following way:
118
119@table @asis
120@item @qcode{"S"}, @nospell{@qcode{"udi"}}
121small: 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"}}
130big: 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"}}
139negative real part: leading block has all eigenvalues in the open left
140half-plane
141
142@item @qcode{"+"}, @nospell{@qcode{"rhp"}}
143non-negative real part: leading block has all eigenvalues in the closed right
144half-plane
145@end table
146
147If a logical vector @var{select} is given instead of a keyword the @code{ordqz}
148function reorders all eigenvalues @code{k} to the left block for which
149@code{select(k)} is true.
150
151Note: 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
682OCTAVE_NAMESPACE_END
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.cc:1744
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
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:1025
error
void error(const char *fmt,...)
Definition: error.cc:980
panic_impossible
#define panic_impossible()
Definition: error.h:411
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
ComplexMatrix
class OCTAVE_API ComplexMatrix
Definition: mx-fwd.h:32
octave::math::isinteger
bool isinteger(double x)
Definition: lo-mappers.h:225
octave::NONE
@ NONE
Definition: gl-render.cc:85
octave::err_nonconformant
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
Definition: lo-array-errwarn.cc:71
octave_stdout
#define octave_stdout
Definition: pager.h:314
abs
static T abs(T x)
Definition: pr-output.cc:1678