octave-lyh: libinterp/corefcn/qz.cc annotate

annotate libinterp/corefcn/qz.cc @ 17501:278ef6bd821d

isprop.m: Overhaul function. * scripts/plot/isprop.m: Add note that input may also be an array of object handles. Use variable names in error() messages. Avoid unnecessary catch block in try/catch statement. Add %!error input validation tests.

author	Rik <rik@octave.org>
date	Wed, 25 Sep 2013 11:57:10 -0700
parents	bc924baa2c4e
children

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

Mercurial > hg > octave-lyh

annotate libinterp/corefcn/qz.cc @ 17501:278ef6bd821d