Mercurial > octave
annotate libinterp/corefcn/qz.cc @ 27923:bd51beb6205e
update formatting of copyright notices
* Use <https://octave.org/copyright/> instead of
<https://octave.org/COPYRIGHT.html/>.
* For consistency with other comments in the Octave sources, use
C++-style comments for copyright blocks in C and C++ files.
* Use delimiters above and below copyright blocks that are appropriate
for the language used in the file.
* Eliminate extra spacing inside copyright blocks.
* lex.ll (looks_like_copyright): Also allow newlines and carriage
returns before the word "Copyright".
* scripts/mk-doc.pl (gethelp): Also skip empty comment lines.
* bp-table.cc, type.m: Adjust tests.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Wed, 08 Jan 2020 11:59:41 -0500 |
| parents | 1891570abac8 |
| children | d9d028b479ac 0a5b15007766 |
| 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 // |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
3 // Copyright (C) 1998-2020 The Octave Project Developers |
|
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 | 25 |
| 26 // Generalized eigenvalue balancing via LAPACK | |
| 3911 | 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 | 30 |
| 31 #undef DEBUG | |
| 32 #undef DEBUG_SORT | |
| 33 #undef DEBUG_EIG | |
| 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 | 38 |
| 23430 | 39 #include <cctype> |
| 40 #include <cmath> | |
| 4051 | 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 | 44 #endif |
| 3183 | 45 |
| 4153 | 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 | 49 #include "quit.h" |
| 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 | 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 | 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 | 58 #endif |
| 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 | 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 |
|
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
100 $$ 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
|
101 $$ 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
|
102 $$ AA = Q^T AZ, BB = Q^T BZ $$ |
|
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
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 |
|
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
109 @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
|
110 @var{W}' * @var{A} = diag (@var{lambda}) * @var{W}' * @var{B} |
|
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
111 @var{AA} = @var{Q} * @var{A} * @var{Z}, @var{BB} = @var{Q} * @var{B} * @var{Z} |
|
21966
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
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}. |
|
25734
c7095a755185
New function "ordeig" (patch #9670)
Sébastien Villemot <sebastien@debian.org>
parents:
25570
diff
changeset
|
168 @seealso{eig, ordeig, balance, lu, chol, hess, qr, qzhess, schur, svd} |
|
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 | 170 { |
| 23430 | 171 int nargin = args.length (); |
| 3183 | 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 | 176 #endif |
| 3183 | 177 |
| 3185 | 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 | 182 error ("qz: invalid number of output arguments for form 3 call"); |
| 3183 | 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 | 186 #endif |
| 3183 | 187 |
| 10311 | 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 | 192 |
|
24116
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
193 if (nargin == 3) |
| 3185 | 194 { |
| 23430 | 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 | 201 |
| 23430 | 202 std::string valid_opts = "NSB+-"; |
| 203 | |
| 204 if (valid_opts.find_first_of (ord_job) == std::string::npos) | |
| 205 error ("qz: invalid order option '%c'", ord_job); | |
| 3185 | 206 |
| 207 // overflow constant required by dlag2 | |
| 4552 | 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 | 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 | 216 #endif |
| 3183 | 217 |
| 10311 | 218 // Some machines (e.g., DEC alpha) get safmin = 0; |
| 219 // for these, use eps instead to avoid problems in dlag2. | |
| 3185 | 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 | 225 #endif |
| 3183 | 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 | 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 | 235 #endif |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
236 } |
| 3183 | 237 } |
| 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 | 241 #endif |
| 242 | |
| 23430 | 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 | 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 | 250 #endif |
| 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 | 253 { |
| 23430 | 254 warn_empty_arg ("qz: A"); |
| 3185 | 255 return octave_value_list (2, Matrix ()); |
| 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 | 259 |
| 23430 | 260 // Matrix A: get value. |
| 3183 | 261 Matrix aa; |
| 262 ComplexMatrix caa; | |
| 3185 | 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 | 265 caa = args(0).complex_matrix_value (); |
| 3185 | 266 else |
| 3183 | 267 aa = args(0).matrix_value (); |
| 3185 | 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 | 271 #endif |
| 3183 | 272 |
| 10311 | 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 | 279 |
| 3183 | 280 Matrix bb; |
| 281 ComplexMatrix cbb; | |
| 3185 | 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 | 284 cbb = args(1).complex_matrix_value (); |
| 285 else | |
| 286 bb = args(1).matrix_value (); | |
| 3185 | 287 |
| 23430 | 288 // Both matrices loaded, now check whether to calculate complex or real val. |
| 3185 | 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 | 291 |
| 3185 | 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 | 294 |
| 23430 | 295 // First, declare variables used in both the real and complex cases. |
| 296 // FIXME: There are a lot of excess variables declared. | |
| 297 // Probably a better way to handle this. | |
| 298 Matrix QQ (nn,nn), ZZ (nn,nn), VR (nn,nn), VL (nn,nn); | |
| 299 RowVector alphar (nn), alphai (nn), betar (nn); | |
| 300 ComplexRowVector xalpha (nn), xbeta (nn); | |
| 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 | 303 char comp_q = (nargout >= 3 ? 'V' : 'N'); |
| 304 char comp_z = ((nargout >= 4 || nargin == 3)? 'V' : 'N'); | |
| 3183 | 305 |
| 23430 | 306 // Initialize Q, Z to identity matrix if either is needed |
| 307 if (comp_q == 'V' || comp_z == 'V') | |
| 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 | 311 std::fill_n (QQptr, QQ.numel (), 0.0); |
| 312 std::fill_n (ZZptr, ZZ.numel (), 0.0); | |
| 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 | 315 QQ(i,i) = 1.0; |
| 316 ZZ(i,i) = 1.0; | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
317 } |
| 23430 | 318 } |
| 3183 | 319 |
| 10311 | 320 // Always perform permutation balancing. |
| 4552 | 321 const char bal_job = 'P'; |
| 10311 | 322 RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn); |
| 3183 | 323 |
| 3185 | 324 if (complex_case) |
| 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 | 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 | 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 | 334 cbb = ComplexMatrix (bb); |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
335 |
| 23430 | 336 if (comp_q == 'V') |
| 10311 | 337 CQ = ComplexMatrix (QQ); |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
338 |
| 23430 | 339 if (comp_z == 'V') |
| 10311 | 340 CZ = ComplexMatrix (ZZ); |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
341 |
| 10311 | 342 F77_XFCN (zggbal, ZGGBAL, |
| 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 | 346 nn, ilo, ihi, lscale.fortran_vec (), |
| 347 rscale.fortran_vec (), work.fortran_vec (), info | |
| 348 F77_CHAR_ARG_LEN (1))); | |
| 3185 | 349 } |
| 3183 | 350 else |
| 3185 | 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 | 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 | 355 #endif |
| 3183 | 356 |
| 3185 | 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 | 363 } |
| 3183 | 364 |
| 23430 | 365 // Only permutation balance above is done. Skip scaling balance. |
| 366 | |
| 367 #if 0 | |
| 3183 | 368 // Since we just want the balancing matrices, we can use dggbal |
| 10311 | 369 // for both the real and complex cases; left first |
| 3185 | 370 |
| 23430 | 371 if (comp_q == 'V') |
| 3185 | 372 { |
| 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 | 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 | 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 | 384 #endif |
|
25103
078b795c5219
maint: style check C++ ahead of 4.4 release.
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
385 } |
| 3183 | 386 |
| 387 // then right | |
| 23430 | 388 if (comp_z == 'V') |
| 3185 | 389 { |
| 4552 | 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 | 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 | 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 | 401 #endif |
| 10311 | 402 } |
| 403 #endif | |
| 3183 | 404 |
| 23430 | 405 char qz_job = (nargout < 2 ? 'E' : 'S'); |
| 3185 | 406 |
| 3183 | 407 if (complex_case) |
| 3185 | 408 { |
| 10311 | 409 // Complex case. |
| 410 | |
| 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 | 413 // The R matrix of QR decomposition for cbb. |
| 414 cbb = cbqr.R (); | |
| 415 // (Q*)caa for following work. | |
| 416 caa = (cbqr.Q ().hermitian ()) * caa; | |
| 417 CQ = CQ * cbqr.Q (); | |
| 418 | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
419 F77_XFCN (zgghrd, ZGGHRD, |
| 23430 | 420 (F77_CONST_CHAR_ARG2 (&comp_q, 1), |
| 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 | 428 |
| 23430 | 429 ComplexRowVector cwork (nn); |
| 10311 | 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 | 433 F77_CONST_CHAR_ARG2 (&comp_q, 1), |
| 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 | 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 | 447 |
| 23430 | 448 if (comp_q == 'V') |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
449 { |
| 10311 | 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 | 459 |
| 23430 | 460 if (comp_z == 'V') |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
461 { |
| 23430 | 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 | 471 |
| 472 } | |
| 10311 | 473 else |
| 3185 | 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 | 477 #endif |
| 3183 | 478 |
| 10311 | 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 | 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 | 485 #endif |
| 3183 | 486 |
| 3185 | 487 bb = bqr.R (); |
| 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 | 491 #endif |
| 3183 | 492 |
| 10311 | 493 aa = (bqr.Q ()).transpose () * aa; |
| 3185 | 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 | 498 |
| 23430 | 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 | 501 #endif |
| 3183 | 502 |
| 23430 | 503 if (comp_q == 'V') |
| 10311 | 504 QQ = QQ * bqr.Q (); |
| 3183 | 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 | 508 #endif |
| 3183 | 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 | 513 #endif |
| 3183 | 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 | 516 F77_XFCN (dgghrd, DGGHRD, |
| 23430 | 517 (F77_CONST_CHAR_ARG2 (&comp_q, 1), |
| 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 | 524 |
| 23430 | 525 // Check if just computing generalized eigenvalues, |
| 526 // or if we're actually computing the decomposition. | |
| 3183 | 527 |
| 10311 | 528 // Reduce to generalized Schur form. |
| 3185 | 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 | 531 F77_CONST_CHAR_ARG2 (&comp_q, 1), |
| 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 | 540 |
| 23430 | 541 if (comp_q == 'V') |
| 10311 | 542 { |
| 543 F77_XFCN (dggbak, DGGBAK, | |
| 544 (F77_CONST_CHAR_ARG2 (&bal_job, 1), | |
| 545 F77_CONST_CHAR_ARG2 ("L", 1), | |
| 546 nn, ilo, ihi, lscale.data (), rscale.data (), | |
| 547 nn, QQ.fortran_vec (), nn, info | |
| 548 F77_CHAR_ARG_LEN (1) | |
| 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 | 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 | 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 | 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 | 574 } |
| 3183 | 575 |
| 10311 | 576 // Order the QZ decomposition? |
| 23430 | 577 if (ord_job != 'N') |
| 3183 | 578 { |
| 3185 | 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 | 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 | 586 #endif |
| 3183 | 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 | 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 | 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 | 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 | 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 | 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 | 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 | 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 | 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 | 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 | 659 #endif |
| 3183 | 660 } |
| 3185 | 661 |
| 23430 | 662 // Compute the generalized eigenvalues as well? |
| 3183 | 663 ComplexColumnVector gev; |
| 3185 | 664 |
| 665 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) | |
| 3183 | 666 { |
| 3185 | 667 if (complex_case) |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
668 { |
| 23430 | 669 ComplexColumnVector tmp (nn); |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
670 |
| 23430 | 671 for (F77_INT i = 0; i < nn; i++) |
| 672 tmp(i) = xalpha(i) / xbeta(i); | |
| 10311 | 673 |
| 674 gev = tmp; | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
675 } |
| 3185 | 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 | 680 #endif |
| 3183 | 681 |
| 10311 | 682 // Return finite generalized eigenvalues. |
| 23430 | 683 ComplexColumnVector tmp (nn); |
| 3185 | 684 |
| 23430 | 685 F77_INT cnt = 0; |
| 686 for (F77_INT i = 0; i < nn; i++) | |
| 687 if (betar(i) != 0) | |
| 688 tmp(cnt++) = Complex (alphar(i), alphai(i)) / betar(i); | |
| 3185 | 689 |
| 23430 | 690 tmp.resize (cnt); // Trim vector to number of return values |
| 10311 | 691 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
692 gev = tmp; |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
693 } |
| 3183 | 694 } |
| 695 | |
| 10311 | 696 // Right, left eigenvector matrices. |
| 3185 | 697 if (nargout >= 5) |
| 3183 | 698 { |
| 10311 | 699 // Which side to compute? |
| 700 char side = (nargout == 5 ? 'R' : 'B'); | |
| 701 // Compute all of them and backtransform | |
| 23430 | 702 char howmany = 'B'; |
| 10311 | 703 // 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
|
704 F77_INT *select = nullptr; |
| 3185 | 705 |
| 706 if (complex_case) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
707 { |
| 10311 | 708 CVL = CQ; |
| 709 CVR = CZ; | |
| 710 ComplexRowVector cwork2 (2 * nn); | |
| 711 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
|
712 F77_INT m; |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
713 |
| 10311 | 714 F77_XFCN (ztgevc, ZTGEVC, |
| 715 (F77_CONST_CHAR_ARG2 (&side, 1), | |
| 23430 | 716 F77_CONST_CHAR_ARG2 (&howmany, 1), |
|
22407
34ce5be04942
maint: Style check C++ code in libinterp/.
Rik <rik@octave.org>
parents:
22372
diff
changeset
|
717 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
|
718 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), |
|
34ce5be04942
maint: Style check C++ code in libinterp/.
Rik <rik@octave.org>
parents:
22372
diff
changeset
|
719 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
|
720 F77_DBLE_CMPLX_ARG (CVR.fortran_vec ()), nn, nn, |
| 23430 | 721 m, F77_DBLE_CMPLX_ARG (cwork2.fortran_vec ()), |
| 722 rwork2.fortran_vec (), info | |
| 10311 | 723 F77_CHAR_ARG_LEN (1) |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
724 F77_CHAR_ARG_LEN (1))); |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
725 } |
| 3185 | 726 else |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
727 { |
|
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
728 #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
|
729 octave_stdout << "qz: computing generalized eigenvectors" << std::endl; |
| 3185 | 730 #endif |
| 731 | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
732 VL = QQ; |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
733 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
|
734 F77_INT m; |
| 10311 | 735 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
736 F77_XFCN (dtgevc, DTGEVC, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
737 (F77_CONST_CHAR_ARG2 (&side, 1), |
| 23430 | 738 F77_CONST_CHAR_ARG2 (&howmany, 1), |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
739 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
740 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
|
741 m, work.fortran_vec (), info |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
742 F77_CHAR_ARG_LEN (1) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
743 F77_CHAR_ARG_LEN (1))); |
| 3185 | 744 |
| 10311 | 745 // Now construct the complex form of VV, WW. |
| 23430 | 746 F77_INT j = 0; |
| 3185 | 747 |
| 23430 | 748 while (j < nn) |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
749 { |
|
22860
0b1e25cc4457
eliminate use of OCTAVE_QUIT macro in C++ sources
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
750 octave_quit (); |
| 4153 | 751 |
| 10311 | 752 // See if real or complex eigenvalue. |
| 753 | |
| 754 // Column increment; assume complex eigenvalue. | |
| 755 int cinc = 2; | |
| 3185 | 756 |
| 23430 | 757 if (j == (nn-1)) |
| 10311 | 758 // Single column. |
| 759 cinc = 1; | |
| 23430 | 760 else if (aa(j+1,j) == 0) |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
761 cinc = 1; |
| 3185 | 762 |
| 10311 | 763 // Now copy the eigenvector (s) to CVR, CVL. |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
764 if (cinc == 1) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
765 { |
| 23430 | 766 for (F77_INT i = 0; i < nn; i++) |
| 767 CVR(i,j) = VR(i,j); | |
| 3185 | 768 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
769 if (side == 'B') |
| 23430 | 770 for (F77_INT i = 0; i < nn; i++) |
| 771 CVL(i,j) = VL(i,j); | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
772 } |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
773 else |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
774 { |
| 10311 | 775 // Double column; complex vector. |
| 3185 | 776 |
| 23430 | 777 for (F77_INT i = 0; i < nn; i++) |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
778 { |
| 23430 | 779 CVR(i,j) = Complex (VR(i,j), VR(i,j+1)); |
| 780 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
|
781 } |
| 3183 | 782 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
783 if (side == 'B') |
| 23430 | 784 for (F77_INT i = 0; i < nn; i++) |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
785 { |
| 23430 | 786 CVL(i,j) = Complex (VL(i,j), VL(i,j+1)); |
| 787 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
|
788 } |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
789 } |
| 3185 | 790 |
| 10311 | 791 // Advance to next eigenvectors (if any). |
| 23430 | 792 j += cinc; |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
793 } |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
794 } |
| 3183 | 795 } |
| 3185 | 796 |
| 23430 | 797 octave_value_list retval (nargout); |
| 798 | |
| 3185 | 799 switch (nargout) |
| 800 { | |
| 801 case 7: | |
| 802 retval(6) = gev; | |
|
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
803 OCTAVE_FALLTHROUGH; |
| 3185 | 804 |
| 10311 | 805 case 6: |
| 23430 | 806 // Return left eigenvectors. |
| 3185 | 807 retval(5) = CVL; |
|
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
808 OCTAVE_FALLTHROUGH; |
| 3185 | 809 |
| 10311 | 810 case 5: |
| 23430 | 811 // Return right eigenvectors. |
| 3185 | 812 retval(4) = CVR; |
|
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
813 OCTAVE_FALLTHROUGH; |
| 3185 | 814 |
| 815 case 4: | |
| 816 if (nargin == 3) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
817 { |
|
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
818 #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
|
819 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
|
820 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
|
821 octave_stdout << std::endl; |
| 3185 | 822 #endif |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
823 retval(3) = gev; |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
824 } |
| 3185 | 825 else |
| 10311 | 826 { |
| 827 if (complex_case) | |
| 828 retval(3) = CZ; | |
| 829 else | |
| 830 retval(3) = ZZ; | |
| 831 } | |
|
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
832 OCTAVE_FALLTHROUGH; |
| 3185 | 833 |
| 834 case 3: | |
| 835 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
|
836 { |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
837 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
|
838 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
|
839 else |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
840 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
|
841 } |
| 3185 | 842 else |
| 10311 | 843 { |
| 844 if (complex_case) | |
| 845 retval(2) = CQ.hermitian (); | |
| 846 else | |
| 847 retval(2) = QQ.transpose (); | |
| 848 } | |
|
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
849 OCTAVE_FALLTHROUGH; |
| 10311 | 850 |
| 3185 | 851 case 2: |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
852 { |
| 10311 | 853 if (complex_case) |
| 854 { | |
|
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
855 #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
|
856 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
|
857 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
|
858 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
|
859 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
|
860 octave_stdout << std::endl; |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
861 #endif |
| 10311 | 862 retval(1) = cbb; |
| 863 retval(0) = caa; | |
| 864 } | |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
865 else |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
866 { |
|
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
867 #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
|
868 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
|
869 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
|
870 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
|
871 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
|
872 octave_stdout << std::endl; |
| 3185 | 873 #endif |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
874 retval(1) = bb; |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
875 retval(0) = aa; |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
876 } |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
877 } |
| 10311 | 878 break; |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
879 |
| 3185 | 880 case 1: |
| 881 case 0: | |
|
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
882 #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
|
883 octave_stdout << "qz: retval(0) = gev = " << gev << std::endl; |
| 3185 | 884 #endif |
| 885 retval(0) = gev; | |
| 886 break; | |
| 887 | |
| 888 default: | |
|
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
889 error ("qz: too many return arguments"); |
| 3185 | 890 break; |
|
25103
078b795c5219
maint: style check C++ ahead of 4.4 release.
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
891 } |
| 3183 | 892 |
|
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
893 #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
|
894 octave_stdout << "qz: exiting (at long last)" << std::endl; |
| 3185 | 895 #endif |
| 3183 | 896 |
| 897 return retval; | |
| 898 } | |
|
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
899 |
|
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 %!shared a, b, c |
|
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
902 %! a = [1 2; 0 3]; |
|
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
903 %! b = [1 0; 0 0]; |
|
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
904 %! c = [0 1; 0 0]; |
|
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
905 %!assert (qz (a,b), 1) |
|
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
906 %!assert (isempty (qz (a,c))) |
|
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
907 |
| 23430 | 908 ## 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
|
909 %!test |
|
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
910 %! a = [ 10 1 2; |
|
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
911 %! 1 2 -1; |
|
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
912 %! 1 1 2]; |
|
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
913 %! 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
|
914 %! [aa, bb, q, z, v, w, lambda] = qz (a, b); |
|
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
915 %! sz = length (lambda); |
|
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
916 %! observed = (b * v * diag ([lambda;0])) (:, 1:sz); |
|
18678
6113e0c6920b
maint: Clean up extra spaces before/after parentheses.
Rik <rik@octave.org>
parents:
17787
diff
changeset
|
917 %! assert ((a*v)(:, 1:sz), observed, norm (observed) * 1e-14); |
|
13088
4ffea87ad71b
codesprint: Fix tolerance for qz.cc tests.
Ben Abbott <bpabbott@mac.com>
parents:
13074
diff
changeset
|
918 %! observed = (diag ([lambda;0]) * w' * b) (1:sz, :); |
|
18678
6113e0c6920b
maint: Clean up extra spaces before/after parentheses.
Rik <rik@octave.org>
parents:
17787
diff
changeset
|
919 %! assert ((w'*a)(1:sz, :) , observed, norm (observed) * 1e-13); |
|
13088
4ffea87ad71b
codesprint: Fix tolerance for qz.cc tests.
Ben Abbott <bpabbott@mac.com>
parents:
13074
diff
changeset
|
920 %! 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
|
921 %! assert (q * b * z, bb, norm (bb) * 1e-14); |
|
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
922 |
|
14876
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
923 %!test |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
924 %! 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
|
925 %! 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
|
926 %! [AA, BB, Q, Z1] = qz (A, B); |
| 23430 | 927 %! [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
|
928 %! 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
|
929 |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
930 %!test |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
931 %! 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
|
932 %! 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
|
933 %! -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
|
934 %! -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
|
935 %! 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
|
936 %! 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
|
937 %! 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
|
938 %! 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
|
939 %! [AA, BB, Z, lambda] = qz (A, B, "+"); |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
940 %! assert (all (real (lambda(1:3)) >= 0)) |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
941 %! assert (real (lambda(4) < 0)) |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
942 %! [AA, BB, Z, lambda] = qz (A, B, "-"); |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
943 %! assert (real (lambda(1) < 0)) |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
944 %! assert (all (real (lambda(2:4)) >= 0)) |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
945 %! [AA, BB, Z, lambda] = qz (A, B, "B"); |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
946 %! assert (all (abs (lambda(1:3)) >= 1)) |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
947 %! assert (abs (lambda(4) < 1)) |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
948 %! [AA, BB, Z, lambda] = qz (A, B, "S"); |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
949 %! assert (abs (lambda(1) < 1)) |
|
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
950 %! 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
|
951 */ |