octave: libinterp/corefcn/qz.cc annotate

annotate libinterp/corefcn/qz.cc @ 23582:0cc2011d800e

maint: Deprecate is_real_type and replace with isreal. * ov.h (is_real_type): Use OCTAVE_DEPRECATED macro around function. * ov.h (isreal): New function. * QtHandlesUtils.cc, __luinc__.cc, bsxfun.cc, cellfun.cc, conv2.cc, data.cc, det.cc, eig.cc, ellipj.cc, fft.cc, fft2.cc, fftn.cc, find.cc, graphics.cc, graphics.in.h, gsvd.cc, hess.cc, inv.cc, lu.cc, mex.cc, oct-stream.cc, pinv.cc, psi.cc, qz.cc, schur.cc, svd.cc, typecast.cc, chol.cc, dmperm.cc, qr.cc, symbfact.cc, symrcm.cc, ov-base-int.h, ov-base.h, ov-bool-mat.h, ov-bool-sparse.h, ov-bool.h, ov-ch-mat.h, ov-float.h, ov-flt-re-diag.h, ov-flt-re-mat.h, ov-java.cc, ov-lazy-idx.h, ov-perm.h, ov-range.h, ov-re-diag.h, ov-re-mat.h, ov-re-sparse.h, ov-scalar.h, pt-tm-const.cc: Replace instances of is_real_type with isreal.

author	Rik <rik@octave.org>
date	Tue, 13 Jun 2017 07:53:53 -0700
parents	c3075ae020e1
children	bd77ab816e43

