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