Mercurial > hg > octave-lyh

comparison src/DLD-FUNCTIONS/qz.cc @ 3185:9580887dd160

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1998-09-26 02:45:55 by jwe]
author jwe
date Sat, 26 Sep 1998 02:45:59 +0000
parents 5edc539c2f80
children 3d3eca53ecca
comparison
equal deleted inserted replaced
3184:3988763ec9d3 3185:9580887dd160
45 #include "help.h" 45 #include "help.h"
46 #include "oct-obj.h" 46 #include "oct-obj.h"
47 #include "oct-map.h" 47 #include "oct-map.h"
48 #include "ov.h" 48 #include "ov.h"
49 #include "pager.h" 49 #include "pager.h"
50 #if defined(DEBUG) || defined(DEBUG_SORT) 50 #if defined (DEBUG) || defined (DEBUG_SORT)
51 #include "pr-output.h" 51 #include "pr-output.h"
52 #endif 52 #endif
53 #include "symtab.h" 53 #include "symtab.h"
54 #include "utils.h" 54 #include "utils.h"
55 #include "variables.h" 55 #include "variables.h"
56 56
57 typedef int (*sort_function) (const int& LSIZE, const double& ALPHA, 57 typedef int (*sort_function) (const int& LSIZE, const double& ALPHA,
58 const double& BETA, const double& S, const double& P); 58 const double& BETA, const double& S,
59 const double& P);
59 60
60 extern "C" 61 extern "C"
61 { 62 {
62 int F77_FCN( dggbal, DGGBAL) (const char* JOB, const int& N, 63 int F77_FCN (dggbal, DGGBAL) (const char* JOB, const int& N,
63 double* A, const int& LDA, double* B, const int& LDB, 64 double* A, const int& LDA, double* B,
64 int& ILO, int & IHI, double* LSCALE, 65 const int& LDB, int& ILO, int& IHI,
65 double* RSCALE, double* WORK, int& INFO, long ); 66 double* LSCALE, double* RSCALE,
66 67 double* WORK, int& INFO, long);
67 int F77_FCN( dggbak, DGGBAK) (const char* JOB, const char* SIDE, 68
68 const int& N, const int& ILO, const int& IHI, 69 int F77_FCN (dggbak, DGGBAK) (const char* JOB, const char* SIDE,
69 double* LSCALE, double* RSCALE, int& M, 70 const int& N, const int& ILO,
70 double* V, const int& LDV, int& INFO, long, long); 71 const int& IHI, double* LSCALE,
71 72 double* RSCALE, int& M, double* V,
72 int F77_FCN( dgghrd, DGGHRD) ( const char* COMPQ, const char* COMPZ, 73 const int& LDV, int& INFO, long, long);
73 const int& N, const int& ILO, const int& IHI, double* A, const int& LDA, 74
74 double* B, const int& LDB, double* Q, const int& LDQ, double* Z, 75 int F77_FCN (dgghrd, DGGHRD) (const char* COMPQ, const char* COMPZ,
75 const int& LDZ, int& INFO, const long, const long); 76 const int& N, const int& ILO,
76 77 const int& IHI, double* A,
77 int F77_FCN( dhgeqz, DHGEQZ) ( const char* JOB, const char* COMPQ, 78 const int& LDA, double* B,
78 const char* COMPZ, const int& N, const int& ILO, const int& IHI, 79 const int& LDB, double* Q,
79 double* A, const int& LDA, double* B, const int& LDB, 80 const int& LDQ, double* Z,
80 double* ALPHAR, double* ALPHAI, double* BETA, double* Q, 81 const int& LDZ, int& INFO, long, long);
81 const int& LDQ, double* Z, const int& LDZ, double* WORK, 82
82 const int& LWORK, int& INFO, const long, const long, const long ); 83 int F77_FCN (dhgeqz, DHGEQZ) (const char* JOB, const char* COMPQ,
83 84 const char* COMPZ, const int& N,
84 int F77_FCN( dlag2, DLAG2) ( double* A, const int& LDA, double* B, 85 const int& ILO, const int& IHI,
85 const int& LDB, const double& SAFMIN, double& SCALE1, double& SCALE2, 86 double* A, const int& LDA, double* B,
86 double& WR1, double& WR2, double& WI ); 87 const int& LDB, double* ALPHAR,
88 double* ALPHAI, double* BETA, double* Q,
89 const int& LDQ, double* Z,
90 const int& LDZ, double* WORK,
91 const int& LWORK, int& INFO,
92 long, long, long);
93
94 int F77_FCN (dlag2, DLAG2) (double* A, const int& LDA, double* B,
95 const int& LDB, const double& SAFMIN,
96 double& SCALE1, double& SCALE2,
97 double& WR1, double& WR2, double& WI);
87 98
88 // Van Dooren's code (netlib.org: toms/590) for reordering 99 // Van Dooren's code (netlib.org: toms/590) for reordering
89 // GEP. Only processes Z, not Q. 100 // GEP. Only processes Z, not Q.
90 int F77_FCN( dsubsp, DSUBSP) ( const int& NMAX, const int& N, double* A, 101 int F77_FCN (dsubsp, DSUBSP) (const int& NMAX, const int& N, double* A,
91 double* B, double* Z, sort_function, const double& EPS, 102 double* B, double* Z, sort_function,
92 int& NDIM, int& FAIL, int* IND); 103 const double& EPS, int& NDIM, int& FAIL,
104 int* IND);
93 105
94 // documentation for DTGEVC incorrectly states that VR, VL are 106 // documentation for DTGEVC incorrectly states that VR, VL are
95 // complex*16; they are declared in DTGEVC as double precision 107 // complex*16; they are declared in DTGEVC as double precision
96 // (probably a cut and paste problem fro ZTGEVC) 108 // (probably a cut and paste problem fro ZTGEVC)
97 int F77_FCN( dtgevc, DTGEVC) ( const char* SIDE, const char* HOWMNY, 109 int F77_FCN (dtgevc, DTGEVC) (const char* SIDE, const char* HOWMNY,
98 int* SELECT, const int& N, double* A, const int& LDA, double* B, 110 int* SELECT, const int& N, double* A,
99 const int& LDB, double* VL, const int& LDVL, double* VR, 111 const int& LDA, double* B,
100 const int& LDVR, const int& MM, int& M, double* WORK, int& INFO, 112 const int& LDB, double* VL,
101 long, long ); 113 const int& LDVL, double* VR,
102 114 const int& LDVR, const int& MM,
103 int F77_FCN ( xdlamch, XDLAMCH) (const char* cmach, double& retval, long); 115 int& M, double* WORK, int& INFO,
104 int F77_FCN ( xdlange, XDLANGE) (const char*, const int&, 116 long, long);
117
118 int F77_FCN (xdlamch, XDLAMCH) (const char* cmach, double& retval, long);
119
120 int F77_FCN (xdlange, XDLANGE) (const char*, const int&,
105 const int&, const double*, 121 const int&, const double*,
106 const int&, double*, double&); 122 const int&, double*, double&);
107 } 123 }
108 124
109 // fcrhp, fin, fout, folhp: 125 // fcrhp, fin, fout, folhp:
112 // fin: |lambda| < 1 128 // fin: |lambda| < 1
113 // fout: |lambda| >= 1 129 // fout: |lambda| >= 1
114 // fcrhp: real(lambda) >= 0 130 // fcrhp: real(lambda) >= 0
115 // folhp: real(lambda) < 0 131 // folhp: real(lambda) < 0
116 132
117 static int fcrhp(const int& lsize, const double& alpha, 133 static int
118 const double& beta, const double& s, const double& p) 134 fcrhp (const int& lsize, const double& alpha,
135 const double& beta, const double& s, const double&)
119 { 136 {
120 if(lsize == 1) 137 if (lsize == 1)
121 return (alpha*beta >= 0 ? 1 : -1); 138 return (alpha*beta >= 0 ? 1 : -1);
122 else 139 else
123 return (s >= 0 ? 1 : -1); 140 return (s >= 0 ? 1 : -1);
124 } 141 }
125 static int fin(const int& lsize, const double& alpha, 142
126 const double& beta, const double& s, const double& p) 143 static int
144 fin (const int& lsize, const double& alpha,
145 const double& beta, const double&, const double& p)
127 { 146 {
128 int retval; 147 int retval;
129 if(lsize == 1) 148
130 retval = (fabs(alpha) < fabs(beta) ? 1 : -1); 149 if (lsize == 1)
131 else 150 retval = (fabs (alpha) < fabs (beta) ? 1 : -1);
132 retval = (fabs(p) < 1 ? 1 : -1); 151 else
133 152 retval = (fabs (p) < 1 ? 1 : -1);
134 #ifdef DEBUG 153
154 #ifdef DEBUG
135 cout << "qz: fin: retval=" << retval << endl; 155 cout << "qz: fin: retval=" << retval << endl;
136 #endif 156 #endif
157
137 return retval; 158 return retval;
138 } 159 }
139 static int folhp(const int& lsize, const double& alpha, 160
140 const double& beta, const double& s, const double& p) 161 static int
162 folhp (const int& lsize, const double& alpha,
163 const double& beta, const double& s, const double&)
141 { 164 {
142 if(lsize == 1) 165 if (lsize == 1)
143 return (alpha*beta < 0 ? 1 : -1); 166 return (alpha*beta < 0 ? 1 : -1);
144 else 167 else
145 return (s < 0 ? 1 : -1); 168 return (s < 0 ? 1 : -1);
146 } 169 }
147 static int fout(const int& lsize, const double& alpha, 170
148 const double& beta, const double& s, const double& p) 171 static int
172 fout (const int& lsize, const double& alpha,
173 const double& beta, const double&, const double& p)
149 { 174 {
150 if(lsize == 1) 175 if (lsize == 1)
151 return (fabs(alpha) >= fabs(beta) ? 1 : -1); 176 return (fabs (alpha) >= fabs (beta) ? 1 : -1);
152 else 177 else
153 return (fabs(p) >= 1 ? 1 : -1); 178 return (fabs (p) >= 1 ? 1 : -1);
154 } 179 }
155 180
156 DEFUN_DLD (qz, args, nargout, 181 DEFUN_DLD (qz, args, nargout,
157 "Usage: lambda = qz(A,B) form [1]\n\ 182 "Usage:\n\
158 [AA,BB,Q,Z{,V,W,lambda}] = qz(A,B) form [2]\n\ 183
159 [AA,BB,Z{,lambda}] = qz(A,B,opt) form [3]\n\ 184 lambda = qz (A, B) form [1]\n\
185 [AA, BB, Q, Z {, V, W, lambda}] = qz (A, B) form [2]\n\
186 [AA, BB, Z{, lambda}] = qz (A, B, opt) form [3]\n\
187 \n\
160 Generalized eigenvalue problem A x = s B x \n\ 188 Generalized eigenvalue problem A x = s B x \n\
161 189 \n\
162 Form [1]: Computes the generalized eigenvalues lambda of (A - sB).\n\ 190 Form [1]: Computes the generalized eigenvalues lambda of (A - sB).\n\
191 \n\
163 Form [2]: Computes qz decomposition, generalized eigenvectors, and \n\ 192 Form [2]: Computes qz decomposition, generalized eigenvectors, and \n\
164 generalized eigenvalues of (A - sB)\n\ 193 generalized eigenvalues of (A - sB)\n\
165 A V = B V diag(lambda)\n\ 194 A V = B V diag (lambda)\n\
166 W' A = diag(lambda) W' B\n\ 195 W' A = diag (lambda) W' B\n\
167 AA = Q'*A*Z, BB = Q'*B*Z with Q, Z orthogonal (unitary)= I\n\ 196 AA = Q'*A*Z, BB = Q'*B*Z with Q, Z orthogonal (unitary)= I\n\
197 \n\
168 Form [3]: As in form [2], but allows ordering of generalized eigenpairs\n\ 198 Form [3]: As in form [2], but allows ordering of generalized eigenpairs\n\
169 for (e.g.) solution of discrete time algebraic Riccati equations.\n\ 199 for (e.g.) solution of discrete time algebraic Riccati equations.\n\
170 Form 3 is not available for complex matrices and does not compute\n\ 200 Form 3 is not available for complex matrices and does not compute\n\
171 the generalized eigenvectors V, W, nor the orthogonal matrix Q.\n\ 201 the generalized eigenvectors V, W, nor the orthogonal matrix Q.\n\
172 \n\ 202 \n\
173 opt: for ordering eigenvalues of the GEP pencil. The leading block\n\ 203 opt: for ordering eigenvalues of the GEP pencil. The leading block\n\
174 of the revised pencil contains all eigenvalues that satisfy:\n\ 204 of the revised pencil contains all eigenvalues that satisfy:\n\
205 \n\
175 \"N\" = unordered (default) \n\ 206 \"N\" = unordered (default) \n\
176 \"S\" = small: leading block has all |lambda| <=1 \n\ 207 \"S\" = small: leading block has all |lambda| <=1 \n\
177 \"B\" = big: leading block has all |lambda >= 1 \n\ 208 \"B\" = big: leading block has all |lambda >= 1 \n\
178 \"-\" = negative real part: leading block has all eigenvalues\n\ 209 \"-\" = negative real part: leading block has all eigenvalues\n\
179 in the open left half-plant\n\ 210 in the open left half-plant\n\
181 in the closed right half-plane\n\ 212 in the closed right half-plane\n\
182 \n\ 213 \n\
183 Note: Permutation balancing is performed, but not scaling (see balance)\n\ 214 Note: Permutation balancing is performed, but not scaling (see balance)\n\
184 Order of output arguments was selected for compatibility with MATLAB\n\ 215 Order of output arguments was selected for compatibility with MATLAB\n\
185 \n\ 216 \n\
186 See also: balance, dare, eig, schur\n") 217 See also: balance, dare, eig, schur")
187 { 218 {
188 octave_value_list retval; 219 octave_value_list retval;
189 int nargin = args.length (); 220 int nargin = args.length ();
190 221
191 #ifdef DEBUG 222 #ifdef DEBUG
192 cout << "qz: nargin = " << nargin << ", nargout = " << nargout << endl; 223 cout << "qz: nargin = " << nargin << ", nargout = " << nargout << endl;
193 #endif 224 #endif
194 225
195 if (nargin < 2 || nargin > 3 || nargout > 7 ) 226 if (nargin < 2 || nargin > 3 || nargout > 7)
196 { 227 {
197 print_usage ("qz"); 228 print_usage ("qz");
198 return retval; 229 return retval;
199 } 230 }
200 else if(nargin == 3 && (nargout < 3 || nargout > 4)) 231 else if (nargin == 3 && (nargout < 3 || nargout > 4))
201 { 232 {
202 error("qz: Illegal number of output arguments for form [3] call"); 233 error ("qz: Illegal number of output arguments for form [3] call");
203 } 234 return retval;
204 235 }
205 #ifdef DEBUG 236
237 #ifdef DEBUG
206 cout << "qz: determine ordering option" << endl; 238 cout << "qz: determine ordering option" << endl;
207 #endif 239 #endif
208 240
209 // Determine ordering option 241 // Determine ordering option
210 string ord_job; 242 string ord_job;
211 static double safmin; 243 static double safmin;
212 if(nargin == 2) 244
245 if (nargin == 2)
213 ord_job = "N"; 246 ord_job = "N";
214 else if( !args(2).is_string() ) 247 else if (!args(2).is_string ())
215 error("qz: argument 3 must be a string"); 248 {
216 else 249 error ("qz: argument 3 must be a string");
217 { 250 return retval;
218 ord_job = args(2).string_value(); 251 }
219 if(ord_job[0] != 'N' && ord_job[0] != 'S' && ord_job[0] != 'B' 252 else
220 && ord_job[0] != '+' && ord_job[0] != '-') 253 {
221 error("qz: illegal order option"); 254 ord_job = args(2).string_value ();
222 255
223 // overflow constant required by dlag2 256 if (ord_job[0] != 'N'
224 F77_XFCN ( xdlamch, XDLAMCH, ("S", safmin, 1L)); 257 && ord_job[0] != 'S'
225 258 && ord_job[0] != 'B'
226 #ifdef DEBUG_EIG 259 && ord_job[0] != '+'
227 cout << "qz: initial value of safmin=" << setiosflags(ios::scientific) 260 && ord_job[0] != '-')
228 << safmin << endl; 261 {
229 #endif 262 error ("qz: illegal order option");
230 263 return retval;
231 // some machines (e.g., DEC alpha) get safmin = 0; 264 }
232 // for these, use eps instead to avoid problems in dlag2 265
233 if(safmin == 0) 266 // overflow constant required by dlag2
234 { 267 F77_FCN (xdlamch, XDLAMCH) ("S", safmin, 1L);
235 #ifdef DEBUG_EIG 268
236 cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << endl; 269 #ifdef DEBUG_EIG
237 #endif 270 cout << "qz: initial value of safmin=" << setiosflags (ios::scientific)
238 271 << safmin << endl;
239 F77_XFCN ( xdlamch, XDLAMCH, ("E", safmin, 1L)); 272 #endif
240 273
241 #ifdef DEBUG_EIG 274 // some machines (e.g., DEC alpha) get safmin = 0;
242 cout << "qz: safmin set to " << setiosflags(ios::scientific) 275 // for these, use eps instead to avoid problems in dlag2
243 << safmin << endl; 276 if (safmin == 0)
244 #endif 277 {
245 } 278 #ifdef DEBUG_EIG
246 } 279 cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << endl;
247 280 #endif
248 #ifdef DEBUG 281
282 F77_FCN (xdlamch, XDLAMCH) ("E", safmin, 1L);
283
284 #ifdef DEBUG_EIG
285 cout << "qz: safmin set to " << setiosflags (ios::scientific)
286 << safmin << endl;
287 #endif
288 }
289 }
290
291 #ifdef DEBUG
249 cout << "qz: check argument 1" << endl; 292 cout << "qz: check argument 1" << endl;
250 #endif 293 #endif
251 294
252 // Argument 1: check if it's o.k. dimensioned 295 // Argument 1: check if it's o.k. dimensioned
253 int nn = args(0).rows(); 296 int nn = args(0).rows ();
254 297
255 #ifdef DEBUG 298 #ifdef DEBUG
256 cout << "argument 1 dimensions: (" << nn << "," << args(0).columns() << ")" 299 cout << "argument 1 dimensions: (" << nn << "," << args(0).columns () << ")"
257 << endl; 300 << endl;
258 #endif 301 #endif
259 int arg_is_empty = empty_arg ("qz", nn, args(0).columns()); 302
303 int arg_is_empty = empty_arg ("qz", nn, args(0).columns ());
304
260 if (arg_is_empty < 0) 305 if (arg_is_empty < 0)
261 { 306 {
262 gripe_empty_arg("qz: parameter 1",0); 307 gripe_empty_arg ("qz: parameter 1", 0);
263 return retval; 308 return retval;
264 } 309 }
265 else if (arg_is_empty > 0) 310 else if (arg_is_empty > 0)
266 { 311 {
267 gripe_empty_arg("qz: parameter 1; continuing",0); 312 gripe_empty_arg ("qz: parameter 1; continuing", 0);
268 return octave_value_list (2, Matrix ()); 313 return octave_value_list (2, Matrix ());
269 } 314 }
270 else if (args(0).columns() != nn) 315 else if (args(0).columns () != nn)
271 { 316 {
272 gripe_square_matrix_required ("qz"); 317 gripe_square_matrix_required ("qz");
273 return retval; 318 return retval;
274 } 319 }
275 320
276 // Argument 1: dimensions look good; get the value 321 // Argument 1: dimensions look good; get the value
277 Matrix aa; 322 Matrix aa;
278 ComplexMatrix caa; 323 ComplexMatrix caa;
279 if (args(0).is_complex_type ()) 324
325 if (args(0).is_complex_type ())
280 caa = args(0).complex_matrix_value (); 326 caa = args(0).complex_matrix_value ();
281 else 327 else
282 aa = args(0).matrix_value (); 328 aa = args(0).matrix_value ();
283 if (error_state) 329
330 if (error_state)
284 return retval; 331 return retval;
285 332
286 #ifdef DEBUG 333 #ifdef DEBUG
287 cout << "qz: check argument 2" << endl; 334 cout << "qz: check argument 2" << endl;
288 #endif 335 #endif
289 336
290 // Extract argument 2 (bb, or cbb if complex) 337 // Extract argument 2 (bb, or cbb if complex)
291 if( (nn != args(1).columns()) || (nn != args(1).rows() )) 338 if ((nn != args(1).columns ()) || (nn != args(1).rows ()))
292 { 339 {
293 gripe_nonconformant (); 340 gripe_nonconformant ();
294 return retval; 341 return retval;
295 } 342 }
343
296 Matrix bb; 344 Matrix bb;
297 ComplexMatrix cbb; 345 ComplexMatrix cbb;
298 if (args(1).is_complex_type ()) 346
347 if (args(1).is_complex_type ())
299 cbb = args(1).complex_matrix_value (); 348 cbb = args(1).complex_matrix_value ();
300 else 349 else
301 bb = args(1).matrix_value (); 350 bb = args(1).matrix_value ();
302 if (error_state) 351
352 if (error_state)
303 return retval; 353 return retval;
304 354
305 // Both matrices loaded, now let's check what kind of arithmetic: 355 // Both matrices loaded, now let's check what kind of arithmetic:
306 //declared static to avoid compiler warnings about long jumps, vforks. 356 //declared static to avoid compiler warnings about long jumps, vforks.
307 static int complex_case 357
308 = (args(0).is_complex_type() || args(1).is_complex_type()); 358 static int complex_case
309 359 = (args(0).is_complex_type () || args(1).is_complex_type ());
310 if(nargin == 3 && complex_case) 360
311 error("qz: cannot re-order complex qz decomposition."); 361 if (nargin == 3 && complex_case)
362 {
363 error ("qz: cannot re-order complex qz decomposition.");
364 return retval;
365 }
312 366
313 // first, declare variables used in both the real and complex case 367 // first, declare variables used in both the real and complex case
314 Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn); 368 Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn);
315 RowVector alphar(nn), alphai(nn), betar(nn); 369 RowVector alphar(nn), alphai(nn), betar(nn);
316 370
317 ComplexMatrix CQ(nn,nn), CZ(nn,nn),CVR(nn,nn),CVL(nn,nn); 371 ComplexMatrix CQ(nn,nn), CZ(nn,nn), CVR(nn,nn), CVL(nn,nn);
318 int ilo, ihi, info; 372 int ilo, ihi, info;
319 char compq = (nargout >= 3 ? 'V' : 'N'), 373 char compq = (nargout >= 3 ? 'V' : 'N');
320 compz = (nargout >= 4 ? 'V' : 'N'); 374 char compz = (nargout >= 4 ? 'V' : 'N');
321 375
322 // initialize Q,Z to identity if we need either of them 376 // initialize Q, Z to identity if we need either of them
323 if(compq == 'V' || compz == 'V') 377 if (compq == 'V' || compz == 'V')
324 for(int ii=0; ii < nn ; ii++) 378 for (int ii = 0; ii < nn; ii++)
325 for( int jj=0; jj < nn ; jj++) 379 for (int jj = 0; jj < nn; jj++)
326 QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0); 380 QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0);
327 381
328 // always perform permutation balancing 382 // always perform permutation balancing
329 char bal_job = 'P'; 383 char bal_job = 'P';
330 RowVector lscale(nn), rscale(nn), work(6*nn); 384 RowVector lscale(nn), rscale(nn), work(6*nn);
331 385
332 if(complex_case) 386 if (complex_case)
333 error("Complex case not implemented yet"); 387 {
388 error ("Complex case not implemented yet");
389 return retval;
390 }
334 else 391 else
335 { 392 {
336 #ifdef DEBUG 393 #ifdef DEBUG
337 if(compq == 'V') 394 if (compq == 'V')
338 cout << "qz: performing balancing; QQ=" << endl << QQ << endl; 395 cout << "qz: performing balancing; QQ=" << endl << QQ << endl;
339 #endif 396 #endif
340 397
341 F77_XFCN( dggbal, DGGBAL, (&bal_job, nn, aa.fortran_vec(), 398 F77_XFCN (dggbal, DGGBAL,
342 nn, bb.fortran_vec() , nn, ilo, ihi, lscale.fortran_vec(), 399 (&bal_job, nn, aa.fortran_vec(), nn, bb.fortran_vec(),
343 rscale.fortran_vec(), work.fortran_vec(), info , 1L)); 400 nn, ilo, ihi, lscale.fortran_vec(),
344 if(f77_exception_encountered) 401 rscale.fortran_vec(), work.fortran_vec(), info, 1L));
345 (*current_liboctave_error_handler) ("unrecoverable error in qz(bal)"); 402
346 } 403 if (f77_exception_encountered)
404 {
405 error ("unrecoverable error in qz (bal)");
406 return retval;
407 }
408 }
347 409
348 // Since we just want the balancing matrices, we can use dggbal 410 // Since we just want the balancing matrices, we can use dggbal
349 // for both the real and complex cases; 411 // for both the real and complex cases;
350 // left first 412 // left first
351 if(compq == 'V') 413
352 { 414 if (compq == 'V')
353 F77_XFCN( dggbak, DGGBAK, (&bal_job, "L", 415 {
354 nn, ilo, ihi, lscale.fortran_vec(), 416 F77_XFCN (dggbak, DGGBAK,
355 rscale.fortran_vec(), nn, QQ.fortran_vec(), 417 (&bal_job, "L", nn, ilo, ihi, lscale.fortran_vec(),
356 nn, info, 1L, 1L)); 418 rscale.fortran_vec(), nn, QQ.fortran_vec(),
357 419 nn, info, 1L, 1L));
358 #ifdef DEBUG 420
359 if(compq == 'V') cout << "qz: balancing done; QQ=" << endl << QQ << endl; 421 #ifdef DEBUG
360 #endif 422 if (compq == 'V')
361 423 cout << "qz: balancing done; QQ=" << endl << QQ << endl;
362 if(f77_exception_encountered) 424 #endif
363 (*current_liboctave_error_handler) ("unrecoverable error in qz(bal-L)"); 425
426 if (f77_exception_encountered)
427 {
428 error ("unrecoverable error in qz (bal-L)");
429 return retval;
430 }
364 } 431 }
365 432
366 // then right 433 // then right
367 if(compz == 'V') 434 if (compz == 'V')
368 { 435 {
369 F77_XFCN(dggbak, DGGBAK, (&bal_job, "R", 436 F77_XFCN (dggbak, DGGBAK, (&bal_job, "R",
370 nn, ilo, ihi, lscale.fortran_vec(), 437 nn, ilo, ihi, lscale.fortran_vec(),
371 rscale.fortran_vec(), nn, ZZ.fortran_vec(), 438 rscale.fortran_vec(), nn, ZZ.fortran_vec(),
372 nn, info, 1L, 1L)); 439 nn, info, 1L, 1L));
373 440
374 #ifdef DEBUG 441 #ifdef DEBUG
375 if(compz == 'V') cout << "qz: balancing done; ZZ=" << endl << ZZ << endl; 442 if (compz == 'V')
376 #endif 443 cout << "qz: balancing done; ZZ=" << endl << ZZ << endl;
377 444 #endif
378 if(f77_exception_encountered) 445
379 (*current_liboctave_error_handler) ("unrecoverable error in qz(bal-R)"); 446 if (f77_exception_encountered)
380 } 447 {
448 error ("unrecoverable error in qz (bal-R)");
449 return retval;
450 }
451 }
381 452
382 static char qz_job; 453 static char qz_job;
383 qz_job = (nargout < 2 ? 'E' : 'S'); 454 qz_job = (nargout < 2 ? 'E' : 'S');
455
384 if (complex_case) 456 if (complex_case)
385 { 457 {
386 // complex case 458 // complex case
387 if (args(0).is_real_type ()) caa = aa; 459 if (args(0).is_real_type ())
388 if (args(1).is_real_type ()) cbb = bb; 460 caa = aa;
389 if(compq == 'V') CQ = QQ; 461
390 if(compz == 'V') CZ = ZZ; 462 if (args(1).is_real_type ())
391 error("complex case not done yet"); 463 cbb = bb;
392 } 464
465 if (compq == 'V')
466 CQ = QQ;
467
468 if (compz == 'V')
469 CZ = ZZ;
470
471 error ("complex case not done yet");
472 return retval;
473 }
393 else // real matrices case 474 else // real matrices case
394 { 475 {
395 #ifdef DEBUG 476 #ifdef DEBUG
396 cout << "qz: peforming qr decomposition of bb" << endl; 477 cout << "qz: peforming qr decomposition of bb" << endl;
397 #endif 478 #endif
398 479
399 // compute the QR factorization of bb 480 // compute the QR factorization of bb
400 QR bqr(bb); 481 QR bqr (bb);
401 482
402 #ifdef DEBUG 483 #ifdef DEBUG
403 cout << "qz: qr(bb) done; now peforming qz decomposition" << endl; 484 cout << "qz: qr (bb) done; now peforming qz decomposition" << endl;
404 #endif 485 #endif
405 486
406 bb = bqr.R(); 487 bb = bqr.R ();
407 #ifdef DEBUG 488
408 cout << "qz: extracted bb" << endl; 489 #ifdef DEBUG
409 #endif 490 cout << "qz: extracted bb" << endl;
410 491 #endif
411 aa = (bqr.Q()).transpose()*aa; 492
412 #ifdef DEBUG 493 aa = (bqr.Q ()).transpose ()*aa;
413 cout << "qz: updated aa " << endl; 494
414 cout << "bqr.Q () = " << endl << bqr.Q () << endl; 495 #ifdef DEBUG
415 if(compq == 'V') cout << "QQ =" << QQ << endl; 496 cout << "qz: updated aa " << endl;
416 #endif 497 cout << "bqr.Q () = " << endl << bqr.Q () << endl;
417 498
418 if(compq == 'V') QQ = QQ*bqr.Q(); 499 if (compq == 'V')
419 500 cout << "QQ =" << QQ << endl;
420 #ifdef DEBUG 501 #endif
421 cout << "qz: precursors done..." << endl; 502
422 #endif 503 if (compq == 'V')
423 504 QQ = QQ*bqr.Q ();
424 505
425 #ifdef DEBUG 506 #ifdef DEBUG
426 cout << "qz: compq = " << compq << ", compz = " << compz << endl; 507 cout << "qz: precursors done..." << endl;
427 #endif 508 #endif
428 509
429 // reduce to generalized hessenberg form 510 #ifdef DEBUG
430 F77_XFCN( dgghrd, DGGHRD, (&compq, &compz, nn, ilo, ihi, aa.fortran_vec(), 511 cout << "qz: compq = " << compq << ", compz = " << compz << endl;
431 nn, bb.fortran_vec(), nn, QQ.fortran_vec(), nn, ZZ.fortran_vec(), 512 #endif
432 nn, info,1L,1L)); 513
433 if(f77_exception_encountered) 514 // reduce to generalized hessenberg form
434 (*current_liboctave_error_handler) ("unrecoverable error in qz(dgghrd)"); 515 F77_XFCN (dgghrd, DGGHRD,
435 516 (&compq, &compz, nn, ilo, ihi, aa.fortran_vec(),
436 // check if just computing generalized eigenvalues or if we're 517 nn, bb.fortran_vec(), nn, QQ.fortran_vec(), nn,
437 // actually computing the decomposition 518 ZZ.fortran_vec(), nn, info, 1L, 1L));
438 519
439 // reduce to generalized Schur form 520 if (f77_exception_encountered)
440 F77_XFCN( dhgeqz, DHGEQZ, ( &qz_job, &compq, &compz, nn, ilo, ihi, 521 {
441 aa.fortran_vec(), nn, bb.fortran_vec(), nn, alphar.fortran_vec(), 522 error ("unrecoverable error in qz (dgghrd)");
442 alphai.fortran_vec(), betar.fortran_vec(), QQ.fortran_vec(), 523 return retval;
443 nn, ZZ.fortran_vec(), nn, work.fortran_vec(), nn, info, 1L, 1L, 1L)); 524 }
444 if(f77_exception_encountered) 525
445 (*current_liboctave_error_handler) ("unrecoverable error in qz(dhgeqz)"); 526 // check if just computing generalized eigenvalues or if we're
446 527 // actually computing the decomposition
447 } 528
529 // reduce to generalized Schur form
530 F77_XFCN (dhgeqz, DHGEQZ,
531 (&qz_job, &compq, &compz, nn, ilo, ihi,
532 aa.fortran_vec(), nn, bb.fortran_vec(), nn,
533 alphar.fortran_vec(), alphai.fortran_vec(),
534 betar.fortran_vec(), QQ.fortran_vec(), nn,
535 ZZ.fortran_vec(), nn, work.fortran_vec(), nn, info,
536 1L, 1L, 1L));
537
538 if (f77_exception_encountered)
539 {
540 error ("unrecoverable error in qz (dhgeqz)");
541 return retval;
542 }
543 }
448 544
449 // order the QZ decomposition? 545 // order the QZ decomposition?
450 if(ord_job[0] != 'N') 546 if (ord_job[0] != 'N')
451 { 547 {
452 if(complex_case) // probably not needed, but better be safe 548 if (complex_case)
453 error("qz: cannot re-order complex qz decomposition."); 549 {
454 550 // probably not needed, but better be safe
455 else 551 error ("qz: cannot re-order complex qz decomposition.");
456 { 552 return retval;
457 #ifdef DEBUG_SORT 553 }
458 cout << "qz: ordering eigenvalues: ord_job = " << ord_job[0] << endl; 554 else
459 #endif 555 {
460 556 #ifdef DEBUG_SORT
461 // declared static to avoid vfork/long jump compiler complaints 557 cout << "qz: ordering eigenvalues: ord_job = " << ord_job[0] << endl;
462 static sort_function sort_test; 558 #endif
463 sort_test = NULL; 559
464 560 // declared static to avoid vfork/long jump compiler complaints
465 switch(ord_job[0]) 561 static sort_function sort_test;
466 { 562 sort_test = NULL;
467 case 'S': 563
468 sort_test = &fin; 564 switch (ord_job[0])
469 break; 565 {
470 case 'B': 566 case 'S':
471 sort_test = &fout; 567 sort_test = &fin;
472 break; 568 break;
473 case '+': 569
474 sort_test = &fcrhp; 570 case 'B':
475 break; 571 sort_test = &fout;
476 case '-': 572 break;
477 sort_test = &folhp; 573
478 break; 574 case '+':
479 default: // this should never happen 575 sort_test = &fcrhp;
480 error("qz: illegal order option"); 576 break;
577
578 case '-':
579 sort_test = &folhp;
580 break;
581
582 default:
583 // illegal order option (should never happen, since we
584 // checked the options at the top).
585 panic_impossible ();
586 break;
481 } 587 }
482 588
483 int ndim, fail, ind[nn]; 589 int ndim, fail, ind[nn];
484 double inf_norm; 590 double inf_norm;
485 F77_XFCN (xdlange, XDLANGE, ("I", nn, nn, aa.fortran_vec (), nn, 591
486 work.fortran_vec (), inf_norm)); 592 F77_XFCN (xdlange, XDLANGE,
487 593 ("I", nn, nn, aa.fortran_vec (), nn,
488 double eps = DBL_EPSILON*inf_norm*nn; 594 work.fortran_vec (), inf_norm));
489 595
490 #ifdef DEBUG_SORT 596 double eps = DBL_EPSILON*inf_norm*nn;
491 cout << "qz: calling dsubsp: aa=" << endl; 597
492 octave_print_internal(cout,aa,0); 598 #ifdef DEBUG_SORT
493 cout << endl << "bb=" << endl; 599 cout << "qz: calling dsubsp: aa=" << endl;
494 octave_print_internal(cout,bb,0); 600 octave_print_internal (cout, aa, 0);
495 if(compz == 'V') 601 cout << endl << "bb=" << endl;
496 { 602 octave_print_internal (cout, bb, 0);
497 cout << endl << "ZZ=" << endl; 603 if (compz == 'V')
498 octave_print_internal(cout,ZZ,0); 604 {
499 } 605 cout << endl << "ZZ=" << endl;
500 cout << endl; 606 octave_print_internal (cout, ZZ, 0);
501 cout << "alphar = " << endl; 607 }
502 octave_print_internal(cout,(Matrix) alphar,0); 608 cout << endl;
503 cout << endl << "alphai = " << endl; 609 cout << "alphar = " << endl;
504 octave_print_internal(cout,(Matrix) alphai,0); 610 octave_print_internal (cout, (Matrix) alphar, 0);
505 cout << endl << "beta = " << endl; 611 cout << endl << "alphai = " << endl;
506 octave_print_internal(cout,(Matrix) betar,0); 612 octave_print_internal (cout, (Matrix) alphai, 0);
507 cout << endl; 613 cout << endl << "beta = " << endl;
508 #endif 614 octave_print_internal (cout, (Matrix) betar, 0);
509 615 cout << endl;
510 F77_XFCN( dsubsp, DSUBSP, (nn,nn,aa.fortran_vec(), bb.fortran_vec(), 616 #endif
511 ZZ.fortran_vec(), sort_test, eps, ndim, fail, ind)); 617
512 618 F77_XFCN (dsubsp, DSUBSP,
513 #ifdef DEBUG 619 (nn, nn, aa.fortran_vec(), bb.fortran_vec(),
514 cout << "qz: back from dsubsp: aa=" << endl; 620 ZZ.fortran_vec(), sort_test, eps, ndim, fail, ind));
515 octave_print_internal(cout,aa,0); 621
516 cout << endl << "bb=" << endl; 622 #ifdef DEBUG
517 octave_print_internal(cout,bb,0); 623 cout << "qz: back from dsubsp: aa=" << endl;
518 if(compz == 'V') 624 octave_print_internal (cout, aa, 0);
519 { 625 cout << endl << "bb=" << endl;
520 cout << endl << "ZZ=" << endl; 626 octave_print_internal (cout, bb, 0);
521 octave_print_internal(cout,ZZ,0); 627 if (compz == 'V')
522 } 628 {
523 cout << endl; 629 cout << endl << "ZZ=" << endl;
524 #endif 630 octave_print_internal (cout, ZZ, 0);
525 631 }
526 // manually update alphar, alphai, betar 632 cout << endl;
527 static int jj; 633 #endif
528 jj=0; 634
529 while(jj < nn) 635 // manually update alphar, alphai, betar
530 { 636 static int jj;
531 #ifdef DEBUG_EIG 637
532 cout << "computing gen eig #" << jj << endl; 638 jj=0;
533 #endif 639 while (jj < nn)
534 640 {
535 static int zcnt; // number of zeros in this block 641 #ifdef DEBUG_EIG
536 if(jj == (nn-1)) 642 cout << "computing gen eig #" << jj << endl;
537 zcnt = 1; 643 #endif
538 else if(aa(jj+1,jj) == 0) 644
539 zcnt = 1; 645 static int zcnt; // number of zeros in this block
540 else zcnt = 2; 646
541 647 if (jj == (nn-1))
542 if(zcnt == 1) // real zero 648 zcnt = 1;
543 { 649 else if (aa(jj+1,jj) == 0)
544 #ifdef DEBUG_EIG 650 zcnt = 1;
545 cout << " single gen eig:" << endl; 651 else zcnt = 2;
546 cout << " alphar(" << jj << ") = " << aa(jj,jj) << endl; 652
547 cout << " betar( " << jj << ") = " << bb(jj,jj) << endl; 653 if (zcnt == 1) // real zero
548 cout << " alphai(" << jj << ") = 0" << endl; 654 {
549 #endif 655 #ifdef DEBUG_EIG
550 656 cout << " single gen eig:" << endl;
551 alphar(jj) = aa(jj,jj); 657 cout << " alphar(" << jj << ") = " << aa(jj,jj) << endl;
552 alphai(jj) = 0; 658 cout << " betar( " << jj << ") = " << bb(jj,jj) << endl;
553 betar(jj) = bb(jj,jj); 659 cout << " alphai(" << jj << ") = 0" << endl;
554 } 660 #endif
555 else // complex conjugate pair 661
556 { 662 alphar(jj) = aa(jj,jj);
557 #ifdef DEBUG_EIG 663 alphai(jj) = 0;
558 cout << "qz: calling dlag2:" << endl; 664 betar(jj) = bb(jj,jj);
559 cout << "safmin=" << setiosflags(ios::scientific) << safmin << endl; 665 }
560 for(int idr = jj ; idr <= jj+1 ; idr++) 666 else
561 { 667 {
562 for(int idc = jj ; idc <= jj+1 ; idc++) 668 // complex conjugate pair
563 { 669 #ifdef DEBUG_EIG
564 cout << "aa(" << idr << "," << idc << ")=" 670 cout << "qz: calling dlag2:" << endl;
565 << aa(idr,idc) << endl; 671 cout << "safmin="
566 cout << "bb(" << idr << "," << idc << ")=" 672 << setiosflags (ios::scientific) << safmin << endl;
567 << bb(idr,idc) << endl; 673
568 } 674 for (int idr = jj; idr <= jj+1; idr++)
569 } 675 {
570 #endif 676 for (int idc = jj; idc <= jj+1; idc++)
571 double scale1, scale2, wr1, wr2, wi; 677 {
572 F77_XFCN( dlag2, DLAG2, ( &aa(jj,jj), nn, &bb(jj,jj), nn, safmin, 678 cout << "aa(" << idr << "," << idc << ")="
573 scale1, scale2, wr1, wr2, wi)); 679 << aa(idr,idc) << endl;
574 680 cout << "bb(" << idr << "," << idc << ")="
575 #ifdef DEBUG_EIG 681 << bb(idr,idc) << endl;
576 cout << "dlag2 returns: scale1=" << scale1 682 }
577 << "\tscale2=" << scale2 << endl 683 }
578 << "\twr1=" << wr1 << "\twr2=" << wr2 684 #endif
579 << "\twi=" << wi << endl; 685
580 #endif 686 double scale1, scale2, wr1, wr2, wi;
581 // just to be safe, check if it's a real pair 687 F77_XFCN (dlag2, DLAG2,
582 if(wi == 0) 688 (&aa(jj,jj), nn, &bb(jj,jj), nn, safmin,
583 { 689 scale1, scale2, wr1, wr2, wi));
584 alphar(jj) = wr1; 690
585 alphai(jj) = 0; 691 #ifdef DEBUG_EIG
586 betar(jj) = scale1; 692 cout << "dlag2 returns: scale1=" << scale1
587 alphar(jj+1) = wr2; 693 << "\tscale2=" << scale2 << endl
588 alphai(jj+1) = 0; 694 << "\twr1=" << wr1 << "\twr2=" << wr2
589 betar(jj+1) = scale2; 695 << "\twi=" << wi << endl;
590 } 696 #endif
591 else 697
592 { 698 // just to be safe, check if it's a real pair
593 alphar(jj) = alphar(jj+1)=wr1; 699 if (wi == 0)
594 alphai(jj) = -(alphai(jj+1) = wi); 700 {
595 betar(jj) = betar(jj+1) = scale1; 701 alphar(jj) = wr1;
596 } 702 alphai(jj) = 0;
597 } 703 betar(jj) = scale1;
598 704 alphar(jj+1) = wr2;
599 jj += zcnt; // advance past this block 705 alphai(jj+1) = 0;
600 } 706 betar(jj+1) = scale2;
601 707 }
602 #ifdef DEBUG_SORT 708 else
603 cout << "qz: back from dsubsp: aa=" << endl; 709 {
604 octave_print_internal(cout,aa,0); 710 alphar(jj) = alphar(jj+1)=wr1;
605 cout << endl << "bb=" << endl; 711 alphai(jj) = -(alphai(jj+1) = wi);
606 octave_print_internal(cout,bb,0); 712 betar(jj) = betar(jj+1) = scale1;
607 if(compz == 'V') 713 }
608 { 714 }
609 cout << endl << "ZZ=" << endl; 715
610 octave_print_internal(cout,ZZ,0); 716 // advance past this block
611 } 717 jj += zcnt;
612 cout << endl << "qz: ndim=" << ndim << endl << "fail=" << fail << endl; 718 }
613 cout << "alphar = " << endl; 719
614 octave_print_internal(cout,(Matrix) alphar,0); 720 #ifdef DEBUG_SORT
615 cout << endl << "alphai = " << endl; 721 cout << "qz: back from dsubsp: aa=" << endl;
616 octave_print_internal(cout,(Matrix) alphai,0); 722 octave_print_internal (cout, aa, 0);
617 cout << endl << "beta = " << endl; 723 cout << endl << "bb=" << endl;
618 octave_print_internal(cout,(Matrix) betar,0); 724 octave_print_internal (cout, bb, 0);
619 cout << endl; 725
620 #endif 726 if (compz == 'V')
621 } 727 {
622 } 728 cout << endl << "ZZ=" << endl;
623 729 octave_print_internal (cout, ZZ, 0);
730 }
731 cout << endl << "qz: ndim=" << ndim << endl
732 << "fail=" << fail << endl;
733 cout << "alphar = " << endl;
734 octave_print_internal (cout, (Matrix) alphar, 0);
735 cout << endl << "alphai = " << endl;
736 octave_print_internal (cout, (Matrix) alphai, 0);
737 cout << endl << "beta = " << endl;
738 octave_print_internal (cout, (Matrix) betar, 0);
739 cout << endl;
740 #endif
741 }
742 }
743
624 // compute generalized eigenvalues? 744 // compute generalized eigenvalues?
625 ComplexColumnVector gev; 745 ComplexColumnVector gev;
626 if(nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) 746
627 { 747 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4))
628 if(complex_case) 748 {
629 error("complex case not yet implemented"); 749 if (complex_case)
630 else 750 {
631 { 751 error ("complex case not yet implemented");
632 #ifdef DEBUG 752 return retval;
633 cout << "qz: computing generalized eigenvalues" << endl; 753 }
634 #endif 754 else
635 755 {
636 // return finite generalized eigenvalues 756 #ifdef DEBUG
637 int ii, cnt = 0; 757 cout << "qz: computing generalized eigenvalues" << endl;
638 for( ii=0 ; ii < nn ; ii++) 758 #endif
639 if(betar(ii) != 0) 759
640 cnt++; 760 // return finite generalized eigenvalues
641 ComplexColumnVector tmp(cnt); 761 int cnt = 0;
642 for( ii=0 ; ii < nn ; ii++) 762
643 if(betar(ii) != 0) 763 for (int ii = 0; ii < nn; ii++)
644 tmp(ii) = Complex(alphar(ii), alphai(ii))/betar(ii); 764 if (betar(ii) != 0)
645 gev = tmp; 765 cnt++;
646 } 766
767 ComplexColumnVector tmp(cnt);
768
769 for (int ii = 0; ii < nn; ii++)
770 if (betar(ii) != 0)
771 tmp(ii) = Complex(alphar(ii), alphai(ii))/betar(ii);
772 gev = tmp;
773 }
774 }
775
776 // right, left eigenvector matrices
777 if (nargout >= 5)
778 {
779 char side = (nargout == 5 ? 'R' : 'B'); // which side to compute?
780 char howmny = 'B'; // compute all of them and backtransform
781 int *select = NULL; // dummy pointer; select is not used.
782 int m;
783
784 if (complex_case)
785 {
786 error ("complex type not yet implemented");
787 return retval;
788 }
789 else
790 {
791 #ifdef DEBUG
792 cout << "qz: computing generalized eigenvectors" << endl;
793 #endif
794
795 VL = QQ;
796 VR = ZZ;
797
798 F77_XFCN (dtgevc, DTGEVC,
799 (&side, &howmny, select, nn, aa.fortran_vec(),
800 nn, bb.fortran_vec(), nn, VL.fortran_vec(), nn,
801 VR.fortran_vec(), nn, nn, m, work.fortran_vec(),
802 info, 1L, 1L));
803
804 if (f77_exception_encountered)
805 {
806 error ("unrecoverable error in qz (dtgevc)");
807 return retval;
808 }
809
810 // now construct the complex form of VV, WW
811 int jj = 0;
812
813 while (jj < nn)
814 {
815 // see if real or complex eigenvalue
816 int cinc = 2; // column increment; assume complex eigenvalue
817
818 if (jj == (nn-1))
819 cinc = 1; // single column
820 else if (aa(jj+1,jj) == 0)
821 cinc = 1;
822
823 // now copy the eigenvector (s) to CVR, CVL
824 if (cinc == 1)
825 {
826 for (int ii = 0; ii < nn; ii++)
827 CVR(ii,jj) = VR(ii,jj);
828
829 if (side == 'B')
830 for (int ii = 0; ii < nn; ii++)
831 CVL(ii,jj) = VL(ii,jj);
832 }
833 else
834 {
835 // double column; complex vector
836
837 for (int ii = 0; ii < nn; ii++)
838 {
839 CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1));
840 CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1));
841 }
842
843 if (side == 'B')
844 for (int ii = 0; ii < nn; ii++)
845 {
846 CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1));
847 CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1));
848 }
849 }
850
851 // advance to next eigenvectors (if any)
852 jj += cinc;
853 }
854 }
855 }
856
857 switch (nargout)
858 {
859 case 7:
860 retval(6) = gev;
861
862 case 6: // return eigenvectors
863 retval(5) = CVL;
864
865 case 5: // return eigenvectors
866 retval(4) = CVR;
867
868 case 4:
869 if (nargin == 3)
870 {
871 #ifdef DEBUG
872 cout << "qz: sort: retval(3) = gev = " << endl;
873 octave_print_internal (cout, gev);
874 cout << endl;
875 #endif
876 retval(3) = gev;
877 }
878 else
879 retval(3) = ZZ;
880
881 case 3:
882 if (nargin == 3)
883 retval(2) = ZZ;
884 else
885 retval(2) = QQ;
886
887 case 2:
888 #ifdef DEBUG
889 cout << "qz: retval (1) = bb = " << endl;
890 octave_print_internal (cout, bb, 0);
891 cout << endl << "qz: retval(0) = aa = " <<endl;
892 octave_print_internal (cout, aa, 0);
893 cout << endl;
894 #endif
895 retval(1) = bb;
896 retval(0) = aa;
897 break;
898
899 case 1:
900 case 0:
901 #ifdef DEBUG
902 cout << "qz: retval(0) = gev = " << gev << endl;
903 #endif
904 retval(0) = gev;
905 break;
906
907 default:
908 error ("qz: too many return arguments.");
909 break;
647 } 910 }
648 911
649 // right, left eigenvector matrices 912 #ifdef DEBUG
650 if(nargout >= 5)
651 {
652 char side = (nargout == 5 ? 'R' : 'B'), // which side to compute?
653 howmny = 'B'; // compute all of them and backtransform
654 int *select = NULL; // dummy pointer; select is not used.
655 int m;
656
657 if(complex_case)
658 error("complex type not yet implemented");
659 else
660 {
661 #ifdef DEBUG
662 cout << "qz: computing generalized eigenvectors" << endl;
663 #endif
664
665 VL = QQ;
666 VR = ZZ;
667
668 F77_XFCN( dtgevc, DTGEVC, ( &side, &howmny, select, nn, aa.fortran_vec(),
669 nn, bb.fortran_vec(), nn, VL.fortran_vec(), nn, VR.fortran_vec(),
670 nn, nn, m, work.fortran_vec(), info, 1L, 1L ));
671 if(f77_exception_encountered)
672 (*current_liboctave_error_handler)
673 ("unrecoverable error in qz(dtgevc)");
674
675 // now construct the complex form of VV, WW
676 int jj = 0;
677 while(jj < nn)
678 {
679 // see if real or complex eigenvalue
680 int cinc = 2; // column increment; assume complex eigenvalue
681 if(jj == (nn-1))
682 cinc = 1; // single column
683 else if(aa(jj+1,jj) == 0)
684 cinc = 1;
685
686 // now copy the eigenvector(s) to CVR, CVL
687 if(cinc == 1)
688 {
689 int ii;
690 for(ii = 0; ii < nn ; ii++)
691 CVR(ii,jj) = VR(ii,jj);
692 if(side == 'B')
693 for(ii = 0; ii < nn ; ii++)
694 CVL(ii,jj) = VL(ii,jj);
695 }
696 else // double column; complex vector
697 {
698 int ii;
699 for(ii = 0; ii < nn ; ii++)
700 {
701 CVR(ii,jj) = Complex(VR(ii,jj),VR(ii,jj+1));
702 CVR(ii,jj+1) = Complex(VR(ii,jj),-VR(ii,jj+1));
703 }
704 if(side == 'B')
705 for(ii = 0; ii < nn ; ii++)
706 {
707 CVL(ii,jj) = Complex(VL(ii,jj),VL(ii,jj+1));
708 CVL(ii,jj+1) = Complex(VL(ii,jj),-VL(ii,jj+1));
709 }
710 }
711 jj += cinc; // advance to next eigenvectors (if any)
712 }
713 }
714 }
715
716 switch(nargout)
717 {
718 case 7:
719 retval(6) = gev;
720 case 6: // return eigenvectors
721 retval(5) = CVL;
722 case 5: // return eigenvectors
723 retval(4) = CVR;
724 case 4:
725 if(nargin == 3)
726 {
727 #ifdef DEBUG
728 cout << "qz: sort: retval(3) = gev = " << endl;
729 octave_print_internal(cout,gev);
730 cout << endl;
731 #endif
732 retval(3) = gev;
733 }
734 else retval(3) = ZZ;
735 case 3:
736 if(nargin == 3)
737 retval(2) = ZZ;
738 else
739 retval(2) = QQ;
740 case 2:
741 #ifdef DEBUG
742 cout << "qz: retval(1) = bb = " << endl;
743 octave_print_internal(cout,bb,0);
744 cout << endl << "qz: retval(0) = aa = " <<endl;
745 octave_print_internal(cout,aa,0);
746 cout << endl;
747 #endif
748 retval(1) = bb;
749 retval(0) = aa;
750 break;
751 case 1:
752 case 0:
753 #ifdef DEBUG
754 cout << "qz: retval(0) = gev = " << gev << endl;
755 #endif
756 retval(0) = gev;
757 break;
758 default:
759 error("qz: too many return arguments. Sorry. ");
760 }
761
762 #ifdef DEBUG
763 cout << "qz: exiting (at long last)" << endl; 913 cout << "qz: exiting (at long last)" << endl;
764 #endif 914 #endif
765 915
766 return retval; 916 return retval;
767 } 917 }
768 918
769 /* 919 /*
770 ;;; Local Variables: *** 920 ;;; Local Variables: ***