Mercurial > octave
annotate src/DLD-FUNCTIONS/qz.cc @ 9758:09da0bd91412
Periodic grammar check of Octave documentation files to ensure common format
Use double spaces after a period at start of new sentence
Do not use @iftex blocks where @tex block is sufficient
author | Rik <rdrider0-list@yahoo.com> |
---|---|
date | Fri, 23 Oct 2009 09:28:14 -0700 |
parents | 923c7cb7f13f |
children | 2c279308f6ab |
rev | line source |
---|---|
3183 | 1 /* |
2 | |
8920 | 3 Copyright (C) 1998, 1999, 2000, 2002, 2003, 2004, 2005, 2006, 2007, |
4 2008, 2009 A. S. Hodel | |
3183 | 5 |
6 This file is part of Octave. | |
7 | |
8 Octave is free software; you can redistribute it and/or modify it | |
9 under the terms of the GNU General Public License as published by the | |
7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
11 option) any later version. | |
3183 | 12 |
13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
16 for more details. | |
17 | |
18 You should have received a copy of the GNU General Public License | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
3183 | 21 |
22 */ | |
23 | |
24 // Generalized eigenvalue balancing via LAPACK | |
3911 | 25 |
26 // Author: A. S. Hodel <scotte@eng.auburn.edu> | |
3183 | 27 |
28 #undef DEBUG | |
29 #undef DEBUG_SORT | |
30 #undef DEBUG_EIG | |
31 | |
32 #include "config.h" | |
33 | |
34 #include <cfloat> | |
4051 | 35 |
3523 | 36 #include <iostream> |
4051 | 37 #include <iomanip> |
3183 | 38 |
39 #include "CmplxQRP.h" | |
40 #include "dbleQR.h" | |
4153 | 41 #include "f77-fcn.h" |
7231 | 42 #include "lo-math.h" |
4153 | 43 #include "quit.h" |
44 | |
3183 | 45 #include "defun-dld.h" |
46 #include "error.h" | |
47 #include "gripes.h" | |
48 #include "oct-obj.h" | |
49 #include "oct-map.h" | |
50 #include "ov.h" | |
51 #include "pager.h" | |
3185 | 52 #if defined (DEBUG) || defined (DEBUG_SORT) |
3183 | 53 #include "pr-output.h" |
54 #endif | |
55 #include "symtab.h" | |
56 #include "utils.h" | |
57 #include "variables.h" | |
58 | |
5275 | 59 typedef octave_idx_type (*sort_function) (const octave_idx_type& LSIZE, const double& ALPHA, |
3185 | 60 const double& BETA, const double& S, |
61 const double& P); | |
3183 | 62 |
63 extern "C" | |
64 { | |
4552 | 65 F77_RET_T |
66 F77_FUNC (dggbal, DGGBAL) (F77_CONST_CHAR_ARG_DECL, | |
5275 | 67 const octave_idx_type& N, double* A, const octave_idx_type& LDA, |
68 double* B, const octave_idx_type& LDB, octave_idx_type& ILO, | |
69 octave_idx_type& IHI, double* LSCALE, double* RSCALE, | |
70 double* WORK, octave_idx_type& INFO | |
4552 | 71 F77_CHAR_ARG_LEN_DECL); |
3183 | 72 |
4552 | 73 F77_RET_T |
74 F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL, | |
75 F77_CONST_CHAR_ARG_DECL, | |
5275 | 76 const octave_idx_type& N, const octave_idx_type& ILO, |
77 const octave_idx_type& IHI, const double* LSCALE, | |
78 const double* RSCALE, octave_idx_type& M, double* V, | |
79 const octave_idx_type& LDV, octave_idx_type& INFO | |
4552 | 80 F77_CHAR_ARG_LEN_DECL |
81 F77_CHAR_ARG_LEN_DECL); | |
3183 | 82 |
4552 | 83 F77_RET_T |
84 F77_FUNC (dgghrd, DGGHRD) (F77_CONST_CHAR_ARG_DECL, | |
85 F77_CONST_CHAR_ARG_DECL, | |
5275 | 86 const octave_idx_type& N, const octave_idx_type& ILO, |
87 const octave_idx_type& IHI, double* A, | |
88 const octave_idx_type& LDA, double* B, | |
89 const octave_idx_type& LDB, double* Q, | |
90 const octave_idx_type& LDQ, double* Z, | |
91 const octave_idx_type& LDZ, octave_idx_type& INFO | |
4552 | 92 F77_CHAR_ARG_LEN_DECL |
93 F77_CHAR_ARG_LEN_DECL); | |
3183 | 94 |
4552 | 95 F77_RET_T |
96 F77_FUNC (dhgeqz, DHGEQZ) (F77_CONST_CHAR_ARG_DECL, | |
97 F77_CONST_CHAR_ARG_DECL, | |
98 F77_CONST_CHAR_ARG_DECL, | |
5275 | 99 const octave_idx_type& N, const octave_idx_type& ILO, const octave_idx_type& IHI, |
100 double* A, const octave_idx_type& LDA, double* B, | |
101 const octave_idx_type& LDB, double* ALPHAR, | |
4552 | 102 double* ALPHAI, double* BETA, double* Q, |
5275 | 103 const octave_idx_type& LDQ, double* Z, |
104 const octave_idx_type& LDZ, double* WORK, | |
105 const octave_idx_type& LWORK, octave_idx_type& INFO | |
4552 | 106 F77_CHAR_ARG_LEN_DECL |
107 F77_CHAR_ARG_LEN_DECL | |
108 F77_CHAR_ARG_LEN_DECL); | |
3183 | 109 |
4552 | 110 F77_RET_T |
5275 | 111 F77_FUNC (dlag2, DLAG2) (const double* A, const octave_idx_type& LDA, const double* B, |
112 const octave_idx_type& LDB, const double& SAFMIN, | |
4552 | 113 double& SCALE1, double& SCALE2, |
114 double& WR1, double& WR2, double& WI); | |
3183 | 115 |
116 // Van Dooren's code (netlib.org: toms/590) for reordering | |
117 // GEP. Only processes Z, not Q. | |
4552 | 118 F77_RET_T |
5275 | 119 F77_FUNC (dsubsp, DSUBSP) (const octave_idx_type& NMAX, const octave_idx_type& N, double* A, |
4552 | 120 double* B, double* Z, sort_function, |
5275 | 121 const double& EPS, octave_idx_type& NDIM, octave_idx_type& FAIL, |
122 octave_idx_type* IND); | |
3183 | 123 |
124 // documentation for DTGEVC incorrectly states that VR, VL are | |
125 // complex*16; they are declared in DTGEVC as double precision | |
126 // (probably a cut and paste problem fro ZTGEVC) | |
4552 | 127 F77_RET_T |
128 F77_FUNC (dtgevc, DTGEVC) (F77_CONST_CHAR_ARG_DECL, | |
129 F77_CONST_CHAR_ARG_DECL, | |
5275 | 130 octave_idx_type* SELECT, const octave_idx_type& N, double* A, |
131 const octave_idx_type& LDA, double* B, | |
132 const octave_idx_type& LDB, double* VL, | |
133 const octave_idx_type& LDVL, double* VR, | |
134 const octave_idx_type& LDVR, const octave_idx_type& MM, | |
135 octave_idx_type& M, double* WORK, octave_idx_type& INFO | |
4552 | 136 F77_CHAR_ARG_LEN_DECL |
137 F77_CHAR_ARG_LEN_DECL); | |
3183 | 138 |
4552 | 139 F77_RET_T |
140 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG_DECL, | |
141 double& retval | |
142 F77_CHAR_ARG_LEN_DECL); | |
3185 | 143 |
4552 | 144 F77_RET_T |
145 F77_FUNC (xdlange, XDLANGE) (F77_CONST_CHAR_ARG_DECL, | |
5275 | 146 const octave_idx_type&, const octave_idx_type&, const double*, |
147 const octave_idx_type&, double*, double& | |
4552 | 148 F77_CHAR_ARG_LEN_DECL); |
3183 | 149 } |
150 | |
151 // fcrhp, fin, fout, folhp: | |
152 // routines for ordering of generalized eigenvalues | |
153 // return 1 if test is passed, 0 otherwise | |
154 // fin: |lambda| < 1 | |
155 // fout: |lambda| >= 1 | |
156 // fcrhp: real(lambda) >= 0 | |
157 // folhp: real(lambda) < 0 | |
158 | |
5275 | 159 static octave_idx_type |
160 fcrhp (const octave_idx_type& lsize, const double& alpha, | |
3185 | 161 const double& beta, const double& s, const double&) |
3183 | 162 { |
3185 | 163 if (lsize == 1) |
3183 | 164 return (alpha*beta >= 0 ? 1 : -1); |
3185 | 165 else |
3183 | 166 return (s >= 0 ? 1 : -1); |
167 } | |
3185 | 168 |
5275 | 169 static octave_idx_type |
170 fin (const octave_idx_type& lsize, const double& alpha, | |
3185 | 171 const double& beta, const double&, const double& p) |
3183 | 172 { |
5275 | 173 octave_idx_type retval; |
3183 | 174 |
3185 | 175 if (lsize == 1) |
176 retval = (fabs (alpha) < fabs (beta) ? 1 : -1); | |
177 else | |
178 retval = (fabs (p) < 1 ? 1 : -1); | |
179 | |
180 #ifdef DEBUG | |
3538 | 181 std::cout << "qz: fin: retval=" << retval << std::endl; |
3185 | 182 #endif |
183 | |
3183 | 184 return retval; |
185 } | |
3185 | 186 |
5275 | 187 static octave_idx_type |
188 folhp (const octave_idx_type& lsize, const double& alpha, | |
3185 | 189 const double& beta, const double& s, const double&) |
3183 | 190 { |
3185 | 191 if (lsize == 1) |
3183 | 192 return (alpha*beta < 0 ? 1 : -1); |
3185 | 193 else |
3183 | 194 return (s < 0 ? 1 : -1); |
195 } | |
3185 | 196 |
5275 | 197 static octave_idx_type |
198 fout (const octave_idx_type& lsize, const double& alpha, | |
3185 | 199 const double& beta, const double&, const double& p) |
3183 | 200 { |
3185 | 201 if (lsize == 1) |
202 return (fabs (alpha) >= fabs (beta) ? 1 : -1); | |
203 else | |
204 return (fabs (p) >= 1 ? 1 : -1); | |
3183 | 205 } |
206 | |
207 DEFUN_DLD (qz, args, nargout, | |
3372 | 208 "-*- texinfo -*-\n\ |
209 @deftypefn {Loadable Function} {@var{lambda} =} qz (@var{a}, @var{b})\n\ | |
210 Generalized eigenvalue problem @math{A x = s B x},\n\ | |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8927
diff
changeset
|
211 @var{QZ} decomposition. There are three ways to call this function:\n\ |
3372 | 212 @enumerate\n\ |
213 @item @code{lambda = qz(A,B)}\n\ | |
3185 | 214 \n\ |
5016 | 215 Computes the generalized eigenvalues\n\ |
216 @tex\n\ | |
217 $\\lambda$\n\ | |
218 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
219 @ifnottex\n\ |
5016 | 220 @var{lambda}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
221 @end ifnottex\n\ |
5016 | 222 of @math{(A - s B)}.\n\ |
3500 | 223 @item @code{[AA, BB, Q, Z, V, W, lambda] = qz (A, B)}\n\ |
3185 | 224 \n\ |
3372 | 225 Computes qz decomposition, generalized eigenvectors, and \n\ |
5481 | 226 generalized eigenvalues of @math{(A - sB)}\n\ |
5016 | 227 @tex\n\ |
228 $$ AV = BV{ \\rm diag }(\\lambda) $$\n\ | |
229 $$ W^T A = { \\rm diag }(\\lambda)W^T B $$\n\ | |
230 $$ AA = Q^T AZ, BB = Q^T BZ $$\n\ | |
231 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
232 @ifnottex\n\ |
3372 | 233 @example\n\ |
234 @group\n\ | |
5481 | 235 \n\ |
236 A*V = B*V*diag(lambda)\n\ | |
237 W'*A = diag(lambda)*W'*B\n\ | |
238 AA = Q'*A*Z, BB = Q'*B*Z\n\ | |
239 \n\ | |
3372 | 240 @end group\n\ |
241 @end example\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
242 @end ifnottex\n\ |
5016 | 243 with @var{Q} and @var{Z} orthogonal (unitary)= @var{I}\n\ |
3185 | 244 \n\ |
5016 | 245 @item @code{[AA,BB,Z@{, lambda@}] = qz(A,B,opt)}\n\ |
3185 | 246 \n\ |
3372 | 247 As in form [2], but allows ordering of generalized eigenpairs\n\ |
5481 | 248 for (e.g.) solution of discrete time algebraic Riccati equations.\n\ |
249 Form 3 is not available for complex matrices, and does not compute\n\ | |
250 the generalized eigenvectors @var{V}, @var{W}, nor the orthogonal matrix @var{Q}.\n\ | |
3372 | 251 @table @var\n\ |
252 @item opt\n\ | |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8927
diff
changeset
|
253 for ordering eigenvalues of the GEP pencil. The leading block\n\ |
5481 | 254 of the revised pencil contains all eigenvalues that satisfy:\n\ |
3372 | 255 @table @code\n\ |
256 @item \"N\"\n\ | |
5481 | 257 = unordered (default) \n\ |
3183 | 258 \n\ |
3372 | 259 @item \"S\"\n\ |
260 = small: leading block has all |lambda| <=1 \n\ | |
3185 | 261 \n\ |
3372 | 262 @item \"B\"\n\ |
5481 | 263 = big: leading block has all |lambda| >= 1 \n\ |
3372 | 264 \n\ |
265 @item \"-\"\n\ | |
5481 | 266 = negative real part: leading block has all eigenvalues\n\ |
267 in the open left half-plane\n\ | |
3372 | 268 \n\ |
269 @item \"+\"\n\ | |
9209
923c7cb7f13f
Simplify TeXinfo files by eliminating redundant @iftex followed by @tex construction.
Rik <rdrider0-list@yahoo.com>
parents:
9064
diff
changeset
|
270 = non-negative real part: leading block has all eigenvalues\n\ |
5481 | 271 in the closed right half-plane\n\ |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8927
diff
changeset
|
272 @end table\n\ |
3372 | 273 @end table\n\ |
274 @end enumerate\n\ | |
3183 | 275 \n\ |
3372 | 276 Note: qz performs permutation balancing, but not scaling (see balance).\n\ |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8927
diff
changeset
|
277 Order of output arguments was selected for compatibility with @sc{matlab}\n\ |
3183 | 278 \n\ |
8286
6f2d95255911
fix @seealso references to point to existing anchors
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
7520
diff
changeset
|
279 @seealso{balance, eig, schur}\n\ |
3372 | 280 @end deftypefn") |
3183 | 281 { |
282 octave_value_list retval; | |
283 int nargin = args.length (); | |
284 | |
3185 | 285 #ifdef DEBUG |
3538 | 286 std::cout << "qz: nargin = " << nargin << ", nargout = " << nargout << std::endl; |
3185 | 287 #endif |
3183 | 288 |
3185 | 289 if (nargin < 2 || nargin > 3 || nargout > 7) |
290 { | |
5823 | 291 print_usage (); |
3185 | 292 return retval; |
293 } | |
294 else if (nargin == 3 && (nargout < 3 || nargout > 4)) | |
295 { | |
3427 | 296 error ("qz: invalid number of output arguments for form [3] call"); |
3185 | 297 return retval; |
298 } | |
3183 | 299 |
3185 | 300 #ifdef DEBUG |
3538 | 301 std::cout << "qz: determine ordering option" << std::endl; |
3185 | 302 #endif |
3183 | 303 |
304 // Determine ordering option | |
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
305 volatile char ord_job = 0; |
3183 | 306 static double safmin; |
3185 | 307 |
308 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
|
309 ord_job = 'N'; |
3185 | 310 else if (!args(2).is_string ()) |
311 { | |
312 error ("qz: argument 3 must be a string"); | |
313 return retval; | |
314 } | |
315 else | |
316 { | |
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
317 std::string tmp = args(2).string_value (); |
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
318 |
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
319 if (! tmp.empty ()) |
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
320 ord_job = tmp[0]; |
3183 | 321 |
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
322 if (! (ord_job == 'N' || ord_job == 'n' |
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
323 || ord_job == 'S' || ord_job == 's' |
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
324 || ord_job == 'B' || ord_job == 'b' |
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
325 || ord_job == '+' || ord_job == '-')) |
3185 | 326 { |
3427 | 327 error ("qz: invalid order option"); |
3185 | 328 return retval; |
329 } | |
330 | |
331 // overflow constant required by dlag2 | |
4552 | 332 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1), |
333 safmin | |
334 F77_CHAR_ARG_LEN (1)); | |
3183 | 335 |
3185 | 336 #ifdef DEBUG_EIG |
3538 | 337 std::cout << "qz: initial value of safmin=" << setiosflags (std::ios::scientific) |
338 << safmin << std::endl; | |
3185 | 339 #endif |
3183 | 340 |
3185 | 341 // some machines (e.g., DEC alpha) get safmin = 0; |
342 // for these, use eps instead to avoid problems in dlag2 | |
343 if (safmin == 0) | |
344 { | |
345 #ifdef DEBUG_EIG | |
3538 | 346 std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl; |
3185 | 347 #endif |
3183 | 348 |
4552 | 349 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("E", 1), |
350 safmin | |
351 F77_CHAR_ARG_LEN (1)); | |
3185 | 352 |
353 #ifdef DEBUG_EIG | |
3538 | 354 std::cout << "qz: safmin set to " << setiosflags (std::ios::scientific) |
355 << safmin << std::endl; | |
3185 | 356 #endif |
357 } | |
3183 | 358 } |
359 | |
3185 | 360 #ifdef DEBUG |
3538 | 361 std::cout << "qz: check argument 1" << std::endl; |
3185 | 362 #endif |
3183 | 363 |
364 // Argument 1: check if it's o.k. dimensioned | |
5275 | 365 octave_idx_type nn = args(0).rows (); |
3185 | 366 |
367 #ifdef DEBUG | |
3531 | 368 std::cout << "argument 1 dimensions: (" << nn << "," << args(0).columns () << ")" |
3538 | 369 << std::endl; |
3185 | 370 #endif |
371 | |
372 int arg_is_empty = empty_arg ("qz", nn, args(0).columns ()); | |
373 | |
3183 | 374 if (arg_is_empty < 0) |
3185 | 375 { |
376 gripe_empty_arg ("qz: parameter 1", 0); | |
377 return retval; | |
378 } | |
3183 | 379 else if (arg_is_empty > 0) |
3185 | 380 { |
381 gripe_empty_arg ("qz: parameter 1; continuing", 0); | |
382 return octave_value_list (2, Matrix ()); | |
383 } | |
384 else if (args(0).columns () != nn) | |
385 { | |
386 gripe_square_matrix_required ("qz"); | |
387 return retval; | |
388 } | |
3183 | 389 |
390 // Argument 1: dimensions look good; get the value | |
391 Matrix aa; | |
392 ComplexMatrix caa; | |
3185 | 393 |
394 if (args(0).is_complex_type ()) | |
3183 | 395 caa = args(0).complex_matrix_value (); |
3185 | 396 else |
3183 | 397 aa = args(0).matrix_value (); |
3185 | 398 |
399 if (error_state) | |
3183 | 400 return retval; |
401 | |
3185 | 402 #ifdef DEBUG |
3538 | 403 std::cout << "qz: check argument 2" << std::endl; |
3185 | 404 #endif |
3183 | 405 |
406 // Extract argument 2 (bb, or cbb if complex) | |
3185 | 407 if ((nn != args(1).columns ()) || (nn != args(1).rows ())) |
408 { | |
409 gripe_nonconformant (); | |
410 return retval; | |
411 } | |
412 | |
3183 | 413 Matrix bb; |
414 ComplexMatrix cbb; | |
3185 | 415 |
416 if (args(1).is_complex_type ()) | |
3183 | 417 cbb = args(1).complex_matrix_value (); |
418 else | |
419 bb = args(1).matrix_value (); | |
3185 | 420 |
421 if (error_state) | |
3183 | 422 return retval; |
423 | |
424 // Both matrices loaded, now let's check what kind of arithmetic: | |
425 //declared static to avoid compiler warnings about long jumps, vforks. | |
3185 | 426 |
427 static int complex_case | |
428 = (args(0).is_complex_type () || args(1).is_complex_type ()); | |
3183 | 429 |
3185 | 430 if (nargin == 3 && complex_case) |
431 { | |
432 error ("qz: cannot re-order complex qz decomposition."); | |
433 return retval; | |
434 } | |
3183 | 435 |
436 // first, declare variables used in both the real and complex case | |
437 Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn); | |
438 RowVector alphar(nn), alphai(nn), betar(nn); | |
439 | |
3185 | 440 ComplexMatrix CQ(nn,nn), CZ(nn,nn), CVR(nn,nn), CVL(nn,nn); |
5275 | 441 octave_idx_type ilo, ihi, info; |
3185 | 442 char compq = (nargout >= 3 ? 'V' : 'N'); |
443 char compz = (nargout >= 4 ? 'V' : 'N'); | |
3183 | 444 |
3185 | 445 // initialize Q, Z to identity if we need either of them |
446 if (compq == 'V' || compz == 'V') | |
5275 | 447 for (octave_idx_type ii = 0; ii < nn; ii++) |
448 for (octave_idx_type jj = 0; jj < nn; jj++) | |
4153 | 449 { |
450 OCTAVE_QUIT; | |
451 QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0); | |
452 } | |
3183 | 453 |
3185 | 454 // always perform permutation balancing |
4552 | 455 const char bal_job = 'P'; |
3183 | 456 RowVector lscale(nn), rscale(nn), work(6*nn); |
457 | |
3185 | 458 if (complex_case) |
459 { | |
460 error ("Complex case not implemented yet"); | |
461 return retval; | |
462 } | |
3183 | 463 else |
3185 | 464 { |
465 #ifdef DEBUG | |
466 if (compq == 'V') | |
3538 | 467 std::cout << "qz: performing balancing; QQ=" << std::endl << QQ << std::endl; |
3185 | 468 #endif |
3183 | 469 |
3185 | 470 F77_XFCN (dggbal, DGGBAL, |
4552 | 471 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
472 nn, aa.fortran_vec (), nn, bb.fortran_vec (), | |
473 nn, ilo, ihi, lscale.fortran_vec (), | |
474 rscale.fortran_vec (), work.fortran_vec (), info | |
475 F77_CHAR_ARG_LEN (1))); | |
3185 | 476 } |
3183 | 477 |
478 // Since we just want the balancing matrices, we can use dggbal | |
479 // for both the real and complex cases; | |
480 // left first | |
3185 | 481 |
482 if (compq == 'V') | |
483 { | |
484 F77_XFCN (dggbak, DGGBAK, | |
4552 | 485 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
486 F77_CONST_CHAR_ARG2 ("L", 1), | |
4566 | 487 nn, ilo, ihi, lscale.data (), rscale.data (), |
488 nn, QQ.fortran_vec (), nn, info | |
4552 | 489 F77_CHAR_ARG_LEN (1) |
490 F77_CHAR_ARG_LEN (1))); | |
3183 | 491 |
3185 | 492 #ifdef DEBUG |
493 if (compq == 'V') | |
3538 | 494 std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl; |
3185 | 495 #endif |
3183 | 496 } |
497 | |
498 // then right | |
3185 | 499 if (compz == 'V') |
500 { | |
4552 | 501 F77_XFCN (dggbak, DGGBAK, |
502 (F77_CONST_CHAR_ARG2 (&bal_job, 1), | |
503 F77_CONST_CHAR_ARG2 ("R", 1), | |
4566 | 504 nn, ilo, ihi, lscale.data (), rscale.data (), |
505 nn, ZZ.fortran_vec (), nn, info | |
4552 | 506 F77_CHAR_ARG_LEN (1) |
507 F77_CHAR_ARG_LEN (1))); | |
3183 | 508 |
3185 | 509 #ifdef DEBUG |
510 if (compz == 'V') | |
3538 | 511 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl; |
3185 | 512 #endif |
513 } | |
3183 | 514 |
515 static char qz_job; | |
516 qz_job = (nargout < 2 ? 'E' : 'S'); | |
3185 | 517 |
3183 | 518 if (complex_case) |
3185 | 519 { |
520 // complex case | |
521 if (args(0).is_real_type ()) | |
3587 | 522 caa = ComplexMatrix (aa); |
3185 | 523 |
524 if (args(1).is_real_type ()) | |
3587 | 525 cbb = ComplexMatrix (bb); |
3185 | 526 |
527 if (compq == 'V') | |
3587 | 528 CQ = ComplexMatrix (QQ); |
3185 | 529 |
530 if (compz == 'V') | |
3587 | 531 CZ = ComplexMatrix (ZZ); |
3185 | 532 |
533 error ("complex case not done yet"); | |
534 return retval; | |
535 } | |
3183 | 536 else // real matrices case |
3185 | 537 { |
538 #ifdef DEBUG | |
3538 | 539 std::cout << "qz: peforming qr decomposition of bb" << std::endl; |
3185 | 540 #endif |
3183 | 541 |
3185 | 542 // compute the QR factorization of bb |
543 QR bqr (bb); | |
544 | |
545 #ifdef DEBUG | |
3538 | 546 std::cout << "qz: qr (bb) done; now peforming qz decomposition" << std::endl; |
3185 | 547 #endif |
3183 | 548 |
3185 | 549 bb = bqr.R (); |
550 | |
551 #ifdef DEBUG | |
3538 | 552 std::cout << "qz: extracted bb" << std::endl; |
3185 | 553 #endif |
3183 | 554 |
3185 | 555 aa = (bqr.Q ()).transpose ()*aa; |
556 | |
557 #ifdef DEBUG | |
3538 | 558 std::cout << "qz: updated aa " << std::endl; |
559 std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl; | |
3183 | 560 |
3185 | 561 if (compq == 'V') |
3538 | 562 std::cout << "QQ =" << QQ << std::endl; |
3185 | 563 #endif |
3183 | 564 |
3185 | 565 if (compq == 'V') |
566 QQ = QQ*bqr.Q (); | |
3183 | 567 |
3185 | 568 #ifdef DEBUG |
3538 | 569 std::cout << "qz: precursors done..." << std::endl; |
3185 | 570 #endif |
3183 | 571 |
3185 | 572 #ifdef DEBUG |
3538 | 573 std::cout << "qz: compq = " << compq << ", compz = " << compz << std::endl; |
3185 | 574 #endif |
3183 | 575 |
3185 | 576 // reduce to generalized hessenberg form |
577 F77_XFCN (dgghrd, DGGHRD, | |
4552 | 578 (F77_CONST_CHAR_ARG2 (&compq, 1), |
579 F77_CONST_CHAR_ARG2 (&compz, 1), | |
580 nn, ilo, ihi, aa.fortran_vec (), | |
581 nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn, | |
582 ZZ.fortran_vec (), nn, info | |
583 F77_CHAR_ARG_LEN (1) | |
584 F77_CHAR_ARG_LEN (1))); | |
3183 | 585 |
3185 | 586 // check if just computing generalized eigenvalues or if we're |
587 // actually computing the decomposition | |
3183 | 588 |
3185 | 589 // reduce to generalized Schur form |
590 F77_XFCN (dhgeqz, DHGEQZ, | |
4552 | 591 (F77_CONST_CHAR_ARG2 (&qz_job, 1), |
592 F77_CONST_CHAR_ARG2 (&compq, 1), | |
593 F77_CONST_CHAR_ARG2 (&compz, 1), | |
594 nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (), | |
595 nn, alphar.fortran_vec (), alphai.fortran_vec (), | |
596 betar.fortran_vec (), QQ.fortran_vec (), nn, | |
597 ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info | |
598 F77_CHAR_ARG_LEN (1) | |
599 F77_CHAR_ARG_LEN (1) | |
600 F77_CHAR_ARG_LEN (1))); | |
3185 | 601 } |
3183 | 602 |
603 // order the QZ decomposition? | |
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
604 if (! (ord_job == 'N' || ord_job == 'n')) |
3183 | 605 { |
3185 | 606 if (complex_case) |
607 { | |
608 // probably not needed, but better be safe | |
609 error ("qz: cannot re-order complex qz decomposition."); | |
610 return retval; | |
611 } | |
612 else | |
613 { | |
614 #ifdef DEBUG_SORT | |
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
615 std::cout << "qz: ordering eigenvalues: ord_job = " |
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
616 << ord_job << std::endl; |
3185 | 617 #endif |
3183 | 618 |
3185 | 619 // declared static to avoid vfork/long jump compiler complaints |
620 static sort_function sort_test; | |
7520 | 621 sort_test = 0; |
3183 | 622 |
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
623 switch (ord_job) |
3185 | 624 { |
625 case 'S': | |
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
626 case 's': |
3185 | 627 sort_test = &fin; |
628 break; | |
3183 | 629 |
3185 | 630 case 'B': |
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
631 case 'b': |
3185 | 632 sort_test = &fout; |
633 break; | |
3183 | 634 |
3185 | 635 case '+': |
636 sort_test = &fcrhp; | |
637 break; | |
3183 | 638 |
3185 | 639 case '-': |
640 sort_test = &folhp; | |
641 break; | |
642 | |
643 default: | |
3427 | 644 // invalid order option (should never happen, since we |
3185 | 645 // checked the options at the top). |
646 panic_impossible (); | |
647 break; | |
3550 | 648 } |
3183 | 649 |
5275 | 650 octave_idx_type ndim, fail; |
3185 | 651 double inf_norm; |
652 | |
653 F77_XFCN (xdlange, XDLANGE, | |
4552 | 654 (F77_CONST_CHAR_ARG2 ("I", 1), |
4566 | 655 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm |
4552 | 656 F77_CHAR_ARG_LEN (1))); |
3185 | 657 |
658 double eps = DBL_EPSILON*inf_norm*nn; | |
659 | |
660 #ifdef DEBUG_SORT | |
3538 | 661 std::cout << "qz: calling dsubsp: aa=" << std::endl; |
3531 | 662 octave_print_internal (std::cout, aa, 0); |
3538 | 663 std::cout << std::endl << "bb=" << std::endl; |
3531 | 664 octave_print_internal (std::cout, bb, 0); |
3185 | 665 if (compz == 'V') |
666 { | |
3538 | 667 std::cout << std::endl << "ZZ=" << std::endl; |
3531 | 668 octave_print_internal (std::cout, ZZ, 0); |
3185 | 669 } |
3538 | 670 std::cout << std::endl; |
671 std::cout << "alphar = " << std::endl; | |
3531 | 672 octave_print_internal (std::cout, (Matrix) alphar, 0); |
3538 | 673 std::cout << std::endl << "alphai = " << std::endl; |
3531 | 674 octave_print_internal (std::cout, (Matrix) alphai, 0); |
3538 | 675 std::cout << std::endl << "beta = " << std::endl; |
3531 | 676 octave_print_internal (std::cout, (Matrix) betar, 0); |
3538 | 677 std::cout << std::endl; |
3185 | 678 #endif |
679 | |
5275 | 680 Array<octave_idx_type> ind (nn); |
3550 | 681 |
3185 | 682 F77_XFCN (dsubsp, DSUBSP, |
4552 | 683 (nn, nn, aa.fortran_vec (), bb.fortran_vec (), |
684 ZZ.fortran_vec (), sort_test, eps, ndim, fail, | |
3550 | 685 ind.fortran_vec ())); |
3185 | 686 |
687 #ifdef DEBUG | |
3538 | 688 std::cout << "qz: back from dsubsp: aa=" << std::endl; |
3531 | 689 octave_print_internal (std::cout, aa, 0); |
3538 | 690 std::cout << std::endl << "bb=" << std::endl; |
3531 | 691 octave_print_internal (std::cout, bb, 0); |
3185 | 692 if (compz == 'V') |
693 { | |
3538 | 694 std::cout << std::endl << "ZZ=" << std::endl; |
3531 | 695 octave_print_internal (std::cout, ZZ, 0); |
3185 | 696 } |
3538 | 697 std::cout << std::endl; |
3185 | 698 #endif |
699 | |
700 // manually update alphar, alphai, betar | |
701 static int jj; | |
702 | |
703 jj=0; | |
704 while (jj < nn) | |
705 { | |
706 #ifdef DEBUG_EIG | |
3538 | 707 std::cout << "computing gen eig #" << jj << std::endl; |
3185 | 708 #endif |
709 | |
710 static int zcnt; // number of zeros in this block | |
711 | |
712 if (jj == (nn-1)) | |
713 zcnt = 1; | |
714 else if (aa(jj+1,jj) == 0) | |
715 zcnt = 1; | |
716 else zcnt = 2; | |
717 | |
718 if (zcnt == 1) // real zero | |
719 { | |
720 #ifdef DEBUG_EIG | |
3538 | 721 std::cout << " single gen eig:" << std::endl; |
722 std::cout << " alphar(" << jj << ") = " << aa(jj,jj) << std::endl; | |
723 std::cout << " betar( " << jj << ") = " << bb(jj,jj) << std::endl; | |
724 std::cout << " alphai(" << jj << ") = 0" << std::endl; | |
3185 | 725 #endif |
726 | |
727 alphar(jj) = aa(jj,jj); | |
728 alphai(jj) = 0; | |
729 betar(jj) = bb(jj,jj); | |
730 } | |
731 else | |
732 { | |
733 // complex conjugate pair | |
734 #ifdef DEBUG_EIG | |
3538 | 735 std::cout << "qz: calling dlag2:" << std::endl; |
3531 | 736 std::cout << "safmin=" |
3538 | 737 << setiosflags (std::ios::scientific) << safmin << std::endl; |
3185 | 738 |
739 for (int idr = jj; idr <= jj+1; idr++) | |
740 { | |
741 for (int idc = jj; idc <= jj+1; idc++) | |
742 { | |
3531 | 743 std::cout << "aa(" << idr << "," << idc << ")=" |
3538 | 744 << aa(idr,idc) << std::endl; |
3531 | 745 std::cout << "bb(" << idr << "," << idc << ")=" |
3538 | 746 << bb(idr,idc) << std::endl; |
3185 | 747 } |
748 } | |
749 #endif | |
750 | |
5775 | 751 // FIXME -- probably should be using |
4566 | 752 // fortran_vec instead of &aa(jj,jj) here. |
753 | |
3185 | 754 double scale1, scale2, wr1, wr2, wi; |
4629 | 755 const double *aa_ptr = aa.data () + jj*nn+jj; |
756 const double *bb_ptr = bb.data () + jj*nn+jj; | |
3185 | 757 F77_XFCN (dlag2, DLAG2, |
4629 | 758 (aa_ptr, nn, bb_ptr, nn, safmin, |
3185 | 759 scale1, scale2, wr1, wr2, wi)); |
760 | |
761 #ifdef DEBUG_EIG | |
3531 | 762 std::cout << "dlag2 returns: scale1=" << scale1 |
3538 | 763 << "\tscale2=" << scale2 << std::endl |
3185 | 764 << "\twr1=" << wr1 << "\twr2=" << wr2 |
3538 | 765 << "\twi=" << wi << std::endl; |
3185 | 766 #endif |
767 | |
768 // just to be safe, check if it's a real pair | |
769 if (wi == 0) | |
770 { | |
771 alphar(jj) = wr1; | |
772 alphai(jj) = 0; | |
773 betar(jj) = scale1; | |
774 alphar(jj+1) = wr2; | |
775 alphai(jj+1) = 0; | |
776 betar(jj+1) = scale2; | |
777 } | |
778 else | |
779 { | |
780 alphar(jj) = alphar(jj+1)=wr1; | |
781 alphai(jj) = -(alphai(jj+1) = wi); | |
782 betar(jj) = betar(jj+1) = scale1; | |
783 } | |
784 } | |
785 | |
786 // advance past this block | |
787 jj += zcnt; | |
788 } | |
789 | |
790 #ifdef DEBUG_SORT | |
3538 | 791 std::cout << "qz: back from dsubsp: aa=" << std::endl; |
3531 | 792 octave_print_internal (std::cout, aa, 0); |
3538 | 793 std::cout << std::endl << "bb=" << std::endl; |
3531 | 794 octave_print_internal (std::cout, bb, 0); |
3185 | 795 |
796 if (compz == 'V') | |
797 { | |
3538 | 798 std::cout << std::endl << "ZZ=" << std::endl; |
3531 | 799 octave_print_internal (std::cout, ZZ, 0); |
3185 | 800 } |
3538 | 801 std::cout << std::endl << "qz: ndim=" << ndim << std::endl |
802 << "fail=" << fail << std::endl; | |
803 std::cout << "alphar = " << std::endl; | |
3531 | 804 octave_print_internal (std::cout, (Matrix) alphar, 0); |
3538 | 805 std::cout << std::endl << "alphai = " << std::endl; |
3531 | 806 octave_print_internal (std::cout, (Matrix) alphai, 0); |
3538 | 807 std::cout << std::endl << "beta = " << std::endl; |
3531 | 808 octave_print_internal (std::cout, (Matrix) betar, 0); |
3538 | 809 std::cout << std::endl; |
3185 | 810 #endif |
811 } | |
3183 | 812 } |
3185 | 813 |
3183 | 814 // compute generalized eigenvalues? |
815 ComplexColumnVector gev; | |
3185 | 816 |
817 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) | |
3183 | 818 { |
3185 | 819 if (complex_case) |
820 { | |
821 error ("complex case not yet implemented"); | |
822 return retval; | |
823 } | |
824 else | |
825 { | |
826 #ifdef DEBUG | |
3538 | 827 std::cout << "qz: computing generalized eigenvalues" << std::endl; |
3185 | 828 #endif |
3183 | 829 |
3185 | 830 // return finite generalized eigenvalues |
831 int cnt = 0; | |
832 | |
833 for (int ii = 0; ii < nn; ii++) | |
834 if (betar(ii) != 0) | |
835 cnt++; | |
836 | |
837 ComplexColumnVector tmp(cnt); | |
838 | |
3671 | 839 cnt = 0; |
3185 | 840 for (int ii = 0; ii < nn; ii++) |
841 if (betar(ii) != 0) | |
3671 | 842 tmp(cnt++) = Complex(alphar(ii), alphai(ii))/betar(ii); |
3185 | 843 gev = tmp; |
844 } | |
3183 | 845 } |
846 | |
847 // right, left eigenvector matrices | |
3185 | 848 if (nargout >= 5) |
3183 | 849 { |
3185 | 850 char side = (nargout == 5 ? 'R' : 'B'); // which side to compute? |
851 char howmny = 'B'; // compute all of them and backtransform | |
7520 | 852 octave_idx_type *select = 0; // dummy pointer; select is not used. |
5275 | 853 octave_idx_type m; |
3185 | 854 |
855 if (complex_case) | |
856 { | |
857 error ("complex type not yet implemented"); | |
858 return retval; | |
859 } | |
860 else | |
861 { | |
862 #ifdef DEBUG | |
3538 | 863 std::cout << "qz: computing generalized eigenvectors" << std::endl; |
3185 | 864 #endif |
865 | |
866 VL = QQ; | |
867 VR = ZZ; | |
868 | |
869 F77_XFCN (dtgevc, DTGEVC, | |
4552 | 870 (F77_CONST_CHAR_ARG2 (&side, 1), |
871 F77_CONST_CHAR_ARG2 (&howmny, 1), | |
872 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (), | |
873 nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn, | |
874 m, work.fortran_vec (), info | |
875 F77_CHAR_ARG_LEN (1) | |
876 F77_CHAR_ARG_LEN (1))); | |
3185 | 877 |
878 // now construct the complex form of VV, WW | |
879 int jj = 0; | |
880 | |
881 while (jj < nn) | |
882 { | |
4153 | 883 OCTAVE_QUIT; |
884 | |
3185 | 885 // see if real or complex eigenvalue |
886 int cinc = 2; // column increment; assume complex eigenvalue | |
887 | |
888 if (jj == (nn-1)) | |
889 cinc = 1; // single column | |
890 else if (aa(jj+1,jj) == 0) | |
891 cinc = 1; | |
892 | |
893 // now copy the eigenvector (s) to CVR, CVL | |
894 if (cinc == 1) | |
895 { | |
896 for (int ii = 0; ii < nn; ii++) | |
897 CVR(ii,jj) = VR(ii,jj); | |
898 | |
899 if (side == 'B') | |
900 for (int ii = 0; ii < nn; ii++) | |
901 CVL(ii,jj) = VL(ii,jj); | |
902 } | |
903 else | |
904 { | |
905 // double column; complex vector | |
906 | |
907 for (int ii = 0; ii < nn; ii++) | |
908 { | |
909 CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1)); | |
910 CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1)); | |
911 } | |
3183 | 912 |
3185 | 913 if (side == 'B') |
914 for (int ii = 0; ii < nn; ii++) | |
915 { | |
916 CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1)); | |
917 CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1)); | |
918 } | |
919 } | |
920 | |
921 // advance to next eigenvectors (if any) | |
922 jj += cinc; | |
923 } | |
924 } | |
3183 | 925 } |
3185 | 926 |
927 switch (nargout) | |
928 { | |
929 case 7: | |
930 retval(6) = gev; | |
931 | |
932 case 6: // return eigenvectors | |
933 retval(5) = CVL; | |
934 | |
935 case 5: // return eigenvectors | |
936 retval(4) = CVR; | |
937 | |
938 case 4: | |
939 if (nargin == 3) | |
940 { | |
941 #ifdef DEBUG | |
3538 | 942 std::cout << "qz: sort: retval(3) = gev = " << std::endl; |
3531 | 943 octave_print_internal (std::cout, gev); |
3538 | 944 std::cout << std::endl; |
3185 | 945 #endif |
946 retval(3) = gev; | |
947 } | |
948 else | |
949 retval(3) = ZZ; | |
950 | |
951 case 3: | |
952 if (nargin == 3) | |
953 retval(2) = ZZ; | |
954 else | |
955 retval(2) = QQ; | |
956 | |
957 case 2: | |
958 #ifdef DEBUG | |
3538 | 959 std::cout << "qz: retval (1) = bb = " << std::endl; |
3531 | 960 octave_print_internal (std::cout, bb, 0); |
3538 | 961 std::cout << std::endl << "qz: retval(0) = aa = " <<std::endl; |
3531 | 962 octave_print_internal (std::cout, aa, 0); |
3538 | 963 std::cout << std::endl; |
3185 | 964 #endif |
965 retval(1) = bb; | |
966 retval(0) = aa; | |
967 break; | |
968 | |
969 case 1: | |
970 case 0: | |
971 #ifdef DEBUG | |
3538 | 972 std::cout << "qz: retval(0) = gev = " << gev << std::endl; |
3185 | 973 #endif |
974 retval(0) = gev; | |
975 break; | |
976 | |
977 default: | |
978 error ("qz: too many return arguments."); | |
979 break; | |
3183 | 980 } |
981 | |
3185 | 982 #ifdef DEBUG |
3538 | 983 std::cout << "qz: exiting (at long last)" << std::endl; |
3185 | 984 #endif |
3183 | 985 |
986 return retval; | |
987 } | |
988 | |
989 /* | |
990 ;;; Local Variables: *** | |
991 ;;; mode: C++ *** | |
992 ;;; End: *** | |
993 */ |