GNU Octave  3.8.0
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-2013 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 <[email protected]>
26 
27 #undef DEBUG
28 #undef DEBUG_SORT
29 #undef DEBUG_EIG
30 
31 #ifdef HAVE_CONFIG_H
32 #include <config.h>
33 #endif
34 
35 #include <cfloat>
36 
37 #include <iostream>
38 #include <iomanip>
39 
40 #include "CmplxQRP.h"
41 #include "CmplxQR.h"
42 #include "dbleQR.h"
43 #include "f77-fcn.h"
44 #include "lo-math.h"
45 #include "quit.h"
46 
47 #include "defun.h"
48 #include "error.h"
49 #include "gripes.h"
50 #include "oct-obj.h"
51 #include "oct-map.h"
52 #include "ov.h"
53 #include "pager.h"
54 #if defined (DEBUG) || defined (DEBUG_SORT)
55 #include "pr-output.h"
56 #endif
57 #include "symtab.h"
58 #include "utils.h"
59 #include "variables.h"
60 
61 typedef octave_idx_type (*sort_function) (const octave_idx_type& LSIZE,
62  const double& ALPHA,
63  const double& BETA, const double& S,
64  const double& P);
65 
66 extern "C"
67 {
68  F77_RET_T
69  F77_FUNC (dggbal, DGGBAL) (F77_CONST_CHAR_ARG_DECL,
70  const octave_idx_type& N, double* A,
71  const octave_idx_type& LDA, double* B,
72  const octave_idx_type& LDB, octave_idx_type& ILO,
73  octave_idx_type& IHI, double* LSCALE,
74  double* RSCALE, double* WORK,
75  octave_idx_type& INFO
76  F77_CHAR_ARG_LEN_DECL);
77 
78 F77_RET_T
79  F77_FUNC (zggbal, ZGGBAL) (F77_CONST_CHAR_ARG_DECL,
80  const octave_idx_type& N, Complex* A,
81  const octave_idx_type& LDA, Complex* B,
82  const octave_idx_type& LDB, octave_idx_type& ILO,
83  octave_idx_type& IHI, double* LSCALE,
84  double* RSCALE, double* WORK,
85  octave_idx_type& INFO
86  F77_CHAR_ARG_LEN_DECL);
87 
88  F77_RET_T
89  F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL,
90  F77_CONST_CHAR_ARG_DECL,
91  const octave_idx_type& N,
92  const octave_idx_type& ILO,
93  const octave_idx_type& IHI,
94  const double* LSCALE, const double* RSCALE,
95  octave_idx_type& M, double* V,
96  const octave_idx_type& LDV, octave_idx_type& INFO
97  F77_CHAR_ARG_LEN_DECL
98  F77_CHAR_ARG_LEN_DECL);
99 
100 F77_RET_T
101  F77_FUNC (zggbak, ZGGBAK) (F77_CONST_CHAR_ARG_DECL,
102  F77_CONST_CHAR_ARG_DECL,
103  const octave_idx_type& N,
104  const octave_idx_type& ILO,
105  const octave_idx_type& IHI,
106  const double* LSCALE, const double* RSCALE,
107  octave_idx_type& M, Complex* V,
108  const octave_idx_type& LDV, octave_idx_type& INFO
109  F77_CHAR_ARG_LEN_DECL
110  F77_CHAR_ARG_LEN_DECL);
111 
112  F77_RET_T
113  F77_FUNC (dgghrd, DGGHRD) (F77_CONST_CHAR_ARG_DECL,
114  F77_CONST_CHAR_ARG_DECL,
115  const octave_idx_type& N,
116  const octave_idx_type& ILO,
117  const octave_idx_type& IHI, double* A,
118  const octave_idx_type& LDA, double* B,
119  const octave_idx_type& LDB, double* Q,
120  const octave_idx_type& LDQ, double* Z,
121  const octave_idx_type& LDZ, octave_idx_type& INFO
122  F77_CHAR_ARG_LEN_DECL
123  F77_CHAR_ARG_LEN_DECL);
124 
125  F77_RET_T
126  F77_FUNC (zgghrd, ZGGHRD) (F77_CONST_CHAR_ARG_DECL,
127  F77_CONST_CHAR_ARG_DECL,
128  const octave_idx_type& N,
129  const octave_idx_type& ILO,
130  const octave_idx_type& IHI, Complex* A,
131  const octave_idx_type& LDA, Complex* B,
132  const octave_idx_type& LDB, Complex* Q,
133  const octave_idx_type& LDQ, Complex* Z,
134  const octave_idx_type& LDZ, octave_idx_type& INFO
135  F77_CHAR_ARG_LEN_DECL
136  F77_CHAR_ARG_LEN_DECL);
137 
138  F77_RET_T
139  F77_FUNC (dhgeqz, DHGEQZ) (F77_CONST_CHAR_ARG_DECL,
140  F77_CONST_CHAR_ARG_DECL,
141  F77_CONST_CHAR_ARG_DECL,
142  const octave_idx_type& N,
143  const octave_idx_type& ILO,
144  const octave_idx_type& IHI,
145  double* A, const octave_idx_type& LDA, double* B,
146  const octave_idx_type& LDB, double* ALPHAR,
147  double* ALPHAI, double* BETA, double* Q,
148  const octave_idx_type& LDQ, double* Z,
149  const octave_idx_type& LDZ, double* WORK,
150  const octave_idx_type& LWORK,
151  octave_idx_type& INFO
152  F77_CHAR_ARG_LEN_DECL
153  F77_CHAR_ARG_LEN_DECL
154  F77_CHAR_ARG_LEN_DECL);
155 
156 F77_RET_T
157  F77_FUNC (zhgeqz, ZHGEQZ) (F77_CONST_CHAR_ARG_DECL,
158  F77_CONST_CHAR_ARG_DECL,
159  F77_CONST_CHAR_ARG_DECL,
160  const octave_idx_type& N,
161  const octave_idx_type& ILO,
162  const octave_idx_type& IHI,
163  Complex* A, const octave_idx_type& LDA,
164  Complex* B, const octave_idx_type& LDB,
165  Complex* ALPHA, Complex* BETA, Complex* CQ,
166  const octave_idx_type& LDQ,
167  Complex* CZ, const octave_idx_type& LDZ,
168  Complex* WORK, const octave_idx_type& LWORK,
169  double* RWORK, octave_idx_type& INFO
170  F77_CHAR_ARG_LEN_DECL
171  F77_CHAR_ARG_LEN_DECL
172  F77_CHAR_ARG_LEN_DECL);
173 
174  F77_RET_T
175  F77_FUNC (dlag2, DLAG2) (const double* A, const octave_idx_type& LDA,
176  const double* B, const octave_idx_type& LDB,
177  const double& SAFMIN, double& SCALE1,
178  double& SCALE2, double& WR1, double& WR2,
179  double& WI);
180 
181  // Van Dooren's code (netlib.org: toms/590) for reordering
182  // GEP. Only processes Z, not Q.
183  F77_RET_T
184  F77_FUNC (dsubsp, DSUBSP) (const octave_idx_type& NMAX,
185  const octave_idx_type& N, double* A,
186  double* B, double* Z, sort_function,
187  const double& EPS, octave_idx_type& NDIM,
188  octave_idx_type& FAIL, octave_idx_type* IND);
189 
190  // Documentation for DTGEVC incorrectly states that VR, VL are
191  // complex*16; they are declared in DTGEVC as double precision
192  // (probably a cut and paste problem fro ZTGEVC).
193  F77_RET_T
194  F77_FUNC (dtgevc, DTGEVC) (F77_CONST_CHAR_ARG_DECL,
195  F77_CONST_CHAR_ARG_DECL,
196  octave_idx_type* SELECT,
197  const octave_idx_type& N, double* A,
198  const octave_idx_type& LDA, double* B,
199  const octave_idx_type& LDB, double* VL,
200  const octave_idx_type& LDVL, double* VR,
201  const octave_idx_type& LDVR,
202  const octave_idx_type& MM, octave_idx_type& M,
203  double* WORK, octave_idx_type& INFO
204  F77_CHAR_ARG_LEN_DECL
205  F77_CHAR_ARG_LEN_DECL);
206 
207 F77_RET_T
208  F77_FUNC (ztgevc, ZTGEVC) (F77_CONST_CHAR_ARG_DECL,
209  F77_CONST_CHAR_ARG_DECL,
210  octave_idx_type* SELECT,
211  const octave_idx_type& N, const Complex* A,
212  const octave_idx_type& LDA,const Complex* B,
213  const octave_idx_type& LDB, Complex* xVL,
214  const octave_idx_type& LDVL, Complex* xVR,
215  const octave_idx_type& LDVR,
216  const octave_idx_type& MM, octave_idx_type& M,
217  Complex* CWORK, double* RWORK,
218  octave_idx_type& INFO
219  F77_CHAR_ARG_LEN_DECL
220  F77_CHAR_ARG_LEN_DECL);
221 
222  F77_RET_T
223  F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG_DECL,
224  double& retval
225  F77_CHAR_ARG_LEN_DECL);
226 
227  F77_RET_T
228  F77_FUNC (xdlange, XDLANGE) (F77_CONST_CHAR_ARG_DECL,
229  const octave_idx_type&,
230  const octave_idx_type&, const double*,
231  const octave_idx_type&, double*, double&
232  F77_CHAR_ARG_LEN_DECL);
233 }
234 
235 // fcrhp, fin, fout, folhp:
236 // routines for ordering of generalized eigenvalues
237 // return 1 if test is passed, 0 otherwise
238 // fin: |lambda| < 1
239 // fout: |lambda| >= 1
240 // fcrhp: real(lambda) >= 0
241 // folhp: real(lambda) < 0
242 
243 static octave_idx_type
244 fcrhp (const octave_idx_type& lsize, const double& alpha,
245  const double& beta, const double& s, const double&)
246 {
247  if (lsize == 1)
248  return (alpha * beta >= 0 ? 1 : -1);
249  else
250  return (s >= 0 ? 1 : -1);
251 }
252 
253 static octave_idx_type
254 fin (const octave_idx_type& lsize, const double& alpha,
255  const double& beta, const double&, const double& p)
256 {
257  octave_idx_type retval;
258 
259  if (lsize == 1)
260  retval = (fabs (alpha) < fabs (beta) ? 1 : -1);
261  else
262  retval = (fabs (p) < 1 ? 1 : -1);
263 
264 #ifdef DEBUG
265  std::cout << "qz: fin: retval=" << retval << std::endl;
266 #endif
267 
268  return retval;
269 }
270 
271 static octave_idx_type
272 folhp (const octave_idx_type& lsize, const double& alpha,
273  const double& beta, const double& s, const double&)
274 {
275  if (lsize == 1)
276  return (alpha * beta < 0 ? 1 : -1);
277  else
278  return (s < 0 ? 1 : -1);
279 }
280 
281 static octave_idx_type
282 fout (const octave_idx_type& lsize, const double& alpha,
283  const double& beta, const double&, const double& p)
284 {
285  if (lsize == 1)
286  return (fabs (alpha) >= fabs (beta) ? 1 : -1);
287  else
288  return (fabs (p) >= 1 ? 1 : -1);
289 }
290 
291 
292 //FIXME: Matlab does not produce lambda as the first output argument.
293 // Compatibility problem?
294 DEFUN (qz, args, nargout,
295  "-*- texinfo -*-\n\
296 @deftypefn {Built-in Function} {@var{lambda} =} qz (@var{A}, @var{B})\n\
297 @deftypefnx {Built-in Function} {@var{lambda} =} qz (@var{A}, @var{B}, @var{opt})\n\
298 QZ@tie{}decomposition of the generalized eigenvalue problem\n\
299 (@math{A x = s B x}). There are three ways to call this function:\n\
300 @enumerate\n\
301 @item @code{@var{lambda} = qz (@var{A}, @var{B})}\n\
302 \n\
303 Computes the generalized eigenvalues\n\
304 @tex\n\
305 $\\lambda$\n\
306 @end tex\n\
307 @ifnottex\n\
308 @var{lambda}\n\
309 @end ifnottex\n\
310 of @math{(A - s B)}.\n\
311 \n\
312 @item @code{[AA, BB, Q, Z, V, W, @var{lambda}] = qz (@var{A}, @var{B})}\n\
313 \n\
314 Computes QZ@tie{}decomposition, generalized eigenvectors, and\n\
315 generalized eigenvalues of @math{(A - s B)}\n\
316 @tex\n\
317 $$ AV = BV{ \\rm diag }(\\lambda) $$\n\
318 $$ W^T A = { \\rm diag }(\\lambda)W^T B $$\n\
319 $$ AA = Q^T AZ, BB = Q^T BZ $$\n\
320 @end tex\n\
321 @ifnottex\n\
322 \n\
323 @example\n\
324 @group\n\
325 \n\
326 A * V = B * V * diag (@var{lambda})\n\
327 W' * A = diag (@var{lambda}) * W' * B\n\
328 AA = Q * A * Z, BB = Q * B * Z\n\
329 \n\
330 @end group\n\
331 @end example\n\
332 \n\
333 @end ifnottex\n\
334 with @var{Q} and @var{Z} orthogonal (unitary)= @var{I}\n\
335 \n\
336 @item @code{[AA,BB,Z@{, @var{lambda}@}] = qz (@var{A}, @var{B}, @var{opt})}\n\
337 \n\
338 As in form [2], but allows ordering of generalized eigenpairs\n\
339 for (e.g.) solution of discrete time algebraic Riccati equations.\n\
340 Form 3 is not available for complex matrices, and does not compute\n\
341 the generalized eigenvectors @var{V}, @var{W}, nor the orthogonal matrix\n\
342 @var{Q}.\n\
343 \n\
344 @table @var\n\
345 @item opt\n\
346 for ordering eigenvalues of the @nospell{GEP} pencil. The leading block\n\
347 of the revised pencil contains all eigenvalues that satisfy:\n\
348 \n\
349 @table @asis\n\
350 @item @qcode{\"N\"}\n\
351 = unordered (default)\n\
352 \n\
353 @item @qcode{\"S\"}\n\
354 = small: leading block has all |lambda| @leq{} 1\n\
355 \n\
356 @item @qcode{\"B\"}\n\
357 = big: leading block has all |lambda| @geq{} 1\n\
358 \n\
359 @item @qcode{\"-\"}\n\
360 = negative real part: leading block has all eigenvalues\n\
361 in the open left half-plane\n\
362 \n\
363 @item @qcode{\"+\"}\n\
364 = non-negative real part: leading block has all eigenvalues\n\
365 in the closed right half-plane\n\
366 @end table\n\
367 @end table\n\
368 @end enumerate\n\
369 \n\
370 Note: @code{qz} performs permutation balancing, but not scaling\n\
371 (@pxref{XREFbalance}). The order of output arguments was selected for\n\
372 compatibility with @sc{matlab}.\n\
373 @seealso{eig, balance, lu, chol, hess, qr, qzhess, schur, svd}\n\
374 @end deftypefn")
375 {
376  octave_value_list retval;
377  int nargin = args.length ();
378 
379 #ifdef DEBUG
380  std::cout << "qz: nargin = " << nargin
381  << ", nargout = " << nargout << std::endl;
382 #endif
383 
384  if (nargin < 2 || nargin > 3 || nargout > 7)
385  {
386  print_usage ();
387  return retval;
388  }
389  else if (nargin == 3 && (nargout < 3 || nargout > 4))
390  {
391  error ("qz: invalid number of output arguments for form [3] call");
392  return retval;
393  }
394 
395 #ifdef DEBUG
396  std::cout << "qz: determine ordering option" << std::endl;
397 #endif
398 
399  // Determine ordering option.
400  volatile char ord_job = 0;
401  static double safmin;
402 
403  if (nargin == 2)
404  ord_job = 'N';
405  else if (!args(2).is_string ())
406  {
407  error ("qz: OPT must be a string");
408  return retval;
409  }
410  else
411  {
412  std::string tmp = args(2).string_value ();
413 
414  if (! tmp.empty ())
415  ord_job = tmp[0];
416 
417  if (! (ord_job == 'N' || ord_job == 'n'
418  || ord_job == 'S' || ord_job == 's'
419  || ord_job == 'B' || ord_job == 'b'
420  || ord_job == '+' || ord_job == '-'))
421  {
422  error ("qz: invalid order option");
423  return retval;
424  }
425 
426  // overflow constant required by dlag2
427  F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1),
428  safmin
429  F77_CHAR_ARG_LEN (1));
430 
431 #ifdef DEBUG_EIG
432  std::cout << "qz: initial value of safmin="
433  << setiosflags (std::ios::scientific)
434  << safmin << std::endl;
435 #endif
436 
437  // Some machines (e.g., DEC alpha) get safmin = 0;
438  // for these, use eps instead to avoid problems in dlag2.
439  if (safmin == 0)
440  {
441 #ifdef DEBUG_EIG
442  std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl;
443 #endif
444 
445  F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("E", 1),
446  safmin
447  F77_CHAR_ARG_LEN (1));
448 
449 #ifdef DEBUG_EIG
450  std::cout << "qz: safmin set to "
451  << setiosflags (std::ios::scientific)
452  << safmin << std::endl;
453 #endif
454  }
455  }
456 
457 #ifdef DEBUG
458  std::cout << "qz: check argument 1" << std::endl;
459 #endif
460 
461  // Argument 1: check if it's o.k. dimensioned.
462  octave_idx_type nn = args(0).rows ();
463 
464 #ifdef DEBUG
465  std::cout << "argument 1 dimensions: ("
466  << nn << "," << args(0).columns () << ")"
467  << std::endl;
468 #endif
469 
470  int arg_is_empty = empty_arg ("qz", nn, args(0).columns ());
471 
472  if (arg_is_empty < 0)
473  {
474  gripe_empty_arg ("qz: parameter 1", 0);
475  return retval;
476  }
477  else if (arg_is_empty > 0)
478  {
479  gripe_empty_arg ("qz: parameter 1; continuing", 0);
480  return octave_value_list (2, Matrix ());
481  }
482  else if (args(0).columns () != nn)
483  {
484  gripe_square_matrix_required ("qz");
485  return retval;
486  }
487 
488  // Argument 1: dimensions look good; get the value.
489  Matrix aa;
490  ComplexMatrix caa;
491 
492  if (args(0).is_complex_type ())
493  caa = args(0).complex_matrix_value ();
494  else
495  aa = args(0).matrix_value ();
496 
497  if (error_state)
498  return retval;
499 
500 #ifdef DEBUG
501  std::cout << "qz: check argument 2" << std::endl;
502 #endif
503 
504  // Extract argument 2 (bb, or cbb if complex).
505  if ((nn != args(1).columns ()) || (nn != args(1).rows ()))
506  {
507  gripe_nonconformant ();
508  return retval;
509  }
510 
511  Matrix bb;
512  ComplexMatrix cbb;
513 
514  if (args(1).is_complex_type ())
515  cbb = args(1).complex_matrix_value ();
516  else
517  bb = args(1).matrix_value ();
518 
519  if (error_state)
520  return retval;
521 
522  // Both matrices loaded, now let's check what kind of arithmetic:
523  // declared volatile to avoid compiler warnings about long jumps,
524  // vforks.
525 
526  volatile int complex_case
527  = (args(0).is_complex_type () || args(1).is_complex_type ());
528 
529  if (nargin == 3 && complex_case)
530  {
531  error ("qz: cannot re-order complex qz decomposition");
532  return retval;
533  }
534 
535  // First, declare variables used in both the real and complex case.
536  Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn);
537  RowVector alphar(nn), alphai(nn), betar(nn);
538  ComplexRowVector xalpha(nn), xbeta(nn);
539  ComplexMatrix CQ(nn,nn), CZ(nn,nn), CVR(nn,nn), CVL(nn,nn);
540  octave_idx_type ilo, ihi, info;
541  char compq = (nargout >= 3 ? 'V' : 'N');
542  char compz = ((nargout >= 4 || nargin == 3)? 'V' : 'N');
543 
544  // Initialize Q, Z to identity if we need either of them.
545  if (compq == 'V' || compz == 'V')
546  for (octave_idx_type ii = 0; ii < nn; ii++)
547  for (octave_idx_type jj = 0; jj < nn; jj++)
548  {
549  OCTAVE_QUIT;
550  QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0);
551  }
552 
553  // Always perform permutation balancing.
554  const char bal_job = 'P';
555  RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn);
556 
557  if (complex_case)
558  {
559 #ifdef DEBUG
560  if (compq == 'V')
561  std::cout << "qz: performing balancing; CQ=" << std::endl
562  << CQ << std::endl;
563 #endif
564  if (args(0).is_real_type ())
565  caa = ComplexMatrix (aa);
566 
567  if (args(1).is_real_type ())
568  cbb = ComplexMatrix (bb);
569 
570  if (compq == 'V')
571  CQ = ComplexMatrix (QQ);
572 
573  if (compz == 'V')
574  CZ = ComplexMatrix (ZZ);
575 
576  F77_XFCN (zggbal, ZGGBAL,
577  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
578  nn, caa.fortran_vec (), nn, cbb.fortran_vec (),
579  nn, ilo, ihi, lscale.fortran_vec (),
580  rscale.fortran_vec (), work.fortran_vec (), info
581  F77_CHAR_ARG_LEN (1)));
582  }
583  else
584  {
585 #ifdef DEBUG
586  if (compq == 'V')
587  std::cout << "qz: performing balancing; QQ=" << std::endl
588  << QQ << std::endl;
589 #endif
590 
591  F77_XFCN (dggbal, DGGBAL,
592  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
593  nn, aa.fortran_vec (), nn, bb.fortran_vec (),
594  nn, ilo, ihi, lscale.fortran_vec (),
595  rscale.fortran_vec (), work.fortran_vec (), info
596  F77_CHAR_ARG_LEN (1)));
597  }
598 
599  // Since we just want the balancing matrices, we can use dggbal
600  // for both the real and complex cases; left first
601 
602 #if 0
603  if (compq == 'V')
604  {
605  F77_XFCN (dggbak, DGGBAK,
606  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
607  F77_CONST_CHAR_ARG2 ("L", 1),
608  nn, ilo, ihi, lscale.data (), rscale.data (),
609  nn, QQ.fortran_vec (), nn, info
610  F77_CHAR_ARG_LEN (1)
611  F77_CHAR_ARG_LEN (1)));
612 
613 #ifdef DEBUG
614  if (compq == 'V')
615  std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl;
616 #endif
617  }
618 
619  // then right
620  if (compz == 'V')
621  {
622  F77_XFCN (dggbak, DGGBAK,
623  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
624  F77_CONST_CHAR_ARG2 ("R", 1),
625  nn, ilo, ihi, lscale.data (), rscale.data (),
626  nn, ZZ.fortran_vec (), nn, info
627  F77_CHAR_ARG_LEN (1)
628  F77_CHAR_ARG_LEN (1)));
629 
630 #ifdef DEBUG
631  if (compz == 'V')
632  std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
633 #endif
634  }
635 #endif
636 
637  static char qz_job;
638  qz_job = (nargout < 2 ? 'E' : 'S');
639 
640  if (complex_case)
641  {
642  // Complex case.
643 
644  // The QR decomposition of cbb.
645  ComplexQR cbqr (cbb);
646  // The R matrix of QR decomposition for cbb.
647  cbb = cbqr.R ();
648  // (Q*)caa for following work.
649  caa = (cbqr.Q ().hermitian ()) * caa;
650  CQ = CQ * cbqr.Q ();
651 
652  F77_XFCN (zgghrd, ZGGHRD,
653  (F77_CONST_CHAR_ARG2 (&compq, 1),
654  F77_CONST_CHAR_ARG2 (&compz, 1),
655  nn, ilo, ihi, caa.fortran_vec (),
656  nn, cbb.fortran_vec (), nn, CQ.fortran_vec (), nn,
657  CZ.fortran_vec (), nn, info
658  F77_CHAR_ARG_LEN (1)
659  F77_CHAR_ARG_LEN (1)));
660 
661  ComplexRowVector cwork (1 * nn);
662 
663  F77_XFCN (zhgeqz, ZHGEQZ,
664  (F77_CONST_CHAR_ARG2 (&qz_job, 1),
665  F77_CONST_CHAR_ARG2 (&compq, 1),
666  F77_CONST_CHAR_ARG2 (&compz, 1),
667  nn, ilo, ihi,
668  caa.fortran_vec (), nn,
669  cbb.fortran_vec (),nn,
670  xalpha.fortran_vec (), xbeta.fortran_vec (),
671  CQ.fortran_vec (), nn,
672  CZ.fortran_vec (), nn,
673  cwork.fortran_vec (), nn, rwork.fortran_vec (), info
674  F77_CHAR_ARG_LEN (1)
675  F77_CHAR_ARG_LEN (1)
676  F77_CHAR_ARG_LEN (1)));
677 
678  if (compq == 'V')
679  {
680  // Left eigenvector.
681  F77_XFCN (zggbak, ZGGBAK,
682  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
683  F77_CONST_CHAR_ARG2 ("L", 1),
684  nn, ilo, ihi, lscale.data (), rscale.data (),
685  nn, CQ.fortran_vec (), nn, info
686  F77_CHAR_ARG_LEN (1)
687  F77_CHAR_ARG_LEN (1)));
688  }
689 
690  // Right eigenvector.
691  if (compz == 'V')
692  {
693  F77_XFCN (zggbak, ZGGBAK,
694  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
695  F77_CONST_CHAR_ARG2 ("R", 1),
696  nn, ilo, ihi, lscale.data (), rscale.data (),
697  nn, CZ.fortran_vec (), nn, info
698  F77_CHAR_ARG_LEN (1)
699  F77_CHAR_ARG_LEN (1)));
700  }
701 
702  }
703  else
704  {
705 #ifdef DEBUG
706  std::cout << "qz: peforming qr decomposition of bb" << std::endl;
707 #endif
708 
709  // Compute the QR factorization of bb.
710  QR bqr (bb);
711 
712 #ifdef DEBUG
713  std::cout << "qz: qr (bb) done; now peforming qz decomposition"
714  << std::endl;
715 #endif
716 
717  bb = bqr.R ();
718 
719 #ifdef DEBUG
720  std::cout << "qz: extracted bb" << std::endl;
721 #endif
722 
723  aa = (bqr.Q ()).transpose () * aa;
724 
725 #ifdef DEBUG
726  std::cout << "qz: updated aa " << std::endl;
727  std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl;
728 
729  if (compq == 'V')
730  std::cout << "QQ =" << QQ << std::endl;
731 #endif
732 
733  if (compq == 'V')
734  QQ = QQ * bqr.Q ();
735 
736 #ifdef DEBUG
737  std::cout << "qz: precursors done..." << std::endl;
738 #endif
739 
740 #ifdef DEBUG
741  std::cout << "qz: compq = " << compq << ", compz = " << compz
742  << std::endl;
743 #endif
744 
745  // Reduce to generalized hessenberg form.
746  F77_XFCN (dgghrd, DGGHRD,
747  (F77_CONST_CHAR_ARG2 (&compq, 1),
748  F77_CONST_CHAR_ARG2 (&compz, 1),
749  nn, ilo, ihi, aa.fortran_vec (),
750  nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn,
751  ZZ.fortran_vec (), nn, info
752  F77_CHAR_ARG_LEN (1)
753  F77_CHAR_ARG_LEN (1)));
754 
755  // Check if just computing generalized eigenvalues or if we're
756  // actually computing the decomposition.
757 
758  // Reduce to generalized Schur form.
759  F77_XFCN (dhgeqz, DHGEQZ,
760  (F77_CONST_CHAR_ARG2 (&qz_job, 1),
761  F77_CONST_CHAR_ARG2 (&compq, 1),
762  F77_CONST_CHAR_ARG2 (&compz, 1),
763  nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (),
764  nn, alphar.fortran_vec (), alphai.fortran_vec (),
765  betar.fortran_vec (), QQ.fortran_vec (), nn,
766  ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info
767  F77_CHAR_ARG_LEN (1)
768  F77_CHAR_ARG_LEN (1)
769  F77_CHAR_ARG_LEN (1)));
770 
771  if (compq == 'V')
772  {
773  F77_XFCN (dggbak, DGGBAK,
774  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
775  F77_CONST_CHAR_ARG2 ("L", 1),
776  nn, ilo, ihi, lscale.data (), rscale.data (),
777  nn, QQ.fortran_vec (), nn, info
778  F77_CHAR_ARG_LEN (1)
779  F77_CHAR_ARG_LEN (1)));
780 
781 #ifdef DEBUG
782  if (compq == 'V')
783  std::cout << "qz: balancing done; QQ=" << std::endl
784  << QQ << std::endl;
785 #endif
786  }
787 
788  // then right
789  if (compz == 'V')
790  {
791  F77_XFCN (dggbak, DGGBAK,
792  (F77_CONST_CHAR_ARG2 (&bal_job, 1),
793  F77_CONST_CHAR_ARG2 ("R", 1),
794  nn, ilo, ihi, lscale.data (), rscale.data (),
795  nn, ZZ.fortran_vec (), nn, info
796  F77_CHAR_ARG_LEN (1)
797  F77_CHAR_ARG_LEN (1)));
798 
799 #ifdef DEBUG
800  if (compz == 'V')
801  std::cout << "qz: balancing done; ZZ=" << std::endl
802  << ZZ << std::endl;
803 #endif
804  }
805 
806  }
807 
808  // Order the QZ decomposition?
809  if (! (ord_job == 'N' || ord_job == 'n'))
810  {
811  if (complex_case)
812  {
813  // Probably not needed, but better be safe.
814  error ("qz: cannot re-order complex qz decomposition");
815  return retval;
816  }
817  else
818  {
819 #ifdef DEBUG_SORT
820  std::cout << "qz: ordering eigenvalues: ord_job = "
821  << ord_job << std::endl;
822 #endif
823 
824  // Declared static to avoid vfork/long jump compiler complaints.
825  static sort_function sort_test;
826  sort_test = 0;
827 
828  switch (ord_job)
829  {
830  case 'S':
831  case 's':
832  sort_test = &fin;
833  break;
834 
835  case 'B':
836  case 'b':
837  sort_test = &fout;
838  break;
839 
840  case '+':
841  sort_test = &fcrhp;
842  break;
843 
844  case '-':
845  sort_test = &folhp;
846  break;
847 
848  default:
849  // Invalid order option (should never happen, since we
850  // checked the options at the top).
851  panic_impossible ();
852  break;
853  }
854 
855  octave_idx_type ndim, fail;
856  double inf_norm;
857 
858  F77_XFCN (xdlange, XDLANGE,
859  (F77_CONST_CHAR_ARG2 ("I", 1),
860  nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm
861  F77_CHAR_ARG_LEN (1)));
862 
863  double eps = std::numeric_limits<double>::epsilon () * inf_norm * nn;
864 
865 #ifdef DEBUG_SORT
866  std::cout << "qz: calling dsubsp: aa=" << std::endl;
867  octave_print_internal (std::cout, aa, 0);
868  std::cout << std::endl << "bb=" << std::endl;
869  octave_print_internal (std::cout, bb, 0);
870  if (compz == 'V')
871  {
872  std::cout << std::endl << "ZZ=" << std::endl;
873  octave_print_internal (std::cout, ZZ, 0);
874  }
875  std::cout << std::endl;
876  std::cout << "alphar = " << std::endl;
877  octave_print_internal (std::cout, (Matrix) alphar, 0);
878  std::cout << std::endl << "alphai = " << std::endl;
879  octave_print_internal (std::cout, (Matrix) alphai, 0);
880  std::cout << std::endl << "beta = " << std::endl;
881  octave_print_internal (std::cout, (Matrix) betar, 0);
882  std::cout << std::endl;
883 #endif
884 
885  Array<octave_idx_type> ind (dim_vector (nn, 1));
886 
887  F77_XFCN (dsubsp, DSUBSP,
888  (nn, nn, aa.fortran_vec (), bb.fortran_vec (),
889  ZZ.fortran_vec (), sort_test, eps, ndim, fail,
890  ind.fortran_vec ()));
891 
892 #ifdef DEBUG
893  std::cout << "qz: back from dsubsp: aa=" << std::endl;
894  octave_print_internal (std::cout, aa, 0);
895  std::cout << std::endl << "bb=" << std::endl;
896  octave_print_internal (std::cout, bb, 0);
897  if (compz == 'V')
898  {
899  std::cout << std::endl << "ZZ=" << std::endl;
900  octave_print_internal (std::cout, ZZ, 0);
901  }
902  std::cout << std::endl;
903 #endif
904 
905  // Manually update alphar, alphai, betar.
906  static int jj;
907 
908  jj = 0;
909  while (jj < nn)
910  {
911 #ifdef DEBUG_EIG
912  std::cout << "computing gen eig #" << jj << std::endl;
913 #endif
914 
915  // Number of zeros in this block.
916  static int zcnt;
917 
918  if (jj == (nn-1))
919  zcnt = 1;
920  else if (aa(jj+1,jj) == 0)
921  zcnt = 1;
922  else zcnt = 2;
923 
924  if (zcnt == 1)
925  {
926  // Real zero.
927 #ifdef DEBUG_EIG
928  std::cout << " single gen eig:" << std::endl;
929  std::cout << " alphar(" << jj << ") = " << aa(jj,jj)
930  << std::endl;
931  std::cout << " betar( " << jj << ") = " << bb(jj,jj)
932  << std::endl;
933  std::cout << " alphai(" << jj << ") = 0" << std::endl;
934 #endif
935 
936  alphar(jj) = aa(jj,jj);
937  alphai(jj) = 0;
938  betar(jj) = bb(jj,jj);
939  }
940  else
941  {
942  // Complex conjugate pair.
943 #ifdef DEBUG_EIG
944  std::cout << "qz: calling dlag2:" << std::endl;
945  std::cout << "safmin="
946  << setiosflags (std::ios::scientific)
947  << safmin << std::endl;
948 
949  for (int idr = jj; idr <= jj+1; idr++)
950  {
951  for (int idc = jj; idc <= jj+1; idc++)
952  {
953  std::cout << "aa(" << idr << "," << idc << ")="
954  << aa(idr,idc) << std::endl;
955  std::cout << "bb(" << idr << "," << idc << ")="
956  << bb(idr,idc) << std::endl;
957  }
958  }
959 #endif
960 
961  // FIXME: probably should be using
962  // fortran_vec instead of &aa(jj,jj) here.
963 
964  double scale1, scale2, wr1, wr2, wi;
965  const double *aa_ptr = aa.data () + jj * nn + jj;
966  const double *bb_ptr = bb.data () + jj * nn + jj;
967  F77_XFCN (dlag2, DLAG2,
968  (aa_ptr, nn, bb_ptr, nn, safmin,
969  scale1, scale2, wr1, wr2, wi));
970 
971 #ifdef DEBUG_EIG
972  std::cout << "dlag2 returns: scale1=" << scale1
973  << "\tscale2=" << scale2 << std::endl
974  << "\twr1=" << wr1 << "\twr2=" << wr2
975  << "\twi=" << wi << std::endl;
976 #endif
977 
978  // Just to be safe, check if it's a real pair.
979  if (wi == 0)
980  {
981  alphar(jj) = wr1;
982  alphai(jj) = 0;
983  betar(jj) = scale1;
984  alphar(jj+1) = wr2;
985  alphai(jj+1) = 0;
986  betar(jj+1) = scale2;
987  }
988  else
989  {
990  alphar(jj) = alphar(jj+1) = wr1;
991  alphai(jj) = -(alphai(jj+1) = wi);
992  betar(jj) = betar(jj+1) = scale1;
993  }
994  }
995 
996  // Advance past this block.
997  jj += zcnt;
998  }
999 
1000 #ifdef DEBUG_SORT
1001  std::cout << "qz: back from dsubsp: aa=" << std::endl;
1002  octave_print_internal (std::cout, aa, 0);
1003  std::cout << std::endl << "bb=" << std::endl;
1004  octave_print_internal (std::cout, bb, 0);
1005 
1006  if (compz == 'V')
1007  {
1008  std::cout << std::endl << "ZZ=" << std::endl;
1009  octave_print_internal (std::cout, ZZ, 0);
1010  }
1011  std::cout << std::endl << "qz: ndim=" << ndim << std::endl
1012  << "fail=" << fail << std::endl;
1013  std::cout << "alphar = " << std::endl;
1014  octave_print_internal (std::cout, (Matrix) alphar, 0);
1015  std::cout << std::endl << "alphai = " << std::endl;
1016  octave_print_internal (std::cout, (Matrix) alphai, 0);
1017  std::cout << std::endl << "beta = " << std::endl;
1018  octave_print_internal (std::cout, (Matrix) betar, 0);
1019  std::cout << std::endl;
1020 #endif
1021  }
1022  }
1023 
1024  // Compute generalized eigenvalues?
1025  ComplexColumnVector gev;
1026 
1027  if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4))
1028  {
1029  if (complex_case)
1030  {
1031  int cnt = 0;
1032 
1033  for (int ii = 0; ii < nn; ii++)
1034  cnt++;
1035 
1036  ComplexColumnVector tmp (cnt);
1037 
1038  cnt = 0;
1039  for (int ii = 0; ii < nn; ii++)
1040  tmp(cnt++) = xalpha(ii) / xbeta(ii);
1041 
1042  gev = tmp;
1043  }
1044  else
1045  {
1046 #ifdef DEBUG
1047  std::cout << "qz: computing generalized eigenvalues" << std::endl;
1048 #endif
1049 
1050  // Return finite generalized eigenvalues.
1051  int cnt = 0;
1052 
1053  for (int ii = 0; ii < nn; ii++)
1054  if (betar(ii) != 0)
1055  cnt++;
1056 
1057  ComplexColumnVector tmp (cnt);
1058 
1059  cnt = 0;
1060  for (int ii = 0; ii < nn; ii++)
1061  if (betar(ii) != 0)
1062  tmp(cnt++) = Complex(alphar(ii), alphai(ii))/betar(ii);
1063 
1064  gev = tmp;
1065  }
1066  }
1067 
1068  // Right, left eigenvector matrices.
1069  if (nargout >= 5)
1070  {
1071  // Which side to compute?
1072  char side = (nargout == 5 ? 'R' : 'B');
1073  // Compute all of them and backtransform
1074  char howmny = 'B';
1075  // Dummy pointer; select is not used.
1076  octave_idx_type *select = 0;
1077 
1078  if (complex_case)
1079  {
1080  CVL = CQ;
1081  CVR = CZ;
1082  ComplexRowVector cwork2 (2 * nn);
1083  RowVector rwork2 (8 * nn);
1084  octave_idx_type m;
1085 
1086  F77_XFCN (ztgevc, ZTGEVC,
1087  (F77_CONST_CHAR_ARG2 (&side, 1),
1088  F77_CONST_CHAR_ARG2 (&howmny, 1),
1089  select, nn, caa.fortran_vec (), nn, cbb.fortran_vec (),
1090  nn, CVL.fortran_vec (), nn, CVR.fortran_vec (), nn, nn,
1091  m, cwork2.fortran_vec (), rwork2.fortran_vec (), info
1092  F77_CHAR_ARG_LEN (1)
1093  F77_CHAR_ARG_LEN (1)));
1094  }
1095  else
1096  {
1097 #ifdef DEBUG
1098  std::cout << "qz: computing generalized eigenvectors" << std::endl;
1099 #endif
1100 
1101  VL = QQ;
1102  VR = ZZ;
1103  octave_idx_type m;
1104 
1105  F77_XFCN (dtgevc, DTGEVC,
1106  (F77_CONST_CHAR_ARG2 (&side, 1),
1107  F77_CONST_CHAR_ARG2 (&howmny, 1),
1108  select, nn, aa.fortran_vec (), nn, bb.fortran_vec (),
1109  nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn,
1110  m, work.fortran_vec (), info
1111  F77_CHAR_ARG_LEN (1)
1112  F77_CHAR_ARG_LEN (1)));
1113 
1114  // Now construct the complex form of VV, WW.
1115  int jj = 0;
1116 
1117  while (jj < nn)
1118  {
1119  OCTAVE_QUIT;
1120 
1121  // See if real or complex eigenvalue.
1122 
1123  // Column increment; assume complex eigenvalue.
1124  int cinc = 2;
1125 
1126  if (jj == (nn-1))
1127  // Single column.
1128  cinc = 1;
1129  else if (aa(jj+1,jj) == 0)
1130  cinc = 1;
1131 
1132  // Now copy the eigenvector (s) to CVR, CVL.
1133  if (cinc == 1)
1134  {
1135  for (int ii = 0; ii < nn; ii++)
1136  CVR(ii,jj) = VR(ii,jj);
1137 
1138  if (side == 'B')
1139  for (int ii = 0; ii < nn; ii++)
1140  CVL(ii,jj) = VL(ii,jj);
1141  }
1142  else
1143  {
1144  // Double column; complex vector.
1145 
1146  for (int ii = 0; ii < nn; ii++)
1147  {
1148  CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1));
1149  CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1));
1150  }
1151 
1152  if (side == 'B')
1153  for (int ii = 0; ii < nn; ii++)
1154  {
1155  CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1));
1156  CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1));
1157  }
1158  }
1159 
1160  // Advance to next eigenvectors (if any).
1161  jj += cinc;
1162  }
1163  }
1164  }
1165 
1166  switch (nargout)
1167  {
1168  case 7:
1169  retval(6) = gev;
1170 
1171  case 6:
1172  // Return eigenvectors.
1173  retval(5) = CVL;
1174 
1175  case 5:
1176  // Return eigenvectors.
1177  retval(4) = CVR;
1178 
1179  case 4:
1180  if (nargin == 3)
1181  {
1182 #ifdef DEBUG
1183  std::cout << "qz: sort: retval(3) = gev = " << std::endl;
1184  octave_print_internal (std::cout, gev);
1185  std::cout << std::endl;
1186 #endif
1187  retval(3) = gev;
1188  }
1189  else
1190  {
1191  if (complex_case)
1192  retval(3) = CZ;
1193  else
1194  retval(3) = ZZ;
1195  }
1196 
1197  case 3:
1198  if (nargin == 3)
1199  {
1200  if (complex_case)
1201  retval(2) = CZ;
1202  else
1203  retval(2) = ZZ;
1204  }
1205  else
1206  {
1207  if (complex_case)
1208  retval(2) = CQ.hermitian ();
1209  else
1210  retval(2) = QQ.transpose ();
1211  }
1212 
1213  case 2:
1214  {
1215  if (complex_case)
1216  {
1217 #ifdef DEBUG
1218  std::cout << "qz: retval(1) = cbb = " << std::endl;
1219  octave_print_internal (std::cout, cbb, 0);
1220  std::cout << std::endl << "qz: retval(0) = caa = " <<std::endl;
1221  octave_print_internal (std::cout, caa, 0);
1222  std::cout << std::endl;
1223 #endif
1224  retval(1) = cbb;
1225  retval(0) = caa;
1226  }
1227  else
1228  {
1229 #ifdef DEBUG
1230  std::cout << "qz: retval(1) = bb = " << std::endl;
1231  octave_print_internal (std::cout, bb, 0);
1232  std::cout << std::endl << "qz: retval(0) = aa = " <<std::endl;
1233  octave_print_internal (std::cout, aa, 0);
1234  std::cout << std::endl;
1235 #endif
1236  retval(1) = bb;
1237  retval(0) = aa;
1238  }
1239  }
1240  break;
1241 
1242 
1243  case 1:
1244  case 0:
1245 #ifdef DEBUG
1246  std::cout << "qz: retval(0) = gev = " << gev << std::endl;
1247 #endif
1248  retval(0) = gev;
1249  break;
1250 
1251  default:
1252  error ("qz: too many return arguments");
1253  break;
1254  }
1255 
1256 #ifdef DEBUG
1257  std::cout << "qz: exiting (at long last)" << std::endl;
1258 #endif
1259 
1260  return retval;
1261 }
1262 
1263 /*
1264 %!shared a, b, c
1265 %! a = [1 2; 0 3];
1266 %! b = [1 0; 0 0];
1267 %! c = [0 1; 0 0];
1268 %!assert (qz (a,b), 1)
1269 %!assert (isempty (qz (a,c)))
1270 
1271 ## Exaple 7.7.3 in Golub & Van Loan
1272 %!test
1273 %! a = [ 10 1 2;
1274 %! 1 2 -1;
1275 %! 1 1 2];
1276 %! b = reshape (1:9,3,3);
1277 %! [aa, bb, q, z, v, w, lambda] = qz (a, b);
1278 %! sz = length (lambda);
1279 %! observed = (b * v * diag ([lambda;0])) (:, 1:sz);
1280 %! assert ( (a*v) (:, 1:sz), observed, norm (observed) * 1e-14);
1281 %! observed = (diag ([lambda;0]) * w' * b) (1:sz, :);
1282 %! assert ( (w'*a) (1:sz, :) , observed, norm (observed) * 1e-13);
1283 %! assert (q * a * z, aa, norm (aa) * 1e-14);
1284 %! assert (q * b * z, bb, norm (bb) * 1e-14);
1285 
1286 %!test
1287 %! A = [0, 0, -1, 0; 1, 0, 0, 0; -1, 0, -2, -1; 0, -1, 1, 0];
1288 %! B = [0, 0, 0, 0; 0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1];
1289 %! [AA, BB, Q, Z1] = qz (A, B);
1290 %! [AA, BB, Z2] = qz (A, B, '-');
1291 %! assert (Z1, Z2);
1292 */
Generated on Mon Dec 30 2013 03:04:29 for GNU Octave by   doxygen 1.8.1.2