GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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