GNU Octave  8.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-2023 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 // FIXME: Matlab does not produce lambda as the first output argument.
63 // Compatibility problem?
64 
65 DEFUN (qz, args, nargout,
66  doc: /* -*- texinfo -*-
67 @deftypefn {} {@var{lambda} =} qz (@var{A}, @var{B})
68 @deftypefnx {} {[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}, @var{lambda}] =} qz (@var{A}, @var{B})
69 @deftypefnx {} {[@var{AA}, @var{BB}, @var{Z}] =} qz (@var{A}, @var{B}, @var{opt})
70 @deftypefnx {} {[@var{AA}, @var{BB}, @var{Z}, @var{lambda}] =} qz (@var{A}, @var{B}, @var{opt})
71 Compute the QZ@tie{}decomposition of a generalized eigenvalue problem.
72 
73 The generalized eigenvalue problem is defined as
74 
75 @tex
76 $$A x = \lambda B x$$
77 @end tex
78 @ifnottex
79 
80 @math{A x = @var{lambda} B x}
81 
82 @end ifnottex
83 
84 There are three calling forms of the function:
85 
86 @enumerate
87 @item @code{@var{lambda} = qz (@var{A}, @var{B})}
88 
89 Compute the generalized eigenvalues
90 @tex
91 $\lambda.$
92 @end tex
93 @ifnottex
94 @var{lambda}.
95 @end ifnottex
96 
97 @item @code{[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}, @var{lambda}] = qz (@var{A}, @var{B})}
98 
99 Compute QZ@tie{}decomposition, generalized eigenvectors, and generalized
100 eigenvalues.
101 @tex
102 $$ AA = Q^T AZ, BB = Q^T BZ $$
103 $$ AV = BV{ \rm diag }(\lambda) $$
104 $$ W^T A = { \rm diag }(\lambda)W^T B $$
105 @end tex
106 @ifnottex
107 
108 @example
109 @group
110 
111 @var{AA} = @var{Q} * @var{A} * @var{Z}, @var{BB} = @var{Q} * @var{B} * @var{Z}
112 @var{A} * @var{V} = @var{B} * @var{V} * diag (@var{lambda})
113 @var{W}' * @var{A} = diag (@var{lambda}) * @var{W}' * @var{B}
114 
115 @end group
116 @end example
117 
118 @end ifnottex
119 with @var{Q} and @var{Z} orthogonal (unitary for complex case).
120 
121 @item @code{[@var{AA}, @var{BB}, @var{Z} @{, @var{lambda}@}] = qz (@var{A}, @var{B}, @var{opt})}
122 
123 As in form 2 above, but allows ordering of generalized eigenpairs for, e.g.,
124 solution of discrete time algebraic @nospell{Riccati} equations. Form 3 is not
125 available for complex matrices, and does not compute the generalized
126 eigenvectors @var{V}, @var{W}, nor the orthogonal matrix @var{Q}.
127 
128 @table @var
129 @item opt
130 for ordering eigenvalues of the @nospell{GEP} pencil. The leading block of
131 the revised pencil contains all eigenvalues that satisfy:
132 
133 @table @asis
134 @item @qcode{"N"}
135 unordered (default)
136 
137 @item @qcode{"S"}
138 small: leading block has all
139 @tex
140 $|\lambda| < 1$
141 @end tex
142 @ifnottex
143 |@var{lambda}| < 1
144 @end ifnottex
145 
146 @item @qcode{"B"}
147 big: leading block has all
148 @tex
149 $|\lambda| \geq 1$
150 @end tex
151 @ifnottex
152 |@var{lambda}| @geq{} 1
153 @end ifnottex
154 
155 @item @qcode{"-"}
156 negative real part: leading block has all eigenvalues in the open left
157 half-plane
158 
159 @item @qcode{"+"}
160 non-negative real part: leading block has all eigenvalues in the closed right
161 half-plane
162 @end table
163 @end table
164 @end enumerate
165 
166 Note: @code{qz} performs permutation balancing, but not scaling
167 (@pxref{XREFbalance,,@code{balance}}), which may be lead to less accurate
168 results than @code{eig}. The order of output arguments was selected for
169 compatibility with @sc{matlab}.
170 @seealso{eig, gsvd, balance, chol, hess, lu, qr, qzhess, schur}
171 @end deftypefn */)
172 {
173  int nargin = args.length ();
174 
175 #if defined (DEBUG)
176  octave_stdout << "qz: nargin = " << nargin
177  << ", nargout = " << nargout << std::endl;
178 #endif
179 
180  if (nargin < 2 || nargin > 3 || nargout > 7)
181  print_usage ();
182 
183  if (nargin == 3 && (nargout < 3 || nargout > 4))
184  error ("qz: invalid number of output arguments for form 3 call");
185 
186 #if defined (DEBUG)
187  octave_stdout << "qz: determine ordering option" << std::endl;
188 #endif
189 
190  // Determine ordering option.
191 
192  char ord_job = 'N';
193  double safmin = 0.0;
194 
195  if (nargin == 3)
196  {
197  std::string opt = args(2).xstring_value ("qz: OPT must be a string");
198 
199  if (opt.empty ())
200  error ("qz: OPT must be a non-empty string");
201 
202  ord_job = std::toupper (opt[0]);
203 
204  std::string valid_opts = "NSB+-";
205 
206  if (valid_opts.find_first_of (ord_job) == std::string::npos)
207  error ("qz: invalid order option '%c'", ord_job);
208 
209  // overflow constant required by dlag2
210  F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1),
211  safmin
212  F77_CHAR_ARG_LEN (1));
213 
214 #if defined (DEBUG_EIG)
215  octave_stdout << "qz: initial value of safmin="
216  << setiosflags (std::ios::scientific)
217  << safmin << std::endl;
218 #endif
219 
220  // Some machines (e.g., DEC alpha) get safmin = 0;
221  // for these, use eps instead to avoid problems in dlag2.
222  if (safmin == 0)
223  {
224 #if defined (DEBUG_EIG)
225  octave_stdout << "qz: DANGER WILL ROBINSON: safmin is 0!"
226  << std::endl;
227 #endif
228 
229  F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("E", 1),
230  safmin
231  F77_CHAR_ARG_LEN (1));
232 
233 #if defined (DEBUG_EIG)
234  octave_stdout << "qz: safmin set to "
235  << setiosflags (std::ios::scientific)
236  << safmin << std::endl;
237 #endif
238  }
239  }
240 
241 #if defined (DEBUG)
242  octave_stdout << "qz: check matrix A" << std::endl;
243 #endif
244 
245  // Matrix A: check dimensions.
246  F77_INT nn = to_f77_int (args(0).rows ());
247  F77_INT nc = to_f77_int (args(0).columns ());
248 
249 #if defined (DEBUG)
250  octave_stdout << "Matrix A dimensions: (" << nn << ',' << nc << ')'
251  << std::endl;
252 #endif
253 
254  if (args(0).isempty ())
255  {
256  warn_empty_arg ("qz: A");
257  return octave_value_list (2, Matrix ());
258  }
259  else if (nc != nn)
260  err_square_matrix_required ("qz", "A");
261 
262  // Matrix A: get value.
263  Matrix aa;
264  ComplexMatrix caa;
265 
266  if (args(0).iscomplex ())
267  caa = args(0).complex_matrix_value ();
268  else
269  aa = args(0).matrix_value ();
270 
271 #if defined (DEBUG)
272  octave_stdout << "qz: check matrix B" << std::endl;
273 #endif
274 
275  // Extract argument 2 (bb, or cbb if complex).
276  F77_INT b_nr = to_f77_int (args(1).rows ());
277  F77_INT b_nc = to_f77_int (args(1).columns ());
278 
279  if (nn != b_nc || nn != b_nr)
280  ::err_nonconformant ();
281 
282  Matrix bb;
283  ComplexMatrix cbb;
284 
285  if (args(1).iscomplex ())
286  cbb = args(1).complex_matrix_value ();
287  else
288  bb = args(1).matrix_value ();
289 
290  // Both matrices loaded, now check whether to calculate complex or real val.
291 
292  bool complex_case = (args(0).iscomplex () || args(1).iscomplex ());
293 
294  if (nargin == 3 && complex_case)
295  error ("qz: cannot re-order complex qz decomposition");
296 
297  // First, declare variables used in both the real and complex cases.
298  // FIXME: There are a lot of excess variables declared.
299  // Probably a better way to handle this.
300  Matrix QQ (nn, nn), ZZ (nn, nn), VR (nn, nn), VL (nn, nn);
301  RowVector alphar (nn), alphai (nn), betar (nn);
302  ComplexRowVector xalpha (nn), xbeta (nn);
303  ComplexMatrix CQ (nn, nn), CZ (nn, nn), CVR (nn, nn), CVL (nn, nn);
304  F77_INT ilo, ihi, info;
305  char comp_q = (nargout >= 3 ? 'V' : 'N');
306  char comp_z = ((nargout >= 4 || nargin == 3)? 'V' : 'N');
307 
308  // Initialize Q, Z to identity matrix if either is needed
309  if (comp_q == 'V' || comp_z == 'V')
310  {
311  double *QQptr = QQ.fortran_vec ();
312  double *ZZptr = ZZ.fortran_vec ();
313  std::fill_n (QQptr, QQ.numel (), 0.0);
314  std::fill_n (ZZptr, ZZ.numel (), 0.0);
315  for (F77_INT i = 0; i < nn; i++)
316  {
317  QQ(i, i) = 1.0;
318  ZZ(i, i) = 1.0;
319  }
320  }
321 
322  // Always perform permutation balancing.
323  const char bal_job = 'P';
324  RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn);
325 
326  if (complex_case)
327  {
328 #if defined (DEBUG)
329  if (comp_q == 'V')
330  octave_stdout << "qz: performing balancing; CQ =\n" << CQ << std::endl;
331 #endif
332  if (args(0).isreal ())
333  caa = ComplexMatrix (aa);
334 
335  if (args(1).isreal ())
336  cbb = ComplexMatrix (bb);
337 
338  if (comp_q == 'V')
339  CQ = ComplexMatrix (QQ);
340 
341  if (comp_z == 'V')
342  CZ = ComplexMatrix (ZZ);
343 
344  F77_XFCN (zggbal, ZGGBAL,
345  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
346  nn, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
347  F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()),
348  nn, ilo, ihi, lscale.fortran_vec (),
349  rscale.fortran_vec (), work.fortran_vec (), info
350  F77_CHAR_ARG_LEN (1)));
351  }
352  else
353  {
354 #if defined (DEBUG)
355  if (comp_q == 'V')
356  octave_stdout << "qz: performing balancing; QQ =\n" << QQ << std::endl;
357 #endif
358 
359  F77_XFCN (dggbal, DGGBAL,
360  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
361  nn, aa.fortran_vec (), nn, bb.fortran_vec (),
362  nn, ilo, ihi, lscale.fortran_vec (),
363  rscale.fortran_vec (), work.fortran_vec (), info
364  F77_CHAR_ARG_LEN (1)));
365  }
366 
367  // Only permutation balance above is done. Skip scaling balance.
368 
369 #if 0
370  // Since we just want the balancing matrices, we can use dggbal
371  // for both the real and complex cases; left first
372 
373  if (comp_q == 'V')
374  {
375  F77_XFCN (dggbak, DGGBAK,
376  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
377  F77_CONST_CHAR_ARG2 ("L", 1),
378  nn, ilo, ihi, lscale.data (), rscale.data (),
379  nn, QQ.fortran_vec (), nn, info
380  F77_CHAR_ARG_LEN (1)
381  F77_CHAR_ARG_LEN (1)));
382 
383 #if defined (DEBUG)
384  if (comp_q == 'V')
385  octave_stdout << "qz: balancing done; QQ =\n" << QQ << std::endl;
386 #endif
387  }
388 
389  // then right
390  if (comp_z == 'V')
391  {
392  F77_XFCN (dggbak, DGGBAK,
393  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
394  F77_CONST_CHAR_ARG2 ("R", 1),
395  nn, ilo, ihi, lscale.data (), rscale.data (),
396  nn, ZZ.fortran_vec (), nn, info
397  F77_CHAR_ARG_LEN (1)
398  F77_CHAR_ARG_LEN (1)));
399 
400 #if defined (DEBUG)
401  if (comp_z == 'V')
402  octave_stdout << "qz: balancing done; ZZ=\n" << ZZ << std::endl;
403 #endif
404  }
405 #endif
406 
407  char qz_job = (nargout < 2 ? 'E' : 'S');
408 
409  if (complex_case)
410  {
411  // Complex case.
412 
413  // The QR decomposition of cbb.
414  math::qr<ComplexMatrix> cbqr (cbb);
415  // The R matrix of QR decomposition for cbb.
416  cbb = cbqr.R ();
417  // (Q*)caa for following work.
418  caa = (cbqr.Q ().hermitian ()) * caa;
419  CQ = CQ * cbqr.Q ();
420 
421  F77_XFCN (zgghrd, ZGGHRD,
422  (F77_CONST_CHAR_ARG2 (&comp_q, 1),
423  F77_CONST_CHAR_ARG2 (&comp_z, 1),
424  nn, ilo, ihi, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()),
425  nn, F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), nn,
426  F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn,
427  F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn, info
428  F77_CHAR_ARG_LEN (1)
429  F77_CHAR_ARG_LEN (1)));
430 
431  ComplexRowVector cwork (nn);
432 
433  F77_XFCN (zhgeqz, ZHGEQZ,
434  (F77_CONST_CHAR_ARG2 (&qz_job, 1),
435  F77_CONST_CHAR_ARG2 (&comp_q, 1),
436  F77_CONST_CHAR_ARG2 (&comp_z, 1),
437  nn, ilo, ihi,
438  F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
439  F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), nn,
440  F77_DBLE_CMPLX_ARG (xalpha.fortran_vec ()),
441  F77_DBLE_CMPLX_ARG (xbeta.fortran_vec ()),
442  F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn,
443  F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn,
444  F77_DBLE_CMPLX_ARG (cwork.fortran_vec ()), nn,
445  rwork.fortran_vec (), info
446  F77_CHAR_ARG_LEN (1)
447  F77_CHAR_ARG_LEN (1)
448  F77_CHAR_ARG_LEN (1)));
449 
450  if (comp_q == 'V')
451  {
452  // Left eigenvector.
453  F77_XFCN (zggbak, ZGGBAK,
454  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
455  F77_CONST_CHAR_ARG2 ("L", 1),
456  nn, ilo, ihi, lscale.data (), rscale.data (),
457  nn, F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn, info
458  F77_CHAR_ARG_LEN (1)
459  F77_CHAR_ARG_LEN (1)));
460  }
461 
462  if (comp_z == 'V')
463  {
464  // Right eigenvector.
465  F77_XFCN (zggbak, ZGGBAK,
466  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
467  F77_CONST_CHAR_ARG2 ("R", 1),
468  nn, ilo, ihi, lscale.data (), rscale.data (),
469  nn, F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn, info
470  F77_CHAR_ARG_LEN (1)
471  F77_CHAR_ARG_LEN (1)));
472  }
473 
474  }
475  else
476  {
477 #if defined (DEBUG)
478  octave_stdout << "qz: performing qr decomposition of bb" << std::endl;
479 #endif
480 
481  // Compute the QR factorization of bb.
482  math::qr<Matrix> bqr (bb);
483 
484 #if defined (DEBUG)
485  octave_stdout << "qz: qr (bb) done; now performing qz decomposition"
486  << std::endl;
487 #endif
488 
489  bb = bqr.R ();
490 
491 #if defined (DEBUG)
492  octave_stdout << "qz: extracted bb" << std::endl;
493 #endif
494 
495  aa = (bqr.Q ()).transpose () * aa;
496 
497 #if defined (DEBUG)
498  octave_stdout << "qz: updated aa " << std::endl;
499  octave_stdout << "bqr.Q () =\n" << bqr.Q () << std::endl;
500 
501  if (comp_q == 'V')
502  octave_stdout << "QQ =" << QQ << std::endl;
503 #endif
504 
505  if (comp_q == 'V')
506  QQ = QQ * bqr.Q ();
507 
508 #if defined (DEBUG)
509  octave_stdout << "qz: precursors done..." << std::endl;
510 #endif
511 
512 #if defined (DEBUG)
513  octave_stdout << "qz: comp_q = " << comp_q << ", comp_z = " << comp_z
514  << std::endl;
515 #endif
516 
517  // Reduce to generalized Hessenberg form.
518  F77_XFCN (dgghrd, DGGHRD,
519  (F77_CONST_CHAR_ARG2 (&comp_q, 1),
520  F77_CONST_CHAR_ARG2 (&comp_z, 1),
521  nn, ilo, ihi, aa.fortran_vec (),
522  nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn,
523  ZZ.fortran_vec (), nn, info
524  F77_CHAR_ARG_LEN (1)
525  F77_CHAR_ARG_LEN (1)));
526 
527  // Check if just computing generalized eigenvalues,
528  // or if we're actually computing the decomposition.
529 
530  // Reduce to generalized Schur form.
531  F77_XFCN (dhgeqz, DHGEQZ,
532  (F77_CONST_CHAR_ARG2 (&qz_job, 1),
533  F77_CONST_CHAR_ARG2 (&comp_q, 1),
534  F77_CONST_CHAR_ARG2 (&comp_z, 1),
535  nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (),
536  nn, alphar.fortran_vec (), alphai.fortran_vec (),
537  betar.fortran_vec (), QQ.fortran_vec (), nn,
538  ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info
539  F77_CHAR_ARG_LEN (1)
540  F77_CHAR_ARG_LEN (1)
541  F77_CHAR_ARG_LEN (1)));
542 
543  if (comp_q == 'V')
544  {
545  F77_XFCN (dggbak, DGGBAK,
546  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
547  F77_CONST_CHAR_ARG2 ("L", 1),
548  nn, ilo, ihi, lscale.data (), rscale.data (),
549  nn, QQ.fortran_vec (), nn, info
550  F77_CHAR_ARG_LEN (1)
551  F77_CHAR_ARG_LEN (1)));
552 
553 #if defined (DEBUG)
554  if (comp_q == 'V')
555  octave_stdout << "qz: balancing done; QQ=\n" << QQ << std::endl;
556 #endif
557  }
558 
559  // then right
560  if (comp_z == 'V')
561  {
562  F77_XFCN (dggbak, DGGBAK,
563  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
564  F77_CONST_CHAR_ARG2 ("R", 1),
565  nn, ilo, ihi, lscale.data (), rscale.data (),
566  nn, ZZ.fortran_vec (), nn, info
567  F77_CHAR_ARG_LEN (1)
568  F77_CHAR_ARG_LEN (1)));
569 
570 #if defined (DEBUG)
571  if (comp_z == 'V')
572  octave_stdout << "qz: balancing done; ZZ=\n" << ZZ << std::endl;
573 #endif
574  }
575 
576  }
577 
578  // Order the QZ decomposition?
579  if (ord_job != 'N')
580  {
581  if (complex_case)
582  // Probably not needed, but better be safe.
583  error ("qz: cannot re-order complex QZ decomposition");
584 
585 #if defined (DEBUG_SORT)
586  octave_stdout << "qz: ordering eigenvalues: ord_job = "
587  << ord_job << std::endl;
588 #endif
589 
590  Array<F77_LOGICAL> select (dim_vector (nn, 1));
591 
592  for (int j = 0; j < nn; j++)
593  {
594  switch (ord_job)
595  {
596  case 'S':
597  select(j) = alphar(j)*alphar(j) + alphai(j)*alphai(j) < betar(j)*betar(j);
598  break;
599 
600  case 'B':
601  select(j) = alphar(j)*alphar(j) + alphai(j)*alphai(j) >= betar(j)*betar(j);
602  break;
603 
604  case '+':
605  select(j) = alphar(j) * betar(j) >= 0;
606  break;
607 
608  case '-':
609  select(j) = alphar(j) * betar(j) < 0;
610  break;
611 
612  default:
613  // Invalid order option
614  // (should never happen since options were checked at the top).
615  panic_impossible ();
616  break;
617  }
618  }
619 
620  F77_LOGICAL wantq = 0, wantz = (comp_z == 'V');
621  F77_INT ijob = 0, mm, lrwork3 = 4*nn+16, liwork = nn;
622  F77_DBLE pl, pr;
623  RowVector rwork3(lrwork3);
624  Array<F77_INT> iwork (dim_vector (liwork, 1));
625 
626  F77_XFCN (dtgsen, DTGSEN,
627  (ijob, wantq, wantz,
628  select.fortran_vec (), nn,
629  aa.fortran_vec (), nn,
630  bb.fortran_vec (), nn,
631  alphar.fortran_vec (),
632  alphai.fortran_vec (),
633  betar.fortran_vec (),
634  nullptr, nn,
635  ZZ.fortran_vec (), nn,
636  mm,
637  pl, pr,
638  nullptr,
639  rwork3.fortran_vec (), lrwork3,
640  iwork.fortran_vec (), liwork,
641  info));
642 
643 #if defined (DEBUG_SORT)
644  octave_stdout << "qz: back from dtgsen: aa =\n";
645  octave_print_internal (octave_stdout, aa);
646  octave_stdout << "\nbb =\n";
647  octave_print_internal (octave_stdout, bb);
648  if (comp_z == 'V')
649  {
650  octave_stdout << "\nZZ =\n";
651  octave_print_internal (octave_stdout, ZZ);
652  }
653  octave_stdout << "\nqz: info=" << info;
654  octave_stdout << "\nalphar =\n";
655  octave_print_internal (octave_stdout, Matrix (alphar));
656  octave_stdout << "\nalphai =\n";
657  octave_print_internal (octave_stdout, Matrix (alphai));
658  octave_stdout << "\nbeta =\n";
659  octave_print_internal (octave_stdout, Matrix (betar));
660  octave_stdout << std::endl;
661 #endif
662  }
663 
664  // Compute the generalized eigenvalues as well?
665  ComplexColumnVector gev;
666 
667  if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4))
668  {
669  if (complex_case)
670  {
671  ComplexColumnVector tmp (nn);
672 
673  for (F77_INT i = 0; i < nn; i++)
674  tmp(i) = xalpha(i) / xbeta(i);
675 
676  gev = tmp;
677  }
678  else
679  {
680 #if defined (DEBUG)
681  octave_stdout << "qz: computing generalized eigenvalues" << std::endl;
682 #endif
683 
684  // Return finite generalized eigenvalues.
685  ComplexColumnVector tmp (nn);
686 
687  for (F77_INT i = 0; i < nn; i++)
688  {
689  if (betar(i) != 0)
690  tmp(i) = Complex (alphar(i), alphai(i)) / betar(i);
691  else
692  tmp(i) = numeric_limits<double>::Inf ();
693  }
694 
695  gev = tmp;
696  }
697  }
698 
699  // Right, left eigenvector matrices.
700  if (nargout >= 5)
701  {
702  // Which side to compute?
703  char side = (nargout == 5 ? 'R' : 'B');
704  // Compute all of them and backtransform
705  char howmany = 'B';
706  // Dummy pointer; select is not used.
707  F77_INT *select = nullptr;
708 
709  if (complex_case)
710  {
711  CVL = CQ;
712  CVR = CZ;
713  ComplexRowVector cwork2 (2 * nn);
714  RowVector rwork2 (8 * nn);
715  F77_INT m;
716 
717  F77_XFCN (ztgevc, ZTGEVC,
718  (F77_CONST_CHAR_ARG2 (&side, 1),
719  F77_CONST_CHAR_ARG2 (&howmany, 1),
720  select, nn, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
721  F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()),
722  nn, F77_DBLE_CMPLX_ARG (CVL.fortran_vec ()), nn,
723  F77_DBLE_CMPLX_ARG (CVR.fortran_vec ()), nn, nn,
724  m, F77_DBLE_CMPLX_ARG (cwork2.fortran_vec ()),
725  rwork2.fortran_vec (), info
726  F77_CHAR_ARG_LEN (1)
727  F77_CHAR_ARG_LEN (1)));
728  }
729  else
730  {
731 #if defined (DEBUG)
732  octave_stdout << "qz: computing generalized eigenvectors" << std::endl;
733 #endif
734 
735  VL = QQ;
736  VR = ZZ;
737  F77_INT m;
738 
739  F77_XFCN (dtgevc, DTGEVC,
740  (F77_CONST_CHAR_ARG2 (&side, 1),
741  F77_CONST_CHAR_ARG2 (&howmany, 1),
742  select, nn, aa.fortran_vec (), nn, bb.fortran_vec (),
743  nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn,
744  m, work.fortran_vec (), info
745  F77_CHAR_ARG_LEN (1)
746  F77_CHAR_ARG_LEN (1)));
747 
748  // Now construct the complex form of VV, WW.
749  F77_INT j = 0;
750 
751  while (j < nn)
752  {
753  octave_quit ();
754 
755  // See if real or complex eigenvalue.
756 
757  // Column increment; assume complex eigenvalue.
758  int cinc = 2;
759 
760  if (j == (nn-1))
761  // Single column.
762  cinc = 1;
763  else if (aa(j+1, j) == 0)
764  cinc = 1;
765 
766  // Now copy the eigenvector (s) to CVR, CVL.
767  if (cinc == 1)
768  {
769  for (F77_INT i = 0; i < nn; i++)
770  CVR(i, j) = VR(i, j);
771 
772  if (side == 'B')
773  for (F77_INT i = 0; i < nn; i++)
774  CVL(i, j) = VL(i, j);
775  }
776  else
777  {
778  // Double column; complex vector.
779 
780  for (F77_INT i = 0; i < nn; i++)
781  {
782  CVR(i, j) = Complex (VR(i, j), VR(i, j+1));
783  CVR(i, j+1) = Complex (VR(i, j), -VR(i, j+1));
784  }
785 
786  if (side == 'B')
787  for (F77_INT i = 0; i < nn; i++)
788  {
789  CVL(i, j) = Complex (VL(i, j), VL(i, j+1));
790  CVL(i, j+1) = Complex (VL(i, j), -VL(i, j+1));
791  }
792  }
793 
794  // Advance to next eigenvectors (if any).
795  j += cinc;
796  }
797  }
798  }
799 
800  octave_value_list retval (nargout);
801 
802  switch (nargout)
803  {
804  case 7:
805  retval(6) = gev;
806  OCTAVE_FALLTHROUGH;
807 
808  case 6:
809  // Return left eigenvectors.
810  retval(5) = CVL;
811  OCTAVE_FALLTHROUGH;
812 
813  case 5:
814  // Return right eigenvectors.
815  retval(4) = CVR;
816  OCTAVE_FALLTHROUGH;
817 
818  case 4:
819  if (nargin == 3)
820  {
821 #if defined (DEBUG)
822  octave_stdout << "qz: sort: retval(3) = gev =\n";
823  octave_print_internal (octave_stdout, ComplexMatrix (gev));
824  octave_stdout << std::endl;
825 #endif
826  retval(3) = gev;
827  }
828  else
829  {
830  if (complex_case)
831  retval(3) = CZ;
832  else
833  retval(3) = ZZ;
834  }
835  OCTAVE_FALLTHROUGH;
836 
837  case 3:
838  if (nargin == 3)
839  {
840  if (complex_case)
841  retval(2) = CZ;
842  else
843  retval(2) = ZZ;
844  }
845  else
846  {
847  if (complex_case)
848  retval(2) = CQ.hermitian ();
849  else
850  retval(2) = QQ.transpose ();
851  }
852  OCTAVE_FALLTHROUGH;
853 
854  case 2:
855  {
856  if (complex_case)
857  {
858 #if defined (DEBUG)
859  octave_stdout << "qz: retval(1) = cbb =\n";
860  octave_print_internal (octave_stdout, cbb);
861  octave_stdout << "\nqz: retval(0) = caa =\n";
862  octave_print_internal (octave_stdout, caa);
863  octave_stdout << std::endl;
864 #endif
865  retval(1) = cbb;
866  retval(0) = caa;
867  }
868  else
869  {
870 #if defined (DEBUG)
871  octave_stdout << "qz: retval(1) = bb =\n";
872  octave_print_internal (octave_stdout, bb);
873  octave_stdout << "\nqz: retval(0) = aa =\n";
874  octave_print_internal (octave_stdout, aa);
875  octave_stdout << std::endl;
876 #endif
877  retval(1) = bb;
878  retval(0) = aa;
879  }
880  }
881  break;
882 
883  case 1:
884  case 0:
885 #if defined (DEBUG)
886  octave_stdout << "qz: retval(0) = gev = " << gev << std::endl;
887 #endif
888  retval(0) = gev;
889  break;
890 
891  default:
892  error ("qz: too many return arguments");
893  break;
894  }
895 
896 #if defined (DEBUG)
897  octave_stdout << "qz: exiting (at long last)" << std::endl;
898 #endif
899 
900  return retval;
901 }
902 
903 /*
904 %!test
905 %! A = [1 2; 0 3];
906 %! B = [1 0; 0 0];
907 %! C = [0 1; 0 0];
908 %! assert (qz (A,B), [1; Inf]);
909 %! assert (qz (A,C), [Inf; Inf]);
910 
911 ## Example 7.7.3 in Golub & Van Loan
912 %!test
913 %! a = [ 10 1 2;
914 %! 1 2 -1;
915 %! 1 1 2 ];
916 %! b = reshape (1:9, 3,3);
917 %! [aa, bb, q, z, v, w, lambda] = qz (a, b);
918 %! assert (q * a * z, aa, norm (aa) * 1e-14);
919 %! assert (q * b * z, bb, norm (bb) * 1e-14);
920 %! is_finite = abs (lambda) < 1 / eps (max (a(:)));
921 %! observed = (b * v * diag (lambda))(:,is_finite);
922 %! assert (observed, (a*v)(:,is_finite), norm (observed) * 1e-14);
923 %! observed = (diag (lambda) * w' * b)(is_finite,:);
924 %! assert (observed, (w'*a)(is_finite,:), norm (observed) * 1e-13);
925 
926 %!test
927 %! A = [0, 0, -1, 0; 1, 0, 0, 0; -1, 0, -2, -1; 0, -1, 1, 0];
928 %! B = [0, 0, 0, 0; 0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1];
929 %! [AA, BB, Q, Z1] = qz (A, B);
930 %! [AA, BB, Z2] = qz (A, B, "-");
931 %! assert (Z1, Z2);
932 
933 %!test
934 %! A = [ -1.03428 0.24929 0.43205 -0.12860;
935 %! 1.16228 0.27870 2.12954 0.69250;
936 %! -0.51524 -0.34939 -0.77820 2.13721;
937 %! -1.32941 2.11870 0.72005 1.00835 ];
938 %! B = [ 1.407302 -0.632956 -0.360628 0.068534;
939 %! 0.149898 0.298248 0.991777 0.023652;
940 %! 0.169281 -0.405205 -1.775834 1.511730;
941 %! 0.717770 1.291390 -1.766607 -0.531352 ];
942 %! [AA, BB, Z, lambda] = qz (A, B, "+");
943 %! assert (all (real (lambda(1:3)) >= 0));
944 %! assert (real (lambda(4) < 0));
945 %! [AA, BB, Z, lambda] = qz (A, B, "-");
946 %! assert (real (lambda(1) < 0));
947 %! assert (all (real (lambda(2:4)) >= 0));
948 %! [AA, BB, Z, lambda] = qz (A, B, "B");
949 %! assert (all (abs (lambda(1:3)) >= 1));
950 %! assert (abs (lambda(4) < 1));
951 %! [AA, BB, Z, lambda] = qz (A, B, "S");
952 %! assert (abs (lambda(1) < 1));
953 %! assert (all (abs (lambda(2:4)) >= 1));
954 */
955 
956 OCTAVE_END_NAMESPACE(octave)
OCTAVE_END_NAMESPACE
OCTAVE_END_NAMESPACE(octave)
nn
static F77_INT nn
Definition: DASPK.cc:66
Inf
#define Inf
Definition: Faddeeva.cc:260
Array
N Dimensional Array with copy-on-write semantics.
Definition: Array.h:130
Array::fortran_vec
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array-base.cc:1766
Matrix
Definition: dMatrix.h:42
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
octave_value_list
Definition: ovl.h:44
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
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
void error(const char *fmt,...)
Definition: error.cc:979
panic_impossible
#define panic_impossible()
Definition: error.h:508
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
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
ComplexMatrix
class OCTAVE_API ComplexMatrix
Definition: mx-fwd.h:32
m
T octave_idx_type m
Definition: mx-inlines.cc:773
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:33
octave_stdout
#define octave_stdout
Definition: pager.h:314
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:1762
transpose
static void transpose(octave_idx_type N, const octave_idx_type *ridx, const octave_idx_type *cidx, octave_idx_type *ridx2, octave_idx_type *cidx2)
Definition: symrcm.cc:391
xdlamch
subroutine xdlamch(cmach, retval)
Definition: xdlamch.f:2
F77_FUNC
F77_RET_T F77_FUNC(xerbla, XERBLA)(F77_CONST_CHAR_ARG_DEF(s_arg