Mercurial > octave
annotate libinterp/corefcn/qz.cc @ 25106:d7ad543255c5 stable
doc: Shorten very long first sentences of docstrings (bug #53388).
* bsxfun.cc (Fbsxfun), daspk.cc (Fdaspk), dasrt.cc (Fdasrt), dassl.cc (Fdassl),
gsvd.cc (Fgsvd), load-save.cc (Foctave_core_file_limit), qz.cc (Fqz), svd.cc
(Fsvd), sylvester.cc (Fsylvester), utils.cc (Ferrno), bincoeff.m, bessel.m,
krylov.m, expint.m, moment.m: Shorten very long first sentences.
author | Rik <rik@octave.org> |
---|---|
date | Tue, 03 Apr 2018 21:37:29 -0700 |
parents | 078b795c5219 |
children | da2b38b7a077 |
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 |
25 // Author: A. S. Hodel <scotte@eng.auburn.edu> | |
3183 | 26 |
27 #undef DEBUG | |
28 #undef DEBUG_SORT | |
29 #undef DEBUG_EIG | |
30 | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
31 #if defined (HAVE_CONFIG_H) |
21301
40de9f8f23a6
Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents:
21279
diff
changeset
|
32 # include "config.h" |
9786
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
33 #endif |
3183 | 34 |
23430 | 35 #include <cctype> |
36 #include <cmath> | |
4051 | 37 |
23430 | 38 #if defined (DEBUG) || defined (DEBUG_SORT) || defined (DEBUG_EIG) |
39 # include <iostream> | |
40 # if defined (DEBUG_EIG) | |
41 # include <iomanip> | |
42 # endif | |
43 #endif | |
3183 | 44 |
4153 | 45 #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
|
46 #include "lo-lapack-proto.h" |
21279
eb1524b07fe3
better use of templates for qr classes
John W. Eaton <jwe@octave.org>
parents:
21200
diff
changeset
|
47 #include "qr.h" |
4153 | 48 #include "quit.h" |
49 | |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14876
diff
changeset
|
50 #include "defun.h" |
3183 | 51 #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
|
52 #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
|
53 #include "ovl.h" |
3185 | 54 #if defined (DEBUG) || defined (DEBUG_SORT) |
21200
fcac5dbbf9ed
maint: Indent #ifdef blocks in libinterp.
Rik <rik@octave.org>
parents:
21110
diff
changeset
|
55 # include "pr-output.h" |
3183 | 56 #endif |
57 | |
22964
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
58 typedef F77_INT (*sort_function) (const F77_INT& LSIZE, |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
59 const double& ALPHA, const double& BETA, |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
60 const double& S, const double& P); |
3183 | 61 |
62 extern "C" | |
63 { | |
23430 | 64 // Van Dooren's code (netlib.org: toms/590) for reordering GEP. |
65 // Only processes Z, not Q. | |
4552 | 66 F77_RET_T |
22964
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
67 F77_FUNC (dsubsp, DSUBSP) (const F77_INT& NMAX, const F77_INT& N, |
23449
c763214a8260
maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents:
23430
diff
changeset
|
68 F77_DBLE *A, F77_DBLE *B, F77_DBLE *Z, |
22964
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
69 sort_function, const F77_DBLE& EPS, |
23449
c763214a8260
maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents:
23430
diff
changeset
|
70 F77_INT& NDIM, F77_INT& FAIL, F77_INT *IND); |
3183 | 71 } |
72 | |
73 // fcrhp, fin, fout, folhp: | |
21562
6c2fd62db1f7
maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
74 // Routines for ordering of generalized eigenvalues. |
6c2fd62db1f7
maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
75 // Return 1 if test is passed, 0 otherwise. |
6c2fd62db1f7
maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
76 // fin: |lambda| < 1 |
6c2fd62db1f7
maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
77 // fout: |lambda| >= 1 |
6c2fd62db1f7
maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
78 // fcrhp: real(lambda) >= 0 |
6c2fd62db1f7
maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
79 // folhp: real(lambda) < 0 |
3183 | 80 |
22964
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
81 static F77_INT |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
82 fcrhp (const F77_INT& lsize, const double& alpha, const double& beta, |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
83 const double& s, const double&) |
3183 | 84 { |
3185 | 85 if (lsize == 1) |
10311 | 86 return (alpha * beta >= 0 ? 1 : -1); |
3185 | 87 else |
3183 | 88 return (s >= 0 ? 1 : -1); |
89 } | |
3185 | 90 |
22964
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
91 static F77_INT |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
92 fin (const F77_INT& lsize, const double& alpha, const double& beta, |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
93 const double&, const double& p) |
3183 | 94 { |
22964
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
95 F77_INT retval; |
3183 | 96 |
3185 | 97 if (lsize == 1) |
23430 | 98 retval = (std::abs (alpha) < std::abs (beta) ? 1 : -1); |
3185 | 99 else |
23430 | 100 retval = (std::abs (p) < 1 ? 1 : -1); |
3185 | 101 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
102 #if defined (DEBUG) |
3538 | 103 std::cout << "qz: fin: retval=" << retval << std::endl; |
3185 | 104 #endif |
105 | |
3183 | 106 return retval; |
107 } | |
3185 | 108 |
22964
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
109 static F77_INT |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
110 folhp (const F77_INT& lsize, const double& alpha, const double& beta, |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
111 const double& s, const double&) |
3183 | 112 { |
3185 | 113 if (lsize == 1) |
10311 | 114 return (alpha * beta < 0 ? 1 : -1); |
3185 | 115 else |
3183 | 116 return (s < 0 ? 1 : -1); |
117 } | |
3185 | 118 |
22964
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
119 static F77_INT |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
120 fout (const F77_INT& lsize, const double& alpha, const double& beta, |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
121 const double&, const double& p) |
3183 | 122 { |
3185 | 123 if (lsize == 1) |
23430 | 124 return (std::abs (alpha) >= std::abs (beta) ? 1 : -1); |
3185 | 125 else |
23430 | 126 return (std::abs (p) >= 1 ? 1 : -1); |
3183 | 127 } |
128 | |
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
129 //FIXME: Matlab does not produce lambda as the first output argument. |
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
130 // Compatibility problem? |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14876
diff
changeset
|
131 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
|
132 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
|
133 @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
|
134 @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
|
135 @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
|
136 @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
|
137 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
|
138 |
d7ad543255c5
doc: Shorten very long first sentences of docstrings (bug #53388).
Rik <rik@octave.org>
parents:
25103
diff
changeset
|
139 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
|
140 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
141 @tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
142 $$A x = \lambda B x$$ |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
143 @end tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
144 @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
|
145 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
146 @math{A x = @var{lambda} B x} |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
147 |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
148 @end ifnottex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
149 |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
150 There are three calling forms of the function: |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
151 |
21966
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
152 @enumerate |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
153 @item @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
|
154 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
155 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
|
156 @tex |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
157 $\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
|
158 @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
|
159 @ifnottex |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
160 @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
|
161 @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
|
162 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
163 @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
|
164 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
165 Compute QZ@tie{}decomposition, generalized eigenvectors, and generalized |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
166 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
|
167 @tex |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
168 $$ 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
|
169 $$ 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
|
170 $$ 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
|
171 @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
|
172 @ifnottex |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
173 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
174 @example |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
175 @group |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
176 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
177 @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
|
178 @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
|
179 @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
|
180 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
181 @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
|
182 @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
|
183 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
184 @end ifnottex |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
185 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
|
186 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
187 @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
|
188 |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
189 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
|
190 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
|
191 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
|
192 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
|
193 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
194 @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
|
195 @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
|
196 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
|
197 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
|
198 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
199 @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
|
200 @item @qcode{"N"} |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
201 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
|
202 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
203 @item @qcode{"S"} |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
204 small: leading block has all |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
205 @tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
206 $|\lambda| < 1$ |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
207 @end tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
208 @ifnottex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
209 |@var{lambda}| < 1 |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
210 @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
|
211 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
212 @item @qcode{"B"} |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
213 big: leading block has all |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
214 @tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
215 $|\lambda| \geq 1$ |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
216 @end tex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
217 @ifnottex |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
218 |@var{lambda}| @geq{} 1 |
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
219 @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
|
220 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
221 @item @qcode{"-"} |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
222 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
|
223 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
|
224 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
225 @item @qcode{"+"} |
23428
1bc0e610e293
doc: Redo docstring for qz (bug #50846).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
226 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
|
227 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
|
228 @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
|
229 @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
|
230 @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
|
231 |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
232 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
|
233 (@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
|
234 @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
|
235 @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
|
236 @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
|
237 @end deftypefn */) |
3183 | 238 { |
23430 | 239 int nargin = args.length (); |
3183 | 240 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
241 #if defined (DEBUG) |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
242 std::cout << "qz: nargin = " << nargin |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
243 << ", nargout = " << nargout << std::endl; |
3185 | 244 #endif |
3183 | 245 |
3185 | 246 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
|
247 print_usage (); |
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20700
diff
changeset
|
248 |
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20700
diff
changeset
|
249 if (nargin == 3 && (nargout < 3 || nargout > 4)) |
23430 | 250 error ("qz: invalid number of output arguments for form 3 call"); |
3183 | 251 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
252 #if defined (DEBUG) |
3538 | 253 std::cout << "qz: determine ordering option" << std::endl; |
3185 | 254 #endif |
3183 | 255 |
10311 | 256 // Determine ordering option. |
23430 | 257 // declared volatile to avoid compiler warnings about long jumps vforks. |
24116
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
258 volatile char ord_job = 'N'; |
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
259 double safmin = 0.0; |
3185 | 260 |
24116
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
261 if (nargin == 3) |
3185 | 262 { |
23430 | 263 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
|
264 |
24116
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
265 if (opt.empty ()) |
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
266 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
|
267 |
e1020353fff5
qz: initialize optional ordering flag safely
Mike Miller <mtmiller@octave.org>
parents:
23826
diff
changeset
|
268 ord_job = std::toupper (opt[0]); |
3183 | 269 |
23430 | 270 std::string valid_opts = "NSB+-"; |
271 | |
272 if (valid_opts.find_first_of (ord_job) == std::string::npos) | |
273 error ("qz: invalid order option '%c'", ord_job); | |
3185 | 274 |
275 // overflow constant required by dlag2 | |
4552 | 276 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
|
277 safmin |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
278 F77_CHAR_ARG_LEN (1)); |
3183 | 279 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
280 #if defined (DEBUG_EIG) |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
281 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
|
282 << setiosflags (std::ios::scientific) |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
283 << safmin << std::endl; |
3185 | 284 #endif |
3183 | 285 |
10311 | 286 // Some machines (e.g., DEC alpha) get safmin = 0; |
287 // for these, use eps instead to avoid problems in dlag2. | |
3185 | 288 if (safmin == 0) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
289 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
290 #if defined (DEBUG_EIG) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
291 std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl; |
3185 | 292 #endif |
3183 | 293 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
294 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
|
295 safmin |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
296 F77_CHAR_ARG_LEN (1)); |
3185 | 297 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
298 #if defined (DEBUG_EIG) |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
299 std::cout << "qz: safmin set to " |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
300 << setiosflags (std::ios::scientific) |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
301 << safmin << std::endl; |
3185 | 302 #endif |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
303 } |
3183 | 304 } |
305 | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
306 #if defined (DEBUG) |
23430 | 307 std::cout << "qz: check matrix A" << std::endl; |
3185 | 308 #endif |
309 | |
23430 | 310 // Matrix A: check dimensions. |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22964
diff
changeset
|
311 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
|
312 F77_INT nc = octave::to_f77_int (args(0).columns ()); |
20892 | 313 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
314 #if defined (DEBUG) |
23807
336f89b6208b
Use character literals 'c' rather than string literals "c" when possible.
Rik <rik@octave.org>
parents:
23795
diff
changeset
|
315 std::cout << "Matrix A dimensions: (" << nn << ',' << nc << ')' << std::endl; |
3185 | 316 #endif |
317 | |
23577
80c42f4cca13
maint: Deprecate is_empty and replace with isempty.
Rik <rik@octave.org>
parents:
23477
diff
changeset
|
318 if (args(0).isempty ()) |
3185 | 319 { |
23430 | 320 warn_empty_arg ("qz: A"); |
3185 | 321 return octave_value_list (2, Matrix ()); |
322 } | |
22964
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
323 else if (nc != nn) |
21110
3d0d84305600
Use err_square_matrix_required more widely.
Rik <rik@octave.org>
parents:
21100
diff
changeset
|
324 err_square_matrix_required ("qz", "A"); |
3183 | 325 |
23430 | 326 // Matrix A: get value. |
3183 | 327 Matrix aa; |
328 ComplexMatrix caa; | |
3185 | 329 |
23581
c3075ae020e1
maint: Deprecate is_complex_type and replace with iscomplex.
Rik <rik@octave.org>
parents:
23577
diff
changeset
|
330 if (args(0).iscomplex ()) |
3183 | 331 caa = args(0).complex_matrix_value (); |
3185 | 332 else |
3183 | 333 aa = args(0).matrix_value (); |
3185 | 334 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
335 #if defined (DEBUG) |
23430 | 336 std::cout << "qz: check matrix B" << std::endl; |
3185 | 337 #endif |
3183 | 338 |
10311 | 339 // 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
|
340 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
|
341 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
|
342 |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
343 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
|
344 err_nonconformant (); |
3185 | 345 |
3183 | 346 Matrix bb; |
347 ComplexMatrix cbb; | |
3185 | 348 |
23581
c3075ae020e1
maint: Deprecate is_complex_type and replace with iscomplex.
Rik <rik@octave.org>
parents:
23577
diff
changeset
|
349 if (args(1).iscomplex ()) |
3183 | 350 cbb = args(1).complex_matrix_value (); |
351 else | |
352 bb = args(1).matrix_value (); | |
3185 | 353 |
23430 | 354 // Both matrices loaded, now check whether to calculate complex or real val. |
3185 | 355 |
23430 | 356 // declared volatile to avoid compiler warnings about long jumps vforks. |
357 volatile bool complex_case | |
23581
c3075ae020e1
maint: Deprecate is_complex_type and replace with iscomplex.
Rik <rik@octave.org>
parents:
23577
diff
changeset
|
358 = (args(0).iscomplex () || args(1).iscomplex ()); |
10311 | 359 |
3185 | 360 if (nargin == 3 && complex_case) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
361 error ("qz: cannot re-order complex qz decomposition"); |
3183 | 362 |
23430 | 363 // First, declare variables used in both the real and complex cases. |
364 // FIXME: There are a lot of excess variables declared. | |
365 // Probably a better way to handle this. | |
366 Matrix QQ (nn,nn), ZZ (nn,nn), VR (nn,nn), VL (nn,nn); | |
367 RowVector alphar (nn), alphai (nn), betar (nn); | |
368 ComplexRowVector xalpha (nn), xbeta (nn); | |
369 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
|
370 F77_INT ilo, ihi, info; |
23430 | 371 char comp_q = (nargout >= 3 ? 'V' : 'N'); |
372 char comp_z = ((nargout >= 4 || nargin == 3)? 'V' : 'N'); | |
3183 | 373 |
23430 | 374 // Initialize Q, Z to identity matrix if either is needed |
375 if (comp_q == 'V' || comp_z == 'V') | |
376 { | |
23477
3530b956d707
maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents:
23457
diff
changeset
|
377 double *QQptr = QQ.fortran_vec (); |
3530b956d707
maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents:
23457
diff
changeset
|
378 double *ZZptr = ZZ.fortran_vec (); |
23430 | 379 std::fill_n (QQptr, QQ.numel (), 0.0); |
380 std::fill_n (ZZptr, ZZ.numel (), 0.0); | |
381 for (F77_INT i = 0; i < nn; i++) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
382 { |
23430 | 383 QQ(i,i) = 1.0; |
384 ZZ(i,i) = 1.0; | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
385 } |
23430 | 386 } |
3183 | 387 |
10311 | 388 // Always perform permutation balancing. |
4552 | 389 const char bal_job = 'P'; |
10311 | 390 RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn); |
3183 | 391 |
3185 | 392 if (complex_case) |
393 { | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
394 #if defined (DEBUG) |
23430 | 395 if (comp_q == 'V') |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
396 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
|
397 << CQ << std::endl; |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
398 #endif |
23582
0cc2011d800e
maint: Deprecate is_real_type and replace with isreal.
Rik <rik@octave.org>
parents:
23581
diff
changeset
|
399 if (args(0).isreal ()) |
10311 | 400 caa = ComplexMatrix (aa); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
401 |
23582
0cc2011d800e
maint: Deprecate is_real_type and replace with isreal.
Rik <rik@octave.org>
parents:
23581
diff
changeset
|
402 if (args(1).isreal ()) |
10311 | 403 cbb = ComplexMatrix (bb); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
404 |
23430 | 405 if (comp_q == 'V') |
10311 | 406 CQ = ComplexMatrix (QQ); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
407 |
23430 | 408 if (comp_z == 'V') |
10311 | 409 CZ = ComplexMatrix (ZZ); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
410 |
10311 | 411 F77_XFCN (zggbal, ZGGBAL, |
412 (F77_CONST_CHAR_ARG2 (&bal_job, 1), | |
22407
34ce5be04942
maint: Style check C++ code in libinterp/.
Rik <rik@octave.org>
parents:
22372
diff
changeset
|
413 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
|
414 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), |
10311 | 415 nn, ilo, ihi, lscale.fortran_vec (), |
416 rscale.fortran_vec (), work.fortran_vec (), info | |
417 F77_CHAR_ARG_LEN (1))); | |
3185 | 418 } |
3183 | 419 else |
3185 | 420 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
421 #if defined (DEBUG) |
23430 | 422 if (comp_q == 'V') |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
423 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
|
424 << QQ << std::endl; |
3185 | 425 #endif |
3183 | 426 |
3185 | 427 F77_XFCN (dggbal, DGGBAL, |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
428 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
429 nn, aa.fortran_vec (), nn, bb.fortran_vec (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
430 nn, ilo, ihi, lscale.fortran_vec (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
431 rscale.fortran_vec (), work.fortran_vec (), info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
432 F77_CHAR_ARG_LEN (1))); |
3185 | 433 } |
3183 | 434 |
23430 | 435 // Only permutation balance above is done. Skip scaling balance. |
436 | |
437 #if 0 | |
3183 | 438 // Since we just want the balancing matrices, we can use dggbal |
10311 | 439 // for both the real and complex cases; left first |
3185 | 440 |
23430 | 441 if (comp_q == 'V') |
3185 | 442 { |
443 F77_XFCN (dggbak, DGGBAK, | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
444 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
445 F77_CONST_CHAR_ARG2 ("L", 1), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
446 nn, ilo, ihi, lscale.data (), rscale.data (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
447 nn, QQ.fortran_vec (), nn, info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
448 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
449 F77_CHAR_ARG_LEN (1))); |
3183 | 450 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
451 #if defined (DEBUG) |
23430 | 452 if (comp_q == 'V') |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
453 std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl; |
3185 | 454 #endif |
25103
078b795c5219
maint: style check C++ ahead of 4.4 release.
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
455 } |
3183 | 456 |
457 // then right | |
23430 | 458 if (comp_z == 'V') |
3185 | 459 { |
4552 | 460 F77_XFCN (dggbak, DGGBAK, |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
461 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
462 F77_CONST_CHAR_ARG2 ("R", 1), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
463 nn, ilo, ihi, lscale.data (), rscale.data (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
464 nn, ZZ.fortran_vec (), nn, info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
465 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
466 F77_CHAR_ARG_LEN (1))); |
3183 | 467 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
468 #if defined (DEBUG) |
23430 | 469 if (comp_z == 'V') |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
470 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl; |
3185 | 471 #endif |
10311 | 472 } |
473 #endif | |
3183 | 474 |
23430 | 475 char qz_job = (nargout < 2 ? 'E' : 'S'); |
3185 | 476 |
3183 | 477 if (complex_case) |
3185 | 478 { |
10311 | 479 // Complex case. |
480 | |
481 // 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
|
482 octave::math::qr<ComplexMatrix> cbqr (cbb); |
10311 | 483 // The R matrix of QR decomposition for cbb. |
484 cbb = cbqr.R (); | |
485 // (Q*)caa for following work. | |
486 caa = (cbqr.Q ().hermitian ()) * caa; | |
487 CQ = CQ * cbqr.Q (); | |
488 | |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
489 F77_XFCN (zgghrd, ZGGHRD, |
23430 | 490 (F77_CONST_CHAR_ARG2 (&comp_q, 1), |
491 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
|
492 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
|
493 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
|
494 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
|
495 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
|
496 F77_CHAR_ARG_LEN (1) |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
497 F77_CHAR_ARG_LEN (1))); |
10311 | 498 |
23430 | 499 ComplexRowVector cwork (nn); |
10311 | 500 |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
501 F77_XFCN (zhgeqz, ZHGEQZ, |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
502 (F77_CONST_CHAR_ARG2 (&qz_job, 1), |
23430 | 503 F77_CONST_CHAR_ARG2 (&comp_q, 1), |
504 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
|
505 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
|
506 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
|
507 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
|
508 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
|
509 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
|
510 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
|
511 F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn, |
23430 | 512 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
|
513 rwork.fortran_vec (), info |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
514 F77_CHAR_ARG_LEN (1) |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
515 F77_CHAR_ARG_LEN (1) |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
516 F77_CHAR_ARG_LEN (1))); |
3185 | 517 |
23430 | 518 if (comp_q == 'V') |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
519 { |
10311 | 520 // Left eigenvector. |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
521 F77_XFCN (zggbak, ZGGBAK, |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
522 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
523 F77_CONST_CHAR_ARG2 ("L", 1), |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
524 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
|
525 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
|
526 F77_CHAR_ARG_LEN (1) |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
527 F77_CHAR_ARG_LEN (1))); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
528 } |
3185 | 529 |
23430 | 530 if (comp_z == 'V') |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
531 { |
23430 | 532 // Right eigenvector. |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
533 F77_XFCN (zggbak, ZGGBAK, |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
534 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
535 F77_CONST_CHAR_ARG2 ("R", 1), |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
536 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
|
537 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
|
538 F77_CHAR_ARG_LEN (1) |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
539 F77_CHAR_ARG_LEN (1))); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
540 } |
3185 | 541 |
542 } | |
10311 | 543 else |
3185 | 544 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
545 #if defined (DEBUG) |
23430 | 546 std::cout << "qz: performing qr decomposition of bb" << std::endl; |
3185 | 547 #endif |
3183 | 548 |
10311 | 549 // 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
|
550 octave::math::qr<Matrix> bqr (bb); |
3185 | 551 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
552 #if defined (DEBUG) |
23430 | 553 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
|
554 << std::endl; |
3185 | 555 #endif |
3183 | 556 |
3185 | 557 bb = bqr.R (); |
558 | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
559 #if defined (DEBUG) |
3538 | 560 std::cout << "qz: extracted bb" << std::endl; |
3185 | 561 #endif |
3183 | 562 |
10311 | 563 aa = (bqr.Q ()).transpose () * aa; |
3185 | 564 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
565 #if defined (DEBUG) |
3538 | 566 std::cout << "qz: updated aa " << std::endl; |
567 std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl; | |
3183 | 568 |
23430 | 569 if (comp_q == 'V') |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
570 std::cout << "QQ =" << QQ << std::endl; |
3185 | 571 #endif |
3183 | 572 |
23430 | 573 if (comp_q == 'V') |
10311 | 574 QQ = QQ * bqr.Q (); |
3183 | 575 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
576 #if defined (DEBUG) |
3538 | 577 std::cout << "qz: precursors done..." << std::endl; |
3185 | 578 #endif |
3183 | 579 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
580 #if defined (DEBUG) |
23430 | 581 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
|
582 << std::endl; |
3185 | 583 #endif |
3183 | 584 |
21562
6c2fd62db1f7
maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
585 // Reduce to generalized Hessenberg form. |
3185 | 586 F77_XFCN (dgghrd, DGGHRD, |
23430 | 587 (F77_CONST_CHAR_ARG2 (&comp_q, 1), |
588 F77_CONST_CHAR_ARG2 (&comp_z, 1), | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
589 nn, ilo, ihi, aa.fortran_vec (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
590 nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
591 ZZ.fortran_vec (), nn, info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
592 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
593 F77_CHAR_ARG_LEN (1))); |
3183 | 594 |
23430 | 595 // Check if just computing generalized eigenvalues, |
596 // or if we're actually computing the decomposition. | |
3183 | 597 |
10311 | 598 // Reduce to generalized Schur form. |
3185 | 599 F77_XFCN (dhgeqz, DHGEQZ, |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
600 (F77_CONST_CHAR_ARG2 (&qz_job, 1), |
23430 | 601 F77_CONST_CHAR_ARG2 (&comp_q, 1), |
602 F77_CONST_CHAR_ARG2 (&comp_z, 1), | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
603 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
|
604 nn, alphar.fortran_vec (), alphai.fortran_vec (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
605 betar.fortran_vec (), QQ.fortran_vec (), nn, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
606 ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
607 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
608 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
609 F77_CHAR_ARG_LEN (1))); |
10311 | 610 |
23430 | 611 if (comp_q == 'V') |
10311 | 612 { |
613 F77_XFCN (dggbak, DGGBAK, | |
614 (F77_CONST_CHAR_ARG2 (&bal_job, 1), | |
615 F77_CONST_CHAR_ARG2 ("L", 1), | |
616 nn, ilo, ihi, lscale.data (), rscale.data (), | |
617 nn, QQ.fortran_vec (), nn, info | |
618 F77_CHAR_ARG_LEN (1) | |
619 F77_CHAR_ARG_LEN (1))); | |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
620 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
621 #if defined (DEBUG) |
23430 | 622 if (comp_q == 'V') |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
623 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
|
624 << QQ << std::endl; |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
625 #endif |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
626 } |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
627 |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
628 // then right |
23430 | 629 if (comp_z == 'V') |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
630 { |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
631 F77_XFCN (dggbak, DGGBAK, |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
632 (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
|
633 F77_CONST_CHAR_ARG2 ("R", 1), |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
634 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
|
635 nn, ZZ.fortran_vec (), nn, info |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
636 F77_CHAR_ARG_LEN (1) |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
637 F77_CHAR_ARG_LEN (1))); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
638 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
639 #if defined (DEBUG) |
23430 | 640 if (comp_z == 'V') |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
641 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
|
642 << ZZ << std::endl; |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
643 #endif |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
644 } |
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
645 |
3185 | 646 } |
3183 | 647 |
10311 | 648 // Order the QZ decomposition? |
23430 | 649 if (ord_job != 'N') |
3183 | 650 { |
3185 | 651 if (complex_case) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
652 // Probably not needed, but better be safe. |
23430 | 653 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
|
654 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
655 #if defined (DEBUG_SORT) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
656 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
|
657 << ord_job << std::endl; |
3185 | 658 #endif |
3183 | 659 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
660 // Declared static to avoid vfork/long jump compiler complaints. |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
661 static sort_function sort_test; |
23795
980f39c3ab90
Use C++11 nullptr rather than 0 in code (bug #51565).
Rik <rik@octave.org>
parents:
23662
diff
changeset
|
662 sort_test = nullptr; |
3183 | 663 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
664 switch (ord_job) |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
665 { |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
666 case 'S': |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
667 sort_test = &fin; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
668 break; |
3183 | 669 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
670 case 'B': |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
671 sort_test = &fout; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
672 break; |
3183 | 673 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
674 case '+': |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
675 sort_test = &fcrhp; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
676 break; |
3183 | 677 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
678 case '-': |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
679 sort_test = &folhp; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
680 break; |
3185 | 681 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
682 default: |
23430 | 683 // Invalid order option |
684 // (should never happen since options were checked at the top). | |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
685 panic_impossible (); |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
686 break; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
687 } |
3183 | 688 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
689 double inf_norm; |
3185 | 690 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
691 F77_XFCN (xdlange, XDLANGE, |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
692 (F77_CONST_CHAR_ARG2 ("I", 1), |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
693 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
694 F77_CHAR_ARG_LEN (1))); |
3185 | 695 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
696 double eps = std::numeric_limits<double>::epsilon () * inf_norm * nn; |
3185 | 697 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
698 #if defined (DEBUG_SORT) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
699 std::cout << "qz: calling dsubsp: aa=" << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
700 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
|
701 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
|
702 octave_print_internal (std::cout, bb, 0); |
23430 | 703 if (comp_z == 'V') |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
704 { |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
705 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
|
706 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
|
707 } |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
708 std::cout << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
709 std::cout << "alphar = " << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
710 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
|
711 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
|
712 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
|
713 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
|
714 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
|
715 std::cout << std::endl; |
3185 | 716 #endif |
717 | |
22964
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
718 Array<F77_INT> ind (dim_vector (nn, 1)); |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
719 |
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
720 F77_INT ndim, fail; |
3550 | 721 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
722 F77_XFCN (dsubsp, DSUBSP, |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
723 (nn, nn, aa.fortran_vec (), bb.fortran_vec (), |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
724 ZZ.fortran_vec (), sort_test, eps, ndim, fail, |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
725 ind.fortran_vec ())); |
3185 | 726 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
727 #if defined (DEBUG) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
728 std::cout << "qz: back from dsubsp: aa=" << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
729 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
|
730 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
|
731 octave_print_internal (std::cout, bb, 0); |
23430 | 732 if (comp_z == 'V') |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
733 { |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
734 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
|
735 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
|
736 } |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
737 std::cout << std::endl; |
3185 | 738 #endif |
739 | |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
740 // Manually update alphar, alphai, betar. |
23430 | 741 static F77_INT j; |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
742 |
23430 | 743 j = 0; |
744 while (j < nn) | |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
745 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
746 #if defined (DEBUG_EIG) |
23430 | 747 std::cout << "computing gen eig #" << j << std::endl; |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
748 #endif |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
749 |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
750 // Number of zeros in this block. |
22964
0c12642be005
use F77_INT instead of octave_idx_type for libinterp QZ function
John W. Eaton <jwe@octave.org>
parents:
22860
diff
changeset
|
751 static F77_INT zcnt; |
3185 | 752 |
23430 | 753 if (j == (nn-1)) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
754 zcnt = 1; |
23430 | 755 else if (aa(j+1,j) == 0) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
756 zcnt = 1; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
757 else zcnt = 2; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
758 |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
759 if (zcnt == 1) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
760 { |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
761 // Real zero. |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
762 #if defined (DEBUG_EIG) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
763 std::cout << " single gen eig:" << std::endl; |
23430 | 764 std::cout << " alphar(" << j << ") = " << aa(j,j) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
765 << std::endl; |
23430 | 766 std::cout << " betar(" << j << ") = " << bb(j,j) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
767 << std::endl; |
23430 | 768 std::cout << " alphai(" << j << ") = 0" << std::endl; |
3185 | 769 #endif |
770 | |
23430 | 771 alphar(j) = aa(j,j); |
772 alphai(j) = 0; | |
773 betar(j) = bb(j,j); | |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
774 } |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
775 else |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
776 { |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
777 // Complex conjugate pair. |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
778 #if defined (DEBUG_EIG) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
779 std::cout << "qz: calling dlag2:" << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
780 std::cout << "safmin=" |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
781 << setiosflags (std::ios::scientific) |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
782 << safmin << std::endl; |
3185 | 783 |
23430 | 784 for (F77_INT idr = j; idr <= j+1; idr++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
785 { |
23430 | 786 for (F77_INT idc = j; idc <= j+1; idc++) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
787 { |
23807
336f89b6208b
Use character literals 'c' rather than string literals "c" when possible.
Rik <rik@octave.org>
parents:
23795
diff
changeset
|
788 std::cout << "aa(" << idr << ',' << idc << ")=" |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
789 << aa(idr,idc) << std::endl; |
23807
336f89b6208b
Use character literals 'c' rather than string literals "c" when possible.
Rik <rik@octave.org>
parents:
23795
diff
changeset
|
790 std::cout << "bb(" << idr << ',' << idc << ")=" |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
791 << bb(idr,idc) << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
792 } |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
793 } |
3185 | 794 #endif |
795 | |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
796 // FIXME: probably should be using |
23430 | 797 // fortran_vec instead of &aa(j,j) here. |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
798 |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
799 double scale1, scale2, wr1, wr2, wi; |
23430 | 800 const double *aa_ptr = aa.data () + j * nn + j; |
801 const double *bb_ptr = bb.data () + j * nn + j; | |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
802 F77_XFCN (dlag2, DLAG2, |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
803 (aa_ptr, nn, bb_ptr, nn, safmin, |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
804 scale1, scale2, wr1, wr2, wi)); |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
805 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
806 #if defined (DEBUG_EIG) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
807 std::cout << "dlag2 returns: scale1=" << scale1 |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
808 << "\tscale2=" << scale2 << std::endl |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
809 << "\twr1=" << wr1 << "\twr2=" << wr2 |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
810 << "\twi=" << wi << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
811 #endif |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
812 |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
813 // Just to be safe, check if it's a real pair. |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
814 if (wi == 0) |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
815 { |
23430 | 816 alphar(j) = wr1; |
817 alphai(j) = 0; | |
818 betar(j) = scale1; | |
819 alphar(j+1) = wr2; | |
820 alphai(j+1) = 0; | |
821 betar(j+1) = scale2; | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
822 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
823 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
824 { |
23430 | 825 alphar(j) = alphar(j+1) = wr1; |
826 alphai(j) = -(alphai(j+1) = wi); | |
827 betar(j) = betar(j+1) = scale1; | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
828 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
829 } |
3185 | 830 |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
831 // Advance past this block. |
23430 | 832 j += zcnt; |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
833 } |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
834 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
835 #if defined (DEBUG_SORT) |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
836 std::cout << "qz: back from dsubsp: aa=" << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
837 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
|
838 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
|
839 octave_print_internal (std::cout, bb, 0); |
3185 | 840 |
23430 | 841 if (comp_z == 'V') |
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
842 { |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
843 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
|
844 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
|
845 } |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
846 std::cout << std::endl << "qz: ndim=" << ndim << std::endl |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
847 << "fail=" << fail << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
848 std::cout << "alphar = " << std::endl; |
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
849 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
|
850 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
|
851 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
|
852 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
|
853 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
|
854 std::cout << std::endl; |
3185 | 855 #endif |
3183 | 856 } |
3185 | 857 |
23430 | 858 // Compute the generalized eigenvalues as well? |
3183 | 859 ComplexColumnVector gev; |
3185 | 860 |
861 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) | |
3183 | 862 { |
3185 | 863 if (complex_case) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
864 { |
23430 | 865 ComplexColumnVector tmp (nn); |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
866 |
23430 | 867 for (F77_INT i = 0; i < nn; i++) |
868 tmp(i) = xalpha(i) / xbeta(i); | |
10311 | 869 |
870 gev = tmp; | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
871 } |
3185 | 872 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
873 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
874 #if defined (DEBUG) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
875 std::cout << "qz: computing generalized eigenvalues" << std::endl; |
3185 | 876 #endif |
3183 | 877 |
10311 | 878 // Return finite generalized eigenvalues. |
23430 | 879 ComplexColumnVector tmp (nn); |
3185 | 880 |
23430 | 881 F77_INT cnt = 0; |
882 for (F77_INT i = 0; i < nn; i++) | |
883 if (betar(i) != 0) | |
884 tmp(cnt++) = Complex (alphar(i), alphai(i)) / betar(i); | |
3185 | 885 |
23430 | 886 tmp.resize (cnt); // Trim vector to number of return values |
10311 | 887 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
888 gev = tmp; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
889 } |
3183 | 890 } |
891 | |
10311 | 892 // Right, left eigenvector matrices. |
3185 | 893 if (nargout >= 5) |
3183 | 894 { |
10311 | 895 // Which side to compute? |
896 char side = (nargout == 5 ? 'R' : 'B'); | |
897 // Compute all of them and backtransform | |
23430 | 898 char howmany = 'B'; |
10311 | 899 // 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
|
900 F77_INT *select = nullptr; |
3185 | 901 |
902 if (complex_case) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
903 { |
10311 | 904 CVL = CQ; |
905 CVR = CZ; | |
906 ComplexRowVector cwork2 (2 * nn); | |
907 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
|
908 F77_INT m; |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
909 |
10311 | 910 F77_XFCN (ztgevc, ZTGEVC, |
911 (F77_CONST_CHAR_ARG2 (&side, 1), | |
23430 | 912 F77_CONST_CHAR_ARG2 (&howmany, 1), |
22407
34ce5be04942
maint: Style check C++ code in libinterp/.
Rik <rik@octave.org>
parents:
22372
diff
changeset
|
913 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
|
914 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), |
34ce5be04942
maint: Style check C++ code in libinterp/.
Rik <rik@octave.org>
parents:
22372
diff
changeset
|
915 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
|
916 F77_DBLE_CMPLX_ARG (CVR.fortran_vec ()), nn, nn, |
23430 | 917 m, F77_DBLE_CMPLX_ARG (cwork2.fortran_vec ()), |
918 rwork2.fortran_vec (), info | |
10311 | 919 F77_CHAR_ARG_LEN (1) |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
920 F77_CHAR_ARG_LEN (1))); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
921 } |
3185 | 922 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
923 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
924 #if defined (DEBUG) |
21562
6c2fd62db1f7
maint: Eliminate accidental double spaces in code.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
925 std::cout << "qz: computing generalized eigenvectors" << std::endl; |
3185 | 926 #endif |
927 | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
928 VL = QQ; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
929 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
|
930 F77_INT m; |
10311 | 931 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
932 F77_XFCN (dtgevc, DTGEVC, |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
933 (F77_CONST_CHAR_ARG2 (&side, 1), |
23430 | 934 F77_CONST_CHAR_ARG2 (&howmany, 1), |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
935 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (), |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
936 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
|
937 m, work.fortran_vec (), info |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
938 F77_CHAR_ARG_LEN (1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
939 F77_CHAR_ARG_LEN (1))); |
3185 | 940 |
10311 | 941 // Now construct the complex form of VV, WW. |
23430 | 942 F77_INT j = 0; |
3185 | 943 |
23430 | 944 while (j < nn) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
945 { |
22860
0b1e25cc4457
eliminate use of OCTAVE_QUIT macro in C++ sources
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
946 octave_quit (); |
4153 | 947 |
10311 | 948 // See if real or complex eigenvalue. |
949 | |
950 // Column increment; assume complex eigenvalue. | |
951 int cinc = 2; | |
3185 | 952 |
23430 | 953 if (j == (nn-1)) |
10311 | 954 // Single column. |
955 cinc = 1; | |
23430 | 956 else if (aa(j+1,j) == 0) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
957 cinc = 1; |
3185 | 958 |
10311 | 959 // Now copy the eigenvector (s) to CVR, CVL. |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
960 if (cinc == 1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
961 { |
23430 | 962 for (F77_INT i = 0; i < nn; i++) |
963 CVR(i,j) = VR(i,j); | |
3185 | 964 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
965 if (side == 'B') |
23430 | 966 for (F77_INT i = 0; i < nn; i++) |
967 CVL(i,j) = VL(i,j); | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
968 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
969 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
970 { |
10311 | 971 // Double column; complex vector. |
3185 | 972 |
23430 | 973 for (F77_INT i = 0; i < nn; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
974 { |
23430 | 975 CVR(i,j) = Complex (VR(i,j), VR(i,j+1)); |
976 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
|
977 } |
3183 | 978 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
979 if (side == 'B') |
23430 | 980 for (F77_INT i = 0; i < nn; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
981 { |
23430 | 982 CVL(i,j) = Complex (VL(i,j), VL(i,j+1)); |
983 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
|
984 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
985 } |
3185 | 986 |
10311 | 987 // Advance to next eigenvectors (if any). |
23430 | 988 j += cinc; |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
989 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
990 } |
3183 | 991 } |
3185 | 992 |
23430 | 993 octave_value_list retval (nargout); |
994 | |
3185 | 995 switch (nargout) |
996 { | |
997 case 7: | |
998 retval(6) = gev; | |
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
999 OCTAVE_FALLTHROUGH; |
3185 | 1000 |
10311 | 1001 case 6: |
23430 | 1002 // Return left eigenvectors. |
3185 | 1003 retval(5) = CVL; |
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
1004 OCTAVE_FALLTHROUGH; |
3185 | 1005 |
10311 | 1006 case 5: |
23430 | 1007 // Return right eigenvectors. |
3185 | 1008 retval(4) = CVR; |
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
1009 OCTAVE_FALLTHROUGH; |
3185 | 1010 |
1011 case 4: | |
1012 if (nargin == 3) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1013 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
1014 #if defined (DEBUG) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1015 std::cout << "qz: sort: retval(3) = gev = " << std::endl; |
23430 | 1016 octave_print_internal (std::cout, ComplexMatrix (gev)); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1017 std::cout << std::endl; |
3185 | 1018 #endif |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1019 retval(3) = gev; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1020 } |
3185 | 1021 else |
10311 | 1022 { |
1023 if (complex_case) | |
1024 retval(3) = CZ; | |
1025 else | |
1026 retval(3) = ZZ; | |
1027 } | |
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
1028 OCTAVE_FALLTHROUGH; |
3185 | 1029 |
1030 case 3: | |
1031 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
|
1032 { |
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1033 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
|
1034 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
|
1035 else |
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1036 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
|
1037 } |
3185 | 1038 else |
10311 | 1039 { |
1040 if (complex_case) | |
1041 retval(2) = CQ.hermitian (); | |
1042 else | |
1043 retval(2) = QQ.transpose (); | |
1044 } | |
23826
d69021d58a61
avoid fallthrough warnings
John W. Eaton <jwe@octave.org>
parents:
23807
diff
changeset
|
1045 OCTAVE_FALLTHROUGH; |
10311 | 1046 |
3185 | 1047 case 2: |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
1048 { |
10311 | 1049 if (complex_case) |
1050 { | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
1051 #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
|
1052 std::cout << "qz: retval(1) = cbb = " << std::endl; |
10311 | 1053 octave_print_internal (std::cout, cbb, 0); |
1054 std::cout << std::endl << "qz: retval(0) = caa = " <<std::endl; | |
1055 octave_print_internal (std::cout, caa, 0); | |
1056 std::cout << std::endl; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
1057 #endif |
10311 | 1058 retval(1) = cbb; |
1059 retval(0) = caa; | |
1060 } | |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1061 else |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1062 { |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
1063 #if defined (DEBUG) |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1064 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
|
1065 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
|
1066 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
|
1067 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
|
1068 std::cout << std::endl; |
3185 | 1069 #endif |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1070 retval(1) = bb; |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1071 retval(0) = aa; |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1072 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
1073 } |
10311 | 1074 break; |
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
1075 |
3185 | 1076 case 1: |
1077 case 0: | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
1078 #if defined (DEBUG) |
3538 | 1079 std::cout << "qz: retval(0) = gev = " << gev << std::endl; |
3185 | 1080 #endif |
1081 retval(0) = gev; | |
1082 break; | |
1083 | |
1084 default: | |
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
1085 error ("qz: too many return arguments"); |
3185 | 1086 break; |
25103
078b795c5219
maint: style check C++ ahead of 4.4 release.
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
1087 } |
3183 | 1088 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21562
diff
changeset
|
1089 #if defined (DEBUG) |
3538 | 1090 std::cout << "qz: exiting (at long last)" << std::endl; |
3185 | 1091 #endif |
3183 | 1092 |
1093 return retval; | |
1094 } | |
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1095 |
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1096 /* |
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1097 %!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
|
1098 %! 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
|
1099 %! 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
|
1100 %! 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
|
1101 %!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
|
1102 %!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
|
1103 |
23430 | 1104 ## 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
|
1105 %!test |
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1106 %! a = [ 10 1 2; |
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1107 %! 1 2 -1; |
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1108 %! 1 1 2]; |
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
1109 %! 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
|
1110 %! [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
|
1111 %! sz = length (lambda); |
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
1112 %! 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
|
1113 %! 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
|
1114 %! 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
|
1115 %! 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
|
1116 %! 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
|
1117 %! 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
|
1118 |
14876
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1119 %!test |
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1120 %! 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
|
1121 %! 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
|
1122 %! [AA, BB, Q, Z1] = qz (A, B); |
23430 | 1123 %! [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
|
1124 %! assert (Z1, Z2); |
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1125 */ |