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