GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
qz.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1998-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 balancing via LAPACK
27
28// Originally written by A. S. Hodel <scotte@eng.auburn.edu>, but is
29// substantially different with the change to use LAPACK.
30
31#undef DEBUG
32#undef DEBUG_SORT
33#undef DEBUG_EIG
34
35#if defined (HAVE_CONFIG_H)
36# include "config.h"
37#endif
38
39#include <cctype>
40#include <cmath>
41
42#if defined (DEBUG_EIG)
43# include <iomanip>
44#endif
45
46#include "f77-fcn.h"
47#include "lo-lapack-proto.h"
48#include "qr.h"
49#include "quit.h"
50
51#include "defun.h"
52#include "error.h"
53#include "errwarn.h"
54#include "ovl.h"
55#if defined (DEBUG) || defined (DEBUG_SORT)
56# include "pager.h"
57# include "pr-output.h"
58#endif
59
61
62DEFUN (qz, args, nargout,
63 doc: /* -*- texinfo -*-
64@deftypefn {} {[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}] =} qz (@var{A}, @var{B})
65@deftypefnx {} {[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}] =} qz (@var{A}, @var{B}, @var{opt})
66Compute the QZ@tie{}decomposition of a generalized eigenvalue problem.
67
68The generalized eigenvalue problem is defined as
69
70@tex
71$$A x = \lambda B x$$
72@end tex
73@ifnottex
74
75@math{A x = @var{lambda} B x}
76
77@end ifnottex
78
79There are two calling forms of the function:
80
81@enumerate
82@item @code{[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}] = qz (@var{A}, @var{B})}
83
84Compute the complex QZ@tie{}decomposition, generalized eigenvectors, and
85generalized eigenvalues.
86@tex
87$$ AA = Q \cdot A \cdot Z, BB = Q \cdot B \cdot Z $$
88$$ A \cdot V \cdot {\rm diag}(BB) = B \cdot V \cdot {\rm diag}(AA) $$
89$$ {\rm diag}(BB) \cdot W^T \cdot A = {\rm diag}(AA) \cdot W^T \cdot B $$
90@end tex
91@ifnottex
92
93@example
94@group
95
96@var{AA} = @var{Q} * @var{A} * @var{Z}, @var{BB} = @var{Q} * @var{B} * @var{Z}
97@var{A} * @var{V} * diag (diag (@var{BB})) = @var{B} * @var{V} * diag (diag (@var{AA}))
98diag (diag (@var{BB})) * @var{W}' * @var{A} = diag (diag (@var{AA})) * @var{W}' * @var{B}
99
100@end group
101@end example
102
103@end ifnottex
104with @var{AA} and @var{BB} upper triangular, and @var{Q} and @var{Z}
105unitary. The matrices @var{V} and @var{W} respectively contain the right
106and left generalized eigenvectors.
107
108@item @code{[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}] = qz (@var{A}, @var{B}, @var{opt})}
109
110The @var{opt} argument must be equal to either @qcode{"real"} or
111@qcode{"complex"}. If it is equal to @qcode{"complex"}, then this
112calling form is equivalent to the first one with only two input
113arguments.
114
115If @var{opt} is equal to @qcode{"real"}, then the real QZ decomposition
116is computed. In particular, @var{AA} is only guaranteed to be
117quasi-upper triangular with 1-by-1 and 2-by-2 blocks on the diagonal,
118and @var{Q} and @var{Z} are orthogonal. The identities mentioned above
119for right and left generalized eigenvectors are only verified if
120@var{AA} is upper triangular (i.e., when all the generalized eigenvalues
121are real, in which case the real and complex QZ coincide).
122
123@end enumerate
124
125Note: @code{qz} performs permutation balancing, but not scaling
126(@pxref{XREFbalance,,@code{balance}}), which may lead to less accurate
127results than @code{eig}. The order of output arguments was selected for
128compatibility with @sc{matlab}.
129
130For eigenvalues, use the @code{ordeig} function to compute them based on
131the resulting @var{AA} and @var{BB} matrices.
132
133@seealso{eig, gsvd, balance, chol, hess, lu, qr, qzhess, schur, ordeig}
134@end deftypefn */)
135{
136 int nargin = args.length ();
137
138#if defined (DEBUG)
139 octave_stdout << "qz: nargin = " << nargin
140 << ", nargout = " << nargout << std::endl;
141#endif
142
143 if (nargin < 2 || nargin > 3 || nargout > 6)
144 print_usage ();
145
146 bool complex_case = true;
147
148 if (nargin == 3)
149 {
150 std::string opt = args(2).xstring_value ("qz: OPT must be a string");
151
152 if (opt == "real")
153 {
154 if (args(0).iscomplex () || args(1).iscomplex ())
155 {
156 warning ("qz: ignoring 'real' option with complex matrices");
157 complex_case = true;
158 }
159 else
160 complex_case = false;
161 }
162 else if (opt == "complex")
163 complex_case = true;
164 else
165 error ("qz: OPT must be 'real' or 'complex'");
166 }
167
168#if defined (DEBUG)
169 octave_stdout << "qz: check matrix A" << std::endl;
170#endif
171
172 // Matrix A: check dimensions.
173 F77_INT nn = to_f77_int (args(0).rows ());
174 F77_INT nc = to_f77_int (args(0).columns ());
175
176#if defined (DEBUG)
177 octave_stdout << "Matrix A dimensions: (" << nn << ',' << nc << ')'
178 << std::endl;
179#endif
180
181 if (nc != nn)
182 err_square_matrix_required ("qz", "A");
183
184 // Matrix A: get value.
185 Matrix aa;
186 ComplexMatrix caa;
187
188 if (args(0).iscomplex ())
189 caa = args(0).complex_matrix_value ();
190 else
191 aa = args(0).matrix_value ();
192
193#if defined (DEBUG)
194 octave_stdout << "qz: check matrix B" << std::endl;
195#endif
196
197 // Extract argument 2 (bb, or cbb if complex).
198 F77_INT b_nr = to_f77_int (args(1).rows ());
199 F77_INT b_nc = to_f77_int (args(1).columns ());
200
201 if (nn != b_nc || nn != b_nr)
203
204 Matrix bb;
205 ComplexMatrix cbb;
206
207 if (args(1).iscomplex ())
208 cbb = args(1).complex_matrix_value ();
209 else
210 bb = args(1).matrix_value ();
211
212 // First, declare variables used in both the real and complex cases.
213 // FIXME: There are a lot of excess variables declared.
214 // Probably a better way to handle this.
215 Matrix QQ (nn, nn), ZZ (nn, nn), VR (nn, nn), VL (nn, nn);
216 RowVector alphar (nn), alphai (nn), betar (nn);
217 ComplexRowVector xalpha (nn), xbeta (nn);
218 ComplexMatrix CQ (nn, nn), CZ (nn, nn), CVR (nn, nn), CVL (nn, nn);
219 F77_INT ilo, ihi, info;
220 char comp_q = (nargout >= 3 ? 'V' : 'N');
221 char comp_z = (nargout >= 4 ? 'V' : 'N');
222
223 // Initialize Q, Z to identity matrix if either is needed
224 if (comp_q == 'V' || comp_z == 'V')
225 {
226 double *QQptr = QQ.rwdata ();
227 double *ZZptr = ZZ.rwdata ();
228 std::fill_n (QQptr, QQ.numel (), 0.0);
229 std::fill_n (ZZptr, ZZ.numel (), 0.0);
230 for (F77_INT i = 0; i < nn; i++)
231 {
232 QQ(i, i) = 1.0;
233 ZZ(i, i) = 1.0;
234 }
235 }
236
237 // Always perform permutation balancing.
238 const char bal_job = 'P';
239 RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn);
240
241 if (complex_case)
242 {
243#if defined (DEBUG)
244 if (comp_q == 'V')
245 octave_stdout << "qz: performing balancing; CQ =\n" << CQ << std::endl;
246#endif
247 if (args(0).isreal ())
248 caa = ComplexMatrix (aa);
249
250 if (args(1).isreal ())
251 cbb = ComplexMatrix (bb);
252
253 if (comp_q == 'V')
254 CQ = ComplexMatrix (QQ);
255
256 if (comp_z == 'V')
257 CZ = ComplexMatrix (ZZ);
258
259 F77_XFCN (zggbal, ZGGBAL,
260 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
261 nn, F77_DBLE_CMPLX_ARG (caa.rwdata ()), nn,
262 F77_DBLE_CMPLX_ARG (cbb.rwdata ()),
263 nn, ilo, ihi, lscale.rwdata (),
264 rscale.rwdata (), work.rwdata (), info
265 F77_CHAR_ARG_LEN (1)));
266 }
267 else
268 {
269#if defined (DEBUG)
270 if (comp_q == 'V')
271 octave_stdout << "qz: performing balancing; QQ =\n" << QQ << std::endl;
272#endif
273
274 F77_XFCN (dggbal, DGGBAL,
275 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
276 nn, aa.rwdata (), nn, bb.rwdata (),
277 nn, ilo, ihi, lscale.rwdata (),
278 rscale.rwdata (), work.rwdata (), info
279 F77_CHAR_ARG_LEN (1)));
280 }
281
282 // Only permutation balance above is done. Skip scaling balance.
283
284 char qz_job = (nargout < 2 ? 'E' : 'S');
285
286 if (complex_case)
287 {
288 // Complex case.
289
290 // The QR decomposition of cbb.
291 math::qr<ComplexMatrix> cbqr (cbb);
292 // The R matrix of QR decomposition for cbb.
293 cbb = cbqr.R ();
294 // (Q*)caa for following work.
295 caa = (cbqr.Q ().hermitian ()) * caa;
296 CQ = CQ * cbqr.Q ();
297
298 F77_XFCN (zgghrd, ZGGHRD,
299 (F77_CONST_CHAR_ARG2 (&comp_q, 1),
300 F77_CONST_CHAR_ARG2 (&comp_z, 1),
301 nn, ilo, ihi, F77_DBLE_CMPLX_ARG (caa.rwdata ()),
302 nn, F77_DBLE_CMPLX_ARG (cbb.rwdata ()), nn,
303 F77_DBLE_CMPLX_ARG (CQ.rwdata ()), nn,
304 F77_DBLE_CMPLX_ARG (CZ.rwdata ()), nn, info
305 F77_CHAR_ARG_LEN (1)
306 F77_CHAR_ARG_LEN (1)));
307
308 ComplexRowVector cwork (nn);
309
310 F77_XFCN (zhgeqz, ZHGEQZ,
311 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
312 F77_CONST_CHAR_ARG2 (&comp_q, 1),
313 F77_CONST_CHAR_ARG2 (&comp_z, 1),
314 nn, ilo, ihi,
315 F77_DBLE_CMPLX_ARG (caa.rwdata ()), nn,
316 F77_DBLE_CMPLX_ARG (cbb.rwdata ()), nn,
317 F77_DBLE_CMPLX_ARG (xalpha.rwdata ()),
318 F77_DBLE_CMPLX_ARG (xbeta.rwdata ()),
319 F77_DBLE_CMPLX_ARG (CQ.rwdata ()), nn,
320 F77_DBLE_CMPLX_ARG (CZ.rwdata ()), nn,
321 F77_DBLE_CMPLX_ARG (cwork.rwdata ()), nn,
322 rwork.rwdata (), info
323 F77_CHAR_ARG_LEN (1)
324 F77_CHAR_ARG_LEN (1)
325 F77_CHAR_ARG_LEN (1)));
326
327 if (comp_q == 'V')
328 {
329 // Left eigenvector.
330 F77_XFCN (zggbak, ZGGBAK,
331 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
332 F77_CONST_CHAR_ARG2 ("L", 1),
333 nn, ilo, ihi, lscale.data (), rscale.data (),
334 nn, F77_DBLE_CMPLX_ARG (CQ.rwdata ()), nn, info
335 F77_CHAR_ARG_LEN (1)
336 F77_CHAR_ARG_LEN (1)));
337 }
338
339 if (comp_z == 'V')
340 {
341 // Right eigenvector.
342 F77_XFCN (zggbak, ZGGBAK,
343 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
344 F77_CONST_CHAR_ARG2 ("R", 1),
345 nn, ilo, ihi, lscale.data (), rscale.data (),
346 nn, F77_DBLE_CMPLX_ARG (CZ.rwdata ()), nn, info
347 F77_CHAR_ARG_LEN (1)
348 F77_CHAR_ARG_LEN (1)));
349 }
350
351 }
352 else
353 {
354#if defined (DEBUG)
355 octave_stdout << "qz: performing qr decomposition of bb" << std::endl;
356#endif
357
358 // Compute the QR factorization of bb.
359 math::qr<Matrix> bqr (bb);
360
361#if defined (DEBUG)
362 octave_stdout << "qz: qr (bb) done; now performing qz decomposition"
363 << std::endl;
364#endif
365
366 bb = bqr.R ();
367
368#if defined (DEBUG)
369 octave_stdout << "qz: extracted bb" << std::endl;
370#endif
371
372 aa = (bqr.Q ()).transpose () * aa;
373
374#if defined (DEBUG)
375 octave_stdout << "qz: updated aa " << std::endl;
376 octave_stdout << "bqr.Q () =\n" << bqr.Q () << std::endl;
377
378 if (comp_q == 'V')
379 octave_stdout << "QQ =" << QQ << std::endl;
380#endif
381
382 if (comp_q == 'V')
383 QQ = QQ * bqr.Q ();
384
385#if defined (DEBUG)
386 octave_stdout << "qz: precursors done..." << std::endl;
387#endif
388
389#if defined (DEBUG)
390 octave_stdout << "qz: comp_q = " << comp_q << ", comp_z = " << comp_z
391 << std::endl;
392#endif
393
394 // Reduce to generalized Hessenberg form.
395 F77_XFCN (dgghrd, DGGHRD,
396 (F77_CONST_CHAR_ARG2 (&comp_q, 1),
397 F77_CONST_CHAR_ARG2 (&comp_z, 1),
398 nn, ilo, ihi, aa.rwdata (),
399 nn, bb.rwdata (), nn, QQ.rwdata (), nn,
400 ZZ.rwdata (), nn, info
401 F77_CHAR_ARG_LEN (1)
402 F77_CHAR_ARG_LEN (1)));
403
404 // Check if just computing generalized eigenvalues,
405 // or if we're actually computing the decomposition.
406
407 // Reduce to generalized Schur form.
408 F77_XFCN (dhgeqz, DHGEQZ,
409 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
410 F77_CONST_CHAR_ARG2 (&comp_q, 1),
411 F77_CONST_CHAR_ARG2 (&comp_z, 1),
412 nn, ilo, ihi, aa.rwdata (), nn, bb.rwdata (),
413 nn, alphar.rwdata (), alphai.rwdata (),
414 betar.rwdata (), QQ.rwdata (), nn,
415 ZZ.rwdata (), nn, work.rwdata (), nn, info
416 F77_CHAR_ARG_LEN (1)
417 F77_CHAR_ARG_LEN (1)
418 F77_CHAR_ARG_LEN (1)));
419
420 if (comp_q == 'V')
421 {
422 F77_XFCN (dggbak, DGGBAK,
423 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
424 F77_CONST_CHAR_ARG2 ("L", 1),
425 nn, ilo, ihi, lscale.data (), rscale.data (),
426 nn, QQ.rwdata (), nn, info
427 F77_CHAR_ARG_LEN (1)
428 F77_CHAR_ARG_LEN (1)));
429
430#if defined (DEBUG)
431 if (comp_q == 'V')
432 octave_stdout << "qz: balancing done; QQ=\n" << QQ << std::endl;
433#endif
434 }
435
436 // then right
437 if (comp_z == 'V')
438 {
439 F77_XFCN (dggbak, DGGBAK,
440 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
441 F77_CONST_CHAR_ARG2 ("R", 1),
442 nn, ilo, ihi, lscale.data (), rscale.data (),
443 nn, ZZ.rwdata (), nn, info
444 F77_CHAR_ARG_LEN (1)
445 F77_CHAR_ARG_LEN (1)));
446
447#if defined (DEBUG)
448 if (comp_z == 'V')
449 octave_stdout << "qz: balancing done; ZZ=\n" << ZZ << std::endl;
450#endif
451 }
452
453 }
454
455 // Right, left eigenvector matrices.
456 if (nargout >= 5)
457 {
458 // Which side to compute?
459 char side = (nargout == 5 ? 'R' : 'B');
460 // Compute all of them and backtransform
461 char howmany = 'B';
462 // Dummy pointer; select is not used.
463 F77_INT *select = nullptr;
464
465 if (complex_case)
466 {
467 CVL = CQ;
468 CVR = CZ;
469 ComplexRowVector cwork2 (2 * nn);
470 RowVector rwork2 (8 * nn);
471 F77_INT m;
472
473 F77_XFCN (ztgevc, ZTGEVC,
474 (F77_CONST_CHAR_ARG2 (&side, 1),
475 F77_CONST_CHAR_ARG2 (&howmany, 1),
476 select, nn, F77_DBLE_CMPLX_ARG (caa.rwdata ()), nn,
477 F77_DBLE_CMPLX_ARG (cbb.rwdata ()),
478 nn, F77_DBLE_CMPLX_ARG (CVL.rwdata ()), nn,
479 F77_DBLE_CMPLX_ARG (CVR.rwdata ()), nn, nn,
480 m, F77_DBLE_CMPLX_ARG (cwork2.rwdata ()),
481 rwork2.rwdata (), info
482 F77_CHAR_ARG_LEN (1)
483 F77_CHAR_ARG_LEN (1)));
484 }
485 else
486 {
487#if defined (DEBUG)
488 octave_stdout << "qz: computing generalized eigenvectors" << std::endl;
489#endif
490
491 VL = QQ;
492 VR = ZZ;
493 F77_INT m;
494
495 F77_XFCN (dtgevc, DTGEVC,
496 (F77_CONST_CHAR_ARG2 (&side, 1),
497 F77_CONST_CHAR_ARG2 (&howmany, 1),
498 select, nn, aa.rwdata (), nn, bb.rwdata (),
499 nn, VL.rwdata (), nn, VR.rwdata (), nn, nn,
500 m, work.rwdata (), info
501 F77_CHAR_ARG_LEN (1)
502 F77_CHAR_ARG_LEN (1)));
503 }
504 }
505
506 octave_value_list retval (nargout);
507
508 switch (nargout)
509 {
510 case 6:
511 // Return left eigenvectors.
512 if (complex_case)
513 retval(5) = CVL;
514 else
515 retval(5) = VL;
516 OCTAVE_FALLTHROUGH;
517
518 case 5:
519 // Return right eigenvectors.
520 if (complex_case)
521 retval(4) = CVR;
522 else
523 retval(4) = VR;
524 OCTAVE_FALLTHROUGH;
525
526 case 4:
527 if (complex_case)
528 retval(3) = CZ;
529 else
530 retval(3) = ZZ;
531 OCTAVE_FALLTHROUGH;
532
533 case 3:
534 if (complex_case)
535 retval(2) = CQ.hermitian ();
536 else
537 retval(2) = QQ.transpose ();
538 OCTAVE_FALLTHROUGH;
539
540 case 2:
541 case 1:
542 case 0:
543 {
544 if (complex_case)
545 {
546#if defined (DEBUG)
547 octave_stdout << "qz: retval(1) = cbb =\n";
549 octave_stdout << "\nqz: retval(0) = caa =\n";
551 octave_stdout << std::endl;
552#endif
553 retval(1) = cbb;
554 retval(0) = caa;
555 }
556 else
557 {
558#if defined (DEBUG)
559 octave_stdout << "qz: retval(1) = bb =\n";
561 octave_stdout << "\nqz: retval(0) = aa =\n";
563 octave_stdout << std::endl;
564#endif
565 retval(1) = bb;
566 retval(0) = aa;
567 }
568 }
569 break;
570
571 }
572
573 // FIXME: The API for qz changed in version 9.
574 // These warnings can be removed in Octave version 11.
575 if (nargout == 1)
576 {
577 warning_with_id ("Octave:qz:single-arg-out",
578 "qz: requesting a single output argument no longer returns eigenvalues since version 9");
579 disable_warning ("Octave:qz:single-arg-out");
580 }
581
582 if (nargin == 2 && args(0).isreal () && args(1).isreal ()
583 && retval(0).iscomplex ())
584 {
585 warning_with_id ("Octave:qz:complex-default",
586 "qz: returns the complex QZ by default on real matrices since version 9");
587 disable_warning ("Octave:qz:complex-default");
588 }
589
590#if defined (DEBUG)
591 octave_stdout << "qz: exiting (at long last)" << std::endl;
592#endif
593
594 return retval;
595}
596
597/*
598## Example 7.7.3 in Golub & Van Loan (3rd ed.; not present in 4th ed.)
599## The generalized eigenvalues are real, so the real and complex QZ coincide.
600%!test
601%! a = [ 10 1 2;
602%! 1 2 -1;
603%! 1 1 2 ];
604%! b = reshape (1:9, 3,3);
605%! [aa, bb, q, z, v, w] = qz (a, b);
606%! assert (isreal (aa) && isreal (bb) && isreal (q) && isreal (z) ...
607%! && isreal (v) && isreal (w));
608%! assert (istriu (aa) && istriu (bb));
609%! assert (q * q', eye(3), 1e-15);
610%! assert (z * z', eye(3), 1e-15);
611%! assert (q * a * z, aa, norm (aa) * 1e-14);
612%! assert (q * b * z, bb, norm (bb) * 1e-14);
613%! assert (a * v * diag (diag (bb)), b * v * diag (diag (aa)), 1e-13);
614%! assert (diag (diag (bb)) * w' * a, diag (diag (aa)) * w' * b, 1e-13);
615
616## An example with real matrices where some generalized eigenvalues are
617## complex, so the real and complex QZ differ.
618%!test
619%! A = [ -1.03428 0.24929 0.43205 -0.12860;
620%! 1.16228 0.27870 2.12954 0.69250;
621%! -0.51524 -0.34939 -0.77820 2.13721;
622%! -1.32941 2.11870 0.72005 1.00835 ];
623%! B = [ 1.407302 -0.632956 -0.360628 0.068534;
624%! 0.149898 0.298248 0.991777 0.023652;
625%! 0.169281 -0.405205 -1.775834 1.511730;
626%! 0.717770 1.291390 -1.766607 -0.531352 ];
627%! [CAA, CBB, CQ, CZ, CV, CW] = qz (A, B, 'complex');
628%! assert (iscomplex (CAA) && iscomplex (CBB) && iscomplex (CQ) ...
629%! && iscomplex (CZ) && iscomplex (CV) && iscomplex (CW));
630%! assert (istriu (CAA) && istriu (CBB));
631%! assert (CQ * CQ', eye(4), 1e-14);
632%! assert (CZ * CZ', eye(4), 1e-14);
633%! assert (CQ * A * CZ, CAA, norm (CAA) * 1e-14);
634%! assert (CQ * B * CZ, CBB, norm (CBB) * 1e-14);
635%! assert (A * CV * diag (diag (CBB)), B * CV * diag (diag (CAA)), 1e-13);
636%! assert (diag (diag (CBB)) * CW' * A, diag (diag (CAA)) * CW' * B, 1e-13);
637%! [AA, BB, Q, Z, V, W] = qz (A, B, 'real');
638%! assert (isreal (AA) && isreal (BB) && isreal (Q) && isreal (Z) ...
639%! && isreal (V) && isreal (W));
640%! assert (istriu (BB));
641%! assert (Q * Q', eye(4), 1e-14);
642%! assert (Z * Z', eye(4), 1e-14);
643%! assert (Q * A * Z, AA, norm (AA) * 1e-14);
644%! assert (Q * B * Z, BB, norm (BB) * 1e-14);
645
646## Test input validation
647%!error <Invalid call> qz ()
648%!error <Invalid call> qz (1)
649%!error <Invalid call> qz (1,2,3,4)
650%!error <Invalid call> [r1,r2,r3,r4,r5,r6,r7] = qz (1,2)
651%!error <OPT must be a string> qz (1,2, 3)
652%!error <OPT must be 'real' or 'complex'> qz (1,2, 'foobar')
653%!warning <ignoring 'real' option with complex matrices> qz (2i, 3, 'real');
654%!warning <ignoring 'real' option with complex matrices> qz (2, 3i, 'real');
655
656*/
657
658OCTAVE_END_NAMESPACE(octave)
T * rwdata()
Size of the specified dimension.
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
void print_usage()
Definition defun-int.h:72
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition defun.h:56
void warning(const char *fmt,...)
Definition error.cc:1078
void warning_with_id(const char *id, const char *fmt,...)
Definition error.cc:1093
void error(const char *fmt,...)
Definition error.cc:1003
void disable_warning(const std::string &id)
Definition error.cc:1864
void err_square_matrix_required(const char *fcn, const char *name)
Definition errwarn.cc:122
void err_nonconformant()
Definition errwarn.cc:95
#define F77_DBLE_CMPLX_ARG(x)
Definition f77-fcn.h:316
#define F77_XFCN(f, F, args)
Definition f77-fcn.h:45
octave_f77_int_type F77_INT
Definition f77-fcn.h:306
F77_RET_T F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * VR
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX F77_DBLE_CMPLX F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX * CZ
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX F77_DBLE_CMPLX F77_DBLE_CMPLX * CQ
F77_RET_T F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * VL
#define octave_stdout
Definition pager.h:301
void octave_print_internal(std::ostream &os, const float_display_format &fmt, bool d, bool pr_as_read_syntax)