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