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