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