Mercurial > hg > octave-lyh
annotate src/DLD-FUNCTIONS/qz.cc @ 8517:81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
| author | sh@sh-laptop |
|---|---|
| date | Wed, 14 Jan 2009 20:44:25 -0500 |
| parents | 6f2d95255911 |
| children | eb63fbe60fab |
| rev | line source |
|---|---|
| 3183 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 1998, 1999, 2000, 2002, 2003, 2004, 2005, 2006, 2007 |
| 4 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\ | |
| 5016 | 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 @iftex\n\ | |
| 217 @tex\n\ | |
| 218 $\\lambda$\n\ | |
| 219 @end tex\n\ | |
| 220 @end iftex\n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
221 @ifnottex\n\ |
| 5016 | 222 @var{lambda}\n\ |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
223 @end ifnottex\n\ |
| 5016 | 224 of @math{(A - s B)}.\n\ |
| 3500 | 225 @item @code{[AA, BB, Q, Z, V, W, lambda] = qz (A, B)}\n\ |
| 3185 | 226 \n\ |
| 3372 | 227 Computes qz decomposition, generalized eigenvectors, and \n\ |
| 5481 | 228 generalized eigenvalues of @math{(A - sB)}\n\ |
| 5016 | 229 @iftex\n\ |
| 230 @tex\n\ | |
| 231 $$ AV = BV{ \\rm diag }(\\lambda) $$\n\ | |
| 232 $$ W^T A = { \\rm diag }(\\lambda)W^T B $$\n\ | |
| 233 $$ AA = Q^T AZ, BB = Q^T BZ $$\n\ | |
| 234 @end tex\n\ | |
| 235 @end iftex\n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
236 @ifnottex\n\ |
| 3372 | 237 @example\n\ |
| 238 @group\n\ | |
| 5481 | 239 \n\ |
| 240 A*V = B*V*diag(lambda)\n\ | |
| 241 W'*A = diag(lambda)*W'*B\n\ | |
| 242 AA = Q'*A*Z, BB = Q'*B*Z\n\ | |
| 243 \n\ | |
| 3372 | 244 @end group\n\ |
| 245 @end example\n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
246 @end ifnottex\n\ |
| 5016 | 247 with @var{Q} and @var{Z} orthogonal (unitary)= @var{I}\n\ |
| 3185 | 248 \n\ |
| 5016 | 249 @item @code{[AA,BB,Z@{, lambda@}] = qz(A,B,opt)}\n\ |
| 3185 | 250 \n\ |
| 3372 | 251 As in form [2], but allows ordering of generalized eigenpairs\n\ |
| 5481 | 252 for (e.g.) solution of discrete time algebraic Riccati equations.\n\ |
| 253 Form 3 is not available for complex matrices, and does not compute\n\ | |
| 254 the generalized eigenvectors @var{V}, @var{W}, nor the orthogonal matrix @var{Q}.\n\ | |
| 3372 | 255 @table @var\n\ |
| 256 @item opt\n\ | |
| 5481 | 257 for ordering eigenvalues of the GEP pencil. The leading block\n\ |
| 258 of the revised pencil contains all eigenvalues that satisfy:\n\ | |
| 3372 | 259 @table @code\n\ |
| 260 @item \"N\"\n\ | |
| 5481 | 261 = unordered (default) \n\ |
| 3183 | 262 \n\ |
| 3372 | 263 @item \"S\"\n\ |
| 264 = small: leading block has all |lambda| <=1 \n\ | |
| 3185 | 265 \n\ |
| 3372 | 266 @item \"B\"\n\ |
| 5481 | 267 = big: leading block has all |lambda| >= 1 \n\ |
| 3372 | 268 \n\ |
| 269 @item \"-\"\n\ | |
| 5481 | 270 = negative real part: leading block has all eigenvalues\n\ |
| 271 in the open left half-plane\n\ | |
| 3372 | 272 \n\ |
| 273 @item \"+\"\n\ | |
| 5481 | 274 = nonnegative real part: leading block has all eigenvalues\n\ |
| 275 in the closed right half-plane\n\ | |
| 3372 | 276 @end table\n\ |
| 277 @end table\n\ | |
| 278 @end enumerate\n\ | |
| 3183 | 279 \n\ |
| 3372 | 280 Note: qz performs permutation balancing, but not scaling (see balance).\n\ |
| 5481 | 281 Order of output arguments was selected for compatibility with MATLAB\n\ |
| 3183 | 282 \n\ |
|
8286
6f2d95255911
fix @seealso references to point to existing anchors
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
7520
diff
changeset
|
283 @seealso{balance, eig, schur}\n\ |
| 3372 | 284 @end deftypefn") |
| 3183 | 285 { |
| 286 octave_value_list retval; | |
| 287 int nargin = args.length (); | |
| 288 | |
| 3185 | 289 #ifdef DEBUG |
| 3538 | 290 std::cout << "qz: nargin = " << nargin << ", nargout = " << nargout << std::endl; |
| 3185 | 291 #endif |
| 3183 | 292 |
| 3185 | 293 if (nargin < 2 || nargin > 3 || nargout > 7) |
| 294 { | |
| 5823 | 295 print_usage (); |
| 3185 | 296 return retval; |
| 297 } | |
| 298 else if (nargin == 3 && (nargout < 3 || nargout > 4)) | |
| 299 { | |
| 3427 | 300 error ("qz: invalid number of output arguments for form [3] call"); |
| 3185 | 301 return retval; |
| 302 } | |
| 3183 | 303 |
| 3185 | 304 #ifdef DEBUG |
| 3538 | 305 std::cout << "qz: determine ordering option" << std::endl; |
| 3185 | 306 #endif |
| 3183 | 307 |
| 308 // Determine ordering option | |
| 3523 | 309 std::string ord_job; |
| 3183 | 310 static double safmin; |
| 3185 | 311 |
| 312 if (nargin == 2) | |
| 3183 | 313 ord_job = "N"; |
| 3185 | 314 else if (!args(2).is_string ()) |
| 315 { | |
| 316 error ("qz: argument 3 must be a string"); | |
| 317 return retval; | |
| 318 } | |
| 319 else | |
| 320 { | |
| 321 ord_job = args(2).string_value (); | |
| 3183 | 322 |
| 3185 | 323 if (ord_job[0] != 'N' |
| 324 && ord_job[0] != 'S' | |
| 325 && ord_job[0] != 'B' | |
| 326 && ord_job[0] != '+' | |
| 327 && ord_job[0] != '-') | |
| 328 { | |
| 3427 | 329 error ("qz: invalid order option"); |
| 3185 | 330 return retval; |
| 331 } | |
| 332 | |
| 333 // overflow constant required by dlag2 | |
| 4552 | 334 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1), |
| 335 safmin | |
| 336 F77_CHAR_ARG_LEN (1)); | |
| 3183 | 337 |
| 3185 | 338 #ifdef DEBUG_EIG |
| 3538 | 339 std::cout << "qz: initial value of safmin=" << setiosflags (std::ios::scientific) |
| 340 << safmin << std::endl; | |
| 3185 | 341 #endif |
| 3183 | 342 |
| 3185 | 343 // some machines (e.g., DEC alpha) get safmin = 0; |
| 344 // for these, use eps instead to avoid problems in dlag2 | |
| 345 if (safmin == 0) | |
| 346 { | |
| 347 #ifdef DEBUG_EIG | |
| 3538 | 348 std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl; |
| 3185 | 349 #endif |
| 3183 | 350 |
| 4552 | 351 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("E", 1), |
| 352 safmin | |
| 353 F77_CHAR_ARG_LEN (1)); | |
| 3185 | 354 |
| 355 #ifdef DEBUG_EIG | |
| 3538 | 356 std::cout << "qz: safmin set to " << setiosflags (std::ios::scientific) |
| 357 << safmin << std::endl; | |
| 3185 | 358 #endif |
| 359 } | |
| 3183 | 360 } |
| 361 | |
| 3185 | 362 #ifdef DEBUG |
| 3538 | 363 std::cout << "qz: check argument 1" << std::endl; |
| 3185 | 364 #endif |
| 3183 | 365 |
| 366 // Argument 1: check if it's o.k. dimensioned | |
| 5275 | 367 octave_idx_type nn = args(0).rows (); |
| 3185 | 368 |
| 369 #ifdef DEBUG | |
| 3531 | 370 std::cout << "argument 1 dimensions: (" << nn << "," << args(0).columns () << ")" |
| 3538 | 371 << std::endl; |
| 3185 | 372 #endif |
| 373 | |
| 374 int arg_is_empty = empty_arg ("qz", nn, args(0).columns ()); | |
| 375 | |
| 3183 | 376 if (arg_is_empty < 0) |
| 3185 | 377 { |
| 378 gripe_empty_arg ("qz: parameter 1", 0); | |
| 379 return retval; | |
| 380 } | |
| 3183 | 381 else if (arg_is_empty > 0) |
| 3185 | 382 { |
| 383 gripe_empty_arg ("qz: parameter 1; continuing", 0); | |
| 384 return octave_value_list (2, Matrix ()); | |
| 385 } | |
| 386 else if (args(0).columns () != nn) | |
| 387 { | |
| 388 gripe_square_matrix_required ("qz"); | |
| 389 return retval; | |
| 390 } | |
| 3183 | 391 |
| 392 // Argument 1: dimensions look good; get the value | |
| 393 Matrix aa; | |
| 394 ComplexMatrix caa; | |
| 3185 | 395 |
| 396 if (args(0).is_complex_type ()) | |
| 3183 | 397 caa = args(0).complex_matrix_value (); |
| 3185 | 398 else |
| 3183 | 399 aa = args(0).matrix_value (); |
| 3185 | 400 |
| 401 if (error_state) | |
| 3183 | 402 return retval; |
| 403 | |
| 3185 | 404 #ifdef DEBUG |
| 3538 | 405 std::cout << "qz: check argument 2" << std::endl; |
| 3185 | 406 #endif |
| 3183 | 407 |
| 408 // Extract argument 2 (bb, or cbb if complex) | |
| 3185 | 409 if ((nn != args(1).columns ()) || (nn != args(1).rows ())) |
| 410 { | |
| 411 gripe_nonconformant (); | |
| 412 return retval; | |
| 413 } | |
| 414 | |
| 3183 | 415 Matrix bb; |
| 416 ComplexMatrix cbb; | |
| 3185 | 417 |
| 418 if (args(1).is_complex_type ()) | |
| 3183 | 419 cbb = args(1).complex_matrix_value (); |
| 420 else | |
| 421 bb = args(1).matrix_value (); | |
| 3185 | 422 |
| 423 if (error_state) | |
| 3183 | 424 return retval; |
| 425 | |
| 426 // Both matrices loaded, now let's check what kind of arithmetic: | |
| 427 //declared static to avoid compiler warnings about long jumps, vforks. | |
| 3185 | 428 |
| 429 static int complex_case | |
| 430 = (args(0).is_complex_type () || args(1).is_complex_type ()); | |
| 3183 | 431 |
| 3185 | 432 if (nargin == 3 && complex_case) |
| 433 { | |
| 434 error ("qz: cannot re-order complex qz decomposition."); | |
| 435 return retval; | |
| 436 } | |
| 3183 | 437 |
| 438 // first, declare variables used in both the real and complex case | |
| 439 Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn); | |
| 440 RowVector alphar(nn), alphai(nn), betar(nn); | |
| 441 | |
| 3185 | 442 ComplexMatrix CQ(nn,nn), CZ(nn,nn), CVR(nn,nn), CVL(nn,nn); |
| 5275 | 443 octave_idx_type ilo, ihi, info; |
| 3185 | 444 char compq = (nargout >= 3 ? 'V' : 'N'); |
| 445 char compz = (nargout >= 4 ? 'V' : 'N'); | |
| 3183 | 446 |
| 3185 | 447 // initialize Q, Z to identity if we need either of them |
| 448 if (compq == 'V' || compz == 'V') | |
| 5275 | 449 for (octave_idx_type ii = 0; ii < nn; ii++) |
| 450 for (octave_idx_type jj = 0; jj < nn; jj++) | |
| 4153 | 451 { |
| 452 OCTAVE_QUIT; | |
| 453 QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0); | |
| 454 } | |
| 3183 | 455 |
| 3185 | 456 // always perform permutation balancing |
| 4552 | 457 const char bal_job = 'P'; |
| 3183 | 458 RowVector lscale(nn), rscale(nn), work(6*nn); |
| 459 | |
| 3185 | 460 if (complex_case) |
| 461 { | |
| 462 error ("Complex case not implemented yet"); | |
| 463 return retval; | |
| 464 } | |
| 3183 | 465 else |
| 3185 | 466 { |
| 467 #ifdef DEBUG | |
| 468 if (compq == 'V') | |
| 3538 | 469 std::cout << "qz: performing balancing; QQ=" << std::endl << QQ << std::endl; |
| 3185 | 470 #endif |
| 3183 | 471 |
| 3185 | 472 F77_XFCN (dggbal, DGGBAL, |
| 4552 | 473 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
| 474 nn, aa.fortran_vec (), nn, bb.fortran_vec (), | |
| 475 nn, ilo, ihi, lscale.fortran_vec (), | |
| 476 rscale.fortran_vec (), work.fortran_vec (), info | |
| 477 F77_CHAR_ARG_LEN (1))); | |
| 3185 | 478 } |
| 3183 | 479 |
| 480 // Since we just want the balancing matrices, we can use dggbal | |
| 481 // for both the real and complex cases; | |
| 482 // left first | |
| 3185 | 483 |
| 484 if (compq == 'V') | |
| 485 { | |
| 486 F77_XFCN (dggbak, DGGBAK, | |
| 4552 | 487 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
| 488 F77_CONST_CHAR_ARG2 ("L", 1), | |
| 4566 | 489 nn, ilo, ihi, lscale.data (), rscale.data (), |
| 490 nn, QQ.fortran_vec (), nn, info | |
| 4552 | 491 F77_CHAR_ARG_LEN (1) |
| 492 F77_CHAR_ARG_LEN (1))); | |
| 3183 | 493 |
| 3185 | 494 #ifdef DEBUG |
| 495 if (compq == 'V') | |
| 3538 | 496 std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl; |
| 3185 | 497 #endif |
| 3183 | 498 } |
| 499 | |
| 500 // then right | |
| 3185 | 501 if (compz == 'V') |
| 502 { | |
| 4552 | 503 F77_XFCN (dggbak, DGGBAK, |
| 504 (F77_CONST_CHAR_ARG2 (&bal_job, 1), | |
| 505 F77_CONST_CHAR_ARG2 ("R", 1), | |
| 4566 | 506 nn, ilo, ihi, lscale.data (), rscale.data (), |
| 507 nn, ZZ.fortran_vec (), nn, info | |
| 4552 | 508 F77_CHAR_ARG_LEN (1) |
| 509 F77_CHAR_ARG_LEN (1))); | |
| 3183 | 510 |
| 3185 | 511 #ifdef DEBUG |
| 512 if (compz == 'V') | |
| 3538 | 513 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl; |
| 3185 | 514 #endif |
| 515 } | |
| 3183 | 516 |
| 517 static char qz_job; | |
| 518 qz_job = (nargout < 2 ? 'E' : 'S'); | |
| 3185 | 519 |
| 3183 | 520 if (complex_case) |
| 3185 | 521 { |
| 522 // complex case | |
| 523 if (args(0).is_real_type ()) | |
| 3587 | 524 caa = ComplexMatrix (aa); |
| 3185 | 525 |
| 526 if (args(1).is_real_type ()) | |
| 3587 | 527 cbb = ComplexMatrix (bb); |
| 3185 | 528 |
| 529 if (compq == 'V') | |
| 3587 | 530 CQ = ComplexMatrix (QQ); |
| 3185 | 531 |
| 532 if (compz == 'V') | |
| 3587 | 533 CZ = ComplexMatrix (ZZ); |
| 3185 | 534 |
| 535 error ("complex case not done yet"); | |
| 536 return retval; | |
| 537 } | |
| 3183 | 538 else // real matrices case |
| 3185 | 539 { |
| 540 #ifdef DEBUG | |
| 3538 | 541 std::cout << "qz: peforming qr decomposition of bb" << std::endl; |
| 3185 | 542 #endif |
| 3183 | 543 |
| 3185 | 544 // compute the QR factorization of bb |
| 545 QR bqr (bb); | |
| 546 | |
| 547 #ifdef DEBUG | |
| 3538 | 548 std::cout << "qz: qr (bb) done; now peforming qz decomposition" << std::endl; |
| 3185 | 549 #endif |
| 3183 | 550 |
| 3185 | 551 bb = bqr.R (); |
| 552 | |
| 553 #ifdef DEBUG | |
| 3538 | 554 std::cout << "qz: extracted bb" << std::endl; |
| 3185 | 555 #endif |
| 3183 | 556 |
| 3185 | 557 aa = (bqr.Q ()).transpose ()*aa; |
| 558 | |
| 559 #ifdef DEBUG | |
| 3538 | 560 std::cout << "qz: updated aa " << std::endl; |
| 561 std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl; | |
| 3183 | 562 |
| 3185 | 563 if (compq == 'V') |
| 3538 | 564 std::cout << "QQ =" << QQ << std::endl; |
| 3185 | 565 #endif |
| 3183 | 566 |
| 3185 | 567 if (compq == 'V') |
| 568 QQ = QQ*bqr.Q (); | |
| 3183 | 569 |
| 3185 | 570 #ifdef DEBUG |
| 3538 | 571 std::cout << "qz: precursors done..." << std::endl; |
| 3185 | 572 #endif |
| 3183 | 573 |
| 3185 | 574 #ifdef DEBUG |
| 3538 | 575 std::cout << "qz: compq = " << compq << ", compz = " << compz << std::endl; |
| 3185 | 576 #endif |
| 3183 | 577 |
| 3185 | 578 // reduce to generalized hessenberg form |
| 579 F77_XFCN (dgghrd, DGGHRD, | |
| 4552 | 580 (F77_CONST_CHAR_ARG2 (&compq, 1), |
| 581 F77_CONST_CHAR_ARG2 (&compz, 1), | |
| 582 nn, ilo, ihi, aa.fortran_vec (), | |
| 583 nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn, | |
| 584 ZZ.fortran_vec (), nn, info | |
| 585 F77_CHAR_ARG_LEN (1) | |
| 586 F77_CHAR_ARG_LEN (1))); | |
| 3183 | 587 |
| 3185 | 588 // check if just computing generalized eigenvalues or if we're |
| 589 // actually computing the decomposition | |
| 3183 | 590 |
| 3185 | 591 // reduce to generalized Schur form |
| 592 F77_XFCN (dhgeqz, DHGEQZ, | |
| 4552 | 593 (F77_CONST_CHAR_ARG2 (&qz_job, 1), |
| 594 F77_CONST_CHAR_ARG2 (&compq, 1), | |
| 595 F77_CONST_CHAR_ARG2 (&compz, 1), | |
| 596 nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (), | |
| 597 nn, alphar.fortran_vec (), alphai.fortran_vec (), | |
| 598 betar.fortran_vec (), QQ.fortran_vec (), nn, | |
| 599 ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info | |
| 600 F77_CHAR_ARG_LEN (1) | |
| 601 F77_CHAR_ARG_LEN (1) | |
| 602 F77_CHAR_ARG_LEN (1))); | |
| 3185 | 603 } |
| 3183 | 604 |
| 605 // order the QZ decomposition? | |
| 3185 | 606 if (ord_job[0] != 'N') |
| 3183 | 607 { |
| 3185 | 608 if (complex_case) |
| 609 { | |
| 610 // probably not needed, but better be safe | |
| 611 error ("qz: cannot re-order complex qz decomposition."); | |
| 612 return retval; | |
| 613 } | |
| 614 else | |
| 615 { | |
| 616 #ifdef DEBUG_SORT | |
| 3538 | 617 std::cout << "qz: ordering eigenvalues: ord_job = " << ord_job[0] << std::endl; |
| 3185 | 618 #endif |
| 3183 | 619 |
| 3185 | 620 // declared static to avoid vfork/long jump compiler complaints |
| 621 static sort_function sort_test; | |
| 7520 | 622 sort_test = 0; |
| 3183 | 623 |
| 3185 | 624 switch (ord_job[0]) |
| 625 { | |
| 626 case 'S': | |
| 627 sort_test = &fin; | |
| 628 break; | |
| 3183 | 629 |
| 3185 | 630 case 'B': |
| 631 sort_test = &fout; | |
| 632 break; | |
| 3183 | 633 |
| 3185 | 634 case '+': |
| 635 sort_test = &fcrhp; | |
| 636 break; | |
| 3183 | 637 |
| 3185 | 638 case '-': |
| 639 sort_test = &folhp; | |
| 640 break; | |
| 641 | |
| 642 default: | |
| 3427 | 643 // invalid order option (should never happen, since we |
| 3185 | 644 // checked the options at the top). |
| 645 panic_impossible (); | |
| 646 break; | |
| 3550 | 647 } |
| 3183 | 648 |
| 5275 | 649 octave_idx_type ndim, fail; |
| 3185 | 650 double inf_norm; |
| 651 | |
| 652 F77_XFCN (xdlange, XDLANGE, | |
| 4552 | 653 (F77_CONST_CHAR_ARG2 ("I", 1), |
| 4566 | 654 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm |
| 4552 | 655 F77_CHAR_ARG_LEN (1))); |
| 3185 | 656 |
| 657 double eps = DBL_EPSILON*inf_norm*nn; | |
| 658 | |
| 659 #ifdef DEBUG_SORT | |
| 3538 | 660 std::cout << "qz: calling dsubsp: aa=" << std::endl; |
| 3531 | 661 octave_print_internal (std::cout, aa, 0); |
| 3538 | 662 std::cout << std::endl << "bb=" << std::endl; |
| 3531 | 663 octave_print_internal (std::cout, bb, 0); |
| 3185 | 664 if (compz == 'V') |
| 665 { | |
| 3538 | 666 std::cout << std::endl << "ZZ=" << std::endl; |
| 3531 | 667 octave_print_internal (std::cout, ZZ, 0); |
| 3185 | 668 } |
| 3538 | 669 std::cout << std::endl; |
| 670 std::cout << "alphar = " << std::endl; | |
| 3531 | 671 octave_print_internal (std::cout, (Matrix) alphar, 0); |
| 3538 | 672 std::cout << std::endl << "alphai = " << std::endl; |
| 3531 | 673 octave_print_internal (std::cout, (Matrix) alphai, 0); |
| 3538 | 674 std::cout << std::endl << "beta = " << std::endl; |
| 3531 | 675 octave_print_internal (std::cout, (Matrix) betar, 0); |
| 3538 | 676 std::cout << std::endl; |
| 3185 | 677 #endif |
| 678 | |
| 5275 | 679 Array<octave_idx_type> ind (nn); |
| 3550 | 680 |
| 3185 | 681 F77_XFCN (dsubsp, DSUBSP, |
| 4552 | 682 (nn, nn, aa.fortran_vec (), bb.fortran_vec (), |
| 683 ZZ.fortran_vec (), sort_test, eps, ndim, fail, | |
| 3550 | 684 ind.fortran_vec ())); |
| 3185 | 685 |
| 686 #ifdef DEBUG | |
| 3538 | 687 std::cout << "qz: back from dsubsp: aa=" << std::endl; |
| 3531 | 688 octave_print_internal (std::cout, aa, 0); |
| 3538 | 689 std::cout << std::endl << "bb=" << std::endl; |
| 3531 | 690 octave_print_internal (std::cout, bb, 0); |
| 3185 | 691 if (compz == 'V') |
| 692 { | |
| 3538 | 693 std::cout << std::endl << "ZZ=" << std::endl; |
| 3531 | 694 octave_print_internal (std::cout, ZZ, 0); |
| 3185 | 695 } |
| 3538 | 696 std::cout << std::endl; |
| 3185 | 697 #endif |
| 698 | |
| 699 // manually update alphar, alphai, betar | |
| 700 static int jj; | |
| 701 | |
| 702 jj=0; | |
| 703 while (jj < nn) | |
| 704 { | |
| 705 #ifdef DEBUG_EIG | |
| 3538 | 706 std::cout << "computing gen eig #" << jj << std::endl; |
| 3185 | 707 #endif |
| 708 | |
| 709 static int zcnt; // number of zeros in this block | |
| 710 | |
| 711 if (jj == (nn-1)) | |
| 712 zcnt = 1; | |
| 713 else if (aa(jj+1,jj) == 0) | |
| 714 zcnt = 1; | |
| 715 else zcnt = 2; | |
| 716 | |
| 717 if (zcnt == 1) // real zero | |
| 718 { | |
| 719 #ifdef DEBUG_EIG | |
| 3538 | 720 std::cout << " single gen eig:" << std::endl; |
| 721 std::cout << " alphar(" << jj << ") = " << aa(jj,jj) << std::endl; | |
| 722 std::cout << " betar( " << jj << ") = " << bb(jj,jj) << std::endl; | |
| 723 std::cout << " alphai(" << jj << ") = 0" << std::endl; | |
| 3185 | 724 #endif |
| 725 | |
| 726 alphar(jj) = aa(jj,jj); | |
| 727 alphai(jj) = 0; | |
| 728 betar(jj) = bb(jj,jj); | |
| 729 } | |
| 730 else | |
| 731 { | |
| 732 // complex conjugate pair | |
| 733 #ifdef DEBUG_EIG | |
| 3538 | 734 std::cout << "qz: calling dlag2:" << std::endl; |
| 3531 | 735 std::cout << "safmin=" |
| 3538 | 736 << setiosflags (std::ios::scientific) << safmin << std::endl; |
| 3185 | 737 |
| 738 for (int idr = jj; idr <= jj+1; idr++) | |
| 739 { | |
| 740 for (int idc = jj; idc <= jj+1; idc++) | |
| 741 { | |
| 3531 | 742 std::cout << "aa(" << idr << "," << idc << ")=" |
| 3538 | 743 << aa(idr,idc) << std::endl; |
| 3531 | 744 std::cout << "bb(" << idr << "," << idc << ")=" |
| 3538 | 745 << bb(idr,idc) << std::endl; |
| 3185 | 746 } |
| 747 } | |
| 748 #endif | |
| 749 | |
| 5775 | 750 // FIXME -- probably should be using |
| 4566 | 751 // fortran_vec instead of &aa(jj,jj) here. |
| 752 | |
| 3185 | 753 double scale1, scale2, wr1, wr2, wi; |
| 4629 | 754 const double *aa_ptr = aa.data () + jj*nn+jj; |
| 755 const double *bb_ptr = bb.data () + jj*nn+jj; | |
| 3185 | 756 F77_XFCN (dlag2, DLAG2, |
| 4629 | 757 (aa_ptr, nn, bb_ptr, nn, safmin, |
| 3185 | 758 scale1, scale2, wr1, wr2, wi)); |
| 759 | |
| 760 #ifdef DEBUG_EIG | |
| 3531 | 761 std::cout << "dlag2 returns: scale1=" << scale1 |
| 3538 | 762 << "\tscale2=" << scale2 << std::endl |
| 3185 | 763 << "\twr1=" << wr1 << "\twr2=" << wr2 |
| 3538 | 764 << "\twi=" << wi << std::endl; |
| 3185 | 765 #endif |
| 766 | |
| 767 // just to be safe, check if it's a real pair | |
| 768 if (wi == 0) | |
| 769 { | |
| 770 alphar(jj) = wr1; | |
| 771 alphai(jj) = 0; | |
| 772 betar(jj) = scale1; | |
| 773 alphar(jj+1) = wr2; | |
| 774 alphai(jj+1) = 0; | |
| 775 betar(jj+1) = scale2; | |
| 776 } | |
| 777 else | |
| 778 { | |
| 779 alphar(jj) = alphar(jj+1)=wr1; | |
| 780 alphai(jj) = -(alphai(jj+1) = wi); | |
| 781 betar(jj) = betar(jj+1) = scale1; | |
| 782 } | |
| 783 } | |
| 784 | |
| 785 // advance past this block | |
| 786 jj += zcnt; | |
| 787 } | |
| 788 | |
| 789 #ifdef DEBUG_SORT | |
| 3538 | 790 std::cout << "qz: back from dsubsp: aa=" << std::endl; |
| 3531 | 791 octave_print_internal (std::cout, aa, 0); |
| 3538 | 792 std::cout << std::endl << "bb=" << std::endl; |
| 3531 | 793 octave_print_internal (std::cout, bb, 0); |
| 3185 | 794 |
| 795 if (compz == 'V') | |
| 796 { | |
| 3538 | 797 std::cout << std::endl << "ZZ=" << std::endl; |
| 3531 | 798 octave_print_internal (std::cout, ZZ, 0); |
| 3185 | 799 } |
| 3538 | 800 std::cout << std::endl << "qz: ndim=" << ndim << std::endl |
| 801 << "fail=" << fail << std::endl; | |
| 802 std::cout << "alphar = " << std::endl; | |
| 3531 | 803 octave_print_internal (std::cout, (Matrix) alphar, 0); |
| 3538 | 804 std::cout << std::endl << "alphai = " << std::endl; |
| 3531 | 805 octave_print_internal (std::cout, (Matrix) alphai, 0); |
| 3538 | 806 std::cout << std::endl << "beta = " << std::endl; |
| 3531 | 807 octave_print_internal (std::cout, (Matrix) betar, 0); |
| 3538 | 808 std::cout << std::endl; |
| 3185 | 809 #endif |
| 810 } | |
| 3183 | 811 } |
| 3185 | 812 |
| 3183 | 813 // compute generalized eigenvalues? |
| 814 ComplexColumnVector gev; | |
| 3185 | 815 |
| 816 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) | |
| 3183 | 817 { |
| 3185 | 818 if (complex_case) |
| 819 { | |
| 820 error ("complex case not yet implemented"); | |
| 821 return retval; | |
| 822 } | |
| 823 else | |
| 824 { | |
| 825 #ifdef DEBUG | |
| 3538 | 826 std::cout << "qz: computing generalized eigenvalues" << std::endl; |
| 3185 | 827 #endif |
| 3183 | 828 |
| 3185 | 829 // return finite generalized eigenvalues |
| 830 int cnt = 0; | |
| 831 | |
| 832 for (int ii = 0; ii < nn; ii++) | |
| 833 if (betar(ii) != 0) | |
| 834 cnt++; | |
| 835 | |
| 836 ComplexColumnVector tmp(cnt); | |
| 837 | |
| 3671 | 838 cnt = 0; |
| 3185 | 839 for (int ii = 0; ii < nn; ii++) |
| 840 if (betar(ii) != 0) | |
| 3671 | 841 tmp(cnt++) = Complex(alphar(ii), alphai(ii))/betar(ii); |
| 3185 | 842 gev = tmp; |
| 843 } | |
| 3183 | 844 } |
| 845 | |
| 846 // right, left eigenvector matrices | |
| 3185 | 847 if (nargout >= 5) |
| 3183 | 848 { |
| 3185 | 849 char side = (nargout == 5 ? 'R' : 'B'); // which side to compute? |
| 850 char howmny = 'B'; // compute all of them and backtransform | |
| 7520 | 851 octave_idx_type *select = 0; // dummy pointer; select is not used. |
| 5275 | 852 octave_idx_type m; |
| 3185 | 853 |
| 854 if (complex_case) | |
| 855 { | |
| 856 error ("complex type not yet implemented"); | |
| 857 return retval; | |
| 858 } | |
| 859 else | |
| 860 { | |
| 861 #ifdef DEBUG | |
| 3538 | 862 std::cout << "qz: computing generalized eigenvectors" << std::endl; |
| 3185 | 863 #endif |
| 864 | |
| 865 VL = QQ; | |
| 866 VR = ZZ; | |
| 867 | |
| 868 F77_XFCN (dtgevc, DTGEVC, | |
| 4552 | 869 (F77_CONST_CHAR_ARG2 (&side, 1), |
| 870 F77_CONST_CHAR_ARG2 (&howmny, 1), | |
| 871 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (), | |
| 872 nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn, | |
| 873 m, work.fortran_vec (), info | |
| 874 F77_CHAR_ARG_LEN (1) | |
| 875 F77_CHAR_ARG_LEN (1))); | |
| 3185 | 876 |
| 877 // now construct the complex form of VV, WW | |
| 878 int jj = 0; | |
| 879 | |
| 880 while (jj < nn) | |
| 881 { | |
| 4153 | 882 OCTAVE_QUIT; |
| 883 | |
| 3185 | 884 // see if real or complex eigenvalue |
| 885 int cinc = 2; // column increment; assume complex eigenvalue | |
| 886 | |
| 887 if (jj == (nn-1)) | |
| 888 cinc = 1; // single column | |
| 889 else if (aa(jj+1,jj) == 0) | |
| 890 cinc = 1; | |
| 891 | |
| 892 // now copy the eigenvector (s) to CVR, CVL | |
| 893 if (cinc == 1) | |
| 894 { | |
| 895 for (int ii = 0; ii < nn; ii++) | |
| 896 CVR(ii,jj) = VR(ii,jj); | |
| 897 | |
| 898 if (side == 'B') | |
| 899 for (int ii = 0; ii < nn; ii++) | |
| 900 CVL(ii,jj) = VL(ii,jj); | |
| 901 } | |
| 902 else | |
| 903 { | |
| 904 // double column; complex vector | |
| 905 | |
| 906 for (int ii = 0; ii < nn; ii++) | |
| 907 { | |
| 908 CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1)); | |
| 909 CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1)); | |
| 910 } | |
| 3183 | 911 |
| 3185 | 912 if (side == 'B') |
| 913 for (int ii = 0; ii < nn; ii++) | |
| 914 { | |
| 915 CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1)); | |
| 916 CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1)); | |
| 917 } | |
| 918 } | |
| 919 | |
| 920 // advance to next eigenvectors (if any) | |
| 921 jj += cinc; | |
| 922 } | |
| 923 } | |
| 3183 | 924 } |
| 3185 | 925 |
| 926 switch (nargout) | |
| 927 { | |
| 928 case 7: | |
| 929 retval(6) = gev; | |
| 930 | |
| 931 case 6: // return eigenvectors | |
| 932 retval(5) = CVL; | |
| 933 | |
| 934 case 5: // return eigenvectors | |
| 935 retval(4) = CVR; | |
| 936 | |
| 937 case 4: | |
| 938 if (nargin == 3) | |
| 939 { | |
| 940 #ifdef DEBUG | |
| 3538 | 941 std::cout << "qz: sort: retval(3) = gev = " << std::endl; |
| 3531 | 942 octave_print_internal (std::cout, gev); |
| 3538 | 943 std::cout << std::endl; |
| 3185 | 944 #endif |
| 945 retval(3) = gev; | |
| 946 } | |
| 947 else | |
| 948 retval(3) = ZZ; | |
| 949 | |
| 950 case 3: | |
| 951 if (nargin == 3) | |
| 952 retval(2) = ZZ; | |
| 953 else | |
| 954 retval(2) = QQ; | |
| 955 | |
| 956 case 2: | |
| 957 #ifdef DEBUG | |
| 3538 | 958 std::cout << "qz: retval (1) = bb = " << std::endl; |
| 3531 | 959 octave_print_internal (std::cout, bb, 0); |
| 3538 | 960 std::cout << std::endl << "qz: retval(0) = aa = " <<std::endl; |
| 3531 | 961 octave_print_internal (std::cout, aa, 0); |
| 3538 | 962 std::cout << std::endl; |
| 3185 | 963 #endif |
| 964 retval(1) = bb; | |
| 965 retval(0) = aa; | |
| 966 break; | |
| 967 | |
| 968 case 1: | |
| 969 case 0: | |
| 970 #ifdef DEBUG | |
| 3538 | 971 std::cout << "qz: retval(0) = gev = " << gev << std::endl; |
| 3185 | 972 #endif |
| 973 retval(0) = gev; | |
| 974 break; | |
| 975 | |
| 976 default: | |
| 977 error ("qz: too many return arguments."); | |
| 978 break; | |
| 3183 | 979 } |
| 980 | |
| 3185 | 981 #ifdef DEBUG |
| 3538 | 982 std::cout << "qz: exiting (at long last)" << std::endl; |
| 3185 | 983 #endif |
| 3183 | 984 |
| 985 return retval; | |
| 986 } | |
| 987 | |
| 988 /* | |
| 989 ;;; Local Variables: *** | |
| 990 ;;; mode: C++ *** | |
| 991 ;;; End: *** | |
| 992 */ |