GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
ordqz.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 2020-2025 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_BEGIN_NAMESPACE(octave)
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 case NONE:
333 error ("ordqz: invalid select mode NONE");
334 break;
335
336 // We should have handled all possible enum values above.
337 // Rely on compiler diagnostics to warn if we haven't.
338 // For example, GCC's -Wswitch option, enabled by -Wall,
339 // will provide a warning.
340 }
341 }
342
343 F77_LOGICAL wantq, wantz;
344 wantq = 1, wantz = 1;
345 F77_INT ijob, mm, lrwork3, liwork, info;
346 ijob = 0, lrwork3 = 1, liwork = 1;
347 F77_DBLE pl, pr;
348 ComplexRowVector work3 (lrwork3);
349 Array<F77_INT> iwork (dim_vector (liwork, 1));
350
351 F77_XFCN (ztgsen, ZTGSEN,
352 (ijob, wantq, wantz,
353 select.rwdata (), nn,
354 F77_DBLE_CMPLX_ARG (caa.rwdata ()), nn,
355 F77_DBLE_CMPLX_ARG (cbb.rwdata ()), nn,
356 F77_DBLE_CMPLX_ARG (alpha.rwdata ()),
357 F77_DBLE_CMPLX_ARG (beta.rwdata ()),
358 F77_DBLE_CMPLX_ARG (cqq.rwdata ()), nn,
359 F77_DBLE_CMPLX_ARG (czz.rwdata ()), nn,
360 mm,
361 pl, pr,
362 nullptr,
363 F77_DBLE_CMPLX_ARG (work3.rwdata ()), lrwork3,
364 iwork.rwdata (), liwork,
365 info));
366 if (info != 0)
367 error_with_id ("Octave:ordqz:ztgsen_failed",
368 "ordqz: failed to reorder eigenvalues");
369 }
370 else
371 {
372 // Extract eigenvalues
373 RowVector alphar (dim_vector (nn, 1));
374 RowVector alphai (dim_vector (nn, 1));
375 RowVector beta (dim_vector (nn, 1));
376
377 octave_idx_type k;
378
379 k = 0;
380 while (k < nn)
381 {
382#if defined (DEBUG)
383 octave_stdout << "ordqz: k = " << k << " nn = " << nn << " \n";
384#endif
385 if ((k < nn-1 && aa(k+1, k) == 0.0) || k == nn-1)
386 {
387 alphar(k) = aa(k, k);
388 alphai(k) = 0.0;
389 beta(k) = bb(k, k);
390 k++;
391 }
392 else
393 {
394 double ar[2], ai[2], b[2], work[4];
395 char qz_job = 'E';
396 char comp_q = 'N';
397 char comp_z = 'N';
398 F77_INT nl = 2;
399 F77_INT ilo = 1;
400 F77_INT ihi = 2;
401 F77_INT lwork = 4;
402 F77_INT info = 0;
403 double *aa_vec = aa.rwdata ();
404 double *bb_vec = bb.rwdata ();
405
406 F77_XFCN (dhgeqz, DHGEQZ,
407 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
408 F77_CONST_CHAR_ARG2 (&comp_q, 1),
409 F77_CONST_CHAR_ARG2 (&comp_z, 1),
410 nl, ilo, ihi,
411 &aa_vec[k+k*nn], nn,
412 &bb_vec[k+k*nn], nn,
413 ar, ai, b,
414 nullptr, nn,
415 nullptr, nn, work, lwork, info
416 F77_CHAR_ARG_LEN (1)
417 F77_CHAR_ARG_LEN (1)
418 F77_CHAR_ARG_LEN (1)));
419 if (info != 0)
420 error("ordqz: failed to extract eigenvalues");
421
422 alphar(k) = ar[0];
423 alphar(k+1) = ar[1];
424 alphai(k) = ai[0];
425 alphai(k+1) = ai[1];
426 beta(k) = b[0];
427 beta(k+1) = b[1];
428
429 k += 2;
430 }
431
432 }
433
434 for (k = 0; k < nn; k++)
435 {
436 switch (select_mode)
437 {
438 case LHP:
439 select(k) = alphar(k) * beta(k) < 0;
440 break;
441 case RHP:
442 select(k) = alphar(k) * beta(k) > 0;
443 break;
444 case UDI:
445 select(k) = alphar(k)*alphar(k)
446 + alphai(k)*alphai(k) < beta(k)*beta(k);
447 break;
448 case UDO:
449 select(k) = alphar(k)*alphar(k)
450 + alphai(k)*alphai(k) > beta(k)*beta(k);
451 break;
452 case VEC:
453 if (select_array(k) != 0.0)
454 select(k) = true;
455 else
456 select(k) = false;
457 break;
458 case NONE:
459 error ("ordqz: invalid select mode NONE");
460 break;
461
462 // We should have handled all possible enum values above.
463 // Rely on compiler diagnostics to warn if we haven't.
464 // For example, GCC's -Wswitch option, enabled by -Wall,
465 // will provide a warning.
466 }
467 }
468
469 F77_LOGICAL wantq, wantz;
470 wantq = 1, wantz = 1;
471 F77_INT ijob, mm, lrwork3, liwork, info;
472 ijob = 0, lrwork3 = 4*nn+16, liwork = nn;
473 F77_DBLE pl, pr;
474 RowVector rwork3 (lrwork3);
475 Array<F77_INT> iwork (dim_vector (liwork, 1));
476
477 F77_XFCN (dtgsen, DTGSEN,
478 (ijob, wantq, wantz,
479 select.rwdata (), nn,
480 aa.rwdata (), nn,
481 bb.rwdata (), nn,
482 alphar.rwdata (),
483 alphai.rwdata (),
484 beta.rwdata (),
485 qq.rwdata (), nn,
486 zz.rwdata (), nn,
487 mm,
488 pl, pr,
489 nullptr,
490 rwork3.rwdata (), lrwork3,
491 iwork.rwdata (), liwork,
492 info));
493 if (info != 0)
494 error("ordqz: failed to reorder eigenvalues");
495 }
496
497 octave_value_list retval (nargout);
498 switch (nargout)
499 {
500 case 4:
501 if (complex_case)
502 retval(3) = czz;
503 else
504 retval(3) = zz;
505 OCTAVE_FALLTHROUGH;
506 case 3:
507 if (complex_case)
508 retval(2) = cqq.hermitian();
509 else
510 retval(2) = qq.transpose();
511 OCTAVE_FALLTHROUGH;
512 case 2:
513 if (complex_case)
514 retval(1) = cbb;
515 else
516 retval(1) = bb;
517 OCTAVE_FALLTHROUGH;
518 case 1:
519 if (complex_case)
520 retval(0) = caa;
521 else
522 retval(0) = aa;
523 break;
524 case 0:
525 if (complex_case)
526 retval(0) = caa;
527 else
528 retval(0) = aa;
529 break;
530 }
531
532 return retval;
533}
534
535
536/*
537%!shared A, B, AA, BB, QQ, ZZ, AC, BC, AAC, BBC, QQC, ZZC, select, selectc
538%! A = [ -1.03428 0.24929 0.43205 -0.12860;
539%! 1.16228 0.27870 2.12954 0.69250;
540%! -0.51524 -0.34939 -0.77820 2.13721;
541%! -1.32941 2.11870 0.72005 1.00835 ];
542%! B = [ 1.407302 -0.632956 -0.360628 0.068534;
543%! 0.149898 0.298248 0.991777 0.023652;
544%! 0.169281 -0.405205 -1.775834 1.511730;
545%! 0.717770 1.291390 -1.766607 -0.531352 ];
546%! AC = [ 0.4577 + 0.7199i 0.1476 + 0.6946i 0.6202 + 0.2092i 0.7559 + 0.2759i;
547%! 0.5868 + 0.7275i 0.9174 + 0.8781i 0.6741 + 0.1985i 0.4320 + 0.7023i;
548%! 0.2408 + 0.6359i 0.2959 + 0.8501i 0.3904 + 0.5613i 0.5000 + 0.1428i;
549%! 0.8177 + 0.8581i 0.2583 + 0.8970i 0.7706 + 0.5451i 0.1068 + 0.1650i];
550%! BC = [ 0.089898 + 0.209257i 0.157769 + 0.311387i 0.018926 + 0.622517i 0.058825 + 0.374647i;
551%! 0.009367 + 0.098211i 0.736087 + 0.095797i 0.973192 + 0.583765i 0.434018 + 0.461909i;
552%! 0.880784 + 0.868215i 0.032839 + 0.569461i 0.873437 + 0.266081i 0.739426 + 0.362017i;
553%! 0.121649 + 0.115111i 0.426695 + 0.492222i 0.247670 + 0.034414i 0.771629 + 0.078153i];
554%! [AA, BB, QQ, ZZ] = qz (A, B);
555%! [AAC, BBC, QQC, ZZC] = qz (AC, BC);
556%! select = [0 0 1 1];
557%! selectc = [0 0 0 1];
558
559%!test
560%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "rhp");
561%! assert (all (real (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 0));
562%! assert (all (real (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 0));
563%! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
564%! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
565
566%!test
567%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "+");
568%! assert (all (real (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 0));
569%! assert (all (real (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 0));
570
571%!test
572%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "lhp");
573%! assert (all (real (eig (AAX(2:4,2:4), BBX(2:4,2:4))) >= 0));
574%! assert (all (real (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 0));
575%! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
576%! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
577
578%!test
579%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "-");
580%! assert (all (real (eig (AAX(2:4,2:4), BBX(2:4,2:4))) >= 0));
581%! assert (all (real (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 0));
582
583%!test
584%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "udi");
585%! assert (all (abs (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 1));
586%! assert (all (abs (eig (AAX(2:4,2:4), BBX(2:4,2:4))) > 1));
587%! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
588%! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
589
590%!test
591%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "S");
592%! assert (all (abs (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 1));
593%! assert (all (abs (eig (AAX(2:4,2:4), BBX(2:4,2:4))) > 1));
594
595%!test
596%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "udo");
597%! assert (all (abs (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 1));
598%! assert (all (abs (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 1));
599%! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
600%! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
601
602%!test
603%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "B");
604%! assert (all (abs (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 1));
605%! assert (all (abs (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 1));
606
607%!test
608%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, select);
609%! assert (all (iscomplex (eig (AAX(1:2,1:2), BBX(1:2,1:2)))));
610%! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
611%! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
612
613%!test
614%! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "rhp");
615%! assert (all (real (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) >= 0));
616%! assert (all (real (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) < 0));
617%! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
618%! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
619
620%!test
621%! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "lhp");
622%! assert (all (real (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) < 0));
623%! assert (all (real (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) >= 0));
624%! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
625%! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
626
627%!test
628%! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "udi");
629%! assert (all (abs (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) < 1));
630%! assert (all (abs (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) >= 1));
631%! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
632%! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
633
634%!test
635%! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "udo");
636%! assert (all (abs (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) >= 1));
637%! assert (all (abs (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) < 1));
638%! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
639%! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
640
641%!test
642%! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, selectc);
643%! ev = abs (eig (AACX(1:1,1:1), BBCX(1:1,1:1)));
644%! assert(ev > 0.6 && ev < 0.7);
645%! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
646%! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
647
648%!test
649%! A = toeplitz ([1,2,3,4]);
650%! [B, A] = qr (A);
651%! B = B';
652%! [AA, BB, Q, Z] = qz (A, B);
653%! [AAS, BBS, QS, ZS] = ordqz (AA, BB, Q, Z, "lhp");
654%! E2 = ordeig (AAS, BBS);
655%! ECOMP = [-3.414213562373092; -1.099019513592784;
656%! -0.5857864376269046; 9.099019513592784];
657%! assert (norm (ECOMP - E2, "Inf"), 0, 1e-8);
658
659## Test input validation
660%!error <Invalid call> ordqz ()
661%!error <Invalid call> ordqz (eye (2))
662%!error <Invalid call> ordqz (eye (2), eye (2))
663%!error <Invalid call> ordqz (eye (2), eye (2), eye (2))
664%!error <Invalid call> ordqz (eye (2), eye (2), eye (2), eye (2))
665%!error id=Octave:ordqz:unknown-keyword
666%! ordqz (eye (2), eye (2), eye (2), eye (2), "foobar");
667%!error id=Octave:ordqz:select-not-vector
668%! ordqz (eye (2), eye (2), eye (2), eye (2), eye (2));
669%!error id=Octave:ordqz:unknown-arg
670%! ordqz (eye (2), eye (2), eye (2), eye (2), {"foobar"});
671%!error id=Octave:ordqz:nargout
672%! [a,b,c,d,e] = ordqz (eye (2), eye (2), eye (2), eye (2), "udi");
673%!warning <A: argument is empty matrix> ordqz ([], [], [], [], "udi");
674%!error <A must be a square matrix> ordqz (ones (1,2), [], [], [], "udi");
675%!error <nonconformant matrices>
676%! ordqz (eye (3), eye (2), eye (2), eye (2), "udi");
677%!error <nonconformant matrices>
678%! ordqz (eye (2), eye (3), eye (2), eye (2), "udi");
679%!error <nonconformant matrices>
680%! ordqz (eye (2), eye (2), eye (3), eye (2), "udi");
681%!error <nonconformant matrices>
682%! ordqz (eye (2), eye (2), eye (2), eye (3), "udi");
683%!error <SELECT vector .* wrong number of elements>
684%! ordqz (eye (2), eye (2), eye (2), eye (2), ones (1,5));
685%!error <quasi upper triangular matrices are not allowed with complex data>
686%! AA = zeros (2, 2);
687%! AA(2,1) = i;
688%! ordqz (AA, eye (2), eye (2), eye (2), "udi");
689
690*/
691
692OCTAVE_END_NAMESPACE(octave)
Array
N Dimensional Array with copy-on-write semantics.
Definition Array.h:130
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1787
ComplexMatrix::hermitian
ComplexMatrix hermitian() const
Definition CMatrix.h:168
Matrix::transpose
Matrix transpose() const
Definition dMatrix.h:138
dim_vector
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:90
octave_idx_type
octave_value_list
Definition ovl.h:42
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()
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:1048
error
void error(const char *fmt,...)
Definition error.cc:1003
err_square_matrix_required
void err_square_matrix_required(const char *fcn, const char *name)
Definition errwarn.cc:122
err_nonconformant
void err_nonconformant()
Definition errwarn.cc:95
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
isinteger
bool isinteger(double x)
Definition lo-mappers.h:225
octave_stdout
#define octave_stdout
Definition pager.h:301
abs
template int8_t abs(int8_t)