Mercurial > octave
annotate libinterp/corefcn/qz.cc @ 25570:da2b38b7a077 stable
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
* qz.cc: Use LAPACK's DTGSEN to compute and update generalized
eigenvalues. Selection of eigenvalue cluster is now done by filling a
vector of booleans, rather than through a selection function. New
tests for the four variants of the ordered QZ.
* lo-lapack-proto.h (DTGSEN): New prototype.
* liboctave/external/ordered-qz: Delete directory.
* liboctave/external/module.mk: Update.
author | Sébastien Villemot <sebastien@debian.org> |
---|---|
date | Sun, 13 May 2018 12:37:37 +0200 |
parents | d7ad543255c5 |
children | c7095a755185 |
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 |
23430 | 39 #if defined (DEBUG) || defined (DEBUG_SORT) || defined (DEBUG_EIG) |
40 # include <iostream> | |
41 # if defined (DEBUG_EIG) | |
42 # include <iomanip> | |
43 # endif | |
44 #endif | |
3183 | 45 |
4153 | 46 #include "f77-fcn.h" |
22322
93b3cdd36854
move most f77 function decls to separate header files
John W. Eaton <jwe@octave.org>
parents:
22317
diff
changeset
|
47 #include "lo-lapack-proto.h" |
21279
eb1524b07fe3
better use of templates for qr classes
John W. Eaton <jwe@octave.org>
parents:
21200
diff
changeset
|
48 #include "qr.h" |
4153 | 49 #include "quit.h" |
50 | |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14876
diff
changeset
|
51 #include "defun.h" |
3183 | 52 #include "error.h" |
21100
e39e05d90788
Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents:
21024
diff
changeset
|
53 #include "errwarn.h" |
20940
48b2ad5ee801
maint: Rename oct-obj.[cc|h] to ovl.[cc|h] for clarity.
Rik <rik@octave.org>
parents:
20892
diff
changeset
|
54 #include "ovl.h" |
3185 | 55 #if defined (DEBUG) || defined (DEBUG_SORT) |
21200
fcac5dbbf9ed
maint: Indent #ifdef blocks in libinterp.
Rik <rik@octave.org>
parents:
21110
diff
changeset
|
56 # include "pr-output.h" |
3183 | 57 #endif |
58 | |
25570
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
59 // 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
|
60 // Compatibility problem? |
3185 | 61 |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14876
diff
changeset
|
62 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
|
63 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
|
64 @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
|
65 @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
|
66 @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
|
67 @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
|
68 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
|
69 |
d7ad543255c5
doc: Shorten very long first sentences of docstrings (bug #53388).
Rik <rik@octave.org>
parents:
25103
diff
changeset
|
70 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
|
71 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
72 @tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
73 $$A x = \lambda B x$$ |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
74 @end tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
75 @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
|
76 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
77 @math{A x = @var{lambda} B x} |
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 @end ifnottex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
80 |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
81 There are three calling forms of the function: |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
82 |
21966
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
83 @enumerate |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
84 @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
|
85 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
86 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
|
87 @tex |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
88 $\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
|
89 @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
|
90 @ifnottex |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
91 @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
|
92 @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
|
93 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
94 @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
|
95 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
96 Compute QZ@tie{}decomposition, generalized eigenvectors, and generalized |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
97 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
|
98 @tex |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
99 $$ 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
|
100 $$ 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
|
101 $$ 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
|
102 @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
|
103 @ifnottex |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
104 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
105 @example |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
106 @group |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
107 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
108 @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
|
109 @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
|
110 @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
|
111 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
112 @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
|
113 @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
|
114 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
115 @end ifnottex |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
116 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
|
117 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
118 @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
|
119 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
120 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
|
121 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
|
122 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
|
123 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
|
124 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
125 @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
|
126 @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
|
127 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
|
128 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
|
129 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
130 @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
|
131 @item @qcode{"N"} |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
132 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
|
133 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
134 @item @qcode{"S"} |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
135 small: leading block has all |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
136 @tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
137 $|\lambda| < 1$ |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
138 @end tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
139 @ifnottex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
140 |@var{lambda}| < 1 |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
141 @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
|
142 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
143 @item @qcode{"B"} |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
144 big: leading block has all |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
145 @tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
146 $|\lambda| \geq 1$ |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
147 @end tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
148 @ifnottex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
149 |@var{lambda}| @geq{} 1 |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
150 @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
|
151 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
152 @item @qcode{"-"} |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
153 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
|
154 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
|
155 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
156 @item @qcode{"+"} |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
157 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
|
158 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
|
159 @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
|
160 @end table |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
161 @end enumerate |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
162 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
163 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
|
164 (@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
|
165 @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
|
166 @sc{matlab}. |
21966
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
167 @seealso{eig, balance, lu, chol, hess, qr, qzhess, schur, svd} |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
168 @end deftypefn */) |
3183 | 169 { |
23430 | 170 int nargin = args.length (); |
3183 | 171 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
172 #if defined (DEBUG) |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
173 std::cout << "qz: nargin = " << nargin |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
174 << ", nargout = " << nargout << std::endl; |
3185 | 175 #endif |
3183 | 176 |
3185 | 177 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
|
178 print_usage (); |
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20700
diff
changeset
|
179 |
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20700
diff
changeset
|
180 if (nargin == 3 && (nargout < 3 || nargout > 4)) |
23430 | 181 error ("qz: invalid number of output arguments for form 3 call"); |
3183 | 182 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
183 #if defined (DEBUG) |
3538 | 184 std::cout << "qz: determine ordering option" << std::endl; |
3185 | 185 #endif |
3183 | 186 |
10311 | 187 // 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
|
188 |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
189 char ord_job = 'N'; |
24116
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
190 double safmin = 0.0; |
3185 | 191 |
24116
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
192 if (nargin == 3) |
3185 | 193 { |
23430 | 194 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
|
195 |
24116
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
196 if (opt.empty ()) |
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
197 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
|
198 |
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
199 ord_job = std::toupper (opt[0]); |
3183 | 200 |
23430 | 201 std::string valid_opts = "NSB+-"; |
202 | |
203 if (valid_opts.find_first_of (ord_job) == std::string::npos) | |
204 error ("qz: invalid order option '%c'", ord_job); | |
3185 | 205 |
206 // overflow constant required by dlag2 | |
4552 | 207 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
|
208 safmin |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
209 F77_CHAR_ARG_LEN (1)); |
3183 | 210 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
211 #if defined (DEBUG_EIG) |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
212 std::cout << "qz: initial value of safmin=" |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
213 << setiosflags (std::ios::scientific) |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
214 << safmin << std::endl; |
3185 | 215 #endif |
3183 | 216 |
10311 | 217 // Some machines (e.g., DEC alpha) get safmin = 0; |
218 // for these, use eps instead to avoid problems in dlag2. | |
3185 | 219 if (safmin == 0) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
220 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
221 #if defined (DEBUG_EIG) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
222 std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl; |
3185 | 223 #endif |
3183 | 224 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
225 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
|
226 safmin |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
227 F77_CHAR_ARG_LEN (1)); |
3185 | 228 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
229 #if defined (DEBUG_EIG) |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
230 std::cout << "qz: safmin set to " |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
231 << setiosflags (std::ios::scientific) |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
232 << safmin << std::endl; |
3185 | 233 #endif |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
234 } |
3183 | 235 } |
236 | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
237 #if defined (DEBUG) |
23430 | 238 std::cout << "qz: check matrix A" << std::endl; |
3185 | 239 #endif |
240 | |
23430 | 241 // Matrix A: check dimensions. |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22964
diff
changeset
|
242 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
|
243 F77_INT nc = octave::to_f77_int (args(0).columns ()); |
20892 | 244 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
245 #if defined (DEBUG) |
23807
336f89b6208b
Use character literals 'c' rather than string literals "c" when possible.
Rik <rik@octave.org>
parents:
23795
diff
changeset
|
246 std::cout << "Matrix A dimensions: (" << nn << ',' << nc << ')' << 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) |
23430 | 267 std::cout << "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') |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
325 std::cout << "qz: performing balancing; CQ=" << std::endl |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
326 << CQ << std::endl; |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
327 #endif |
23582
0cc2011d800e
maint: Deprecate is_real_type and replace with isreal.
Rik <rik@octave.org>
parents:
23581
diff
changeset
|
328 if (args(0).isreal ()) |
10311 | 329 caa = ComplexMatrix (aa); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
330 |
23582
0cc2011d800e
maint: Deprecate is_real_type and replace with isreal.
Rik <rik@octave.org>
parents:
23581
diff
changeset
|
331 if (args(1).isreal ()) |
10311 | 332 cbb = ComplexMatrix (bb); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
333 |
23430 | 334 if (comp_q == 'V') |
10311 | 335 CQ = ComplexMatrix (QQ); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
336 |
23430 | 337 if (comp_z == 'V') |
10311 | 338 CZ = ComplexMatrix (ZZ); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
339 |
10311 | 340 F77_XFCN (zggbal, ZGGBAL, |
341 (F77_CONST_CHAR_ARG2 (&bal_job, 1), | |
22407
34ce5be04942
maint: Style check C++ code in libinterp/.
Rik <rik@octave.org>
parents:
22372
diff
changeset
|
342 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
|
343 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), |
10311 | 344 nn, ilo, ihi, lscale.fortran_vec (), |
345 rscale.fortran_vec (), work.fortran_vec (), info | |
346 F77_CHAR_ARG_LEN (1))); | |
3185 | 347 } |
3183 | 348 else |
3185 | 349 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
350 #if defined (DEBUG) |
23430 | 351 if (comp_q == 'V') |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
352 std::cout << "qz: performing balancing; QQ=" << std::endl |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
353 << QQ << std::endl; |
3185 | 354 #endif |
3183 | 355 |
3185 | 356 F77_XFCN (dggbal, DGGBAL, |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
357 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
358 nn, aa.fortran_vec (), nn, bb.fortran_vec (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
359 nn, ilo, ihi, lscale.fortran_vec (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
360 rscale.fortran_vec (), work.fortran_vec (), info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
361 F77_CHAR_ARG_LEN (1))); |
3185 | 362 } |
3183 | 363 |
23430 | 364 // Only permutation balance above is done. Skip scaling balance. |
365 | |
366 #if 0 | |
3183 | 367 // Since we just want the balancing matrices, we can use dggbal |
10311 | 368 // for both the real and complex cases; left first |
3185 | 369 |
23430 | 370 if (comp_q == 'V') |
3185 | 371 { |
372 F77_XFCN (dggbak, DGGBAK, | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
373 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
374 F77_CONST_CHAR_ARG2 ("L", 1), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
375 nn, ilo, ihi, lscale.data (), rscale.data (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
376 nn, QQ.fortran_vec (), nn, info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
377 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
378 F77_CHAR_ARG_LEN (1))); |
3183 | 379 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
380 #if defined (DEBUG) |
23430 | 381 if (comp_q == 'V') |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
382 std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl; |
3185 | 383 #endif |
25103
078b795c5219
maint: style check C++ ahead of 4.4 release.
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
384 } |
3183 | 385 |
386 // then right | |
23430 | 387 if (comp_z == 'V') |
3185 | 388 { |
4552 | 389 F77_XFCN (dggbak, DGGBAK, |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
390 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
391 F77_CONST_CHAR_ARG2 ("R", 1), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
392 nn, ilo, ihi, lscale.data (), rscale.data (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
393 nn, ZZ.fortran_vec (), nn, info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
394 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
395 F77_CHAR_ARG_LEN (1))); |
3183 | 396 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
397 #if defined (DEBUG) |
23430 | 398 if (comp_z == 'V') |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
399 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl; |
3185 | 400 #endif |
10311 | 401 } |
402 #endif | |
3183 | 403 |
23430 | 404 char qz_job = (nargout < 2 ? 'E' : 'S'); |
3185 | 405 |
3183 | 406 if (complex_case) |
3185 | 407 { |
10311 | 408 // Complex case. |
409 | |
410 // 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
|
411 octave::math::qr<ComplexMatrix> cbqr (cbb); |
10311 | 412 // The R matrix of QR decomposition for cbb. |
413 cbb = cbqr.R (); | |
414 // (Q*)caa for following work. | |
415 caa = (cbqr.Q ().hermitian ()) * caa; | |
416 CQ = CQ * cbqr.Q (); | |
417 | |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
418 F77_XFCN (zgghrd, ZGGHRD, |
23430 | 419 (F77_CONST_CHAR_ARG2 (&comp_q, 1), |
420 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
|
421 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
|
422 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
|
423 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
|
424 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
|
425 F77_CHAR_ARG_LEN (1) |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
426 F77_CHAR_ARG_LEN (1))); |
10311 | 427 |
23430 | 428 ComplexRowVector cwork (nn); |
10311 | 429 |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
430 F77_XFCN (zhgeqz, ZHGEQZ, |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
431 (F77_CONST_CHAR_ARG2 (&qz_job, 1), |
23430 | 432 F77_CONST_CHAR_ARG2 (&comp_q, 1), |
433 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
|
434 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
|
435 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
|
436 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
|
437 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
|
438 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
|
439 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
|
440 F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn, |
23430 | 441 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
|
442 rwork.fortran_vec (), info |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
443 F77_CHAR_ARG_LEN (1) |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
444 F77_CHAR_ARG_LEN (1) |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
445 F77_CHAR_ARG_LEN (1))); |
3185 | 446 |
23430 | 447 if (comp_q == 'V') |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
448 { |
10311 | 449 // Left eigenvector. |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
450 F77_XFCN (zggbak, ZGGBAK, |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
451 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
452 F77_CONST_CHAR_ARG2 ("L", 1), |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
453 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
|
454 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
|
455 F77_CHAR_ARG_LEN (1) |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
456 F77_CHAR_ARG_LEN (1))); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
457 } |
3185 | 458 |
23430 | 459 if (comp_z == 'V') |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
460 { |
23430 | 461 // Right eigenvector. |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
462 F77_XFCN (zggbak, ZGGBAK, |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
463 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
464 F77_CONST_CHAR_ARG2 ("R", 1), |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
465 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
|
466 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
|
467 F77_CHAR_ARG_LEN (1) |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
468 F77_CHAR_ARG_LEN (1))); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
469 } |
3185 | 470 |
471 } | |
10311 | 472 else |
3185 | 473 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
474 #if defined (DEBUG) |
23430 | 475 std::cout << "qz: performing qr decomposition of bb" << std::endl; |
3185 | 476 #endif |
3183 | 477 |
10311 | 478 // 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
|
479 octave::math::qr<Matrix> bqr (bb); |
3185 | 480 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
481 #if defined (DEBUG) |
23430 | 482 std::cout << "qz: qr (bb) done; now performing qz decomposition" |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
483 << std::endl; |
3185 | 484 #endif |
3183 | 485 |
3185 | 486 bb = bqr.R (); |
487 | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
488 #if defined (DEBUG) |
3538 | 489 std::cout << "qz: extracted bb" << std::endl; |
3185 | 490 #endif |
3183 | 491 |
10311 | 492 aa = (bqr.Q ()).transpose () * aa; |
3185 | 493 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
494 #if defined (DEBUG) |
3538 | 495 std::cout << "qz: updated aa " << std::endl; |
496 std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl; | |
3183 | 497 |
23430 | 498 if (comp_q == 'V') |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
499 std::cout << "QQ =" << QQ << std::endl; |
3185 | 500 #endif |
3183 | 501 |
23430 | 502 if (comp_q == 'V') |
10311 | 503 QQ = QQ * bqr.Q (); |
3183 | 504 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
505 #if defined (DEBUG) |
3538 | 506 std::cout << "qz: precursors done..." << std::endl; |
3185 | 507 #endif |
3183 | 508 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
509 #if defined (DEBUG) |
23430 | 510 std::cout << "qz: comp_q = " << comp_q << ", comp_z = " << comp_z |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
511 << std::endl; |
3185 | 512 #endif |
3183 | 513 |
21562
6c2fd62db1f7
maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
514 // Reduce to generalized Hessenberg form. |
3185 | 515 F77_XFCN (dgghrd, DGGHRD, |
23430 | 516 (F77_CONST_CHAR_ARG2 (&comp_q, 1), |
517 F77_CONST_CHAR_ARG2 (&comp_z, 1), | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
518 nn, ilo, ihi, aa.fortran_vec (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
519 nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
520 ZZ.fortran_vec (), nn, info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
521 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
522 F77_CHAR_ARG_LEN (1))); |
3183 | 523 |
23430 | 524 // Check if just computing generalized eigenvalues, |
525 // or if we're actually computing the decomposition. | |
3183 | 526 |
10311 | 527 // Reduce to generalized Schur form. |
3185 | 528 F77_XFCN (dhgeqz, DHGEQZ, |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
529 (F77_CONST_CHAR_ARG2 (&qz_job, 1), |
23430 | 530 F77_CONST_CHAR_ARG2 (&comp_q, 1), |
531 F77_CONST_CHAR_ARG2 (&comp_z, 1), | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
532 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
|
533 nn, alphar.fortran_vec (), alphai.fortran_vec (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
534 betar.fortran_vec (), QQ.fortran_vec (), nn, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
535 ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
536 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
537 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
538 F77_CHAR_ARG_LEN (1))); |
10311 | 539 |
23430 | 540 if (comp_q == 'V') |
10311 | 541 { |
542 F77_XFCN (dggbak, DGGBAK, | |
543 (F77_CONST_CHAR_ARG2 (&bal_job, 1), | |
544 F77_CONST_CHAR_ARG2 ("L", 1), | |
545 nn, ilo, ihi, lscale.data (), rscale.data (), | |
546 nn, QQ.fortran_vec (), nn, info | |
547 F77_CHAR_ARG_LEN (1) | |
548 F77_CHAR_ARG_LEN (1))); | |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
549 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
550 #if defined (DEBUG) |
23430 | 551 if (comp_q == 'V') |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
552 std::cout << "qz: balancing done; QQ=" << std::endl |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
553 << QQ << std::endl; |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
554 #endif |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
555 } |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
556 |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
557 // then right |
23430 | 558 if (comp_z == 'V') |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
559 { |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
560 F77_XFCN (dggbak, DGGBAK, |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
561 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
562 F77_CONST_CHAR_ARG2 ("R", 1), |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
563 nn, ilo, ihi, lscale.data (), rscale.data (), |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
564 nn, ZZ.fortran_vec (), nn, info |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
565 F77_CHAR_ARG_LEN (1) |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
566 F77_CHAR_ARG_LEN (1))); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
567 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
568 #if defined (DEBUG) |
23430 | 569 if (comp_z == 'V') |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
570 std::cout << "qz: balancing done; ZZ=" << std::endl |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
571 << ZZ << std::endl; |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
572 #endif |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
573 } |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
574 |
3185 | 575 } |
3183 | 576 |
10311 | 577 // Order the QZ decomposition? |
23430 | 578 if (ord_job != 'N') |
3183 | 579 { |
3185 | 580 if (complex_case) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
581 // Probably not needed, but better be safe. |
23430 | 582 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
|
583 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
584 #if defined (DEBUG_SORT) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
585 std::cout << "qz: ordering eigenvalues: ord_job = " |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
586 << ord_job << std::endl; |
3185 | 587 #endif |
3183 | 588 |
25570
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
589 Array<F77_LOGICAL> select (dim_vector (nn, 1)); |
3183 | 590 |
25570
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
591 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
|
592 { |
25570
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
593 switch (ord_job) |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
594 { |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
595 case 'S': |
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 'B': |
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)*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
|
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; |
3183 | 606 |
25570
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
607 case '-': |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
608 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
|
609 break; |
3185 | 610 |
25570
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
611 default: |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
612 // Invalid order option |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
613 // (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
|
614 panic_impossible (); |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
615 break; |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
616 } |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
617 } |
3183 | 618 |
25570
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
619 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
|
620 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
|
621 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
|
622 RowVector rwork3(lrwork3); |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
623 Array<F77_INT> iwork (dim_vector (liwork, 1)); |
3185 | 624 |
25570
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
625 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
|
626 (ijob, wantq, wantz, |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
627 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
|
628 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
|
629 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
|
630 alphar.fortran_vec (), |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
631 alphai.fortran_vec (), |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
632 betar.fortran_vec (), |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
633 nullptr, nn, |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
634 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
|
635 mm, |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
636 pl, pr, |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
637 nullptr, |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
638 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
|
639 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
|
640 info)); |
3185 | 641 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
642 #if defined (DEBUG_SORT) |
25570
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
643 std::cout << "qz: back from dtgsen: aa=" << std::endl; |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
644 octave_print_internal (std::cout, aa, 0); |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
645 std::cout << std::endl << "bb=" << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
646 octave_print_internal (std::cout, bb, 0); |
3185 | 647 |
23430 | 648 if (comp_z == 'V') |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
649 { |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
650 std::cout << std::endl << "ZZ=" << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
651 octave_print_internal (std::cout, ZZ, 0); |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
652 } |
25570
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
653 std::cout << std::endl << "qz: info=" << info << std::endl; |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
654 std::cout << "alphar = " << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
655 octave_print_internal (std::cout, (Matrix) alphar, 0); |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
656 std::cout << std::endl << "alphai = " << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
657 octave_print_internal (std::cout, (Matrix) alphai, 0); |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
658 std::cout << std::endl << "beta = " << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
659 octave_print_internal (std::cout, (Matrix) betar, 0); |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
660 std::cout << std::endl; |
3185 | 661 #endif |
3183 | 662 } |
3185 | 663 |
23430 | 664 // Compute the generalized eigenvalues as well? |
3183 | 665 ComplexColumnVector gev; |
3185 | 666 |
667 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) | |
3183 | 668 { |
3185 | 669 if (complex_case) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
670 { |
23430 | 671 ComplexColumnVector tmp (nn); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
672 |
23430 | 673 for (F77_INT i = 0; i < nn; i++) |
674 tmp(i) = xalpha(i) / xbeta(i); | |
10311 | 675 |
676 gev = tmp; | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
677 } |
3185 | 678 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
679 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
680 #if defined (DEBUG) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
681 std::cout << "qz: computing generalized eigenvalues" << std::endl; |
3185 | 682 #endif |
3183 | 683 |
10311 | 684 // Return finite generalized eigenvalues. |
23430 | 685 ComplexColumnVector tmp (nn); |
3185 | 686 |
23430 | 687 F77_INT cnt = 0; |
688 for (F77_INT i = 0; i < nn; i++) | |
689 if (betar(i) != 0) | |
690 tmp(cnt++) = Complex (alphar(i), alphai(i)) / betar(i); | |
3185 | 691 |
23430 | 692 tmp.resize (cnt); // Trim vector to number of return values |
10311 | 693 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
694 gev = tmp; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
695 } |
3183 | 696 } |
697 | |
10311 | 698 // Right, left eigenvector matrices. |
3185 | 699 if (nargout >= 5) |
3183 | 700 { |
10311 | 701 // Which side to compute? |
702 char side = (nargout == 5 ? 'R' : 'B'); | |
703 // Compute all of them and backtransform | |
23430 | 704 char howmany = 'B'; |
10311 | 705 // 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
|
706 F77_INT *select = nullptr; |
3185 | 707 |
708 if (complex_case) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
709 { |
10311 | 710 CVL = CQ; |
711 CVR = CZ; | |
712 ComplexRowVector cwork2 (2 * nn); | |
713 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
|
714 F77_INT m; |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
715 |
10311 | 716 F77_XFCN (ztgevc, ZTGEVC, |
717 (F77_CONST_CHAR_ARG2 (&side, 1), | |
23430 | 718 F77_CONST_CHAR_ARG2 (&howmany, 1), |
22407
34ce5be04942
maint: Style check C++ code in libinterp/.
Rik <rik@octave.org>
parents:
22372
diff
changeset
|
719 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
|
720 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), |
34ce5be04942
maint: Style check C++ code in libinterp/.
Rik <rik@octave.org>
parents:
22372
diff
changeset
|
721 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
|
722 F77_DBLE_CMPLX_ARG (CVR.fortran_vec ()), nn, nn, |
23430 | 723 m, F77_DBLE_CMPLX_ARG (cwork2.fortran_vec ()), |
724 rwork2.fortran_vec (), info | |
10311 | 725 F77_CHAR_ARG_LEN (1) |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
726 F77_CHAR_ARG_LEN (1))); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
727 } |
3185 | 728 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
729 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
730 #if defined (DEBUG) |
21562
6c2fd62db1f7
maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
731 std::cout << "qz: computing generalized eigenvectors" << std::endl; |
3185 | 732 #endif |
733 | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
734 VL = QQ; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
735 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
|
736 F77_INT m; |
10311 | 737 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
738 F77_XFCN (dtgevc, DTGEVC, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
739 (F77_CONST_CHAR_ARG2 (&side, 1), |
23430 | 740 F77_CONST_CHAR_ARG2 (&howmany, 1), |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
741 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
742 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
|
743 m, work.fortran_vec (), info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
744 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
745 F77_CHAR_ARG_LEN (1))); |
3185 | 746 |
10311 | 747 // Now construct the complex form of VV, WW. |
23430 | 748 F77_INT j = 0; |
3185 | 749 |
23430 | 750 while (j < nn) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
751 { |
22860
0b1e25cc4457
eliminate use of OCTAVE_QUIT macro in C++ sources
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
752 octave_quit (); |
4153 | 753 |
10311 | 754 // See if real or complex eigenvalue. |
755 | |
756 // Column increment; assume complex eigenvalue. | |
757 int cinc = 2; | |
3185 | 758 |
23430 | 759 if (j == (nn-1)) |
10311 | 760 // Single column. |
761 cinc = 1; | |
23430 | 762 else if (aa(j+1,j) == 0) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
763 cinc = 1; |
3185 | 764 |
10311 | 765 // Now copy the eigenvector (s) to CVR, CVL. |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
766 if (cinc == 1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
767 { |
23430 | 768 for (F77_INT i = 0; i < nn; i++) |
769 CVR(i,j) = VR(i,j); | |
3185 | 770 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
771 if (side == 'B') |
23430 | 772 for (F77_INT i = 0; i < nn; i++) |
773 CVL(i,j) = VL(i,j); | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
774 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
775 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
776 { |
10311 | 777 // Double column; complex vector. |
3185 | 778 |
23430 | 779 for (F77_INT i = 0; i < nn; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
780 { |
23430 | 781 CVR(i,j) = Complex (VR(i,j), VR(i,j+1)); |
782 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
|
783 } |
3183 | 784 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
785 if (side == 'B') |
23430 | 786 for (F77_INT i = 0; i < nn; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
787 { |
23430 | 788 CVL(i,j) = Complex (VL(i,j), VL(i,j+1)); |
789 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
|
790 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
791 } |
3185 | 792 |
10311 | 793 // Advance to next eigenvectors (if any). |
23430 | 794 j += cinc; |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
795 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
796 } |
3183 | 797 } |
3185 | 798 |
23430 | 799 octave_value_list retval (nargout); |
800 | |
3185 | 801 switch (nargout) |
802 { | |
803 case 7: | |
804 retval(6) = gev; | |
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
805 OCTAVE_FALLTHROUGH; |
3185 | 806 |
10311 | 807 case 6: |
23430 | 808 // Return left eigenvectors. |
3185 | 809 retval(5) = CVL; |
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
810 OCTAVE_FALLTHROUGH; |
3185 | 811 |
10311 | 812 case 5: |
23430 | 813 // Return right eigenvectors. |
3185 | 814 retval(4) = CVR; |
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
815 OCTAVE_FALLTHROUGH; |
3185 | 816 |
817 case 4: | |
818 if (nargin == 3) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
819 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
820 #if defined (DEBUG) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
821 std::cout << "qz: sort: retval(3) = gev = " << std::endl; |
23430 | 822 octave_print_internal (std::cout, ComplexMatrix (gev)); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
823 std::cout << std::endl; |
3185 | 824 #endif |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
825 retval(3) = gev; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
826 } |
3185 | 827 else |
10311 | 828 { |
829 if (complex_case) | |
830 retval(3) = CZ; | |
831 else | |
832 retval(3) = ZZ; | |
833 } | |
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
834 OCTAVE_FALLTHROUGH; |
3185 | 835 |
836 case 3: | |
837 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
|
838 { |
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
839 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
|
840 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
|
841 else |
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
842 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
|
843 } |
3185 | 844 else |
10311 | 845 { |
846 if (complex_case) | |
847 retval(2) = CQ.hermitian (); | |
848 else | |
849 retval(2) = QQ.transpose (); | |
850 } | |
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
851 OCTAVE_FALLTHROUGH; |
10311 | 852 |
3185 | 853 case 2: |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
854 { |
10311 | 855 if (complex_case) |
856 { | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
857 #if defined (DEBUG) |
14844
5bc9b9cb4362
maint: Use Octave coding conventions for cuddled parenthesis in retval assignments.
Rik <octave@nomad.inbox5.com>
parents:
14501
diff
changeset
|
858 std::cout << "qz: retval(1) = cbb = " << std::endl; |
10311 | 859 octave_print_internal (std::cout, cbb, 0); |
860 std::cout << std::endl << "qz: retval(0) = caa = " <<std::endl; | |
861 octave_print_internal (std::cout, caa, 0); | |
862 std::cout << std::endl; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
863 #endif |
10311 | 864 retval(1) = cbb; |
865 retval(0) = caa; | |
866 } | |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
867 else |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
868 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
869 #if defined (DEBUG) |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
870 std::cout << "qz: retval(1) = bb = " << std::endl; |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
871 octave_print_internal (std::cout, bb, 0); |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
872 std::cout << std::endl << "qz: retval(0) = aa = " <<std::endl; |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
873 octave_print_internal (std::cout, aa, 0); |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
874 std::cout << std::endl; |
3185 | 875 #endif |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
876 retval(1) = bb; |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
877 retval(0) = aa; |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
878 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
879 } |
10311 | 880 break; |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
881 |
3185 | 882 case 1: |
883 case 0: | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
884 #if defined (DEBUG) |
3538 | 885 std::cout << "qz: retval(0) = gev = " << gev << std::endl; |
3185 | 886 #endif |
887 retval(0) = gev; | |
888 break; | |
889 | |
890 default: | |
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
891 error ("qz: too many return arguments"); |
3185 | 892 break; |
25103
078b795c5219
maint: style check C++ ahead of 4.4 release.
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
893 } |
3183 | 894 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
895 #if defined (DEBUG) |
3538 | 896 std::cout << "qz: exiting (at long last)" << std::endl; |
3185 | 897 #endif |
3183 | 898 |
899 return retval; | |
900 } | |
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
901 |
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
902 /* |
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
903 %!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
|
904 %! 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
|
905 %! 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
|
906 %! 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
|
907 %!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
|
908 %!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
|
909 |
23430 | 910 ## 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
|
911 %!test |
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
912 %! a = [ 10 1 2; |
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
913 %! 1 2 -1; |
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
914 %! 1 1 2]; |
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
915 %! 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
|
916 %! [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
|
917 %! sz = length (lambda); |
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
918 %! 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
|
919 %! 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
|
920 %! 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
|
921 %! 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
|
922 %! 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
|
923 %! 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
|
924 |
14876
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
925 %!test |
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
926 %! 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
|
927 %! 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
|
928 %! [AA, BB, Q, Z1] = qz (A, B); |
23430 | 929 %! [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
|
930 %! 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
|
931 |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
932 %!test |
da2b38b7a077
Use LAPACK DTGSEN subrouting for ordered-qz computation (bug #53761)
Sébastien Villemot <sebastien@debian.org>
parents:
25106
diff
changeset
|
933 %! 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
|
934 %! 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
|
935 %! -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
|
936 %! -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
|
937 %! 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
|
938 %! 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
|
939 %! 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
|
940 %! 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
|
941 %! [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
|
942 %! 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
|
943 %! 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
|
944 %! [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
|
945 %! 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
|
946 %! 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
|
947 %! [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
|
948 %! 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
|
949 %! 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
|
950 %! [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
|
951 %! 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
|
952 %! 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
|
953 */ |