GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
qz.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1998-2017 A. S. Hodel
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 // Generalized eigenvalue balancing via LAPACK
24 
25 // Author: A. S. Hodel <scotte@eng.auburn.edu>
26 
27 #undef DEBUG
28 #undef DEBUG_SORT
29 #undef DEBUG_EIG
30 
31 #if defined (HAVE_CONFIG_H)
32 # include "config.h"
33 #endif
34 
35 #include <cfloat>
36 
37 #include <iostream>
38 #include <iomanip>
39 
40 #include "f77-fcn.h"
41 #include "lo-lapack-proto.h"
42 #include "lo-math.h"
43 #include "qr.h"
44 #include "quit.h"
45 
46 #include "defun.h"
47 #include "error.h"
48 #include "errwarn.h"
49 #include "ovl.h"
50 #include "oct-map.h"
51 #include "ov.h"
52 #include "pager.h"
53 #if defined (DEBUG) || defined (DEBUG_SORT)
54 # include "pr-output.h"
55 #endif
56 #include "symtab.h"
57 #include "utils.h"
58 #include "variables.h"
59 
61  const double& ALPHA,
62  const double& BETA, const double& S,
63  const double& P);
64 
65 extern "C"
66 {
67  // Van Dooren's code (netlib.org: toms/590) for reordering
68  // GEP. Only processes Z, not Q.
69  F77_RET_T
70  F77_FUNC (dsubsp, DSUBSP) (const F77_INT& NMAX,
71  const F77_INT& N, F77_DBLE* A,
73  const F77_DBLE& EPS, F77_INT& NDIM,
74  F77_INT& FAIL, F77_INT* IND);
75 }
76 
77 // fcrhp, fin, fout, folhp:
78 // Routines for ordering of generalized eigenvalues.
79 // Return 1 if test is passed, 0 otherwise.
80 // fin: |lambda| < 1
81 // fout: |lambda| >= 1
82 // fcrhp: real(lambda) >= 0
83 // folhp: real(lambda) < 0
84 
85 static octave_idx_type
86 fcrhp (const octave_idx_type& lsize, const double& alpha,
87  const double& beta, const double& s, const double&)
88 {
89  if (lsize == 1)
90  return (alpha * beta >= 0 ? 1 : -1);
91  else
92  return (s >= 0 ? 1 : -1);
93 }
94 
95 static octave_idx_type
96 fin (const octave_idx_type& lsize, const double& alpha,
97  const double& beta, const double&, const double& p)
98 {
100 
101  if (lsize == 1)
102  retval = (fabs (alpha) < fabs (beta) ? 1 : -1);
103  else
104  retval = (fabs (p) < 1 ? 1 : -1);
105 
106 #if defined (DEBUG)
107  std::cout << "qz: fin: retval=" << retval << std::endl;
108 #endif
109 
110  return retval;
111 }
112 
113 static octave_idx_type
114 folhp (const octave_idx_type& lsize, const double& alpha,
115  const double& beta, const double& s, const double&)
116 {
117  if (lsize == 1)
118  return (alpha * beta < 0 ? 1 : -1);
119  else
120  return (s < 0 ? 1 : -1);
121 }
122 
123 static octave_idx_type
124 fout (const octave_idx_type& lsize, const double& alpha,
125  const double& beta, const double&, const double& p)
126 {
127  if (lsize == 1)
128  return (fabs (alpha) >= fabs (beta) ? 1 : -1);
129  else
130  return (fabs (p) >= 1 ? 1 : -1);
131 }
132 
133 //FIXME: Matlab does not produce lambda as the first output argument.
134 // Compatibility problem?
135 DEFUN (qz, args, nargout,
136  doc: /* -*- texinfo -*-
137 @deftypefn {} {@var{lambda} =} qz (@var{A}, @var{B})
138 @deftypefnx {} {@var{lambda} =} qz (@var{A}, @var{B}, @var{opt})
139 QZ@tie{}decomposition of the generalized eigenvalue problem
140 (@math{A x = s B x}).
141 
142 There are three ways to call this function:
143 @enumerate
144 @item @code{@var{lambda} = qz (@var{A}, @var{B})}
145 
146 Computes the generalized eigenvalues
147 @tex
148 $\lambda$
149 @end tex
150 @ifnottex
151 @var{lambda}
152 @end ifnottex
153 of @math{(A - s B)}.
154 
155 @item @code{[AA, BB, Q, Z, V, W, @var{lambda}] = qz (@var{A}, @var{B})}
156 
157 Computes QZ@tie{}decomposition, generalized eigenvectors, and generalized
158 eigenvalues of @math{(A - s B)}
159 @tex
160 $$ AV = BV{ \rm diag }(\lambda) $$
161 $$ W^T A = { \rm diag }(\lambda)W^T B $$
162 $$ AA = Q^T AZ, BB = Q^T BZ $$
163 @end tex
164 @ifnottex
165 
166 @example
167 @group
168 
169 A * V = B * V * diag (@var{lambda})
170 W' * A = diag (@var{lambda}) * W' * B
171 AA = Q * A * Z, BB = Q * B * Z
172 
173 @end group
174 @end example
175 
176 @end ifnottex
177 with @var{Q} and @var{Z} orthogonal (unitary)= @var{I}
178 
179 @item @code{[AA,BB,Z@{, @var{lambda}@}] = qz (@var{A}, @var{B}, @var{opt})}
180 
181 As in form [2], but allows ordering of generalized eigenpairs for, e.g.,
182 solution of discrete time algebraic Riccati equations. Form 3 is not
183 available for complex matrices, and does not compute the generalized
184 eigenvectors @var{V}, @var{W}, nor the orthogonal matrix @var{Q}.
185 
186 @table @var
187 @item opt
188 for ordering eigenvalues of the @nospell{GEP} pencil. The leading block of
189 the revised pencil contains all eigenvalues that satisfy:
190 
191 @table @asis
192 @item @qcode{"N"}
193 = unordered (default)
194 
195 @item @qcode{"S"}
196 = small: leading block has all |lambda| @leq{} 1
197 
198 @item @qcode{"B"}
199 = big: leading block has all |lambda| @geq{} 1
200 
201 @item @qcode{"-"}
202 = negative real part: leading block has all eigenvalues
203 in the open left half-plane
204 
205 @item @qcode{"+"}
206 = non-negative real part: leading block has all eigenvalues
207 in the closed right half-plane
208 @end table
209 @end table
210 @end enumerate
211 
212 Note: @code{qz} performs permutation balancing, but not scaling
213 (@pxref{XREFbalance}). The order of output arguments was selected for
214 compatibility with @sc{matlab}.
215 @seealso{eig, balance, lu, chol, hess, qr, qzhess, schur, svd}
216 @end deftypefn */)
217 {
218  volatile int nargin = args.length ();
219 
220 #if defined (DEBUG)
221  std::cout << "qz: nargin = " << nargin
222  << ", nargout = " << nargout << std::endl;
223 #endif
224 
225  if (nargin < 2 || nargin > 3 || nargout > 7)
226  print_usage ();
227 
228  if (nargin == 3 && (nargout < 3 || nargout > 4))
229  error ("qz: invalid number of output arguments for form [3] call");
230 
231 #if defined (DEBUG)
232  std::cout << "qz: determine ordering option" << std::endl;
233 #endif
234 
235  // Determine ordering option.
236  volatile char ord_job = 0;
237  static double safmin;
238 
239  if (nargin == 2)
240  ord_job = 'N';
241  else
242  {
243  std::string tmp = args(2).xstring_value ("qz: OPT must be a string");
244 
245  if (! tmp.empty ())
246  ord_job = tmp[0];
247 
248  if (! (ord_job == 'N' || ord_job == 'n'
249  || ord_job == 'S' || ord_job == 's'
250  || ord_job == 'B' || ord_job == 'b'
251  || ord_job == '+' || ord_job == '-'))
252  error ("qz: invalid order option");
253 
254  // overflow constant required by dlag2
255  F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1),
256  safmin
257  F77_CHAR_ARG_LEN (1));
258 
259 #if defined (DEBUG_EIG)
260  std::cout << "qz: initial value of safmin="
261  << setiosflags (std::ios::scientific)
262  << safmin << std::endl;
263 #endif
264 
265  // Some machines (e.g., DEC alpha) get safmin = 0;
266  // for these, use eps instead to avoid problems in dlag2.
267  if (safmin == 0)
268  {
269 #if defined (DEBUG_EIG)
270  std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl;
271 #endif
272 
273  F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("E", 1),
274  safmin
275  F77_CHAR_ARG_LEN (1));
276 
277 #if defined (DEBUG_EIG)
278  std::cout << "qz: safmin set to "
279  << setiosflags (std::ios::scientific)
280  << safmin << std::endl;
281 #endif
282  }
283  }
284 
285 #if defined (DEBUG)
286  std::cout << "qz: check argument 1" << std::endl;
287 #endif
288 
289  // Argument 1: check if it's okay dimensioned.
290  octave_idx_type nn = args(0).rows ();
291 
292 #if defined (DEBUG)
293  std::cout << "argument 1 dimensions: ("
294  << nn << "," << args(0).columns () << ")"
295  << std::endl;
296 #endif
297 
299 
300  if (args(0).is_empty ())
301  {
302  warn_empty_arg ("qz: parameter 1; continuing");
303  return octave_value_list (2, Matrix ());
304  }
305  else if (args(0).columns () != nn)
306  err_square_matrix_required ("qz", "A");
307 
308  // Argument 1: dimensions look good; get the value.
309  Matrix aa;
310  ComplexMatrix caa;
311 
312  if (args(0).is_complex_type ())
313  caa = args(0).complex_matrix_value ();
314  else
315  aa = args(0).matrix_value ();
316 
317 #if defined (DEBUG)
318  std::cout << "qz: check argument 2" << std::endl;
319 #endif
320 
321  // Extract argument 2 (bb, or cbb if complex).
322  if ((nn != args(1).columns ()) || (nn != args(1).rows ()))
324 
325  Matrix bb;
326  ComplexMatrix cbb;
327 
328  if (args(1).is_complex_type ())
329  cbb = args(1).complex_matrix_value ();
330  else
331  bb = args(1).matrix_value ();
332 
333  // Both matrices loaded, now let's check what kind of arithmetic:
334  // declared volatile to avoid compiler warnings about long jumps,
335  // vforks.
336 
337  volatile int complex_case
338  = (args(0).is_complex_type () || args(1).is_complex_type ());
339 
340  if (nargin == 3 && complex_case)
341  error ("qz: cannot re-order complex qz decomposition");
342 
343  // First, declare variables used in both the real and complex case.
344  Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn);
345  RowVector alphar(nn), alphai(nn), betar(nn);
346  ComplexRowVector xalpha(nn), xbeta(nn);
347  ComplexMatrix CQ(nn,nn), CZ(nn,nn), CVR(nn,nn), CVL(nn,nn);
348  octave_idx_type ilo, ihi, info;
349  char compq = (nargout >= 3 ? 'V' : 'N');
350  char compz = ((nargout >= 4 || nargin == 3)? 'V' : 'N');
351 
352  // Initialize Q, Z to identity if we need either of them.
353  if (compq == 'V' || compz == 'V')
354  for (octave_idx_type ii = 0; ii < nn; ii++)
355  for (octave_idx_type jj = 0; jj < nn; jj++)
356  {
357  OCTAVE_QUIT;
358  QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0);
359  }
360 
361  // Always perform permutation balancing.
362  const char bal_job = 'P';
363  RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn);
364 
365  if (complex_case)
366  {
367 #if defined (DEBUG)
368  if (compq == 'V')
369  std::cout << "qz: performing balancing; CQ=" << std::endl
370  << CQ << std::endl;
371 #endif
372  if (args(0).is_real_type ())
373  caa = ComplexMatrix (aa);
374 
375  if (args(1).is_real_type ())
376  cbb = ComplexMatrix (bb);
377 
378  if (compq == 'V')
379  CQ = ComplexMatrix (QQ);
380 
381  if (compz == 'V')
382  CZ = ComplexMatrix (ZZ);
383 
384  F77_XFCN (zggbal, ZGGBAL,
385  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
386  nn, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
388  nn, ilo, ihi, lscale.fortran_vec (),
389  rscale.fortran_vec (), work.fortran_vec (), info
390  F77_CHAR_ARG_LEN (1)));
391  }
392  else
393  {
394 #if defined (DEBUG)
395  if (compq == 'V')
396  std::cout << "qz: performing balancing; QQ=" << std::endl
397  << QQ << std::endl;
398 #endif
399 
400  F77_XFCN (dggbal, DGGBAL,
401  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
402  nn, aa.fortran_vec (), nn, bb.fortran_vec (),
403  nn, ilo, ihi, lscale.fortran_vec (),
404  rscale.fortran_vec (), work.fortran_vec (), info
405  F77_CHAR_ARG_LEN (1)));
406  }
407 
408  // Since we just want the balancing matrices, we can use dggbal
409  // for both the real and complex cases; left first
410 
411 #if 0
412  if (compq == 'V')
413  {
414  F77_XFCN (dggbak, DGGBAK,
415  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
416  F77_CONST_CHAR_ARG2 ("L", 1),
417  nn, ilo, ihi, lscale.data (), rscale.data (),
418  nn, QQ.fortran_vec (), nn, info
419  F77_CHAR_ARG_LEN (1)
420  F77_CHAR_ARG_LEN (1)));
421 
422 #if defined (DEBUG)
423  if (compq == 'V')
424  std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl;
425 #endif
426  }
427 
428  // then right
429  if (compz == 'V')
430  {
431  F77_XFCN (dggbak, DGGBAK,
432  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
433  F77_CONST_CHAR_ARG2 ("R", 1),
434  nn, ilo, ihi, lscale.data (), rscale.data (),
435  nn, ZZ.fortran_vec (), nn, info
436  F77_CHAR_ARG_LEN (1)
437  F77_CHAR_ARG_LEN (1)));
438 
439 #if defined (DEBUG)
440  if (compz == 'V')
441  std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
442 #endif
443  }
444 #endif
445 
446  static char qz_job;
447  qz_job = (nargout < 2 ? 'E' : 'S');
448 
449  if (complex_case)
450  {
451  // Complex case.
452 
453  // The QR decomposition of cbb.
455  // The R matrix of QR decomposition for cbb.
456  cbb = cbqr.R ();
457  // (Q*)caa for following work.
458  caa = (cbqr.Q ().hermitian ()) * caa;
459  CQ = CQ * cbqr.Q ();
460 
461  F77_XFCN (zgghrd, ZGGHRD,
462  (F77_CONST_CHAR_ARG2 (&compq, 1),
463  F77_CONST_CHAR_ARG2 (&compz, 1),
464  nn, ilo, ihi, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()),
465  nn, F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), nn,
466  F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn,
467  F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn, info
468  F77_CHAR_ARG_LEN (1)
469  F77_CHAR_ARG_LEN (1)));
470 
471  ComplexRowVector cwork (1 * nn);
472 
473  F77_XFCN (zhgeqz, ZHGEQZ,
474  (F77_CONST_CHAR_ARG2 (&qz_job, 1),
475  F77_CONST_CHAR_ARG2 (&compq, 1),
476  F77_CONST_CHAR_ARG2 (&compz, 1),
477  nn, ilo, ihi,
478  F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
479  F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), nn,
480  F77_DBLE_CMPLX_ARG (xalpha.fortran_vec ()),
481  F77_DBLE_CMPLX_ARG (xbeta.fortran_vec ()),
482  F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn,
483  F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn,
484  F77_DBLE_CMPLX_ARG (cwork.fortran_vec ()), nn, rwork.fortran_vec (), info
485  F77_CHAR_ARG_LEN (1)
486  F77_CHAR_ARG_LEN (1)
487  F77_CHAR_ARG_LEN (1)));
488 
489  if (compq == 'V')
490  {
491  // Left eigenvector.
492  F77_XFCN (zggbak, ZGGBAK,
493  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
494  F77_CONST_CHAR_ARG2 ("L", 1),
495  nn, ilo, ihi, lscale.data (), rscale.data (),
496  nn, F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn, info
497  F77_CHAR_ARG_LEN (1)
498  F77_CHAR_ARG_LEN (1)));
499  }
500 
501  // Right eigenvector.
502  if (compz == 'V')
503  {
504  F77_XFCN (zggbak, ZGGBAK,
505  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
506  F77_CONST_CHAR_ARG2 ("R", 1),
507  nn, ilo, ihi, lscale.data (), rscale.data (),
508  nn, F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn, info
509  F77_CHAR_ARG_LEN (1)
510  F77_CHAR_ARG_LEN (1)));
511  }
512 
513  }
514  else
515  {
516 #if defined (DEBUG)
517  std::cout << "qz: peforming qr decomposition of bb" << std::endl;
518 #endif
519 
520  // Compute the QR factorization of bb.
521  octave::math::qr<Matrix> bqr (bb);
522 
523 #if defined (DEBUG)
524  std::cout << "qz: qr (bb) done; now peforming qz decomposition"
525  << std::endl;
526 #endif
527 
528  bb = bqr.R ();
529 
530 #if defined (DEBUG)
531  std::cout << "qz: extracted bb" << std::endl;
532 #endif
533 
534  aa = (bqr.Q ()).transpose () * aa;
535 
536 #if defined (DEBUG)
537  std::cout << "qz: updated aa " << std::endl;
538  std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl;
539 
540  if (compq == 'V')
541  std::cout << "QQ =" << QQ << std::endl;
542 #endif
543 
544  if (compq == 'V')
545  QQ = QQ * bqr.Q ();
546 
547 #if defined (DEBUG)
548  std::cout << "qz: precursors done..." << std::endl;
549 #endif
550 
551 #if defined (DEBUG)
552  std::cout << "qz: compq = " << compq << ", compz = " << compz
553  << std::endl;
554 #endif
555 
556  // Reduce to generalized Hessenberg form.
557  F77_XFCN (dgghrd, DGGHRD,
558  (F77_CONST_CHAR_ARG2 (&compq, 1),
559  F77_CONST_CHAR_ARG2 (&compz, 1),
560  nn, ilo, ihi, aa.fortran_vec (),
561  nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn,
562  ZZ.fortran_vec (), nn, info
563  F77_CHAR_ARG_LEN (1)
564  F77_CHAR_ARG_LEN (1)));
565 
566  // Check if just computing generalized eigenvalues or if we're
567  // actually computing the decomposition.
568 
569  // Reduce to generalized Schur form.
570  F77_XFCN (dhgeqz, DHGEQZ,
571  (F77_CONST_CHAR_ARG2 (&qz_job, 1),
572  F77_CONST_CHAR_ARG2 (&compq, 1),
573  F77_CONST_CHAR_ARG2 (&compz, 1),
574  nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (),
575  nn, alphar.fortran_vec (), alphai.fortran_vec (),
576  betar.fortran_vec (), QQ.fortran_vec (), nn,
577  ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info
578  F77_CHAR_ARG_LEN (1)
579  F77_CHAR_ARG_LEN (1)
580  F77_CHAR_ARG_LEN (1)));
581 
582  if (compq == 'V')
583  {
584  F77_XFCN (dggbak, DGGBAK,
585  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
586  F77_CONST_CHAR_ARG2 ("L", 1),
587  nn, ilo, ihi, lscale.data (), rscale.data (),
588  nn, QQ.fortran_vec (), nn, info
589  F77_CHAR_ARG_LEN (1)
590  F77_CHAR_ARG_LEN (1)));
591 
592 #if defined (DEBUG)
593  if (compq == 'V')
594  std::cout << "qz: balancing done; QQ=" << std::endl
595  << QQ << std::endl;
596 #endif
597  }
598 
599  // then right
600  if (compz == 'V')
601  {
602  F77_XFCN (dggbak, DGGBAK,
603  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
604  F77_CONST_CHAR_ARG2 ("R", 1),
605  nn, ilo, ihi, lscale.data (), rscale.data (),
606  nn, ZZ.fortran_vec (), nn, info
607  F77_CHAR_ARG_LEN (1)
608  F77_CHAR_ARG_LEN (1)));
609 
610 #if defined (DEBUG)
611  if (compz == 'V')
612  std::cout << "qz: balancing done; ZZ=" << std::endl
613  << ZZ << std::endl;
614 #endif
615  }
616 
617  }
618 
619  // Order the QZ decomposition?
620  if (! (ord_job == 'N' || ord_job == 'n'))
621  {
622  if (complex_case)
623  // Probably not needed, but better be safe.
624  error ("qz: cannot re-order complex qz decomposition");
625 
626 #if defined (DEBUG_SORT)
627  std::cout << "qz: ordering eigenvalues: ord_job = "
628  << ord_job << std::endl;
629 #endif
630 
631  // Declared static to avoid vfork/long jump compiler complaints.
632  static sort_function sort_test;
633  sort_test = 0;
634 
635  switch (ord_job)
636  {
637  case 'S':
638  case 's':
639  sort_test = &fin;
640  break;
641 
642  case 'B':
643  case 'b':
644  sort_test = &fout;
645  break;
646 
647  case '+':
648  sort_test = &fcrhp;
649  break;
650 
651  case '-':
652  sort_test = &folhp;
653  break;
654 
655  default:
656  // Invalid order option (should never happen, since we
657  // checked the options at the top).
658  panic_impossible ();
659  break;
660  }
661 
662  octave_idx_type ndim, fail;
663  double inf_norm;
664 
665  F77_XFCN (xdlange, XDLANGE,
666  (F77_CONST_CHAR_ARG2 ("I", 1),
667  nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm
668  F77_CHAR_ARG_LEN (1)));
669 
670  double eps = std::numeric_limits<double>::epsilon () * inf_norm * nn;
671 
672 #if defined (DEBUG_SORT)
673  std::cout << "qz: calling dsubsp: aa=" << std::endl;
674  octave_print_internal (std::cout, aa, 0);
675  std::cout << std::endl << "bb=" << std::endl;
676  octave_print_internal (std::cout, bb, 0);
677  if (compz == 'V')
678  {
679  std::cout << std::endl << "ZZ=" << std::endl;
680  octave_print_internal (std::cout, ZZ, 0);
681  }
682  std::cout << std::endl;
683  std::cout << "alphar = " << std::endl;
684  octave_print_internal (std::cout, (Matrix) alphar, 0);
685  std::cout << std::endl << "alphai = " << std::endl;
686  octave_print_internal (std::cout, (Matrix) alphai, 0);
687  std::cout << std::endl << "beta = " << std::endl;
688  octave_print_internal (std::cout, (Matrix) betar, 0);
689  std::cout << std::endl;
690 #endif
691 
693 
694  F77_XFCN (dsubsp, DSUBSP,
695  (nn, nn, aa.fortran_vec (), bb.fortran_vec (),
696  ZZ.fortran_vec (), sort_test, eps, ndim, fail,
697  ind.fortran_vec ()));
698 
699 #if defined (DEBUG)
700  std::cout << "qz: back from dsubsp: aa=" << std::endl;
701  octave_print_internal (std::cout, aa, 0);
702  std::cout << std::endl << "bb=" << std::endl;
703  octave_print_internal (std::cout, bb, 0);
704  if (compz == 'V')
705  {
706  std::cout << std::endl << "ZZ=" << std::endl;
707  octave_print_internal (std::cout, ZZ, 0);
708  }
709  std::cout << std::endl;
710 #endif
711 
712  // Manually update alphar, alphai, betar.
713  static int jj;
714 
715  jj = 0;
716  while (jj < nn)
717  {
718 #if defined (DEBUG_EIG)
719  std::cout << "computing gen eig #" << jj << std::endl;
720 #endif
721 
722  // Number of zeros in this block.
723  static int zcnt;
724 
725  if (jj == (nn-1))
726  zcnt = 1;
727  else if (aa(jj+1,jj) == 0)
728  zcnt = 1;
729  else zcnt = 2;
730 
731  if (zcnt == 1)
732  {
733  // Real zero.
734 #if defined (DEBUG_EIG)
735  std::cout << " single gen eig:" << std::endl;
736  std::cout << " alphar(" << jj << ") = " << aa(jj,jj)
737  << std::endl;
738  std::cout << " betar(" << jj << ") = " << bb(jj,jj)
739  << std::endl;
740  std::cout << " alphai(" << jj << ") = 0" << std::endl;
741 #endif
742 
743  alphar(jj) = aa(jj,jj);
744  alphai(jj) = 0;
745  betar(jj) = bb(jj,jj);
746  }
747  else
748  {
749  // Complex conjugate pair.
750 #if defined (DEBUG_EIG)
751  std::cout << "qz: calling dlag2:" << std::endl;
752  std::cout << "safmin="
753  << setiosflags (std::ios::scientific)
754  << safmin << std::endl;
755 
756  for (int idr = jj; idr <= jj+1; idr++)
757  {
758  for (int idc = jj; idc <= jj+1; idc++)
759  {
760  std::cout << "aa(" << idr << "," << idc << ")="
761  << aa(idr,idc) << std::endl;
762  std::cout << "bb(" << idr << "," << idc << ")="
763  << bb(idr,idc) << std::endl;
764  }
765  }
766 #endif
767 
768  // FIXME: probably should be using
769  // fortran_vec instead of &aa(jj,jj) here.
770 
771  double scale1, scale2, wr1, wr2, wi;
772  const double *aa_ptr = aa.data () + jj * nn + jj;
773  const double *bb_ptr = bb.data () + jj * nn + jj;
774  F77_XFCN (dlag2, DLAG2,
775  (aa_ptr, nn, bb_ptr, nn, safmin,
776  scale1, scale2, wr1, wr2, wi));
777 
778 #if defined (DEBUG_EIG)
779  std::cout << "dlag2 returns: scale1=" << scale1
780  << "\tscale2=" << scale2 << std::endl
781  << "\twr1=" << wr1 << "\twr2=" << wr2
782  << "\twi=" << wi << std::endl;
783 #endif
784 
785  // Just to be safe, check if it's a real pair.
786  if (wi == 0)
787  {
788  alphar(jj) = wr1;
789  alphai(jj) = 0;
790  betar(jj) = scale1;
791  alphar(jj+1) = wr2;
792  alphai(jj+1) = 0;
793  betar(jj+1) = scale2;
794  }
795  else
796  {
797  alphar(jj) = alphar(jj+1) = wr1;
798  alphai(jj) = -(alphai(jj+1) = wi);
799  betar(jj) = betar(jj+1) = scale1;
800  }
801  }
802 
803  // Advance past this block.
804  jj += zcnt;
805  }
806 
807 #if defined (DEBUG_SORT)
808  std::cout << "qz: back from dsubsp: aa=" << std::endl;
809  octave_print_internal (std::cout, aa, 0);
810  std::cout << std::endl << "bb=" << std::endl;
811  octave_print_internal (std::cout, bb, 0);
812 
813  if (compz == 'V')
814  {
815  std::cout << std::endl << "ZZ=" << std::endl;
816  octave_print_internal (std::cout, ZZ, 0);
817  }
818  std::cout << std::endl << "qz: ndim=" << ndim << std::endl
819  << "fail=" << fail << std::endl;
820  std::cout << "alphar = " << std::endl;
821  octave_print_internal (std::cout, (Matrix) alphar, 0);
822  std::cout << std::endl << "alphai = " << std::endl;
823  octave_print_internal (std::cout, (Matrix) alphai, 0);
824  std::cout << std::endl << "beta = " << std::endl;
825  octave_print_internal (std::cout, (Matrix) betar, 0);
826  std::cout << std::endl;
827 #endif
828  }
829 
830  // Compute generalized eigenvalues?
832 
833  if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4))
834  {
835  if (complex_case)
836  {
837  int cnt = 0;
838 
839  for (int ii = 0; ii < nn; ii++)
840  cnt++;
841 
842  ComplexColumnVector tmp (cnt);
843 
844  cnt = 0;
845  for (int ii = 0; ii < nn; ii++)
846  tmp(cnt++) = xalpha(ii) / xbeta(ii);
847 
848  gev = tmp;
849  }
850  else
851  {
852 #if defined (DEBUG)
853  std::cout << "qz: computing generalized eigenvalues" << std::endl;
854 #endif
855 
856  // Return finite generalized eigenvalues.
857  int cnt = 0;
858 
859  for (int ii = 0; ii < nn; ii++)
860  if (betar(ii) != 0)
861  cnt++;
862 
863  ComplexColumnVector tmp (cnt);
864 
865  cnt = 0;
866  for (int ii = 0; ii < nn; ii++)
867  if (betar(ii) != 0)
868  tmp(cnt++) = Complex(alphar(ii), alphai(ii))/betar(ii);
869 
870  gev = tmp;
871  }
872  }
873 
874  // Right, left eigenvector matrices.
875  if (nargout >= 5)
876  {
877  // Which side to compute?
878  char side = (nargout == 5 ? 'R' : 'B');
879  // Compute all of them and backtransform
880  char howmny = 'B';
881  // Dummy pointer; select is not used.
882  octave_idx_type *select = 0;
883 
884  if (complex_case)
885  {
886  CVL = CQ;
887  CVR = CZ;
888  ComplexRowVector cwork2 (2 * nn);
889  RowVector rwork2 (8 * nn);
891 
892  F77_XFCN (ztgevc, ZTGEVC,
893  (F77_CONST_CHAR_ARG2 (&side, 1),
894  F77_CONST_CHAR_ARG2 (&howmny, 1),
895  select, nn, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
897  nn, F77_DBLE_CMPLX_ARG (CVL.fortran_vec ()), nn,
898  F77_DBLE_CMPLX_ARG (CVR.fortran_vec ()), nn, nn,
899  m, F77_DBLE_CMPLX_ARG (cwork2.fortran_vec ()), rwork2.fortran_vec (), info
900  F77_CHAR_ARG_LEN (1)
901  F77_CHAR_ARG_LEN (1)));
902  }
903  else
904  {
905 #if defined (DEBUG)
906  std::cout << "qz: computing generalized eigenvectors" << std::endl;
907 #endif
908 
909  VL = QQ;
910  VR = ZZ;
912 
913  F77_XFCN (dtgevc, DTGEVC,
914  (F77_CONST_CHAR_ARG2 (&side, 1),
915  F77_CONST_CHAR_ARG2 (&howmny, 1),
916  select, nn, aa.fortran_vec (), nn, bb.fortran_vec (),
917  nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn,
918  m, work.fortran_vec (), info
919  F77_CHAR_ARG_LEN (1)
920  F77_CHAR_ARG_LEN (1)));
921 
922  // Now construct the complex form of VV, WW.
923  int jj = 0;
924 
925  while (jj < nn)
926  {
927  OCTAVE_QUIT;
928 
929  // See if real or complex eigenvalue.
930 
931  // Column increment; assume complex eigenvalue.
932  int cinc = 2;
933 
934  if (jj == (nn-1))
935  // Single column.
936  cinc = 1;
937  else if (aa(jj+1,jj) == 0)
938  cinc = 1;
939 
940  // Now copy the eigenvector (s) to CVR, CVL.
941  if (cinc == 1)
942  {
943  for (int ii = 0; ii < nn; ii++)
944  CVR(ii,jj) = VR(ii,jj);
945 
946  if (side == 'B')
947  for (int ii = 0; ii < nn; ii++)
948  CVL(ii,jj) = VL(ii,jj);
949  }
950  else
951  {
952  // Double column; complex vector.
953 
954  for (int ii = 0; ii < nn; ii++)
955  {
956  CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1));
957  CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1));
958  }
959 
960  if (side == 'B')
961  for (int ii = 0; ii < nn; ii++)
962  {
963  CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1));
964  CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1));
965  }
966  }
967 
968  // Advance to next eigenvectors (if any).
969  jj += cinc;
970  }
971  }
972  }
973 
974  switch (nargout)
975  {
976  case 7:
977  retval(6) = gev;
978 
979  case 6:
980  // Return eigenvectors.
981  retval(5) = CVL;
982 
983  case 5:
984  // Return eigenvectors.
985  retval(4) = CVR;
986 
987  case 4:
988  if (nargin == 3)
989  {
990 #if defined (DEBUG)
991  std::cout << "qz: sort: retval(3) = gev = " << std::endl;
992  octave_print_internal (std::cout, gev);
993  std::cout << std::endl;
994 #endif
995  retval(3) = gev;
996  }
997  else
998  {
999  if (complex_case)
1000  retval(3) = CZ;
1001  else
1002  retval(3) = ZZ;
1003  }
1004 
1005  case 3:
1006  if (nargin == 3)
1007  {
1008  if (complex_case)
1009  retval(2) = CZ;
1010  else
1011  retval(2) = ZZ;
1012  }
1013  else
1014  {
1015  if (complex_case)
1016  retval(2) = CQ.hermitian ();
1017  else
1018  retval(2) = QQ.transpose ();
1019  }
1020 
1021  case 2:
1022  {
1023  if (complex_case)
1024  {
1025 #if defined (DEBUG)
1026  std::cout << "qz: retval(1) = cbb = " << std::endl;
1027  octave_print_internal (std::cout, cbb, 0);
1028  std::cout << std::endl << "qz: retval(0) = caa = " <<std::endl;
1029  octave_print_internal (std::cout, caa, 0);
1030  std::cout << std::endl;
1031 #endif
1032  retval(1) = cbb;
1033  retval(0) = caa;
1034  }
1035  else
1036  {
1037 #if defined (DEBUG)
1038  std::cout << "qz: retval(1) = bb = " << std::endl;
1039  octave_print_internal (std::cout, bb, 0);
1040  std::cout << std::endl << "qz: retval(0) = aa = " <<std::endl;
1041  octave_print_internal (std::cout, aa, 0);
1042  std::cout << std::endl;
1043 #endif
1044  retval(1) = bb;
1045  retval(0) = aa;
1046  }
1047  }
1048  break;
1049 
1050  case 1:
1051  case 0:
1052 #if defined (DEBUG)
1053  std::cout << "qz: retval(0) = gev = " << gev << std::endl;
1054 #endif
1055  retval(0) = gev;
1056  break;
1057 
1058  default:
1059  error ("qz: too many return arguments");
1060  break;
1061  }
1062 
1063 #if defined (DEBUG)
1064  std::cout << "qz: exiting (at long last)" << std::endl;
1065 #endif
1066 
1067  return retval;
1068 }
1069 
1070 /*
1071 %!shared a, b, c
1072 %! a = [1 2; 0 3];
1073 %! b = [1 0; 0 0];
1074 %! c = [0 1; 0 0];
1075 %!assert (qz (a,b), 1)
1076 %!assert (isempty (qz (a,c)))
1077 
1078 ## Exaple 7.7.3 in Golub & Van Loan
1079 %!test
1080 %! a = [ 10 1 2;
1081 %! 1 2 -1;
1082 %! 1 1 2];
1083 %! b = reshape (1:9,3,3);
1084 %! [aa, bb, q, z, v, w, lambda] = qz (a, b);
1085 %! sz = length (lambda);
1086 %! observed = (b * v * diag ([lambda;0])) (:, 1:sz);
1087 %! assert ((a*v)(:, 1:sz), observed, norm (observed) * 1e-14);
1088 %! observed = (diag ([lambda;0]) * w' * b) (1:sz, :);
1089 %! assert ((w'*a)(1:sz, :) , observed, norm (observed) * 1e-13);
1090 %! assert (q * a * z, aa, norm (aa) * 1e-14);
1091 %! assert (q * b * z, bb, norm (bb) * 1e-14);
1092 
1093 %!test
1094 %! A = [0, 0, -1, 0; 1, 0, 0, 0; -1, 0, -2, -1; 0, -1, 1, 0];
1095 %! B = [0, 0, 0, 0; 0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1];
1096 %! [AA, BB, Q, Z1] = qz (A, B);
1097 %! [AA, BB, Z2] = qz (A, B, '-');
1098 %! assert (Z1, Z2);
1099 */
void warn_empty_arg(const char *name)
Definition: errwarn.cc:326
octave_idx_type(* sort_function)(const octave_idx_type &LSIZE, const double &ALPHA, const double &BETA, const double &S, const double &P)
Definition: qz.cc:60
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
subroutine xdlamch(cmach, retval)
Definition: xdlamch.f:1
T Q(void) const
Definition: qr.h:78
F77_RET_T const F77_INT & N
Definition: qz.cc:71
ar P
Definition: __luinc__.cc:49
ind
Definition: sub2ind.cc:107
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
F77_RET_T const F77_INT F77_DBLE * A
Definition: qz.cc:71
static octave_idx_type nn
Definition: DASPK.cc:71
F77_RET_T const F77_INT F77_DBLE F77_DBLE F77_DBLE const F77_DBLE F77_INT F77_INT & FAIL
Definition: qz.cc:71
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:345
subroutine xdlange(norm, m, n, a, lda, work, retval)
Definition: xdlange.f:1
F77_RET_T const F77_INT F77_DBLE F77_DBLE * B
Definition: qz.cc:71
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:382
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
void error(const char *fmt,...)
Definition: error.cc:570
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 * ALPHA
F77_RET_T const F77_INT F77_DBLE F77_DBLE F77_DBLE const F77_DBLE & EPS
Definition: qz.cc:71
void err_square_matrix_required(const char *fcn, const char *name)
Definition: errwarn.cc:112
s
Definition: file-io.cc:2682
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:52
JNIEnv void * args
Definition: ov-java.cc:67
F77_RET_T const F77_INT F77_DBLE F77_DBLE F77_DBLE * Z
Definition: qz.cc:71
ComplexMatrix hermitian(void) const
Definition: CMatrix.h:163
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function xample nargout(@histc)
Definition: ov-usr-fcn.cc:935
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
F77_RET_T F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * VL
static octave_idx_type fcrhp(const octave_idx_type &lsize, const double &alpha, const double &beta, const double &s, const double &)
Definition: qz.cc:86
int nargin
Definition: graphics.cc:10115
const T * data(void) const
Definition: Array.h:582
double tmp
Definition: data.cc:6300
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
octave_value retval
Definition: data.cc:6294
#define panic_impossible()
Definition: error.h:40
subroutine dsubsp(NMAX, N, A, B, Z, FTEST, EPS, NDIM, FAIL, IND)
Definition: dsubsp.f:1
F77_RET_T const F77_INT F77_DBLE F77_DBLE F77_DBLE const F77_DBLE F77_INT F77_INT F77_INT * IND
Definition: qz.cc:71
Definition: dMatrix.h:37
#define F77_INT
Definition: f77-fcn.h:335
#define eps(C)
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE F77_DBLE F77_DBLE * BETA
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
static octave_idx_type fout(const octave_idx_type &lsize, const double &alpha, const double &beta, const double &, const double &p)
Definition: qz.cc:124
static octave_idx_type fin(const octave_idx_type &lsize, const double &alpha, const double &beta, const double &, const double &p)
Definition: qz.cc:96
static double wi[256]
Definition: randmtzig.cc:440
p
Definition: lu.cc:138
void octave_print_internal(std::ostream &, char, bool)
Definition: pr-output.cc:1722
static octave_idx_type folhp(const octave_idx_type &lsize, const double &alpha, const double &beta, const double &s, const double &)
Definition: qz.cc:114
T R(void) const
Definition: qr.h:80
F77_RET_T F77_FUNC(dsubsp, DSUBSP)(const F77_INT &NMAX
std::complex< double > Complex
Definition: oct-cmplx.h:31
const T * fortran_vec(void) const
Definition: Array.h:584
F77_RET_T const F77_INT F77_DBLE F77_DBLE F77_DBLE const F77_DBLE F77_INT & NDIM
Definition: qz.cc:71
#define OCTAVE_QUIT
Definition: quit.h:212
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
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
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string
Definition: utils.cc:854
#define F77_DBLE
Definition: f77-fcn.h:331