Mercurial > hg > octave-lyh

comparison src/DLD-FUNCTIONS/qz.cc @ 10311:a217e1d74353

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
untabify qz.cc
author John W. Eaton <jwe@octave.org>
date Thu, 11 Feb 2010 11:57:36 -0500
parents 4e4270ab70d6
children 12884915a8e4
comparison
equal deleted inserted replaced
10310:cd14d826025f 10311:a217e1d74353
65 65
66 extern "C" 66 extern "C"
67 { 67 {
68 F77_RET_T 68 F77_RET_T
69 F77_FUNC (dggbal, DGGBAL) (F77_CONST_CHAR_ARG_DECL, 69 F77_FUNC (dggbal, DGGBAL) (F77_CONST_CHAR_ARG_DECL,
70 const octave_idx_type& N, double* A, const octave_idx_type& LDA, 70 const octave_idx_type& N, double* A,
71 double* B, const octave_idx_type& LDB, octave_idx_type& ILO, 71 const octave_idx_type& LDA, double* B,
72 octave_idx_type& IHI, double* LSCALE, double* RSCALE, 72 const octave_idx_type& LDB, octave_idx_type& ILO,
73 double* WORK, octave_idx_type& INFO 73 octave_idx_type& IHI, double* LSCALE,
74 double* RSCALE, double* WORK,
75 octave_idx_type& INFO
74 F77_CHAR_ARG_LEN_DECL); 76 F77_CHAR_ARG_LEN_DECL);
75 77
76 F77_RET_T 78 F77_RET_T
77 F77_FUNC (zggbal, ZGGBAL) (F77_CONST_CHAR_ARG_DECL, 79 F77_FUNC (zggbal, ZGGBAL) (F77_CONST_CHAR_ARG_DECL,
78 const octave_idx_type& N, Complex* A, const octave_idx_type& LDA, 80 const octave_idx_type& N, Complex* A,
79 Complex* B, const octave_idx_type& LDB, octave_idx_type& ILO, 81 const octave_idx_type& LDA, Complex* B,
80 octave_idx_type& IHI, double* LSCALE, double* RSCALE, 82 const octave_idx_type& LDB, octave_idx_type& ILO,
81 double* WORK, octave_idx_type& INFO 83 octave_idx_type& IHI, double* LSCALE,
82 F77_CHAR_ARG_LEN_DECL); 84 double* RSCALE, double* WORK,
85 octave_idx_type& INFO
86 F77_CHAR_ARG_LEN_DECL);
83 87
84 F77_RET_T 88 F77_RET_T
85 F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL, 89 F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL,
86 F77_CONST_CHAR_ARG_DECL, 90 F77_CONST_CHAR_ARG_DECL,
87 const octave_idx_type& N, const octave_idx_type& ILO, 91 const octave_idx_type& N,
88 const octave_idx_type& IHI, const double* LSCALE, 92 const octave_idx_type& ILO,
89 const double* RSCALE, octave_idx_type& M, double* V, 93 const octave_idx_type& IHI,
94 const double* LSCALE, const double* RSCALE,
95 octave_idx_type& M, double* V,
90 const octave_idx_type& LDV, octave_idx_type& INFO 96 const octave_idx_type& LDV, octave_idx_type& INFO
91 F77_CHAR_ARG_LEN_DECL 97 F77_CHAR_ARG_LEN_DECL
92 F77_CHAR_ARG_LEN_DECL); 98 F77_CHAR_ARG_LEN_DECL);
93 99
94 F77_RET_T 100 F77_RET_T
95 F77_FUNC (zggbak, ZGGBAK) (F77_CONST_CHAR_ARG_DECL, 101 F77_FUNC (zggbak, ZGGBAK) (F77_CONST_CHAR_ARG_DECL,
96 F77_CONST_CHAR_ARG_DECL, 102 F77_CONST_CHAR_ARG_DECL,
97 const octave_idx_type& N, const octave_idx_type& ILO, 103 const octave_idx_type& N,
98 const octave_idx_type& IHI, const double* LSCALE, 104 const octave_idx_type& ILO,
99 const double* RSCALE, octave_idx_type& M, Complex* V, 105 const octave_idx_type& IHI,
106 const double* LSCALE, const double* RSCALE,
107 octave_idx_type& M, Complex* V,
100 const octave_idx_type& LDV, octave_idx_type& INFO 108 const octave_idx_type& LDV, octave_idx_type& INFO
101 F77_CHAR_ARG_LEN_DECL 109 F77_CHAR_ARG_LEN_DECL
102 F77_CHAR_ARG_LEN_DECL); 110 F77_CHAR_ARG_LEN_DECL);
103 111
104 F77_RET_T 112 F77_RET_T
105 F77_FUNC (dgghrd, DGGHRD) (F77_CONST_CHAR_ARG_DECL, 113 F77_FUNC (dgghrd, DGGHRD) (F77_CONST_CHAR_ARG_DECL,
106 F77_CONST_CHAR_ARG_DECL, 114 F77_CONST_CHAR_ARG_DECL,
107 const octave_idx_type& N, const octave_idx_type& ILO, 115 const octave_idx_type& N,
116 const octave_idx_type& ILO,
108 const octave_idx_type& IHI, double* A, 117 const octave_idx_type& IHI, double* A,
109 const octave_idx_type& LDA, double* B, 118 const octave_idx_type& LDA, double* B,
110 const octave_idx_type& LDB, double* Q, 119 const octave_idx_type& LDB, double* Q,
111 const octave_idx_type& LDQ, double* Z, 120 const octave_idx_type& LDQ, double* Z,
112 const octave_idx_type& LDZ, octave_idx_type& INFO 121 const octave_idx_type& LDZ, octave_idx_type& INFO
114 F77_CHAR_ARG_LEN_DECL); 123 F77_CHAR_ARG_LEN_DECL);
115 124
116 F77_RET_T 125 F77_RET_T
117 F77_FUNC (zgghrd, ZGGHRD) (F77_CONST_CHAR_ARG_DECL, 126 F77_FUNC (zgghrd, ZGGHRD) (F77_CONST_CHAR_ARG_DECL,
118 F77_CONST_CHAR_ARG_DECL, 127 F77_CONST_CHAR_ARG_DECL,
119 const octave_idx_type& N, const octave_idx_type& ILO, 128 const octave_idx_type& N,
129 const octave_idx_type& ILO,
120 const octave_idx_type& IHI, Complex* A, 130 const octave_idx_type& IHI, Complex* A,
121 const octave_idx_type& LDA, Complex* B, 131 const octave_idx_type& LDA, Complex* B,
122 const octave_idx_type& LDB, Complex* Q, 132 const octave_idx_type& LDB, Complex* Q,
123 const octave_idx_type& LDQ, Complex* Z, 133 const octave_idx_type& LDQ, Complex* Z,
124 const octave_idx_type& LDZ, octave_idx_type& INFO 134 const octave_idx_type& LDZ, octave_idx_type& INFO
127 137
128 F77_RET_T 138 F77_RET_T
129 F77_FUNC (dhgeqz, DHGEQZ) (F77_CONST_CHAR_ARG_DECL, 139 F77_FUNC (dhgeqz, DHGEQZ) (F77_CONST_CHAR_ARG_DECL,
130 F77_CONST_CHAR_ARG_DECL, 140 F77_CONST_CHAR_ARG_DECL,
131 F77_CONST_CHAR_ARG_DECL, 141 F77_CONST_CHAR_ARG_DECL,
132 const octave_idx_type& N, const octave_idx_type& ILO, const octave_idx_type& IHI, 142 const octave_idx_type& N,
143 const octave_idx_type& ILO,
144 const octave_idx_type& IHI,
133 double* A, const octave_idx_type& LDA, double* B, 145 double* A, const octave_idx_type& LDA, double* B,
134 const octave_idx_type& LDB, double* ALPHAR, 146 const octave_idx_type& LDB, double* ALPHAR,
135 double* ALPHAI, double* BETA, double* Q, 147 double* ALPHAI, double* BETA, double* Q,
136 const octave_idx_type& LDQ, double* Z, 148 const octave_idx_type& LDQ, double* Z,
137 const octave_idx_type& LDZ, double* WORK, 149 const octave_idx_type& LDZ, double* WORK,
138 const octave_idx_type& LWORK, octave_idx_type& INFO 150 const octave_idx_type& LWORK,
151 octave_idx_type& INFO
139 F77_CHAR_ARG_LEN_DECL 152 F77_CHAR_ARG_LEN_DECL
140 F77_CHAR_ARG_LEN_DECL 153 F77_CHAR_ARG_LEN_DECL
141 F77_CHAR_ARG_LEN_DECL); 154 F77_CHAR_ARG_LEN_DECL);
142 155
143 F77_RET_T 156 F77_RET_T
144 F77_FUNC (zhgeqz, ZHGEQZ) (F77_CONST_CHAR_ARG_DECL, 157 F77_FUNC (zhgeqz, ZHGEQZ) (F77_CONST_CHAR_ARG_DECL,
145 F77_CONST_CHAR_ARG_DECL, 158 F77_CONST_CHAR_ARG_DECL,
146 F77_CONST_CHAR_ARG_DECL, 159 F77_CONST_CHAR_ARG_DECL,
147 const octave_idx_type& N, const octave_idx_type& ILO, const octave_idx_type& IHI, 160 const octave_idx_type& N,
148 Complex* A, const octave_idx_type& LDA, Complex* B, 161 const octave_idx_type& ILO,
149 const octave_idx_type& LDB, Complex* ALPHA, Complex* BETA, Complex* CQ, const octave_idx_type& LDQ, 162 const octave_idx_type& IHI,
150 Complex* CZ, const octave_idx_type& LDZ, 163 Complex* A, const octave_idx_type& LDA,
151 Complex* WORK, 164 Complex* B, const octave_idx_type& LDB,
152 const octave_idx_type& LWORK, double* RWORK, 165 Complex* ALPHA, Complex* BETA, Complex* CQ,
153 octave_idx_type& INFO 166 const octave_idx_type& LDQ,
167 Complex* CZ, const octave_idx_type& LDZ,
168 Complex* WORK, const octave_idx_type& LWORK,
169 double* RWORK, octave_idx_type& INFO
154 F77_CHAR_ARG_LEN_DECL 170 F77_CHAR_ARG_LEN_DECL
155 F77_CHAR_ARG_LEN_DECL 171 F77_CHAR_ARG_LEN_DECL
156 F77_CHAR_ARG_LEN_DECL); 172 F77_CHAR_ARG_LEN_DECL);
157 173
158 F77_RET_T 174 F77_RET_T
159 F77_FUNC (dlag2, DLAG2) (const double* A, const octave_idx_type& LDA, const double* B, 175 F77_FUNC (dlag2, DLAG2) (const double* A, const octave_idx_type& LDA,
160 const octave_idx_type& LDB, const double& SAFMIN, 176 const double* B, const octave_idx_type& LDB,
161 double& SCALE1, double& SCALE2, 177 const double& SAFMIN, double& SCALE1,
162 double& WR1, double& WR2, double& WI); 178 double& SCALE2, double& WR1, double& WR2,
179 double& WI);
163 180
164 // Van Dooren's code (netlib.org: toms/590) for reordering 181 // Van Dooren's code (netlib.org: toms/590) for reordering
165 // GEP. Only processes Z, not Q. 182 // GEP. Only processes Z, not Q.
166 F77_RET_T 183 F77_RET_T
167 F77_FUNC (dsubsp, DSUBSP) (const octave_idx_type& NMAX, const octave_idx_type& N, double* A, 184 F77_FUNC (dsubsp, DSUBSP) (const octave_idx_type& NMAX,
185 const octave_idx_type& N, double* A,
168 double* B, double* Z, sort_function, 186 double* B, double* Z, sort_function,
169 const double& EPS, octave_idx_type& NDIM, octave_idx_type& FAIL, 187 const double& EPS, octave_idx_type& NDIM,
170 octave_idx_type* IND); 188 octave_idx_type& FAIL, octave_idx_type* IND);
171 189
172 // documentation for DTGEVC incorrectly states that VR, VL are 190 // Documentation for DTGEVC incorrectly states that VR, VL are
173 // complex*16; they are declared in DTGEVC as double precision 191 // complex*16; they are declared in DTGEVC as double precision
174 // (probably a cut and paste problem fro ZTGEVC) 192 // (probably a cut and paste problem fro ZTGEVC).
175 F77_RET_T 193 F77_RET_T
176 F77_FUNC (dtgevc, DTGEVC) (F77_CONST_CHAR_ARG_DECL, 194 F77_FUNC (dtgevc, DTGEVC) (F77_CONST_CHAR_ARG_DECL,
177 F77_CONST_CHAR_ARG_DECL, 195 F77_CONST_CHAR_ARG_DECL,
178 octave_idx_type* SELECT, const octave_idx_type& N, double* A, 196 octave_idx_type* SELECT,
197 const octave_idx_type& N, double* A,
179 const octave_idx_type& LDA, double* B, 198 const octave_idx_type& LDA, double* B,
180 const octave_idx_type& LDB, double* VL, 199 const octave_idx_type& LDB, double* VL,
181 const octave_idx_type& LDVL, double* VR, 200 const octave_idx_type& LDVL, double* VR,
182 const octave_idx_type& LDVR, const octave_idx_type& MM, 201 const octave_idx_type& LDVR,
183 octave_idx_type& M, double* WORK, octave_idx_type& INFO 202 const octave_idx_type& MM, octave_idx_type& M,
203 double* WORK, octave_idx_type& INFO
184 F77_CHAR_ARG_LEN_DECL 204 F77_CHAR_ARG_LEN_DECL
185 F77_CHAR_ARG_LEN_DECL); 205 F77_CHAR_ARG_LEN_DECL);
186 206
187 F77_RET_T 207 F77_RET_T
188 F77_FUNC (ztgevc, ZTGEVC) (F77_CONST_CHAR_ARG_DECL, 208 F77_FUNC (ztgevc, ZTGEVC) (F77_CONST_CHAR_ARG_DECL,
189 F77_CONST_CHAR_ARG_DECL, 209 F77_CONST_CHAR_ARG_DECL,
190 octave_idx_type* SELECT, const octave_idx_type& N,const Complex* A, 210 octave_idx_type* SELECT,
211 const octave_idx_type& N, const Complex* A,
191 const octave_idx_type& LDA,const Complex* B, 212 const octave_idx_type& LDA,const Complex* B,
192 const octave_idx_type& LDB, Complex* xVL, 213 const octave_idx_type& LDB, Complex* xVL,
193 const octave_idx_type& LDVL, Complex* xVR, 214 const octave_idx_type& LDVL, Complex* xVR,
194 const octave_idx_type& LDVR, const octave_idx_type& MM, 215 const octave_idx_type& LDVR,
195 octave_idx_type& M, Complex* CWORK, double* RWORK, octave_idx_type& INFO 216 const octave_idx_type& MM, octave_idx_type& M,
217 Complex* CWORK, double* RWORK,
218 octave_idx_type& INFO
196 F77_CHAR_ARG_LEN_DECL 219 F77_CHAR_ARG_LEN_DECL
197 F77_CHAR_ARG_LEN_DECL); 220 F77_CHAR_ARG_LEN_DECL);
198 221
199 F77_RET_T 222 F77_RET_T
200 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG_DECL, 223 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG_DECL,
201 double& retval 224 double& retval
202 F77_CHAR_ARG_LEN_DECL); 225 F77_CHAR_ARG_LEN_DECL);
203 226
204 F77_RET_T 227 F77_RET_T
205 F77_FUNC (xdlange, XDLANGE) (F77_CONST_CHAR_ARG_DECL, 228 F77_FUNC (xdlange, XDLANGE) (F77_CONST_CHAR_ARG_DECL,
206 const octave_idx_type&, const octave_idx_type&, const double*, 229 const octave_idx_type&,
230 const octave_idx_type&, const double*,
207 const octave_idx_type&, double*, double& 231 const octave_idx_type&, double*, double&
208 F77_CHAR_ARG_LEN_DECL); 232 F77_CHAR_ARG_LEN_DECL);
209 } 233 }
210 234
211 // fcrhp, fin, fout, folhp: 235 // fcrhp, fin, fout, folhp:
219 static octave_idx_type 243 static octave_idx_type
220 fcrhp (const octave_idx_type& lsize, const double& alpha, 244 fcrhp (const octave_idx_type& lsize, const double& alpha,
221 const double& beta, const double& s, const double&) 245 const double& beta, const double& s, const double&)
222 { 246 {
223 if (lsize == 1) 247 if (lsize == 1)
224 return (alpha*beta >= 0 ? 1 : -1); 248 return (alpha * beta >= 0 ? 1 : -1);
225 else 249 else
226 return (s >= 0 ? 1 : -1); 250 return (s >= 0 ? 1 : -1);
227 } 251 }
228 252
229 static octave_idx_type 253 static octave_idx_type
247 static octave_idx_type 271 static octave_idx_type
248 folhp (const octave_idx_type& lsize, const double& alpha, 272 folhp (const octave_idx_type& lsize, const double& alpha,
249 const double& beta, const double& s, const double&) 273 const double& beta, const double& s, const double&)
250 { 274 {
251 if (lsize == 1) 275 if (lsize == 1)
252 return (alpha*beta < 0 ? 1 : -1); 276 return (alpha * beta < 0 ? 1 : -1);
253 else 277 else
254 return (s < 0 ? 1 : -1); 278 return (s < 0 ? 1 : -1);
255 } 279 }
256 280
257 static octave_idx_type 281 static octave_idx_type
291 @end tex\n\ 315 @end tex\n\
292 @ifnottex\n\ 316 @ifnottex\n\
293 @example\n\ 317 @example\n\
294 @group\n\ 318 @group\n\
295 \n\ 319 \n\
296 A*V = B*V*diag(lambda)\n\ 320 A * V = B * V * diag (lambda)\n\
297 W'*A = diag(lambda)*W'*B\n\ 321 W' * A = diag (lambda) * W' * B\n\
298 AA = Q*A*Z, BB = Q*B*Z\n\ 322 AA = Q * A * Z, BB = Q * B * Z\n\
299 \n\ 323 \n\
300 @end group\n\ 324 @end group\n\
301 @end example\n\ 325 @end example\n\
302 @end ifnottex\n\ 326 @end ifnottex\n\
303 with @var{Q} and @var{Z} orthogonal (unitary)= @var{I}\n\ 327 with @var{Q} and @var{Z} orthogonal (unitary)= @var{I}\n\
359 383
360 #ifdef DEBUG 384 #ifdef DEBUG
361 std::cout << "qz: determine ordering option" << std::endl; 385 std::cout << "qz: determine ordering option" << std::endl;
362 #endif 386 #endif
363 387
364 // Determine ordering option 388 // Determine ordering option.
365 volatile char ord_job = 0; 389 volatile char ord_job = 0;
366 static double safmin; 390 static double safmin;
367 391
368 if (nargin == 2) 392 if (nargin == 2)
369 ord_job = 'N'; 393 ord_job = 'N';
396 #ifdef DEBUG_EIG 420 #ifdef DEBUG_EIG
397 std::cout << "qz: initial value of safmin=" << setiosflags (std::ios::scientific) 421 std::cout << "qz: initial value of safmin=" << setiosflags (std::ios::scientific)
398 << safmin << std::endl; 422 << safmin << std::endl;
399 #endif 423 #endif
400 424
401 // some machines (e.g., DEC alpha) get safmin = 0; 425 // Some machines (e.g., DEC alpha) get safmin = 0;
402 // for these, use eps instead to avoid problems in dlag2 426 // for these, use eps instead to avoid problems in dlag2.
403 if (safmin == 0) 427 if (safmin == 0)
404 { 428 {
405 #ifdef DEBUG_EIG 429 #ifdef DEBUG_EIG
406 std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl; 430 std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl;
407 #endif 431 #endif
419 443
420 #ifdef DEBUG 444 #ifdef DEBUG
421 std::cout << "qz: check argument 1" << std::endl; 445 std::cout << "qz: check argument 1" << std::endl;
422 #endif 446 #endif
423 447
424 // Argument 1: check if it's o.k. dimensioned 448 // Argument 1: check if it's o.k. dimensioned.
425 octave_idx_type nn = args(0).rows (); 449 octave_idx_type nn = args(0).rows ();
426 450
427 #ifdef DEBUG 451 #ifdef DEBUG
428 std::cout << "argument 1 dimensions: (" << nn << "," << args(0).columns () << ")" 452 std::cout << "argument 1 dimensions: (" << nn << "," << args(0).columns () << ")"
429 << std::endl; 453 << std::endl;
445 { 469 {
446 gripe_square_matrix_required ("qz"); 470 gripe_square_matrix_required ("qz");
447 return retval; 471 return retval;
448 } 472 }
449 473
450 // Argument 1: dimensions look good; get the value 474 // Argument 1: dimensions look good; get the value.
451 Matrix aa; 475 Matrix aa;
452 ComplexMatrix caa; 476 ComplexMatrix caa;
453 477
454 if (args(0).is_complex_type ()) 478 if (args(0).is_complex_type ())
455 caa = args(0).complex_matrix_value (); 479 caa = args(0).complex_matrix_value ();
461 485
462 #ifdef DEBUG 486 #ifdef DEBUG
463 std::cout << "qz: check argument 2" << std::endl; 487 std::cout << "qz: check argument 2" << std::endl;
464 #endif 488 #endif
465 489
466 // Extract argument 2 (bb, or cbb if complex) 490 // Extract argument 2 (bb, or cbb if complex).
467 if ((nn != args(1).columns ()) || (nn != args(1).rows ())) 491 if ((nn != args(1).columns ()) || (nn != args(1).rows ()))
468 { 492 {
469 gripe_nonconformant (); 493 gripe_nonconformant ();
470 return retval; 494 return retval;
471 } 495 }
480 504
481 if (error_state) 505 if (error_state)
482 return retval; 506 return retval;
483 507
484 // Both matrices loaded, now let's check what kind of arithmetic: 508 // Both matrices loaded, now let's check what kind of arithmetic:
485 //declared static to avoid compiler warnings about long jumps, vforks. 509 // declared volatile to avoid compiler warnings about long jumps,
486 510 // vforks.
487 int complex_case 511
512 volatile int complex_case
488 = (args(0).is_complex_type () || args(1).is_complex_type ()); 513 = (args(0).is_complex_type () || args(1).is_complex_type ());
489 // static complex_case causing random segfault, so it is removed 514
490 if (nargin == 3 && complex_case) 515 if (nargin == 3 && complex_case)
491 { 516 {
492 error ("qz: cannot re-order complex qz decomposition."); 517 error ("qz: cannot re-order complex qz decomposition.");
493 return retval; 518 return retval;
494 } 519 }
495 520
496 // first, declare variables used in both the real and complex case 521 // First, declare variables used in both the real and complex case.
497 Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn); 522 Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn);
498 RowVector alphar(nn), alphai(nn), betar(nn); 523 RowVector alphar(nn), alphai(nn), betar(nn);
499 ComplexRowVector xalpha(nn), xbeta(nn); 524 ComplexRowVector xalpha(nn), xbeta(nn);
500 ComplexMatrix CQ(nn,nn), CZ(nn,nn), CVR(nn,nn), CVL(nn,nn); 525 ComplexMatrix CQ(nn,nn), CZ(nn,nn), CVR(nn,nn), CVL(nn,nn);
501 octave_idx_type ilo, ihi, info; 526 octave_idx_type ilo, ihi, info;
502 char compq = (nargout >= 3 ? 'V' : 'N'); 527 char compq = (nargout >= 3 ? 'V' : 'N');
503 char compz = (nargout >= 4 ? 'V' : 'N'); 528 char compz = (nargout >= 4 ? 'V' : 'N');
504 529
505 // initialize Q, Z to identity if we need either of them 530 // Initialize Q, Z to identity if we need either of them.
506 if (compq == 'V' || compz == 'V') 531 if (compq == 'V' || compz == 'V')
507 for (octave_idx_type ii = 0; ii < nn; ii++) 532 for (octave_idx_type ii = 0; ii < nn; ii++)
508 for (octave_idx_type jj = 0; jj < nn; jj++) 533 for (octave_idx_type jj = 0; jj < nn; jj++)
509 { 534 {
510 OCTAVE_QUIT; 535 OCTAVE_QUIT;
511 QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0); 536 QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0);
512 } 537 }
513 538
514 // always perform permutation balancing 539 // Always perform permutation balancing.
515 const char bal_job = 'P'; 540 const char bal_job = 'P';
516 RowVector lscale(nn), rscale(nn), work(6*nn), rwork(nn); 541 RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn);
517 542
518 if (complex_case) 543 if (complex_case)
519 { 544 {
520 #ifdef DEBUG 545 #ifdef DEBUG
521 if (compq == 'V') 546 if (compq == 'V')
522 std::cout << "qz: performing balancing; CQ=" << std::endl << CQ << std::endl; 547 std::cout << "qz: performing balancing; CQ=" << std::endl << CQ << std::endl;
523 #endif 548 #endif
524 if (args(0).is_real_type ()) 549 if (args(0).is_real_type ())
525 caa = ComplexMatrix (aa); 550 caa = ComplexMatrix (aa);
526 551
527 if (args(1).is_real_type ()) 552 if (args(1).is_real_type ())
528 cbb = ComplexMatrix (bb); 553 cbb = ComplexMatrix (bb);
529 554
530 if (compq == 'V') 555 if (compq == 'V')
531 CQ = ComplexMatrix (QQ); 556 CQ = ComplexMatrix (QQ);
532 557
533 if (compz == 'V') 558 if (compz == 'V')
534 CZ = ComplexMatrix (ZZ); 559 CZ = ComplexMatrix (ZZ);
535 560
536 F77_XFCN (zggbal, ZGGBAL, 561 F77_XFCN (zggbal, ZGGBAL,
537 (F77_CONST_CHAR_ARG2 (&bal_job, 1), 562 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
538 nn, caa.fortran_vec (), nn, cbb.fortran_vec (), 563 nn, caa.fortran_vec (), nn, cbb.fortran_vec (),
539 nn, ilo, ihi, lscale.fortran_vec (), 564 nn, ilo, ihi, lscale.fortran_vec (),
540 rscale.fortran_vec (), work.fortran_vec (), info 565 rscale.fortran_vec (), work.fortran_vec (), info
541 F77_CHAR_ARG_LEN (1))); 566 F77_CHAR_ARG_LEN (1)));
542 } 567 }
543 else 568 else
544 { 569 {
545 #ifdef DEBUG 570 #ifdef DEBUG
546 if (compq == 'V') 571 if (compq == 'V')
554 rscale.fortran_vec (), work.fortran_vec (), info 579 rscale.fortran_vec (), work.fortran_vec (), info
555 F77_CHAR_ARG_LEN (1))); 580 F77_CHAR_ARG_LEN (1)));
556 } 581 }
557 582
558 // Since we just want the balancing matrices, we can use dggbal 583 // Since we just want the balancing matrices, we can use dggbal
559 // for both the real and complex cases; 584 // for both the real and complex cases; left first
560 // left first 585
561 586 #if 0
562 /* if (compq == 'V') 587 if (compq == 'V')
563 { 588 {
564 F77_XFCN (dggbak, DGGBAK, 589 F77_XFCN (dggbak, DGGBAK,
565 (F77_CONST_CHAR_ARG2 (&bal_job, 1), 590 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
566 F77_CONST_CHAR_ARG2 ("L", 1), 591 F77_CONST_CHAR_ARG2 ("L", 1),
567 nn, ilo, ihi, lscale.data (), rscale.data (), 592 nn, ilo, ihi, lscale.data (), rscale.data (),
588 613
589 #ifdef DEBUG 614 #ifdef DEBUG
590 if (compz == 'V') 615 if (compz == 'V')
591 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl; 616 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
592 #endif 617 #endif
593 } */ 618 }
619 #endif
594 620
595 static char qz_job; 621 static char qz_job;
596 qz_job = (nargout < 2 ? 'E' : 'S'); 622 qz_job = (nargout < 2 ? 'E' : 'S');
597 623
598 if (complex_case) 624 if (complex_case)
599 { 625 {
600 // complex case 626 // Complex case.
601 ComplexQR cbqr (cbb); // declare cbqr as the QR decomposition of cbb 627
602 cbb = cbqr.R (); // the R matrix of QR decomposition for cbb 628 // The QR decomposition of cbb.
603 caa = (cbqr.Q ().hermitian ())*caa; // (Q*)caa for following work 629 ComplexQR cbqr (cbb);
604 //if (compq == 'V') 630 // The R matrix of QR decomposition for cbb.
605 CQ = CQ*cbqr.Q (); 631 cbb = cbqr.R ();
632 // (Q*)caa for following work.
633 caa = (cbqr.Q ().hermitian ()) * caa;
634 CQ = CQ * cbqr.Q ();
635
606 F77_XFCN (zgghrd, ZGGHRD, 636 F77_XFCN (zgghrd, ZGGHRD,
607 (F77_CONST_CHAR_ARG2 (&compq, 1), 637 (F77_CONST_CHAR_ARG2 (&compq, 1),
608 F77_CONST_CHAR_ARG2 (&compz, 1), 638 F77_CONST_CHAR_ARG2 (&compz, 1),
609 nn, ilo, ihi, caa.fortran_vec (), 639 nn, ilo, ihi, caa.fortran_vec (),
610 nn, cbb.fortran_vec (), nn, CQ.fortran_vec (), nn, 640 nn, cbb.fortran_vec (), nn, CQ.fortran_vec (), nn,
611 CZ.fortran_vec (), nn, info 641 CZ.fortran_vec (), nn, info
612 F77_CHAR_ARG_LEN (1) 642 F77_CHAR_ARG_LEN (1)
613 F77_CHAR_ARG_LEN (1))); 643 F77_CHAR_ARG_LEN (1)));
614 ComplexRowVector cwork(1*nn); 644
645 ComplexRowVector cwork (1 * nn);
646
615 F77_XFCN (zhgeqz, ZHGEQZ, 647 F77_XFCN (zhgeqz, ZHGEQZ,
616 (F77_CONST_CHAR_ARG2 (&qz_job, 1), 648 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
617 F77_CONST_CHAR_ARG2 (&compq, 1), 649 F77_CONST_CHAR_ARG2 (&compq, 1),
618 F77_CONST_CHAR_ARG2 (&compz, 1), 650 F77_CONST_CHAR_ARG2 (&compz, 1),
619 nn, ilo, ihi, 651 nn, ilo, ihi,
627 F77_CHAR_ARG_LEN (1) 659 F77_CHAR_ARG_LEN (1)
628 F77_CHAR_ARG_LEN (1))); 660 F77_CHAR_ARG_LEN (1)));
629 661
630 if (compq == 'V') 662 if (compq == 'V')
631 { 663 {
632 // Left eigenvector 664 // Left eigenvector.
633 F77_XFCN (zggbak, ZGGBAK, 665 F77_XFCN (zggbak, ZGGBAK,
634 (F77_CONST_CHAR_ARG2 (&bal_job, 1), 666 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
635 F77_CONST_CHAR_ARG2 ("L", 1), 667 F77_CONST_CHAR_ARG2 ("L", 1),
636 nn, ilo, ihi, lscale.data (), rscale.data (), 668 nn, ilo, ihi, lscale.data (), rscale.data (),
637 nn, CQ.fortran_vec (), nn, info 669 nn, CQ.fortran_vec (), nn, info
638 F77_CHAR_ARG_LEN (1) 670 F77_CHAR_ARG_LEN (1)
639 F77_CHAR_ARG_LEN (1))); 671 F77_CHAR_ARG_LEN (1)));
640 } 672 }
641 673
642 // then right 674 // Right eigenvector.
643 if (compz == 'V') 675 if (compz == 'V')
644 { 676 {
645 F77_XFCN (zggbak, ZGGBAK, 677 F77_XFCN (zggbak, ZGGBAK,
646 (F77_CONST_CHAR_ARG2 (&bal_job, 1), 678 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
647 F77_CONST_CHAR_ARG2 ("R", 1), 679 F77_CONST_CHAR_ARG2 ("R", 1),
650 F77_CHAR_ARG_LEN (1) 682 F77_CHAR_ARG_LEN (1)
651 F77_CHAR_ARG_LEN (1))); 683 F77_CHAR_ARG_LEN (1)));
652 } 684 }
653 685
654 } 686 }
655 else // real matrices case 687 else
656 { 688 {
657 #ifdef DEBUG 689 #ifdef DEBUG
658 std::cout << "qz: peforming qr decomposition of bb" << std::endl; 690 std::cout << "qz: peforming qr decomposition of bb" << std::endl;
659 #endif 691 #endif
660 692
661 // compute the QR factorization of bb 693 // Compute the QR factorization of bb.
662 QR bqr (bb); 694 QR bqr (bb);
663 695
664 #ifdef DEBUG 696 #ifdef DEBUG
665 std::cout << "qz: qr (bb) done; now peforming qz decomposition" << std::endl; 697 std::cout << "qz: qr (bb) done; now peforming qz decomposition" << std::endl;
666 #endif 698 #endif
669 701
670 #ifdef DEBUG 702 #ifdef DEBUG
671 std::cout << "qz: extracted bb" << std::endl; 703 std::cout << "qz: extracted bb" << std::endl;
672 #endif 704 #endif
673 705
674 aa = (bqr.Q ()).transpose ()*aa; 706 aa = (bqr.Q ()).transpose () * aa;
675 707
676 #ifdef DEBUG 708 #ifdef DEBUG
677 std::cout << "qz: updated aa " << std::endl; 709 std::cout << "qz: updated aa " << std::endl;
678 std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl; 710 std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl;
679 711
680 if (compq == 'V') 712 if (compq == 'V')
681 std::cout << "QQ =" << QQ << std::endl; 713 std::cout << "QQ =" << QQ << std::endl;
682 #endif 714 #endif
683 715
684 if (compq == 'V') 716 if (compq == 'V')
685 QQ = QQ*bqr.Q (); 717 QQ = QQ * bqr.Q ();
686 718
687 #ifdef DEBUG 719 #ifdef DEBUG
688 std::cout << "qz: precursors done..." << std::endl; 720 std::cout << "qz: precursors done..." << std::endl;
689 #endif 721 #endif
690 722
691 #ifdef DEBUG 723 #ifdef DEBUG
692 std::cout << "qz: compq = " << compq << ", compz = " << compz << std::endl; 724 std::cout << "qz: compq = " << compq << ", compz = " << compz << std::endl;
693 #endif 725 #endif
694 726
695 // reduce to generalized hessenberg form 727 // Reduce to generalized hessenberg form.
696 F77_XFCN (dgghrd, DGGHRD, 728 F77_XFCN (dgghrd, DGGHRD,
697 (F77_CONST_CHAR_ARG2 (&compq, 1), 729 (F77_CONST_CHAR_ARG2 (&compq, 1),
698 F77_CONST_CHAR_ARG2 (&compz, 1), 730 F77_CONST_CHAR_ARG2 (&compz, 1),
699 nn, ilo, ihi, aa.fortran_vec (), 731 nn, ilo, ihi, aa.fortran_vec (),
700 nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn, 732 nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn,
701 ZZ.fortran_vec (), nn, info 733 ZZ.fortran_vec (), nn, info
702 F77_CHAR_ARG_LEN (1) 734 F77_CHAR_ARG_LEN (1)
703 F77_CHAR_ARG_LEN (1))); 735 F77_CHAR_ARG_LEN (1)));
704 736
705 // check if just computing generalized eigenvalues or if we're 737 // Check if just computing generalized eigenvalues or if we're
706 // actually computing the decomposition 738 // actually computing the decomposition.
707 739
708 // reduce to generalized Schur form 740 // Reduce to generalized Schur form.
709 F77_XFCN (dhgeqz, DHGEQZ, 741 F77_XFCN (dhgeqz, DHGEQZ,
710 (F77_CONST_CHAR_ARG2 (&qz_job, 1), 742 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
711 F77_CONST_CHAR_ARG2 (&compq, 1), 743 F77_CONST_CHAR_ARG2 (&compq, 1),
712 F77_CONST_CHAR_ARG2 (&compz, 1), 744 F77_CONST_CHAR_ARG2 (&compz, 1),
713 nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (), 745 nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (),
715 betar.fortran_vec (), QQ.fortran_vec (), nn, 747 betar.fortran_vec (), QQ.fortran_vec (), nn,
716 ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info 748 ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info
717 F77_CHAR_ARG_LEN (1) 749 F77_CHAR_ARG_LEN (1)
718 F77_CHAR_ARG_LEN (1) 750 F77_CHAR_ARG_LEN (1)
719 F77_CHAR_ARG_LEN (1))); 751 F77_CHAR_ARG_LEN (1)));
752
720 if (compq == 'V') 753 if (compq == 'V')
721 { 754 {
722 F77_XFCN (dggbak, DGGBAK, 755 F77_XFCN (dggbak, DGGBAK,
723 (F77_CONST_CHAR_ARG2 (&bal_job, 1), 756 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
724 F77_CONST_CHAR_ARG2 ("L", 1), 757 F77_CONST_CHAR_ARG2 ("L", 1),
725 nn, ilo, ihi, lscale.data (), rscale.data (), 758 nn, ilo, ihi, lscale.data (), rscale.data (),
726 nn, QQ.fortran_vec (), nn, info 759 nn, QQ.fortran_vec (), nn, info
727 F77_CHAR_ARG_LEN (1) 760 F77_CHAR_ARG_LEN (1)
728 F77_CHAR_ARG_LEN (1))); 761 F77_CHAR_ARG_LEN (1)));
729 762
730 #ifdef DEBUG 763 #ifdef DEBUG
731 if (compq == 'V') 764 if (compq == 'V')
732 std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl; 765 std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl;
733 #endif 766 #endif
734 } 767 }
735 768
736 // then right 769 // then right
737 if (compz == 'V') 770 if (compz == 'V')
738 { 771 {
739 F77_XFCN (dggbak, DGGBAK, 772 F77_XFCN (dggbak, DGGBAK,
740 (F77_CONST_CHAR_ARG2 (&bal_job, 1), 773 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
741 F77_CONST_CHAR_ARG2 ("R", 1), 774 F77_CONST_CHAR_ARG2 ("R", 1),
742 nn, ilo, ihi, lscale.data (), rscale.data (), 775 nn, ilo, ihi, lscale.data (), rscale.data (),
743 nn, ZZ.fortran_vec (), nn, info 776 nn, ZZ.fortran_vec (), nn, info
744 F77_CHAR_ARG_LEN (1) 777 F77_CHAR_ARG_LEN (1)
745 F77_CHAR_ARG_LEN (1))); 778 F77_CHAR_ARG_LEN (1)));
746 779
747 #ifdef DEBUG 780 #ifdef DEBUG
748 if (compz == 'V') 781 if (compz == 'V')
749 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl; 782 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
750 #endif 783 #endif
751 } 784 }
752 785
753 } 786 }
754 787
755 // order the QZ decomposition? 788 // Order the QZ decomposition?
756 if (! (ord_job == 'N' || ord_job == 'n')) 789 if (! (ord_job == 'N' || ord_job == 'n'))
757 { 790 {
758 if (complex_case) 791 if (complex_case)
759 { 792 {
760 // probably not needed, but better be safe 793 // Probably not needed, but better be safe.
761 error ("qz: cannot re-order complex qz decomposition."); 794 error ("qz: cannot re-order complex qz decomposition.");
762 return retval; 795 return retval;
763 } 796 }
764 else 797 else
765 { 798 {
766 #ifdef DEBUG_SORT 799 #ifdef DEBUG_SORT
767 std::cout << "qz: ordering eigenvalues: ord_job = " 800 std::cout << "qz: ordering eigenvalues: ord_job = "
768 << ord_job << std::endl; 801 << ord_job << std::endl;
769 #endif 802 #endif
770 803
771 // declared static to avoid vfork/long jump compiler complaints 804 // Declared static to avoid vfork/long jump compiler complaints.
772 static sort_function sort_test; 805 static sort_function sort_test;
773 sort_test = 0; 806 sort_test = 0;
774 807
775 switch (ord_job) 808 switch (ord_job)
776 { 809 {
791 case '-': 824 case '-':
792 sort_test = &folhp; 825 sort_test = &folhp;
793 break; 826 break;
794 827
795 default: 828 default:
796 // invalid order option (should never happen, since we 829 // Invalid order option (should never happen, since we
797 // checked the options at the top). 830 // checked the options at the top).
798 panic_impossible (); 831 panic_impossible ();
799 break; 832 break;
800 } 833 }
801 834
805 F77_XFCN (xdlange, XDLANGE, 838 F77_XFCN (xdlange, XDLANGE,
806 (F77_CONST_CHAR_ARG2 ("I", 1), 839 (F77_CONST_CHAR_ARG2 ("I", 1),
807 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm 840 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm
808 F77_CHAR_ARG_LEN (1))); 841 F77_CHAR_ARG_LEN (1)));
809 842
810 double eps = DBL_EPSILON*inf_norm*nn; 843 double eps = DBL_EPSILON * inf_norm * nn;
811 844
812 #ifdef DEBUG_SORT 845 #ifdef DEBUG_SORT
813 std::cout << "qz: calling dsubsp: aa=" << std::endl; 846 std::cout << "qz: calling dsubsp: aa=" << std::endl;
814 octave_print_internal (std::cout, aa, 0); 847 octave_print_internal (std::cout, aa, 0);
815 std::cout << std::endl << "bb=" << std::endl; 848 std::cout << std::endl << "bb=" << std::endl;
847 octave_print_internal (std::cout, ZZ, 0); 880 octave_print_internal (std::cout, ZZ, 0);
848 } 881 }
849 std::cout << std::endl; 882 std::cout << std::endl;
850 #endif 883 #endif
851 884
852 // manually update alphar, alphai, betar 885 // Manually update alphar, alphai, betar.
853 static int jj; 886 static int jj;
854 887
855 jj=0; 888 jj = 0;
856 while (jj < nn) 889 while (jj < nn)
857 { 890 {
858 #ifdef DEBUG_EIG 891 #ifdef DEBUG_EIG
859 std::cout << "computing gen eig #" << jj << std::endl; 892 std::cout << "computing gen eig #" << jj << std::endl;
860 #endif 893 #endif
861 894
862 static int zcnt; // number of zeros in this block 895 // Number of zeros in this block.
896 static int zcnt;
863 897
864 if (jj == (nn-1)) 898 if (jj == (nn-1))
865 zcnt = 1; 899 zcnt = 1;
866 else if (aa(jj+1,jj) == 0) 900 else if (aa(jj+1,jj) == 0)
867 zcnt = 1; 901 zcnt = 1;
868 else zcnt = 2; 902 else zcnt = 2;
869 903
870 if (zcnt == 1) // real zero 904 if (zcnt == 1)
871 { 905 {
906 // Real zero.
872 #ifdef DEBUG_EIG 907 #ifdef DEBUG_EIG
873 std::cout << " single gen eig:" << std::endl; 908 std::cout << " single gen eig:" << std::endl;
874 std::cout << " alphar(" << jj << ") = " << aa(jj,jj) << std::endl; 909 std::cout << " alphar(" << jj << ") = " << aa(jj,jj) << std::endl;
875 std::cout << " betar( " << jj << ") = " << bb(jj,jj) << std::endl; 910 std::cout << " betar( " << jj << ") = " << bb(jj,jj) << std::endl;
876 std::cout << " alphai(" << jj << ") = 0" << std::endl; 911 std::cout << " alphai(" << jj << ") = 0" << std::endl;
880 alphai(jj) = 0; 915 alphai(jj) = 0;
881 betar(jj) = bb(jj,jj); 916 betar(jj) = bb(jj,jj);
882 } 917 }
883 else 918 else
884 { 919 {
885 // complex conjugate pair 920 // Complex conjugate pair.
886 #ifdef DEBUG_EIG 921 #ifdef DEBUG_EIG
887 std::cout << "qz: calling dlag2:" << std::endl; 922 std::cout << "qz: calling dlag2:" << std::endl;
888 std::cout << "safmin=" 923 std::cout << "safmin="
889 << setiosflags (std::ios::scientific) << safmin << std::endl; 924 << setiosflags (std::ios::scientific) << safmin << std::endl;
890 925
902 937
903 // FIXME -- probably should be using 938 // FIXME -- probably should be using
904 // fortran_vec instead of &aa(jj,jj) here. 939 // fortran_vec instead of &aa(jj,jj) here.
905 940
906 double scale1, scale2, wr1, wr2, wi; 941 double scale1, scale2, wr1, wr2, wi;
907 const double *aa_ptr = aa.data () + jj*nn+jj; 942 const double *aa_ptr = aa.data () + jj * nn + jj;
908 const double *bb_ptr = bb.data () + jj*nn+jj; 943 const double *bb_ptr = bb.data () + jj * nn + jj;
909 F77_XFCN (dlag2, DLAG2, 944 F77_XFCN (dlag2, DLAG2,
910 (aa_ptr, nn, bb_ptr, nn, safmin, 945 (aa_ptr, nn, bb_ptr, nn, safmin,
911 scale1, scale2, wr1, wr2, wi)); 946 scale1, scale2, wr1, wr2, wi));
912 947
913 #ifdef DEBUG_EIG 948 #ifdef DEBUG_EIG
915 << "\tscale2=" << scale2 << std::endl 950 << "\tscale2=" << scale2 << std::endl
916 << "\twr1=" << wr1 << "\twr2=" << wr2 951 << "\twr1=" << wr1 << "\twr2=" << wr2
917 << "\twi=" << wi << std::endl; 952 << "\twi=" << wi << std::endl;
918 #endif 953 #endif
919 954
920 // just to be safe, check if it's a real pair 955 // Just to be safe, check if it's a real pair.
921 if (wi == 0) 956 if (wi == 0)
922 { 957 {
923 alphar(jj) = wr1; 958 alphar(jj) = wr1;
924 alphai(jj) = 0; 959 alphai(jj) = 0;
925 betar(jj) = scale1; 960 betar(jj) = scale1;
927 alphai(jj+1) = 0; 962 alphai(jj+1) = 0;
928 betar(jj+1) = scale2; 963 betar(jj+1) = scale2;
929 } 964 }
930 else 965 else
931 { 966 {
932 alphar(jj) = alphar(jj+1)=wr1; 967 alphar(jj) = alphar(jj+1) = wr1;
933 alphai(jj) = -(alphai(jj+1) = wi); 968 alphai(jj) = -(alphai(jj+1) = wi);
934 betar(jj) = betar(jj+1) = scale1; 969 betar(jj) = betar(jj+1) = scale1;
935 } 970 }
936 } 971 }
937 972
938 // advance past this block 973 // Advance past this block.
939 jj += zcnt; 974 jj += zcnt;
940 } 975 }
941 976
942 #ifdef DEBUG_SORT 977 #ifdef DEBUG_SORT
943 std::cout << "qz: back from dsubsp: aa=" << std::endl; 978 std::cout << "qz: back from dsubsp: aa=" << std::endl;
961 std::cout << std::endl; 996 std::cout << std::endl;
962 #endif 997 #endif
963 } 998 }
964 } 999 }
965 1000
966 // compute generalized eigenvalues? 1001 // Compute generalized eigenvalues?
967 ComplexColumnVector gev; 1002 ComplexColumnVector gev;
968 1003
969 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) 1004 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4))
970 { 1005 {
971 if (complex_case) 1006 if (complex_case)
972 { 1007 {
973 int cnt = 0; 1008 int cnt = 0;
974 1009
975 for (int ii = 0; ii < nn; ii++) 1010 for (int ii = 0; ii < nn; ii++)
976 // if (cbetar(ii) != 0) 1011 cnt++;
977 cnt++; 1012
978 1013 ComplexColumnVector tmp (cnt);
979 ComplexColumnVector tmp(cnt); 1014
980 1015 cnt = 0;
981 cnt = 0; 1016 for (int ii = 0; ii < nn; ii++)
982 for (int ii = 0; ii < nn; ii++) 1017 tmp(cnt++) = xalpha(ii) / xbeta(ii);
983 // if (cbetar(ii) != 0) 1018
984 tmp(cnt++) = xalpha(ii)/xbeta(ii); 1019 gev = tmp;
985 gev = tmp;
986 } 1020 }
987 else 1021 else
988 { 1022 {
989 #ifdef DEBUG 1023 #ifdef DEBUG
990 std::cout << "qz: computing generalized eigenvalues" << std::endl; 1024 std::cout << "qz: computing generalized eigenvalues" << std::endl;
991 #endif 1025 #endif
992 1026
993 // return finite generalized eigenvalues 1027 // Return finite generalized eigenvalues.
994 int cnt = 0; 1028 int cnt = 0;
995 1029
996 for (int ii = 0; ii < nn; ii++) 1030 for (int ii = 0; ii < nn; ii++)
997 if (betar(ii) != 0) 1031 if (betar(ii) != 0)
998 cnt++; 1032 cnt++;
999 1033
1000 ComplexColumnVector tmp(cnt); 1034 ComplexColumnVector tmp (cnt);
1001 1035
1002 cnt = 0; 1036 cnt = 0;
1003 for (int ii = 0; ii < nn; ii++) 1037 for (int ii = 0; ii < nn; ii++)
1004 if (betar(ii) != 0) 1038 if (betar(ii) != 0)
1005 tmp(cnt++) = Complex(alphar(ii), alphai(ii))/betar(ii); 1039 tmp(cnt++) = Complex(alphar(ii), alphai(ii))/betar(ii);
1040
1006 gev = tmp; 1041 gev = tmp;
1007 } 1042 }
1008 } 1043 }
1009 1044
1010 // right, left eigenvector matrices 1045 // Right, left eigenvector matrices.
1011 if (nargout >= 5) 1046 if (nargout >= 5)
1012 { 1047 {
1013 char side = (nargout == 5 ? 'R' : 'B'); // which side to compute? 1048 // Which side to compute?
1014 char howmny = 'B'; // compute all of them and backtransform 1049 char side = (nargout == 5 ? 'R' : 'B');
1015 octave_idx_type *select = 0; // dummy pointer; select is not used. 1050 // Compute all of them and backtransform
1051 char howmny = 'B';
1052 // Dummy pointer; select is not used.
1053 octave_idx_type *select = 0;
1016 1054
1017 if (complex_case) 1055 if (complex_case)
1018 { 1056 {
1057 CVL = CQ;
1058 CVR = CZ;
1059 ComplexRowVector cwork2 (2 * nn);
1060 RowVector rwork2 (8 * nn);
1019 octave_idx_type m; 1061 octave_idx_type m;
1020 CVL=CQ; 1062
1021 CVR=CZ; 1063 F77_XFCN (ztgevc, ZTGEVC,
1022 ComplexRowVector cwork2(2*nn); 1064 (F77_CONST_CHAR_ARG2 (&side, 1),
1023 RowVector rwork2(8*nn); 1065 F77_CONST_CHAR_ARG2 (&howmny, 1),
1024 1066 select, nn, caa.fortran_vec (), nn, cbb.fortran_vec (),
1025 //octave_idx_type n=nn; 1067 nn, CVL.fortran_vec (), nn, CVR.fortran_vec (), nn, nn,
1026 F77_XFCN (ztgevc, ZTGEVC, 1068 m, cwork2.fortran_vec (), rwork2.fortran_vec (), info
1027 (F77_CONST_CHAR_ARG2 (&side, 1), 1069 F77_CHAR_ARG_LEN (1)
1028 F77_CONST_CHAR_ARG2 (&howmny, 1), 1070 F77_CHAR_ARG_LEN (1)));
1029 select, nn, caa.fortran_vec (), nn, cbb.fortran_vec (),
1030 nn, CVL.fortran_vec (), nn, CVR.fortran_vec (), nn, nn,
1031 m, cwork2.fortran_vec (), rwork2.fortran_vec (), info
1032 F77_CHAR_ARG_LEN (1)
1033 F77_CHAR_ARG_LEN (1)));
1034 } 1071 }
1035 else 1072 else
1036 { 1073 {
1037 #ifdef DEBUG 1074 #ifdef DEBUG
1038 std::cout << "qz: computing generalized eigenvectors" << std::endl; 1075 std::cout << "qz: computing generalized eigenvectors" << std::endl;
1039 #endif 1076 #endif
1040 1077
1041 VL = QQ; 1078 VL = QQ;
1042 VR = ZZ; 1079 VR = ZZ;
1043 octave_idx_type m; 1080 octave_idx_type m;
1081
1044 F77_XFCN (dtgevc, DTGEVC, 1082 F77_XFCN (dtgevc, DTGEVC,
1045 (F77_CONST_CHAR_ARG2 (&side, 1), 1083 (F77_CONST_CHAR_ARG2 (&side, 1),
1046 F77_CONST_CHAR_ARG2 (&howmny, 1), 1084 F77_CONST_CHAR_ARG2 (&howmny, 1),
1047 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (), 1085 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (),
1048 nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn, 1086 nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn,
1049 m, work.fortran_vec (), info 1087 m, work.fortran_vec (), info
1050 F77_CHAR_ARG_LEN (1) 1088 F77_CHAR_ARG_LEN (1)
1051 F77_CHAR_ARG_LEN (1))); 1089 F77_CHAR_ARG_LEN (1)));
1052 1090
1053 // now construct the complex form of VV, WW 1091 // Now construct the complex form of VV, WW.
1054 int jj = 0; 1092 int jj = 0;
1055 1093
1056 while (jj < nn) 1094 while (jj < nn)
1057 { 1095 {
1058 OCTAVE_QUIT; 1096 OCTAVE_QUIT;
1059 1097
1060 // see if real or complex eigenvalue 1098 // See if real or complex eigenvalue.
1061 int cinc = 2; // column increment; assume complex eigenvalue 1099
1100 // Column increment; assume complex eigenvalue.
1101 int cinc = 2;
1062 1102
1063 if (jj == (nn-1)) 1103 if (jj == (nn-1))
1064 cinc = 1; // single column 1104 // Single column.
1105 cinc = 1;
1065 else if (aa(jj+1,jj) == 0) 1106 else if (aa(jj+1,jj) == 0)
1066 cinc = 1; 1107 cinc = 1;
1067 1108
1068 // now copy the eigenvector (s) to CVR, CVL 1109 // Now copy the eigenvector (s) to CVR, CVL.
1069 if (cinc == 1) 1110 if (cinc == 1)
1070 { 1111 {
1071 for (int ii = 0; ii < nn; ii++) 1112 for (int ii = 0; ii < nn; ii++)
1072 CVR(ii,jj) = VR(ii,jj); 1113 CVR(ii,jj) = VR(ii,jj);
1073 1114
1075 for (int ii = 0; ii < nn; ii++) 1116 for (int ii = 0; ii < nn; ii++)
1076 CVL(ii,jj) = VL(ii,jj); 1117 CVL(ii,jj) = VL(ii,jj);
1077 } 1118 }
1078 else 1119 else
1079 { 1120 {
1080 // double column; complex vector 1121 // Double column; complex vector.
1081 1122
1082 for (int ii = 0; ii < nn; ii++) 1123 for (int ii = 0; ii < nn; ii++)
1083 { 1124 {
1084 CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1)); 1125 CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1));
1085 CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1)); 1126 CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1));
1091 CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1)); 1132 CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1));
1092 CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1)); 1133 CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1));
1093 } 1134 }
1094 } 1135 }
1095 1136
1096 // advance to next eigenvectors (if any) 1137 // Advance to next eigenvectors (if any).
1097 jj += cinc; 1138 jj += cinc;
1098 } 1139 }
1099 } 1140 }
1100 } 1141 }
1101 1142
1102 switch (nargout) 1143 switch (nargout)
1103 { 1144 {
1104 case 7: 1145 case 7:
1105 retval(6) = gev; 1146 retval(6) = gev;
1106 1147
1107 case 6: // return eigenvectors 1148 case 6:
1149 // Return eigenvectors.
1108 retval(5) = CVL; 1150 retval(5) = CVL;
1109 1151
1110 case 5: // return eigenvectors 1152 case 5:
1153 // Return eigenvectors.
1111 retval(4) = CVR; 1154 retval(4) = CVR;
1112 1155
1113 case 4: 1156 case 4:
1114 if (nargin == 3) 1157 if (nargin == 3)
1115 { 1158 {
1119 std::cout << std::endl; 1162 std::cout << std::endl;
1120 #endif 1163 #endif
1121 retval(3) = gev; 1164 retval(3) = gev;
1122 } 1165 }
1123 else 1166 else
1124 {if (complex_case) 1167 {
1125 retval(3) = CZ; 1168 if (complex_case)
1126 else 1169 retval(3) = CZ;
1127 retval(3) = ZZ; 1170 else
1128 } 1171 retval(3) = ZZ;
1172 }
1129 1173
1130 case 3: 1174 case 3:
1131 if (nargin == 3) 1175 if (nargin == 3)
1132 retval(2) = CZ; 1176 retval(2) = CZ;
1133 else 1177 else
1134 {if (complex_case) 1178 {
1135 retval(2) = CQ.hermitian(); // compabible with MATLAB output 1179 if (complex_case)
1136 else 1180 retval(2) = CQ.hermitian ();
1137 retval(2) = QQ.transpose(); 1181 else
1138 } 1182 retval(2) = QQ.transpose ();
1183 }
1184
1139 case 2: 1185 case 2:
1140 {if (complex_case)
1141 { 1186 {
1142 #ifdef DEBUG 1187 if (complex_case)
1143 std::cout << "qz: retval (1) = cbb = " << std::endl; 1188 {
1144 octave_print_internal (std::cout, cbb, 0); 1189 #ifdef DEBUG
1145 std::cout << std::endl << "qz: retval(0) = caa = " <<std::endl; 1190 std::cout << "qz: retval (1) = cbb = " << std::endl;
1146 octave_print_internal (std::cout, caa, 0); 1191 octave_print_internal (std::cout, cbb, 0);
1147 std::cout << std::endl; 1192 std::cout << std::endl << "qz: retval(0) = caa = " <<std::endl;
1148 #endif 1193 octave_print_internal (std::cout, caa, 0);
1149 retval(1) = cbb; 1194 std::cout << std::endl;
1150 retval(0) = caa; 1195 #endif
1151 } 1196 retval(1) = cbb;
1197 retval(0) = caa;
1198 }
1152 else 1199 else
1153 { // real case 1200 {
1154 #ifdef DEBUG 1201 #ifdef DEBUG
1155 std::cout << "qz: retval (1) = bb = " << std::endl; 1202 std::cout << "qz: retval (1) = bb = " << std::endl;
1156 octave_print_internal (std::cout, bb, 0); 1203 octave_print_internal (std::cout, bb, 0);
1157 std::cout << std::endl << "qz: retval(0) = aa = " <<std::endl; 1204 std::cout << std::endl << "qz: retval(0) = aa = " <<std::endl;
1158 octave_print_internal (std::cout, aa, 0); 1205 octave_print_internal (std::cout, aa, 0);
1159 std::cout << std::endl; 1206 std::cout << std::endl;
1160 #endif 1207 #endif
1161 retval(1) = bb; 1208 retval(1) = bb;
1162 retval(0) = aa; 1209 retval(0) = aa;
1163 } 1210 }
1164 break;} 1211 }
1212 break;
1165 1213
1166 1214
1167 case 1: 1215 case 1:
1168 case 0: 1216 case 0:
1169 #ifdef DEBUG 1217 #ifdef DEBUG