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