rev	line source
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	1 /*
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	2
23219 3ac9f9ecfae5 maint: Update copyright dates. John W. Eaton <jwe@octave.org> parents: 23083 diff changeset	3 Copyright (C) 1998-2017 A. S. Hodel
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	4
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	5 This file is part of Octave.
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	6
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	7 Octave is free software; you can redistribute it and/or modify it
22755 3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22407 diff changeset	8 under the terms of the GNU General Public License as published by
3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22407 diff changeset	9 the Free Software Foundation; either version 3 of the License, or
3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22407 diff changeset	10 (at your option) any later version.
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	11
22755 3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22407 diff changeset	12 Octave is distributed in the hope that it will be useful, but
3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22407 diff changeset	13 WITHOUT ANY WARRANTY; without even the implied warranty of
3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22407 diff changeset	14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22407 diff changeset	15 GNU General Public License for more details.
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	16
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	17 You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	18 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	19 <http://www.gnu.org/licenses/>.
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	20
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	21 */
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	22
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	23 // Generalized eigenvalue balancing via LAPACK
3911 8389e78e67d4 [project @ 2002-04-28 02:15:38 by jwe] jwe parents: 3887 diff changeset	24
8389e78e67d4 [project @ 2002-04-28 02:15:38 by jwe] jwe parents: 3887 diff changeset	25 // Author: A. S. Hodel <scotte@eng.auburn.edu>
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	26
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	27 #undef DEBUG
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	28 #undef DEBUG_SORT
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	29 #undef DEBUG_EIG
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	30
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	31 #if defined (HAVE_CONFIG_H)
21301 40de9f8f23a6 Use '#include "config.h"' rather than <config.h>. Rik <rik@octave.org> parents: 21279 diff changeset	32 # include "config.h"
9786 2c279308f6ab fix includes in some src/DLD-FUNCTIONS files John W. Eaton <jwe@octave.org> parents: 9758 diff changeset	33 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	34
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	35 #include <cctype>
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	36 #include <cmath>
4051 b79da8779a0e [project @ 2002-08-17 19:38:32 by jwe] jwe parents: 3911 diff changeset	37
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	38 #if defined (DEBUG) \|\| defined (DEBUG_SORT) \|\| defined (DEBUG_EIG)
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	39 # include <iostream>
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	40 # if defined (DEBUG_EIG)
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	41 # include <iomanip>
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	42 # endif
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	43 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	44
4153 6b96ce9f5743 [project @ 2002-11-06 20:38:49 by jwe] jwe parents: 4051 diff changeset	45 #include "f77-fcn.h"
22322 93b3cdd36854 move most f77 function decls to separate header files John W. Eaton <jwe@octave.org> parents: 22317 diff changeset	46 #include "lo-lapack-proto.h"
7231 2eb392d058bb [project @ 2007-11-30 18:53:29 by jwe] jwe parents: 7017 diff changeset	47 #include "lo-math.h"
21279 eb1524b07fe3 better use of templates for qr classes John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	48 #include "qr.h"
4153 6b96ce9f5743 [project @ 2002-11-06 20:38:49 by jwe] jwe parents: 4051 diff changeset	49 #include "quit.h"
6b96ce9f5743 [project @ 2002-11-06 20:38:49 by jwe] jwe parents: 4051 diff changeset	50
15039 e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14876 diff changeset	51 #include "defun.h"
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	52 #include "error.h"
21100 e39e05d90788 Switch gripe_XXX to either err_XXX or warn_XXX naming scheme. Rik <rik@octave.org> parents: 21024 diff changeset	53 #include "errwarn.h"
20940 48b2ad5ee801 maint: Rename oct-obj.[cc\|h] to ovl.[cc\|h] for clarity. Rik <rik@octave.org> parents: 20892 diff changeset	54 #include "ovl.h"
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	55 #if defined (DEBUG) \|\| defined (DEBUG_SORT)
21200 fcac5dbbf9ed maint: Indent #ifdef blocks in libinterp. Rik <rik@octave.org> parents: 21110 diff changeset	56 # include "pr-output.h"
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	57 #endif
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	58
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	59 typedef F77_INT (*sort_function) (const F77_INT& LSIZE,
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	60 const double& ALPHA, const double& BETA,
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	61 const double& S, const double& P);
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	62
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	63 extern "C"
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	64 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	65 // Van Dooren's code (netlib.org: toms/590) for reordering GEP.
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	66 // Only processes Z, not Q.
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4153 diff changeset	67 F77_RET_T
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	68 F77_FUNC (dsubsp, DSUBSP) (const F77_INT& NMAX, const F77_INT& N,
23449 c763214a8260 maint: Use convention 'int x' for naming pointers. Rik <rik@octave.org>* parents: 23430 diff changeset	69 F77_DBLE A, F77_DBLE B, F77_DBLE *Z,
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	70 sort_function, const F77_DBLE& EPS,
23449 c763214a8260 maint: Use convention 'int x' for naming pointers. Rik <rik@octave.org>* parents: 23430 diff changeset	71 F77_INT& NDIM, F77_INT& FAIL, F77_INT *IND);
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	72 }
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	73
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	74 // fcrhp, fin, fout, folhp:
21562 6c2fd62db1f7 maint: Eliminate accidental double spaces in code. Rik <rik@octave.org> parents: 21301 diff changeset	75 // Routines for ordering of generalized eigenvalues.
6c2fd62db1f7 maint: Eliminate accidental double spaces in code. Rik <rik@octave.org> parents: 21301 diff changeset	76 // Return 1 if test is passed, 0 otherwise.
6c2fd62db1f7 maint: Eliminate accidental double spaces in code. Rik <rik@octave.org> parents: 21301 diff changeset	77 // fin: \|lambda\| < 1
6c2fd62db1f7 maint: Eliminate accidental double spaces in code. Rik <rik@octave.org> parents: 21301 diff changeset	78 // fout: \|lambda\| >= 1
6c2fd62db1f7 maint: Eliminate accidental double spaces in code. Rik <rik@octave.org> parents: 21301 diff changeset	79 // fcrhp: real(lambda) >= 0
6c2fd62db1f7 maint: Eliminate accidental double spaces in code. Rik <rik@octave.org> parents: 21301 diff changeset	80 // folhp: real(lambda) < 0
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	81
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	82 static F77_INT
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	83 fcrhp (const F77_INT& lsize, const double& alpha, const double& beta,
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	84 const double& s, const double&)
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	85 {
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	86 if (lsize == 1)
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	87 return (alpha * beta >= 0 ? 1 : -1);
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	88 else
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	89 return (s >= 0 ? 1 : -1);
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	90 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	91
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	92 static F77_INT
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	93 fin (const F77_INT& lsize, const double& alpha, const double& beta,
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	94 const double&, const double& p)
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	95 {
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	96 F77_INT retval;
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	97
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	98 if (lsize == 1)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	99 retval = (std::abs (alpha) < std::abs (beta) ? 1 : -1);
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	100 else
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	101 retval = (std::abs (p) < 1 ? 1 : -1);
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	102
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	103 #if defined (DEBUG)
3538 0ff7323dab8b [project @ 2000-02-02 12:36:25 by jwe] jwe parents: 3531 diff changeset	104 std::cout << "qz: fin: retval=" << retval << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	105 #endif
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	106
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	107 return retval;
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	108 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	109
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	110 static F77_INT
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	111 folhp (const F77_INT& lsize, const double& alpha, const double& beta,
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	112 const double& s, const double&)
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	113 {
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	114 if (lsize == 1)
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	115 return (alpha * beta < 0 ? 1 : -1);
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	116 else
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	117 return (s < 0 ? 1 : -1);
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	118 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	119
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	120 static F77_INT
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	121 fout (const F77_INT& lsize, const double& alpha, const double& beta,
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	122 const double&, const double& p)
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	123 {
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	124 if (lsize == 1)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	125 return (std::abs (alpha) >= std::abs (beta) ? 1 : -1);
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	126 else
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	127 return (std::abs (p) >= 1 ? 1 : -1);
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	128 }
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	129
13074 0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	130 //FIXME: Matlab does not produce lambda as the first output argument.
0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	131 // Compatibility problem?
15039 e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14876 diff changeset	132 DEFUN (qz, args, nargout,
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	133 doc: /* -- texinfo --
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	134 @deftypefn {} {@var{lambda} =} qz (@var{A}, @var{B})
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	135 @deftypefnx {} {[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}, @var{lambda}] =} qz (@var{A}, @var{B})
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	136 @deftypefnx {} {[@var{AA}, @var{BB}, @var{Z}] =} qz (@var{A}, @var{B}, @var{opt})
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	137 @deftypefnx {} {[@var{AA}, @var{BB}, @var{Z}, @var{lambda}] =} qz (@var{A}, @var{B}, @var{opt})
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	138 QZ@tie{}decomposition of the generalized eigenvalue problem
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	139 @tex
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	140 $$A x = \lambda B x$$
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	141 @end tex
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	142 @ifnottex
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	143
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	144 @math{A x = @var{lambda} B x}
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	145
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	146 @end ifnottex
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	147
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	148 There are three calling forms of the function:
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	149
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	150 @enumerate
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	151 @item @code{@var{lambda} = qz (@var{A}, @var{B})}
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	152
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	153 Compute the generalized eigenvalues
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	154 @tex
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	155 $\lambda.$
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	156 @end tex
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	157 @ifnottex
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	158 @var{lambda}.
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	159 @end ifnottex
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	160
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	161 @item @code{[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}, @var{lambda}] = qz (@var{A}, @var{B})}
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	162
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	163 Compute QZ@tie{}decomposition, generalized eigenvectors, and generalized
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	164 eigenvalues.
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	165 @tex
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	166 $$ AV = BV{ \rm diag }(\lambda) $$
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	167 $$ W^T A = { \rm diag }(\lambda)W^T B $$
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	168 $$ AA = Q^T AZ, BB = Q^T BZ $$
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	169 @end tex
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	170 @ifnottex
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	171
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	172 @example
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	173 @group
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	174
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	175 @var{A} * @var{V} = @var{B} * @var{V} * diag (@var{lambda})
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	176 @var{W}' * @var{A} = diag (@var{lambda}) * @var{W}' * @var{B}
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	177 @var{AA} = @var{Q} * @var{A} * @var{Z}, @var{BB} = @var{Q} * @var{B} * @var{Z}
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	178
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	179 @end group
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	180 @end example
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	181
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	182 @end ifnottex
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	183 with @var{Q} and @var{Z} orthogonal (unitary for complex case).
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	184
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	185 @item @code{[@var{AA}, @var{BB}, @var{Z} @{, @var{lambda}@}] = qz (@var{A}, @var{B}, @var{opt})}
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	186
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	187 As in form 2 above, but allows ordering of generalized eigenpairs for, e.g.,
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	188 solution of discrete time algebraic Riccati equations. Form 3 is not available
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	189 for complex matrices, and does not compute the generalized eigenvectors
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	190 @var{V}, @var{W}, nor the orthogonal matrix @var{Q}.
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	191
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	192 @table @var
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	193 @item opt
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	194 for ordering eigenvalues of the @nospell{GEP} pencil. The leading block of
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	195 the revised pencil contains all eigenvalues that satisfy:
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	196
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	197 @table @asis
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	198 @item @qcode{"N"}
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	199 unordered (default)
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	200
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	201 @item @qcode{"S"}
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	202 small: leading block has all
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	203 @tex
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	204 $\|\lambda\| < 1$
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	205 @end tex
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	206 @ifnottex
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	207 \|@var{lambda}\| < 1
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	208 @end ifnottex
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	209
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	210 @item @qcode{"B"}
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	211 big: leading block has all
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	212 @tex
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	213 $\|\lambda\| \geq 1$
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	214 @end tex
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	215 @ifnottex
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	216 \|@var{lambda}\| @geq{} 1
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	217 @end ifnottex
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	218
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	219 @item @qcode{"-"}
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	220 negative real part: leading block has all eigenvalues in the open left
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	221 half-plane
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	222
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	223 @item @qcode{"+"}
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	224 non-negative real part: leading block has all eigenvalues in the closed right
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	225 half-plane
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	226 @end table
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	227 @end table
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	228 @end enumerate
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	229
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	230 Note: @code{qz} performs permutation balancing, but not scaling
23428 1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	231 (@pxref{XREFbalance}), which may be lead to less accurate results than
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	232 @code{eig}. The order of output arguments was selected for compatibility with
1bc0e610e293 doc: Redo docstring for qz (bug #50846). Rik <rik@octave.org> parents: 23219 diff changeset	233 @sc{matlab}.
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	234 @seealso{eig, balance, lu, chol, hess, qr, qzhess, schur, svd}
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21751 diff changeset	235 @end deftypefn */)
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	236 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	237 int nargin = args.length ();
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	238
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	239 #if defined (DEBUG)
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	240 std::cout << "qz: nargin = " << nargin
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	241 << ", nargout = " << nargout << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	242 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	243
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	244 if (nargin < 2 \|\| nargin > 3 \|\| nargout > 7)
20802 8bb38ba1bad6 eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20700 diff changeset	245 print_usage ();
8bb38ba1bad6 eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20700 diff changeset	246
8bb38ba1bad6 eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20700 diff changeset	247 if (nargin == 3 && (nargout < 3 \|\| nargout > 4))
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	248 error ("qz: invalid number of output arguments for form 3 call");
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	249
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	250 #if defined (DEBUG)
3538 0ff7323dab8b [project @ 2000-02-02 12:36:25 by jwe] jwe parents: 3531 diff changeset	251 std::cout << "qz: determine ordering option" << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	252 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	253
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	254 // Determine ordering option.
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	255 // declared volatile to avoid compiler warnings about long jumps vforks.
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	256 volatile char ord_job;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	257 double safmin;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	258
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	259 if (nargin == 2)
8927 d75f4ee0538d qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC John W. Eaton <jwe@octave.org> parents: 8920 diff changeset	260 ord_job = 'N';
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	261 else
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	262 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	263 std::string opt = args(2).xstring_value ("qz: OPT must be a string");
8927 d75f4ee0538d qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC John W. Eaton <jwe@octave.org> parents: 8920 diff changeset	264
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	265 if (! opt.empty ())
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	266 ord_job = std::toupper (opt[0]);
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	267
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	268 std::string valid_opts = "NSB+-";
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	269
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	270 if (valid_opts.find_first_of (ord_job) == std::string::npos)
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	271 error ("qz: invalid order option '%c'", ord_job);
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	272
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	273 // overflow constant required by dlag2
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4153 diff changeset	274 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1),
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	275 safmin
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	276 F77_CHAR_ARG_LEN (1));
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	277
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	278 #if defined (DEBUG_EIG)
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	279 std::cout << "qz: initial value of safmin="
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	280 << setiosflags (std::ios::scientific)
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	281 << safmin << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	282 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	283
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	284 // Some machines (e.g., DEC alpha) get safmin = 0;
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	285 // for these, use eps instead to avoid problems in dlag2.
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	286 if (safmin == 0)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	287 {
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	288 #if defined (DEBUG_EIG)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	289 std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	290 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	291
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	292 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("E", 1),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	293 safmin
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	294 F77_CHAR_ARG_LEN (1));
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	295
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	296 #if defined (DEBUG_EIG)
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	297 std::cout << "qz: safmin set to "
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	298 << setiosflags (std::ios::scientific)
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	299 << safmin << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	300 #endif
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	301 }
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	302 }
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	303
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	304 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	305 std::cout << "qz: check matrix A" << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	306 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	307
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	308 // Matrix A: check dimensions.
22988 cd33c785e80e put to_f77_int inside octave namespace John W. Eaton <jwe@octave.org> parents: 22964 diff changeset	309 F77_INT nn = octave::to_f77_int (args(0).rows ());
cd33c785e80e put to_f77_int inside octave namespace John W. Eaton <jwe@octave.org> parents: 22964 diff changeset	310 F77_INT nc = octave::to_f77_int (args(0).columns ());
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	311
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	312 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	313 std::cout << "Matrix A dimensions: (" << nn << "," << nc << ")" << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	314 #endif
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	315
23577 80c42f4cca13 maint: Deprecate is_empty and replace with isempty. Rik <rik@octave.org> parents: 23477 diff changeset	316 if (args(0).isempty ())
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	317 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	318 warn_empty_arg ("qz: A");
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	319 return octave_value_list (2, Matrix ());
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	320 }
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	321 else if (nc != nn)
21110 3d0d84305600 Use err_square_matrix_required more widely. Rik <rik@octave.org> parents: 21100 diff changeset	322 err_square_matrix_required ("qz", "A");
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	323
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	324 // Matrix A: get value.
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	325 Matrix aa;
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	326 ComplexMatrix caa;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	327
23581 c3075ae020e1 maint: Deprecate is_complex_type and replace with iscomplex. Rik <rik@octave.org> parents: 23577 diff changeset	328 if (args(0).iscomplex ())
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	329 caa = args(0).complex_matrix_value ();
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	330 else
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	331 aa = args(0).matrix_value ();
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	332
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	333 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	334 std::cout << "qz: check matrix B" << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	335 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	336
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	337 // Extract argument 2 (bb, or cbb if complex).
22988 cd33c785e80e put to_f77_int inside octave namespace John W. Eaton <jwe@octave.org> parents: 22964 diff changeset	338 F77_INT b_nr = octave::to_f77_int (args(1).rows ());
cd33c785e80e put to_f77_int inside octave namespace John W. Eaton <jwe@octave.org> parents: 22964 diff changeset	339 F77_INT b_nc = octave::to_f77_int (args(1).columns ());
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	340
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	341 if (nn != b_nc \|\| nn != b_nr)
21100 e39e05d90788 Switch gripe_XXX to either err_XXX or warn_XXX naming scheme. Rik <rik@octave.org> parents: 21024 diff changeset	342 err_nonconformant ();
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	343
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	344 Matrix bb;
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	345 ComplexMatrix cbb;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	346
23581 c3075ae020e1 maint: Deprecate is_complex_type and replace with iscomplex. Rik <rik@octave.org> parents: 23577 diff changeset	347 if (args(1).iscomplex ())
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	348 cbb = args(1).complex_matrix_value ();
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	349 else
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	350 bb = args(1).matrix_value ();
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	351
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	352 // Both matrices loaded, now check whether to calculate complex or real val.
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	353
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	354 // declared volatile to avoid compiler warnings about long jumps vforks.
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	355 volatile bool complex_case
23581 c3075ae020e1 maint: Deprecate is_complex_type and replace with iscomplex. Rik <rik@octave.org> parents: 23577 diff changeset	356 = (args(0).iscomplex () \|\| args(1).iscomplex ());
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	357
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	358 if (nargin == 3 && complex_case)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	359 error ("qz: cannot re-order complex qz decomposition");
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	360
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	361 // First, declare variables used in both the real and complex cases.
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	362 // FIXME: There are a lot of excess variables declared.
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	363 // Probably a better way to handle this.
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	364 Matrix QQ (nn,nn), ZZ (nn,nn), VR (nn,nn), VL (nn,nn);
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	365 RowVector alphar (nn), alphai (nn), betar (nn);
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	366 ComplexRowVector xalpha (nn), xbeta (nn);
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	367 ComplexMatrix CQ (nn,nn), CZ (nn,nn), CVR (nn,nn), CVL (nn,nn);
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	368 F77_INT ilo, ihi, info;
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	369 char comp_q = (nargout >= 3 ? 'V' : 'N');
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	370 char comp_z = ((nargout >= 4 \|\| nargin == 3)? 'V' : 'N');
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	371
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	372 // Initialize Q, Z to identity matrix if either is needed
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	373 if (comp_q == 'V' \|\| comp_z == 'V')
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	374 {
23477 3530b956d707 maint: Strip trailing whitespace from source files. John W. Eaton <jwe@octave.org> parents: 23457 diff changeset	375 double *QQptr = QQ.fortran_vec ();
3530b956d707 maint: Strip trailing whitespace from source files. John W. Eaton <jwe@octave.org> parents: 23457 diff changeset	376 double *ZZptr = ZZ.fortran_vec ();
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	377 std::fill_n (QQptr, QQ.numel (), 0.0);
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	378 std::fill_n (ZZptr, ZZ.numel (), 0.0);
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	379 for (F77_INT i = 0; i < nn; i++)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	380 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	381 QQ(i,i) = 1.0;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	382 ZZ(i,i) = 1.0;
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	383 }
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	384 }
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	385
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	386 // Always perform permutation balancing.
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4153 diff changeset	387 const char bal_job = 'P';
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	388 RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn);
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	389
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	390 if (complex_case)
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	391 {
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	392 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	393 if (comp_q == 'V')
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	394 std::cout << "qz: performing balancing; CQ=" << std::endl
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	395 << CQ << std::endl;
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	396 #endif
23582 0cc2011d800e maint: Deprecate is_real_type and replace with isreal. Rik <rik@octave.org> parents: 23581 diff changeset	397 if (args(0).isreal ())
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	398 caa = ComplexMatrix (aa);
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	399
23582 0cc2011d800e maint: Deprecate is_real_type and replace with isreal. Rik <rik@octave.org> parents: 23581 diff changeset	400 if (args(1).isreal ())
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	401 cbb = ComplexMatrix (bb);
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	402
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	403 if (comp_q == 'V')
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	404 CQ = ComplexMatrix (QQ);
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	405
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	406 if (comp_z == 'V')
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	407 CZ = ComplexMatrix (ZZ);
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	408
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	409 F77_XFCN (zggbal, ZGGBAL,
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	410 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
22407 34ce5be04942 maint: Style check C++ code in libinterp/. Rik <rik@octave.org> parents: 22372 diff changeset	411 nn, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
34ce5be04942 maint: Style check C++ code in libinterp/. Rik <rik@octave.org> parents: 22372 diff changeset	412 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()),
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	413 nn, ilo, ihi, lscale.fortran_vec (),
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	414 rscale.fortran_vec (), work.fortran_vec (), info
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	415 F77_CHAR_ARG_LEN (1)));
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	416 }
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	417 else
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	418 {
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	419 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	420 if (comp_q == 'V')
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	421 std::cout << "qz: performing balancing; QQ=" << std::endl
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	422 << QQ << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	423 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	424
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	425 F77_XFCN (dggbal, DGGBAL,
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	426 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	427 nn, aa.fortran_vec (), nn, bb.fortran_vec (),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	428 nn, ilo, ihi, lscale.fortran_vec (),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	429 rscale.fortran_vec (), work.fortran_vec (), info
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	430 F77_CHAR_ARG_LEN (1)));
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	431 }
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	432
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	433 // Only permutation balance above is done. Skip scaling balance.
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	434
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	435 #if 0
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	436 // Since we just want the balancing matrices, we can use dggbal
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	437 // for both the real and complex cases; left first
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	438
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	439 if (comp_q == 'V')
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	440 {
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	441 F77_XFCN (dggbak, DGGBAK,
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	442 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	443 F77_CONST_CHAR_ARG2 ("L", 1),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	444 nn, ilo, ihi, lscale.data (), rscale.data (),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	445 nn, QQ.fortran_vec (), nn, info
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	446 F77_CHAR_ARG_LEN (1)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	447 F77_CHAR_ARG_LEN (1)));
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	448
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	449 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	450 if (comp_q == 'V')
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	451 std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	452 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	453 }
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	454
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	455 // then right
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	456 if (comp_z == 'V')
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	457 {
4552 6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe] jwe parents: 4153 diff changeset	458 F77_XFCN (dggbak, DGGBAK,
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	459 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	460 F77_CONST_CHAR_ARG2 ("R", 1),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	461 nn, ilo, ihi, lscale.data (), rscale.data (),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	462 nn, ZZ.fortran_vec (), nn, info
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	463 F77_CHAR_ARG_LEN (1)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	464 F77_CHAR_ARG_LEN (1)));
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	465
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	466 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	467 if (comp_z == 'V')
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	468 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	469 #endif
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	470 }
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	471 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	472
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	473 char qz_job = (nargout < 2 ? 'E' : 'S');
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	474
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	475 if (complex_case)
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	476 {
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	477 // Complex case.
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	478
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	479 // The QR decomposition of cbb.
22317 6ca3acf5fad8 move some new numeric classes to namespace octave::math John W. Eaton <jwe@octave.org> parents: 22137 diff changeset	480 octave::math::qr<ComplexMatrix> cbqr (cbb);
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	481 // The R matrix of QR decomposition for cbb.
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	482 cbb = cbqr.R ();
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	483 // (Q*)caa for following work.
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	484 caa = (cbqr.Q ().hermitian ()) * caa;
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	485 CQ = CQ * cbqr.Q ();
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	486
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	487 F77_XFCN (zgghrd, ZGGHRD,
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	488 (F77_CONST_CHAR_ARG2 (&comp_q, 1),
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	489 F77_CONST_CHAR_ARG2 (&comp_z, 1),
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	490 nn, ilo, ihi, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()),
22407 34ce5be04942 maint: Style check C++ code in libinterp/. Rik <rik@octave.org> parents: 22372 diff changeset	491 nn, F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), nn,
34ce5be04942 maint: Style check C++ code in libinterp/. Rik <rik@octave.org> parents: 22372 diff changeset	492 F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn,
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	493 F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn, info
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	494 F77_CHAR_ARG_LEN (1)
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	495 F77_CHAR_ARG_LEN (1)));
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	496
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	497 ComplexRowVector cwork (nn);
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	498
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	499 F77_XFCN (zhgeqz, ZHGEQZ,
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	500 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	501 F77_CONST_CHAR_ARG2 (&comp_q, 1),
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	502 F77_CONST_CHAR_ARG2 (&comp_z, 1),
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	503 nn, ilo, ihi,
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	504 F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	505 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), nn,
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	506 F77_DBLE_CMPLX_ARG (xalpha.fortran_vec ()),
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	507 F77_DBLE_CMPLX_ARG (xbeta.fortran_vec ()),
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	508 F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn,
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	509 F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn,
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	510 F77_DBLE_CMPLX_ARG (cwork.fortran_vec ()), nn,
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	511 rwork.fortran_vec (), info
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	512 F77_CHAR_ARG_LEN (1)
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	513 F77_CHAR_ARG_LEN (1)
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	514 F77_CHAR_ARG_LEN (1)));
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	515
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	516 if (comp_q == 'V')
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	517 {
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	518 // Left eigenvector.
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	519 F77_XFCN (zggbak, ZGGBAK,
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	520 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	521 F77_CONST_CHAR_ARG2 ("L", 1),
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	522 nn, ilo, ihi, lscale.data (), rscale.data (),
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	523 nn, F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn, info
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	524 F77_CHAR_ARG_LEN (1)
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	525 F77_CHAR_ARG_LEN (1)));
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	526 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	527
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	528 if (comp_z == 'V')
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	529 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	530 // Right eigenvector.
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	531 F77_XFCN (zggbak, ZGGBAK,
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	532 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	533 F77_CONST_CHAR_ARG2 ("R", 1),
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	534 nn, ilo, ihi, lscale.data (), rscale.data (),
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	535 nn, F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn, info
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	536 F77_CHAR_ARG_LEN (1)
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	537 F77_CHAR_ARG_LEN (1)));
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	538 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	539
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	540 }
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	541 else
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	542 {
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	543 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	544 std::cout << "qz: performing qr decomposition of bb" << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	545 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	546
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	547 // Compute the QR factorization of bb.
22317 6ca3acf5fad8 move some new numeric classes to namespace octave::math John W. Eaton <jwe@octave.org> parents: 22137 diff changeset	548 octave::math::qr<Matrix> bqr (bb);
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	549
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	550 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	551 std::cout << "qz: qr (bb) done; now performing qz decomposition"
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	552 << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	553 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	554
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	555 bb = bqr.R ();
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	556
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	557 #if defined (DEBUG)
3538 0ff7323dab8b [project @ 2000-02-02 12:36:25 by jwe] jwe parents: 3531 diff changeset	558 std::cout << "qz: extracted bb" << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	559 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	560
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	561 aa = (bqr.Q ()).transpose () * aa;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	562
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	563 #if defined (DEBUG)
3538 0ff7323dab8b [project @ 2000-02-02 12:36:25 by jwe] jwe parents: 3531 diff changeset	564 std::cout << "qz: updated aa " << std::endl;
0ff7323dab8b [project @ 2000-02-02 12:36:25 by jwe] jwe parents: 3531 diff changeset	565 std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl;
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	566
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	567 if (comp_q == 'V')
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	568 std::cout << "QQ =" << QQ << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	569 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	570
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	571 if (comp_q == 'V')
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	572 QQ = QQ * bqr.Q ();
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	573
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	574 #if defined (DEBUG)
3538 0ff7323dab8b [project @ 2000-02-02 12:36:25 by jwe] jwe parents: 3531 diff changeset	575 std::cout << "qz: precursors done..." << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	576 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	577
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	578 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	579 std::cout << "qz: comp_q = " << comp_q << ", comp_z = " << comp_z
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	580 << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	581 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	582
21562 6c2fd62db1f7 maint: Eliminate accidental double spaces in code. Rik <rik@octave.org> parents: 21301 diff changeset	583 // Reduce to generalized Hessenberg form.
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	584 F77_XFCN (dgghrd, DGGHRD,
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	585 (F77_CONST_CHAR_ARG2 (&comp_q, 1),
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	586 F77_CONST_CHAR_ARG2 (&comp_z, 1),
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	587 nn, ilo, ihi, aa.fortran_vec (),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	588 nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn,
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	589 ZZ.fortran_vec (), nn, info
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	590 F77_CHAR_ARG_LEN (1)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	591 F77_CHAR_ARG_LEN (1)));
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	592
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	593 // Check if just computing generalized eigenvalues,
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	594 // or if we're actually computing the decomposition.
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	595
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	596 // Reduce to generalized Schur form.
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	597 F77_XFCN (dhgeqz, DHGEQZ,
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	598 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	599 F77_CONST_CHAR_ARG2 (&comp_q, 1),
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	600 F77_CONST_CHAR_ARG2 (&comp_z, 1),
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	601 nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	602 nn, alphar.fortran_vec (), alphai.fortran_vec (),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	603 betar.fortran_vec (), QQ.fortran_vec (), nn,
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	604 ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	605 F77_CHAR_ARG_LEN (1)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	606 F77_CHAR_ARG_LEN (1)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	607 F77_CHAR_ARG_LEN (1)));
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	608
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	609 if (comp_q == 'V')
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	610 {
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	611 F77_XFCN (dggbak, DGGBAK,
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	612 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	613 F77_CONST_CHAR_ARG2 ("L", 1),
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	614 nn, ilo, ihi, lscale.data (), rscale.data (),
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	615 nn, QQ.fortran_vec (), nn, info
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	616 F77_CHAR_ARG_LEN (1)
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	617 F77_CHAR_ARG_LEN (1)));
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	618
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	619 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	620 if (comp_q == 'V')
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	621 std::cout << "qz: balancing done; QQ=" << std::endl
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	622 << QQ << std::endl;
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	623 #endif
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	624 }
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	625
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	626 // then right
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	627 if (comp_z == 'V')
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	628 {
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	629 F77_XFCN (dggbak, DGGBAK,
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	630 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	631 F77_CONST_CHAR_ARG2 ("R", 1),
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	632 nn, ilo, ihi, lscale.data (), rscale.data (),
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	633 nn, ZZ.fortran_vec (), nn, info
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	634 F77_CHAR_ARG_LEN (1)
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	635 F77_CHAR_ARG_LEN (1)));
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	636
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	637 #if defined (DEBUG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	638 if (comp_z == 'V')
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	639 std::cout << "qz: balancing done; ZZ=" << std::endl
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	640 << ZZ << std::endl;
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	641 #endif
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	642 }
4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	643
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	644 }
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	645
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	646 // Order the QZ decomposition?
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	647 if (ord_job != 'N')
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	648 {
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	649 if (complex_case)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	650 // Probably not needed, but better be safe.
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	651 error ("qz: cannot re-order complex QZ decomposition");
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	652
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	653 #if defined (DEBUG_SORT)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	654 std::cout << "qz: ordering eigenvalues: ord_job = "
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	655 << ord_job << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	656 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	657
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	658 // Declared static to avoid vfork/long jump compiler complaints.
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	659 static sort_function sort_test;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	660 sort_test = 0;
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	661
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	662 switch (ord_job)
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	663 {
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	664 case 'S':
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	665 sort_test = &fin;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	666 break;
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	667
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	668 case 'B':
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	669 sort_test = &fout;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	670 break;
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	671
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	672 case '+':
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	673 sort_test = &fcrhp;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	674 break;
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	675
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	676 case '-':
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	677 sort_test = &folhp;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	678 break;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	679
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	680 default:
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	681 // Invalid order option
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	682 // (should never happen since options were checked at the top).
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	683 panic_impossible ();
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	684 break;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	685 }
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	686
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	687 double inf_norm;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	688
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	689 F77_XFCN (xdlange, XDLANGE,
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	690 (F77_CONST_CHAR_ARG2 ("I", 1),
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	691 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	692 F77_CHAR_ARG_LEN (1)));
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	693
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	694 double eps = std::numeric_limits<double>::epsilon () * inf_norm * nn;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	695
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	696 #if defined (DEBUG_SORT)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	697 std::cout << "qz: calling dsubsp: aa=" << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	698 octave_print_internal (std::cout, aa, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	699 std::cout << std::endl << "bb=" << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	700 octave_print_internal (std::cout, bb, 0);
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	701 if (comp_z == 'V')
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	702 {
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	703 std::cout << std::endl << "ZZ=" << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	704 octave_print_internal (std::cout, ZZ, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	705 }
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	706 std::cout << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	707 std::cout << "alphar = " << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	708 octave_print_internal (std::cout, (Matrix) alphar, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	709 std::cout << std::endl << "alphai = " << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	710 octave_print_internal (std::cout, (Matrix) alphai, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	711 std::cout << std::endl << "beta = " << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	712 octave_print_internal (std::cout, (Matrix) betar, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	713 std::cout << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	714 #endif
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	715
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	716 Array<F77_INT> ind (dim_vector (nn, 1));
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	717
0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	718 F77_INT ndim, fail;
3550 bc492f4a94cb [project @ 2000-02-03 01:36:30 by jwe] jwe parents: 3538 diff changeset	719
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	720 F77_XFCN (dsubsp, DSUBSP,
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	721 (nn, nn, aa.fortran_vec (), bb.fortran_vec (),
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	722 ZZ.fortran_vec (), sort_test, eps, ndim, fail,
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	723 ind.fortran_vec ()));
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	724
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	725 #if defined (DEBUG)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	726 std::cout << "qz: back from dsubsp: aa=" << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	727 octave_print_internal (std::cout, aa, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	728 std::cout << std::endl << "bb=" << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	729 octave_print_internal (std::cout, bb, 0);
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	730 if (comp_z == 'V')
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	731 {
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	732 std::cout << std::endl << "ZZ=" << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	733 octave_print_internal (std::cout, ZZ, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	734 }
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	735 std::cout << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	736 #endif
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	737
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	738 // Manually update alphar, alphai, betar.
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	739 static F77_INT j;
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	740
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	741 j = 0;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	742 while (j < nn)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	743 {
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	744 #if defined (DEBUG_EIG)
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	745 std::cout << "computing gen eig #" << j << std::endl;
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	746 #endif
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	747
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	748 // Number of zeros in this block.
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	749 static F77_INT zcnt;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	750
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	751 if (j == (nn-1))
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	752 zcnt = 1;
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	753 else if (aa(j+1,j) == 0)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	754 zcnt = 1;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	755 else zcnt = 2;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	756
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	757 if (zcnt == 1)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	758 {
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	759 // Real zero.
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	760 #if defined (DEBUG_EIG)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	761 std::cout << " single gen eig:" << std::endl;
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	762 std::cout << " alphar(" << j << ") = " << aa(j,j)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	763 << std::endl;
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	764 std::cout << " betar(" << j << ") = " << bb(j,j)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	765 << std::endl;
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	766 std::cout << " alphai(" << j << ") = 0" << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	767 #endif
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	768
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	769 alphar(j) = aa(j,j);
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	770 alphai(j) = 0;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	771 betar(j) = bb(j,j);
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	772 }
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	773 else
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	774 {
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	775 // Complex conjugate pair.
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	776 #if defined (DEBUG_EIG)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	777 std::cout << "qz: calling dlag2:" << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	778 std::cout << "safmin="
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	779 << setiosflags (std::ios::scientific)
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	780 << safmin << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	781
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	782 for (F77_INT idr = j; idr <= j+1; idr++)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	783 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	784 for (F77_INT idc = j; idc <= j+1; idc++)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	785 {
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	786 std::cout << "aa(" << idr << "," << idc << ")="
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	787 << aa(idr,idc) << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	788 std::cout << "bb(" << idr << "," << idc << ")="
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	789 << bb(idr,idc) << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	790 }
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	791 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	792 #endif
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	793
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	794 // FIXME: probably should be using
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	795 // fortran_vec instead of &aa(j,j) here.
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	796
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	797 double scale1, scale2, wr1, wr2, wi;
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	798 const double aa_ptr = aa.data () + j nn + j;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	799 const double bb_ptr = bb.data () + j nn + j;
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	800 F77_XFCN (dlag2, DLAG2,
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	801 (aa_ptr, nn, bb_ptr, nn, safmin,
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	802 scale1, scale2, wr1, wr2, wi));
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	803
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	804 #if defined (DEBUG_EIG)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	805 std::cout << "dlag2 returns: scale1=" << scale1
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	806 << "\tscale2=" << scale2 << std::endl
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	807 << "\twr1=" << wr1 << "\twr2=" << wr2
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	808 << "\twi=" << wi << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	809 #endif
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	810
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	811 // Just to be safe, check if it's a real pair.
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	812 if (wi == 0)
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	813 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	814 alphar(j) = wr1;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	815 alphai(j) = 0;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	816 betar(j) = scale1;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	817 alphar(j+1) = wr2;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	818 alphai(j+1) = 0;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	819 betar(j+1) = scale2;
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	820 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	821 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	822 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	823 alphar(j) = alphar(j+1) = wr1;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	824 alphai(j) = -(alphai(j+1) = wi);
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	825 betar(j) = betar(j+1) = scale1;
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	826 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	827 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	828
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	829 // Advance past this block.
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	830 j += zcnt;
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	831 }
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	832
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	833 #if defined (DEBUG_SORT)
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	834 std::cout << "qz: back from dsubsp: aa=" << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	835 octave_print_internal (std::cout, aa, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	836 std::cout << std::endl << "bb=" << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	837 octave_print_internal (std::cout, bb, 0);
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	838
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	839 if (comp_z == 'V')
20831 35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	840 {
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	841 std::cout << std::endl << "ZZ=" << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	842 octave_print_internal (std::cout, ZZ, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	843 }
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	844 std::cout << std::endl << "qz: ndim=" << ndim << std::endl
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	845 << "fail=" << fail << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	846 std::cout << "alphar = " << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	847 octave_print_internal (std::cout, (Matrix) alphar, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	848 std::cout << std::endl << "alphai = " << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	849 octave_print_internal (std::cout, (Matrix) alphai, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	850 std::cout << std::endl << "beta = " << std::endl;
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	851 octave_print_internal (std::cout, (Matrix) betar, 0);
35241c4b696c eliminate return statements after calls to error John W. Eaton <jwe@octave.org> parents: 20802 diff changeset	852 std::cout << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	853 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	854 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	855
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	856 // Compute the generalized eigenvalues as well?
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	857 ComplexColumnVector gev;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	858
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	859 if (nargout < 2 \|\| nargout == 7 \|\| (nargin == 3 && nargout == 4))
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	860 {
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	861 if (complex_case)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	862 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	863 ComplexColumnVector tmp (nn);
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	864
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	865 for (F77_INT i = 0; i < nn; i++)
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	866 tmp(i) = xalpha(i) / xbeta(i);
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	867
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	868 gev = tmp;
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	869 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	870 else
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	871 {
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	872 #if defined (DEBUG)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	873 std::cout << "qz: computing generalized eigenvalues" << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	874 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	875
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	876 // Return finite generalized eigenvalues.
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	877 ComplexColumnVector tmp (nn);
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	878
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	879 F77_INT cnt = 0;
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	880 for (F77_INT i = 0; i < nn; i++)
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	881 if (betar(i) != 0)
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	882 tmp(cnt++) = Complex (alphar(i), alphai(i)) / betar(i);
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	883
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	884 tmp.resize (cnt); // Trim vector to number of return values
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	885
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	886 gev = tmp;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	887 }
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	888 }
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	889
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	890 // Right, left eigenvector matrices.
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	891 if (nargout >= 5)
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	892 {
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	893 // Which side to compute?
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	894 char side = (nargout == 5 ? 'R' : 'B');
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	895 // Compute all of them and backtransform
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	896 char howmany = 'B';
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	897 // Dummy pointer; select is not used.
23457 21baad6b35c4 maint: Use C++11 nullptr rather than 0 or NULL when possible. Rik <rik@octave.org> parents: 23449 diff changeset	898 F77_INT *select = nullptr;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	899
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	900 if (complex_case)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	901 {
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	902 CVL = CQ;
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	903 CVR = CZ;
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	904 ComplexRowVector cwork2 (2 * nn);
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	905 RowVector rwork2 (8 * nn);
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	906 F77_INT m;
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	907
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	908 F77_XFCN (ztgevc, ZTGEVC,
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	909 (F77_CONST_CHAR_ARG2 (&side, 1),
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	910 F77_CONST_CHAR_ARG2 (&howmany, 1),
22407 34ce5be04942 maint: Style check C++ code in libinterp/. Rik <rik@octave.org> parents: 22372 diff changeset	911 select, nn, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
34ce5be04942 maint: Style check C++ code in libinterp/. Rik <rik@octave.org> parents: 22372 diff changeset	912 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()),
34ce5be04942 maint: Style check C++ code in libinterp/. Rik <rik@octave.org> parents: 22372 diff changeset	913 nn, F77_DBLE_CMPLX_ARG (CVL.fortran_vec ()), nn,
34ce5be04942 maint: Style check C++ code in libinterp/. Rik <rik@octave.org> parents: 22372 diff changeset	914 F77_DBLE_CMPLX_ARG (CVR.fortran_vec ()), nn, nn,
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	915 m, F77_DBLE_CMPLX_ARG (cwork2.fortran_vec ()),
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	916 rwork2.fortran_vec (), info
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	917 F77_CHAR_ARG_LEN (1)
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	918 F77_CHAR_ARG_LEN (1)));
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	919 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	920 else
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	921 {
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	922 #if defined (DEBUG)
21562 6c2fd62db1f7 maint: Eliminate accidental double spaces in code. Rik <rik@octave.org> parents: 21301 diff changeset	923 std::cout << "qz: computing generalized eigenvectors" << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	924 #endif
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	925
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	926 VL = QQ;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	927 VR = ZZ;
22964 0c12642be005 use F77_INT instead of octave_idx_type for libinterp QZ function John W. Eaton <jwe@octave.org> parents: 22860 diff changeset	928 F77_INT m;
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	929
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	930 F77_XFCN (dtgevc, DTGEVC,
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	931 (F77_CONST_CHAR_ARG2 (&side, 1),
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	932 F77_CONST_CHAR_ARG2 (&howmany, 1),
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	933 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (),
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	934 nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn,
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	935 m, work.fortran_vec (), info
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	936 F77_CHAR_ARG_LEN (1)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	937 F77_CHAR_ARG_LEN (1)));
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	938
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	939 // Now construct the complex form of VV, WW.
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	940 F77_INT j = 0;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	941
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	942 while (j < nn)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	943 {
22860 0b1e25cc4457 eliminate use of OCTAVE_QUIT macro in C++ sources John W. Eaton <jwe@octave.org> parents: 22755 diff changeset	944 octave_quit ();
4153 6b96ce9f5743 [project @ 2002-11-06 20:38:49 by jwe] jwe parents: 4051 diff changeset	945
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	946 // See if real or complex eigenvalue.
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	947
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	948 // Column increment; assume complex eigenvalue.
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	949 int cinc = 2;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	950
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	951 if (j == (nn-1))
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	952 // Single column.
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	953 cinc = 1;
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	954 else if (aa(j+1,j) == 0)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	955 cinc = 1;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	956
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	957 // Now copy the eigenvector (s) to CVR, CVL.
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	958 if (cinc == 1)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	959 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	960 for (F77_INT i = 0; i < nn; i++)
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	961 CVR(i,j) = VR(i,j);
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	962
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	963 if (side == 'B')
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	964 for (F77_INT i = 0; i < nn; i++)
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	965 CVL(i,j) = VL(i,j);
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	966 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	967 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	968 {
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	969 // Double column; complex vector.
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	970
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	971 for (F77_INT i = 0; i < nn; i++)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	972 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	973 CVR(i,j) = Complex (VR(i,j), VR(i,j+1));
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	974 CVR(i,j+1) = Complex (VR(i,j), -VR(i,j+1));
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	975 }
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	976
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	977 if (side == 'B')
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	978 for (F77_INT i = 0; i < nn; i++)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	979 {
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	980 CVL(i,j) = Complex (VL(i,j), VL(i,j+1));
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	981 CVL(i,j+1) = Complex (VL(i,j), -VL(i,j+1));
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	982 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	983 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	984
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	985 // Advance to next eigenvectors (if any).
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	986 j += cinc;
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	987 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	988 }
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	989 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	990
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	991 octave_value_list retval (nargout);
32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	992
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	993 switch (nargout)
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	994 {
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	995 case 7:
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	996 retval(6) = gev;
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	997
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	998 case 6:
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	999 // Return left eigenvectors.
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1000 retval(5) = CVL;
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1001
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1002 case 5:
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	1003 // Return right eigenvectors.
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1004 retval(4) = CVR;
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1005
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1006 case 4:
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1007 if (nargin == 3)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	1008 {
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	1009 #if defined (DEBUG)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	1010 std::cout << "qz: sort: retval(3) = gev = " << std::endl;
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	1011 octave_print_internal (std::cout, ComplexMatrix (gev));
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	1012 std::cout << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1013 #endif
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	1014 retval(3) = gev;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9786 diff changeset	1015 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1016 else
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1017 {
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1018 if (complex_case)
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1019 retval(3) = CZ;
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1020 else
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1021 retval(3) = ZZ;
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1022 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1023
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1024 case 3:
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1025 if (nargin == 3)
14876 54a386f2ac4e correcly compute Z for 3-output call to qz (bug #36728) John W. Eaton <jwe@octave.org> parents: 14844 diff changeset	1026 {
54a386f2ac4e correcly compute Z for 3-output call to qz (bug #36728) John W. Eaton <jwe@octave.org> parents: 14844 diff changeset	1027 if (complex_case)
54a386f2ac4e correcly compute Z for 3-output call to qz (bug #36728) John W. Eaton <jwe@octave.org> parents: 14844 diff changeset	1028 retval(2) = CZ;
54a386f2ac4e correcly compute Z for 3-output call to qz (bug #36728) John W. Eaton <jwe@octave.org> parents: 14844 diff changeset	1029 else
54a386f2ac4e correcly compute Z for 3-output call to qz (bug #36728) John W. Eaton <jwe@octave.org> parents: 14844 diff changeset	1030 retval(2) = ZZ;
54a386f2ac4e correcly compute Z for 3-output call to qz (bug #36728) John W. Eaton <jwe@octave.org> parents: 14844 diff changeset	1031 }
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1032 else
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1033 {
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1034 if (complex_case)
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1035 retval(2) = CQ.hermitian ();
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1036 else
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1037 retval(2) = QQ.transpose ();
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1038 }
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1039
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1040 case 2:
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	1041 {
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1042 if (complex_case)
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1043 {
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	1044 #if defined (DEBUG)
14844 5bc9b9cb4362 maint: Use Octave coding conventions for cuddled parenthesis in retval assignments. Rik <octave@nomad.inbox5.com> parents: 14501 diff changeset	1045 std::cout << "qz: retval(1) = cbb = " << std::endl;
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1046 octave_print_internal (std::cout, cbb, 0);
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1047 std::cout << std::endl << "qz: retval(0) = caa = " <<std::endl;
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1048 octave_print_internal (std::cout, caa, 0);
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1049 std::cout << std::endl;
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	1050 #endif
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1051 retval(1) = cbb;
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1052 retval(0) = caa;
a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1053 }
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	1054 else
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	1055 {
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	1056 #if defined (DEBUG)
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	1057 std::cout << "qz: retval(1) = bb = " << std::endl;
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	1058 octave_print_internal (std::cout, bb, 0);
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	1059 std::cout << std::endl << "qz: retval(0) = aa = " <<std::endl;
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	1060 octave_print_internal (std::cout, aa, 0);
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	1061 std::cout << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1062 #endif
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	1063 retval(1) = bb;
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	1064 retval(0) = aa;
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	1065 }
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	1066 }
10311 a217e1d74353 untabify qz.cc John W. Eaton <jwe@octave.org> parents: 10307 diff changeset	1067 break;
10307 4e4270ab70d6 support complex case in qz Jyh-miin Lin <jyhmiinlin@gmail.com> parents: 10155 diff changeset	1068
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1069 case 1:
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1070 case 0:
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	1071 #if defined (DEBUG)
3538 0ff7323dab8b [project @ 2000-02-02 12:36:25 by jwe] jwe parents: 3531 diff changeset	1072 std::cout << "qz: retval(0) = gev = " << gev << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1073 #endif
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1074 retval(0) = gev;
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1075 break;
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1076
9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1077 default:
11553 01f703952eff Improve docstrings for functions in DLD-FUNCTIONS directory. Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	1078 error ("qz: too many return arguments");
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1079 break;
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	1080 }
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	1081
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21562 diff changeset	1082 #if defined (DEBUG)
3538 0ff7323dab8b [project @ 2000-02-02 12:36:25 by jwe] jwe parents: 3531 diff changeset	1083 std::cout << "qz: exiting (at long last)" << std::endl;
3185 9580887dd160 [project @ 1998-09-26 02:45:55 by jwe] jwe parents: 3183 diff changeset	1084 #endif
3183 5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	1085
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	1086 return retval;
5edc539c2f80 [project @ 1998-09-25 03:37:22 by jwe] jwe parents: diff changeset	1087 }
13074 0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	1088
0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	1089 /*
0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	1090 %!shared a, b, c
14501 60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14360 diff changeset	1091 %! a = [1 2; 0 3];
60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14360 diff changeset	1092 %! b = [1 0; 0 0];
60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14360 diff changeset	1093 %! c = [0 1; 0 0];
60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14360 diff changeset	1094 %!assert (qz (a,b), 1)
60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14360 diff changeset	1095 %!assert (isempty (qz (a,c)))
13074 0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	1096
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	1097 ## Example 7.7.3 in Golub & Van Loan
13074 0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	1098 %!test
0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	1099 %! a = [ 10 1 2;
0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	1100 %! 1 2 -1;
0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	1101 %! 1 1 2];
14501 60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14360 diff changeset	1102 %! b = reshape (1:9,3,3);
13074 0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	1103 %! [aa, bb, q, z, v, w, lambda] = qz (a, b);
14501 60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14360 diff changeset	1104 %! sz = length (lambda);
60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14360 diff changeset	1105 %! observed = (b * v * diag ([lambda;0])) (:, 1:sz);
18678 6113e0c6920b maint: Clean up extra spaces before/after parentheses. Rik <rik@octave.org> parents: 17787 diff changeset	1106 %! assert ((av)(:, 1:sz), observed, norm (observed) 1e-14);
13088 4ffea87ad71b codesprint: Fix tolerance for qz.cc tests. Ben Abbott <bpabbott@mac.com> parents: 13074 diff changeset	1107 %! observed = (diag ([lambda;0]) * w' * b) (1:sz, :);
18678 6113e0c6920b maint: Clean up extra spaces before/after parentheses. Rik <rik@octave.org> parents: 17787 diff changeset	1108 %! assert ((w'a)(1:sz, :) , observed, norm (observed) 1e-13);
13088 4ffea87ad71b codesprint: Fix tolerance for qz.cc tests. Ben Abbott <bpabbott@mac.com> parents: 13074 diff changeset	1109 %! assert (q * a * z, aa, norm (aa) * 1e-14);
4ffea87ad71b codesprint: Fix tolerance for qz.cc tests. Ben Abbott <bpabbott@mac.com> parents: 13074 diff changeset	1110 %! assert (q * b * z, bb, norm (bb) * 1e-14);
13074 0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	1111
14876 54a386f2ac4e correcly compute Z for 3-output call to qz (bug #36728) John W. Eaton <jwe@octave.org> parents: 14844 diff changeset	1112 %!test
54a386f2ac4e correcly compute Z for 3-output call to qz (bug #36728) John W. Eaton <jwe@octave.org> parents: 14844 diff changeset	1113 %! A = [0, 0, -1, 0; 1, 0, 0, 0; -1, 0, -2, -1; 0, -1, 1, 0];
54a386f2ac4e correcly compute Z for 3-output call to qz (bug #36728) John W. Eaton <jwe@octave.org> parents: 14844 diff changeset	1114 %! B = [0, 0, 0, 0; 0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1];
54a386f2ac4e correcly compute Z for 3-output call to qz (bug #36728) John W. Eaton <jwe@octave.org> parents: 14844 diff changeset	1115 %! [AA, BB, Q, Z1] = qz (A, B);
23430 32d532b8e7d0 qz.cc: Overhaul qz function. Rik <rik@octave.org> parents: 23429 diff changeset	1116 %! [AA, BB, Z2] = qz (A, B, "-");
14876 54a386f2ac4e correcly compute Z for 3-output call to qz (bug #36728) John W. Eaton <jwe@octave.org> parents: 14844 diff changeset	1117 %! assert (Z1, Z2);
13074 0b789a03bde1 codesprint: Add 3 tests for qz.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12642 diff changeset	1118 */

Mercurial > octave

annotate libinterp/corefcn/qz.cc @ 23582:0cc2011d800e