octave: libinterp/corefcn/qz.cc annotate

annotate libinterp/corefcn/qz.cc @ 29556:d75aa2bf4915 stable

qz.cc: return correct number of eigenvalues (bug #60357). * qz.cc (Fqz): Change docstring to emphasize that QZ factorization is the main output, and eigenvalues are secondary. Don't truncate number of eigenvalues to the number returned by LAPACK. For eigenvalues that are not defined, return Inf as the eig() function does. (grafted from 06f4bb661d966b4ffd504a29685bb265b6350395)

author	Rik <rik@octave.org>
date	Fri, 23 Apr 2021 08:04:48 -0700
parents	0a5b15007766
children	50aa945836c3

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

Mercurial > octave

annotate libinterp/corefcn/qz.cc @ 29556:d75aa2bf4915 stable