1 /* ./src_f77/ctgevc.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static complex c_b1 = {0.f,0.f};
11 static complex c_b2 = {1.f,0.f};
12 static integer c__1 = 1;
13
ctgevc_(char * side,char * howmny,logical * select,integer * n,complex * a,integer * lda,complex * b,integer * ldb,complex * vl,integer * ldvl,complex * vr,integer * ldvr,integer * mm,integer * m,complex * work,real * rwork,integer * info,ftnlen side_len,ftnlen howmny_len)14 /* Subroutine */ int ctgevc_(char *side, char *howmny, logical *select,
15 integer *n, complex *a, integer *lda, complex *b, integer *ldb,
16 complex *vl, integer *ldvl, complex *vr, integer *ldvr, integer *mm,
17 integer *m, complex *work, real *rwork, integer *info, ftnlen
18 side_len, ftnlen howmny_len)
19 {
20 /* System generated locals */
21 integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1,
22 vr_offset, i__1, i__2, i__3, i__4, i__5;
23 real r__1, r__2, r__3, r__4, r__5, r__6;
24 complex q__1, q__2, q__3, q__4;
25
26 /* Builtin functions */
27 double r_imag(complex *);
28 void r_cnjg(complex *, complex *);
29
30 /* Local variables */
31 static complex d__;
32 static integer i__, j;
33 static complex ca, cb;
34 static integer je, im, jr;
35 static real big;
36 static logical lsa, lsb;
37 static real ulp;
38 static complex sum;
39 static integer ibeg, ieig, iend;
40 static real dmin__;
41 static integer isrc;
42 static real temp;
43 static complex suma, sumb;
44 static real xmax, scale;
45 static logical ilall;
46 static integer iside;
47 static real sbeta;
48 extern logical lsame_(char *, char *, ftnlen, ftnlen);
49 extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
50 , complex *, integer *, complex *, integer *, complex *, complex *
51 , integer *, ftnlen);
52 static real small;
53 static logical compl;
54 static real anorm, bnorm;
55 static logical compr, ilbbad;
56 static real acoefa, bcoefa, acoeff;
57 static complex bcoeff;
58 static logical ilback;
59 extern /* Subroutine */ int slabad_(real *, real *);
60 static real ascale, bscale;
61 extern /* Complex */ VOID cladiv_(complex *, complex *, complex *);
62 extern doublereal slamch_(char *, ftnlen);
63 static complex salpha;
64 static real safmin;
65 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
66 static real bignum;
67 static logical ilcomp;
68 static integer ihwmny;
69
70
71 /* -- LAPACK routine (version 3.0) -- */
72 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
73 /* Courant Institute, Argonne National Lab, and Rice University */
74 /* June 30, 1999 */
75
76 /* .. Scalar Arguments .. */
77 /* .. */
78 /* .. Array Arguments .. */
79 /* .. */
80
81
82 /* Purpose */
83 /* ======= */
84
85 /* CTGEVC computes some or all of the right and/or left generalized */
86 /* eigenvectors of a pair of complex upper triangular matrices (A,B). */
87
88 /* The right generalized eigenvector x and the left generalized */
89 /* eigenvector y of (A,B) corresponding to a generalized eigenvalue */
90 /* w are defined by: */
91
92 /* (A - wB) * x = 0 and y**H * (A - wB) = 0 */
93
94 /* where y**H denotes the conjugate tranpose of y. */
95
96 /* If an eigenvalue w is determined by zero diagonal elements of both A */
97 /* and B, a unit vector is returned as the corresponding eigenvector. */
98
99 /* If all eigenvectors are requested, the routine may either return */
100 /* the matrices X and/or Y of right or left eigenvectors of (A,B), or */
101 /* the products Z*X and/or Q*Y, where Z and Q are input unitary */
102 /* matrices. If (A,B) was obtained from the generalized Schur */
103 /* factorization of an original pair of matrices */
104 /* (A0,B0) = (Q*A*Z**H,Q*B*Z**H), */
105 /* then Z*X and Q*Y are the matrices of right or left eigenvectors of */
106 /* A. */
107
108 /* Arguments */
109 /* ========= */
110
111 /* SIDE (input) CHARACTER*1 */
112 /* = 'R': compute right eigenvectors only; */
113 /* = 'L': compute left eigenvectors only; */
114 /* = 'B': compute both right and left eigenvectors. */
115
116 /* HOWMNY (input) CHARACTER*1 */
117 /* = 'A': compute all right and/or left eigenvectors; */
118 /* = 'B': compute all right and/or left eigenvectors, and */
119 /* backtransform them using the input matrices supplied */
120 /* in VR and/or VL; */
121 /* = 'S': compute selected right and/or left eigenvectors, */
122 /* specified by the logical array SELECT. */
123
124 /* SELECT (input) LOGICAL array, dimension (N) */
125 /* If HOWMNY='S', SELECT specifies the eigenvectors to be */
126 /* computed. */
127 /* If HOWMNY='A' or 'B', SELECT is not referenced. */
128 /* To select the eigenvector corresponding to the j-th */
129 /* eigenvalue, SELECT(j) must be set to .TRUE.. */
130
131 /* N (input) INTEGER */
132 /* The order of the matrices A and B. N >= 0. */
133
134 /* A (input) COMPLEX array, dimension (LDA,N) */
135 /* The upper triangular matrix A. */
136
137 /* LDA (input) INTEGER */
138 /* The leading dimension of array A. LDA >= max(1,N). */
139
140 /* B (input) COMPLEX array, dimension (LDB,N) */
141 /* The upper triangular matrix B. B must have real diagonal */
142 /* elements. */
143
144 /* LDB (input) INTEGER */
145 /* The leading dimension of array B. LDB >= max(1,N). */
146
147 /* VL (input/output) COMPLEX array, dimension (LDVL,MM) */
148 /* On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must */
149 /* contain an N-by-N matrix Q (usually the unitary matrix Q */
150 /* of left Schur vectors returned by CHGEQZ). */
151 /* On exit, if SIDE = 'L' or 'B', VL contains: */
152 /* if HOWMNY = 'A', the matrix Y of left eigenvectors of (A,B); */
153 /* if HOWMNY = 'B', the matrix Q*Y; */
154 /* if HOWMNY = 'S', the left eigenvectors of (A,B) specified by */
155 /* SELECT, stored consecutively in the columns of */
156 /* VL, in the same order as their eigenvalues. */
157 /* If SIDE = 'R', VL is not referenced. */
158
159 /* LDVL (input) INTEGER */
160 /* The leading dimension of array VL. */
161 /* LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. */
162
163 /* VR (input/output) COMPLEX array, dimension (LDVR,MM) */
164 /* On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must */
165 /* contain an N-by-N matrix Q (usually the unitary matrix Z */
166 /* of right Schur vectors returned by CHGEQZ). */
167 /* On exit, if SIDE = 'R' or 'B', VR contains: */
168 /* if HOWMNY = 'A', the matrix X of right eigenvectors of (A,B); */
169 /* if HOWMNY = 'B', the matrix Z*X; */
170 /* if HOWMNY = 'S', the right eigenvectors of (A,B) specified by */
171 /* SELECT, stored consecutively in the columns of */
172 /* VR, in the same order as their eigenvalues. */
173 /* If SIDE = 'L', VR is not referenced. */
174
175 /* LDVR (input) INTEGER */
176 /* The leading dimension of the array VR. */
177 /* LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. */
178
179 /* MM (input) INTEGER */
180 /* The number of columns in the arrays VL and/or VR. MM >= M. */
181
182 /* M (output) INTEGER */
183 /* The number of columns in the arrays VL and/or VR actually */
184 /* used to store the eigenvectors. If HOWMNY = 'A' or 'B', M */
185 /* is set to N. Each selected eigenvector occupies one column. */
186
187 /* WORK (workspace) COMPLEX array, dimension (2*N) */
188
189 /* RWORK (workspace) REAL array, dimension (2*N) */
190
191 /* INFO (output) INTEGER */
192 /* = 0: successful exit. */
193 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
194
195 /* ===================================================================== */
196
197 /* .. Parameters .. */
198 /* .. */
199 /* .. Local Scalars .. */
200 /* .. */
201 /* .. External Functions .. */
202 /* .. */
203 /* .. External Subroutines .. */
204 /* .. */
205 /* .. Intrinsic Functions .. */
206 /* .. */
207 /* .. Statement Functions .. */
208 /* .. */
209 /* .. Statement Function definitions .. */
210 /* .. */
211 /* .. Executable Statements .. */
212
213 /* Decode and Test the input parameters */
214
215 /* Parameter adjustments */
216 --select;
217 a_dim1 = *lda;
218 a_offset = 1 + a_dim1;
219 a -= a_offset;
220 b_dim1 = *ldb;
221 b_offset = 1 + b_dim1;
222 b -= b_offset;
223 vl_dim1 = *ldvl;
224 vl_offset = 1 + vl_dim1;
225 vl -= vl_offset;
226 vr_dim1 = *ldvr;
227 vr_offset = 1 + vr_dim1;
228 vr -= vr_offset;
229 --work;
230 --rwork;
231
232 /* Function Body */
233 if (lsame_(howmny, "A", (ftnlen)1, (ftnlen)1)) {
234 ihwmny = 1;
235 ilall = TRUE_;
236 ilback = FALSE_;
237 } else if (lsame_(howmny, "S", (ftnlen)1, (ftnlen)1)) {
238 ihwmny = 2;
239 ilall = FALSE_;
240 ilback = FALSE_;
241 } else if (lsame_(howmny, "B", (ftnlen)1, (ftnlen)1) || lsame_(howmny,
242 "T", (ftnlen)1, (ftnlen)1)) {
243 ihwmny = 3;
244 ilall = TRUE_;
245 ilback = TRUE_;
246 } else {
247 ihwmny = -1;
248 }
249
250 if (lsame_(side, "R", (ftnlen)1, (ftnlen)1)) {
251 iside = 1;
252 compl = FALSE_;
253 compr = TRUE_;
254 } else if (lsame_(side, "L", (ftnlen)1, (ftnlen)1)) {
255 iside = 2;
256 compl = TRUE_;
257 compr = FALSE_;
258 } else if (lsame_(side, "B", (ftnlen)1, (ftnlen)1)) {
259 iside = 3;
260 compl = TRUE_;
261 compr = TRUE_;
262 } else {
263 iside = -1;
264 }
265
266 *info = 0;
267 if (iside < 0) {
268 *info = -1;
269 } else if (ihwmny < 0) {
270 *info = -2;
271 } else if (*n < 0) {
272 *info = -4;
273 } else if (*lda < max(1,*n)) {
274 *info = -6;
275 } else if (*ldb < max(1,*n)) {
276 *info = -8;
277 }
278 if (*info != 0) {
279 i__1 = -(*info);
280 xerbla_("CTGEVC", &i__1, (ftnlen)6);
281 return 0;
282 }
283
284 /* Count the number of eigenvectors */
285
286 if (! ilall) {
287 im = 0;
288 i__1 = *n;
289 for (j = 1; j <= i__1; ++j) {
290 if (select[j]) {
291 ++im;
292 }
293 /* L10: */
294 }
295 } else {
296 im = *n;
297 }
298
299 /* Check diagonal of B */
300
301 ilbbad = FALSE_;
302 i__1 = *n;
303 for (j = 1; j <= i__1; ++j) {
304 if (r_imag(&b[j + j * b_dim1]) != 0.f) {
305 ilbbad = TRUE_;
306 }
307 /* L20: */
308 }
309
310 if (ilbbad) {
311 *info = -7;
312 } else if (compl && *ldvl < *n || *ldvl < 1) {
313 *info = -10;
314 } else if (compr && *ldvr < *n || *ldvr < 1) {
315 *info = -12;
316 } else if (*mm < im) {
317 *info = -13;
318 }
319 if (*info != 0) {
320 i__1 = -(*info);
321 xerbla_("CTGEVC", &i__1, (ftnlen)6);
322 return 0;
323 }
324
325 /* Quick return if possible */
326
327 *m = im;
328 if (*n == 0) {
329 return 0;
330 }
331
332 /* Machine Constants */
333
334 safmin = slamch_("Safe minimum", (ftnlen)12);
335 big = 1.f / safmin;
336 slabad_(&safmin, &big);
337 ulp = slamch_("Epsilon", (ftnlen)7) * slamch_("Base", (ftnlen)4);
338 small = safmin * *n / ulp;
339 big = 1.f / small;
340 bignum = 1.f / (safmin * *n);
341
342 /* Compute the 1-norm of each column of the strictly upper triangular */
343 /* part of A and B to check for possible overflow in the triangular */
344 /* solver. */
345
346 i__1 = a_dim1 + 1;
347 anorm = (r__1 = a[i__1].r, dabs(r__1)) + (r__2 = r_imag(&a[a_dim1 + 1]),
348 dabs(r__2));
349 i__1 = b_dim1 + 1;
350 bnorm = (r__1 = b[i__1].r, dabs(r__1)) + (r__2 = r_imag(&b[b_dim1 + 1]),
351 dabs(r__2));
352 rwork[1] = 0.f;
353 rwork[*n + 1] = 0.f;
354 i__1 = *n;
355 for (j = 2; j <= i__1; ++j) {
356 rwork[j] = 0.f;
357 rwork[*n + j] = 0.f;
358 i__2 = j - 1;
359 for (i__ = 1; i__ <= i__2; ++i__) {
360 i__3 = i__ + j * a_dim1;
361 rwork[j] += (r__1 = a[i__3].r, dabs(r__1)) + (r__2 = r_imag(&a[
362 i__ + j * a_dim1]), dabs(r__2));
363 i__3 = i__ + j * b_dim1;
364 rwork[*n + j] += (r__1 = b[i__3].r, dabs(r__1)) + (r__2 = r_imag(&
365 b[i__ + j * b_dim1]), dabs(r__2));
366 /* L30: */
367 }
368 /* Computing MAX */
369 i__2 = j + j * a_dim1;
370 r__3 = anorm, r__4 = rwork[j] + ((r__1 = a[i__2].r, dabs(r__1)) + (
371 r__2 = r_imag(&a[j + j * a_dim1]), dabs(r__2)));
372 anorm = dmax(r__3,r__4);
373 /* Computing MAX */
374 i__2 = j + j * b_dim1;
375 r__3 = bnorm, r__4 = rwork[*n + j] + ((r__1 = b[i__2].r, dabs(r__1))
376 + (r__2 = r_imag(&b[j + j * b_dim1]), dabs(r__2)));
377 bnorm = dmax(r__3,r__4);
378 /* L40: */
379 }
380
381 ascale = 1.f / dmax(anorm,safmin);
382 bscale = 1.f / dmax(bnorm,safmin);
383
384 /* Left eigenvectors */
385
386 if (compl) {
387 ieig = 0;
388
389 /* Main loop over eigenvalues */
390
391 i__1 = *n;
392 for (je = 1; je <= i__1; ++je) {
393 if (ilall) {
394 ilcomp = TRUE_;
395 } else {
396 ilcomp = select[je];
397 }
398 if (ilcomp) {
399 ++ieig;
400
401 i__2 = je + je * a_dim1;
402 i__3 = je + je * b_dim1;
403 if ((r__2 = a[i__2].r, dabs(r__2)) + (r__3 = r_imag(&a[je +
404 je * a_dim1]), dabs(r__3)) <= safmin && (r__1 = b[
405 i__3].r, dabs(r__1)) <= safmin) {
406
407 /* Singular matrix pencil -- return unit eigenvector */
408
409 i__2 = *n;
410 for (jr = 1; jr <= i__2; ++jr) {
411 i__3 = jr + ieig * vl_dim1;
412 vl[i__3].r = 0.f, vl[i__3].i = 0.f;
413 /* L50: */
414 }
415 i__2 = ieig + ieig * vl_dim1;
416 vl[i__2].r = 1.f, vl[i__2].i = 0.f;
417 goto L140;
418 }
419
420 /* Non-singular eigenvalue: */
421 /* Compute coefficients a and b in */
422 /* H */
423 /* y ( a A - b B ) = 0 */
424
425 /* Computing MAX */
426 i__2 = je + je * a_dim1;
427 i__3 = je + je * b_dim1;
428 r__4 = ((r__2 = a[i__2].r, dabs(r__2)) + (r__3 = r_imag(&a[je
429 + je * a_dim1]), dabs(r__3))) * ascale, r__5 = (r__1 =
430 b[i__3].r, dabs(r__1)) * bscale, r__4 = max(r__4,
431 r__5);
432 temp = 1.f / dmax(r__4,safmin);
433 i__2 = je + je * a_dim1;
434 q__2.r = temp * a[i__2].r, q__2.i = temp * a[i__2].i;
435 q__1.r = ascale * q__2.r, q__1.i = ascale * q__2.i;
436 salpha.r = q__1.r, salpha.i = q__1.i;
437 i__2 = je + je * b_dim1;
438 sbeta = temp * b[i__2].r * bscale;
439 acoeff = sbeta * ascale;
440 q__1.r = bscale * salpha.r, q__1.i = bscale * salpha.i;
441 bcoeff.r = q__1.r, bcoeff.i = q__1.i;
442
443 /* Scale to avoid underflow */
444
445 lsa = dabs(sbeta) >= safmin && dabs(acoeff) < small;
446 lsb = (r__1 = salpha.r, dabs(r__1)) + (r__2 = r_imag(&salpha),
447 dabs(r__2)) >= safmin && (r__3 = bcoeff.r, dabs(r__3)
448 ) + (r__4 = r_imag(&bcoeff), dabs(r__4)) < small;
449
450 scale = 1.f;
451 if (lsa) {
452 scale = small / dabs(sbeta) * dmin(anorm,big);
453 }
454 if (lsb) {
455 /* Computing MAX */
456 r__3 = scale, r__4 = small / ((r__1 = salpha.r, dabs(r__1)
457 ) + (r__2 = r_imag(&salpha), dabs(r__2))) * dmin(
458 bnorm,big);
459 scale = dmax(r__3,r__4);
460 }
461 if (lsa || lsb) {
462 /* Computing MIN */
463 /* Computing MAX */
464 r__5 = 1.f, r__6 = dabs(acoeff), r__5 = max(r__5,r__6),
465 r__6 = (r__1 = bcoeff.r, dabs(r__1)) + (r__2 =
466 r_imag(&bcoeff), dabs(r__2));
467 r__3 = scale, r__4 = 1.f / (safmin * dmax(r__5,r__6));
468 scale = dmin(r__3,r__4);
469 if (lsa) {
470 acoeff = ascale * (scale * sbeta);
471 } else {
472 acoeff = scale * acoeff;
473 }
474 if (lsb) {
475 q__2.r = scale * salpha.r, q__2.i = scale * salpha.i;
476 q__1.r = bscale * q__2.r, q__1.i = bscale * q__2.i;
477 bcoeff.r = q__1.r, bcoeff.i = q__1.i;
478 } else {
479 q__1.r = scale * bcoeff.r, q__1.i = scale * bcoeff.i;
480 bcoeff.r = q__1.r, bcoeff.i = q__1.i;
481 }
482 }
483
484 acoefa = dabs(acoeff);
485 bcoefa = (r__1 = bcoeff.r, dabs(r__1)) + (r__2 = r_imag(&
486 bcoeff), dabs(r__2));
487 xmax = 1.f;
488 i__2 = *n;
489 for (jr = 1; jr <= i__2; ++jr) {
490 i__3 = jr;
491 work[i__3].r = 0.f, work[i__3].i = 0.f;
492 /* L60: */
493 }
494 i__2 = je;
495 work[i__2].r = 1.f, work[i__2].i = 0.f;
496 /* Computing MAX */
497 r__1 = ulp * acoefa * anorm, r__2 = ulp * bcoefa * bnorm,
498 r__1 = max(r__1,r__2);
499 dmin__ = dmax(r__1,safmin);
500
501 /* H */
502 /* Triangular solve of (a A - b B) y = 0 */
503
504 /* H */
505 /* (rowwise in (a A - b B) , or columnwise in a A - b B) */
506
507 i__2 = *n;
508 for (j = je + 1; j <= i__2; ++j) {
509
510 /* Compute */
511 /* j-1 */
512 /* SUM = sum conjg( a*A(k,j) - b*B(k,j) )*x(k) */
513 /* k=je */
514 /* (Scale if necessary) */
515
516 temp = 1.f / xmax;
517 if (acoefa * rwork[j] + bcoefa * rwork[*n + j] > bignum *
518 temp) {
519 i__3 = j - 1;
520 for (jr = je; jr <= i__3; ++jr) {
521 i__4 = jr;
522 i__5 = jr;
523 q__1.r = temp * work[i__5].r, q__1.i = temp *
524 work[i__5].i;
525 work[i__4].r = q__1.r, work[i__4].i = q__1.i;
526 /* L70: */
527 }
528 xmax = 1.f;
529 }
530 suma.r = 0.f, suma.i = 0.f;
531 sumb.r = 0.f, sumb.i = 0.f;
532
533 i__3 = j - 1;
534 for (jr = je; jr <= i__3; ++jr) {
535 r_cnjg(&q__3, &a[jr + j * a_dim1]);
536 i__4 = jr;
537 q__2.r = q__3.r * work[i__4].r - q__3.i * work[i__4]
538 .i, q__2.i = q__3.r * work[i__4].i + q__3.i *
539 work[i__4].r;
540 q__1.r = suma.r + q__2.r, q__1.i = suma.i + q__2.i;
541 suma.r = q__1.r, suma.i = q__1.i;
542 r_cnjg(&q__3, &b[jr + j * b_dim1]);
543 i__4 = jr;
544 q__2.r = q__3.r * work[i__4].r - q__3.i * work[i__4]
545 .i, q__2.i = q__3.r * work[i__4].i + q__3.i *
546 work[i__4].r;
547 q__1.r = sumb.r + q__2.r, q__1.i = sumb.i + q__2.i;
548 sumb.r = q__1.r, sumb.i = q__1.i;
549 /* L80: */
550 }
551 q__2.r = acoeff * suma.r, q__2.i = acoeff * suma.i;
552 r_cnjg(&q__4, &bcoeff);
553 q__3.r = q__4.r * sumb.r - q__4.i * sumb.i, q__3.i =
554 q__4.r * sumb.i + q__4.i * sumb.r;
555 q__1.r = q__2.r - q__3.r, q__1.i = q__2.i - q__3.i;
556 sum.r = q__1.r, sum.i = q__1.i;
557
558 /* Form x(j) = - SUM / conjg( a*A(j,j) - b*B(j,j) ) */
559
560 /* with scaling and perturbation of the denominator */
561
562 i__3 = j + j * a_dim1;
563 q__3.r = acoeff * a[i__3].r, q__3.i = acoeff * a[i__3].i;
564 i__4 = j + j * b_dim1;
565 q__4.r = bcoeff.r * b[i__4].r - bcoeff.i * b[i__4].i,
566 q__4.i = bcoeff.r * b[i__4].i + bcoeff.i * b[i__4]
567 .r;
568 q__2.r = q__3.r - q__4.r, q__2.i = q__3.i - q__4.i;
569 r_cnjg(&q__1, &q__2);
570 d__.r = q__1.r, d__.i = q__1.i;
571 if ((r__1 = d__.r, dabs(r__1)) + (r__2 = r_imag(&d__),
572 dabs(r__2)) <= dmin__) {
573 q__1.r = dmin__, q__1.i = 0.f;
574 d__.r = q__1.r, d__.i = q__1.i;
575 }
576
577 if ((r__1 = d__.r, dabs(r__1)) + (r__2 = r_imag(&d__),
578 dabs(r__2)) < 1.f) {
579 if ((r__1 = sum.r, dabs(r__1)) + (r__2 = r_imag(&sum),
580 dabs(r__2)) >= bignum * ((r__3 = d__.r, dabs(
581 r__3)) + (r__4 = r_imag(&d__), dabs(r__4)))) {
582 temp = 1.f / ((r__1 = sum.r, dabs(r__1)) + (r__2 =
583 r_imag(&sum), dabs(r__2)));
584 i__3 = j - 1;
585 for (jr = je; jr <= i__3; ++jr) {
586 i__4 = jr;
587 i__5 = jr;
588 q__1.r = temp * work[i__5].r, q__1.i = temp *
589 work[i__5].i;
590 work[i__4].r = q__1.r, work[i__4].i = q__1.i;
591 /* L90: */
592 }
593 xmax = temp * xmax;
594 q__1.r = temp * sum.r, q__1.i = temp * sum.i;
595 sum.r = q__1.r, sum.i = q__1.i;
596 }
597 }
598 i__3 = j;
599 q__2.r = -sum.r, q__2.i = -sum.i;
600 cladiv_(&q__1, &q__2, &d__);
601 work[i__3].r = q__1.r, work[i__3].i = q__1.i;
602 /* Computing MAX */
603 i__3 = j;
604 r__3 = xmax, r__4 = (r__1 = work[i__3].r, dabs(r__1)) + (
605 r__2 = r_imag(&work[j]), dabs(r__2));
606 xmax = dmax(r__3,r__4);
607 /* L100: */
608 }
609
610 /* Back transform eigenvector if HOWMNY='B'. */
611
612 if (ilback) {
613 i__2 = *n + 1 - je;
614 cgemv_("N", n, &i__2, &c_b2, &vl[je * vl_dim1 + 1], ldvl,
615 &work[je], &c__1, &c_b1, &work[*n + 1], &c__1, (
616 ftnlen)1);
617 isrc = 2;
618 ibeg = 1;
619 } else {
620 isrc = 1;
621 ibeg = je;
622 }
623
624 /* Copy and scale eigenvector into column of VL */
625
626 xmax = 0.f;
627 i__2 = *n;
628 for (jr = ibeg; jr <= i__2; ++jr) {
629 /* Computing MAX */
630 i__3 = (isrc - 1) * *n + jr;
631 r__3 = xmax, r__4 = (r__1 = work[i__3].r, dabs(r__1)) + (
632 r__2 = r_imag(&work[(isrc - 1) * *n + jr]), dabs(
633 r__2));
634 xmax = dmax(r__3,r__4);
635 /* L110: */
636 }
637
638 if (xmax > safmin) {
639 temp = 1.f / xmax;
640 i__2 = *n;
641 for (jr = ibeg; jr <= i__2; ++jr) {
642 i__3 = jr + ieig * vl_dim1;
643 i__4 = (isrc - 1) * *n + jr;
644 q__1.r = temp * work[i__4].r, q__1.i = temp * work[
645 i__4].i;
646 vl[i__3].r = q__1.r, vl[i__3].i = q__1.i;
647 /* L120: */
648 }
649 } else {
650 ibeg = *n + 1;
651 }
652
653 i__2 = ibeg - 1;
654 for (jr = 1; jr <= i__2; ++jr) {
655 i__3 = jr + ieig * vl_dim1;
656 vl[i__3].r = 0.f, vl[i__3].i = 0.f;
657 /* L130: */
658 }
659
660 }
661 L140:
662 ;
663 }
664 }
665
666 /* Right eigenvectors */
667
668 if (compr) {
669 ieig = im + 1;
670
671 /* Main loop over eigenvalues */
672
673 for (je = *n; je >= 1; --je) {
674 if (ilall) {
675 ilcomp = TRUE_;
676 } else {
677 ilcomp = select[je];
678 }
679 if (ilcomp) {
680 --ieig;
681
682 i__1 = je + je * a_dim1;
683 i__2 = je + je * b_dim1;
684 if ((r__2 = a[i__1].r, dabs(r__2)) + (r__3 = r_imag(&a[je +
685 je * a_dim1]), dabs(r__3)) <= safmin && (r__1 = b[
686 i__2].r, dabs(r__1)) <= safmin) {
687
688 /* Singular matrix pencil -- return unit eigenvector */
689
690 i__1 = *n;
691 for (jr = 1; jr <= i__1; ++jr) {
692 i__2 = jr + ieig * vr_dim1;
693 vr[i__2].r = 0.f, vr[i__2].i = 0.f;
694 /* L150: */
695 }
696 i__1 = ieig + ieig * vr_dim1;
697 vr[i__1].r = 1.f, vr[i__1].i = 0.f;
698 goto L250;
699 }
700
701 /* Non-singular eigenvalue: */
702 /* Compute coefficients a and b in */
703
704 /* ( a A - b B ) x = 0 */
705
706 /* Computing MAX */
707 i__1 = je + je * a_dim1;
708 i__2 = je + je * b_dim1;
709 r__4 = ((r__2 = a[i__1].r, dabs(r__2)) + (r__3 = r_imag(&a[je
710 + je * a_dim1]), dabs(r__3))) * ascale, r__5 = (r__1 =
711 b[i__2].r, dabs(r__1)) * bscale, r__4 = max(r__4,
712 r__5);
713 temp = 1.f / dmax(r__4,safmin);
714 i__1 = je + je * a_dim1;
715 q__2.r = temp * a[i__1].r, q__2.i = temp * a[i__1].i;
716 q__1.r = ascale * q__2.r, q__1.i = ascale * q__2.i;
717 salpha.r = q__1.r, salpha.i = q__1.i;
718 i__1 = je + je * b_dim1;
719 sbeta = temp * b[i__1].r * bscale;
720 acoeff = sbeta * ascale;
721 q__1.r = bscale * salpha.r, q__1.i = bscale * salpha.i;
722 bcoeff.r = q__1.r, bcoeff.i = q__1.i;
723
724 /* Scale to avoid underflow */
725
726 lsa = dabs(sbeta) >= safmin && dabs(acoeff) < small;
727 lsb = (r__1 = salpha.r, dabs(r__1)) + (r__2 = r_imag(&salpha),
728 dabs(r__2)) >= safmin && (r__3 = bcoeff.r, dabs(r__3)
729 ) + (r__4 = r_imag(&bcoeff), dabs(r__4)) < small;
730
731 scale = 1.f;
732 if (lsa) {
733 scale = small / dabs(sbeta) * dmin(anorm,big);
734 }
735 if (lsb) {
736 /* Computing MAX */
737 r__3 = scale, r__4 = small / ((r__1 = salpha.r, dabs(r__1)
738 ) + (r__2 = r_imag(&salpha), dabs(r__2))) * dmin(
739 bnorm,big);
740 scale = dmax(r__3,r__4);
741 }
742 if (lsa || lsb) {
743 /* Computing MIN */
744 /* Computing MAX */
745 r__5 = 1.f, r__6 = dabs(acoeff), r__5 = max(r__5,r__6),
746 r__6 = (r__1 = bcoeff.r, dabs(r__1)) + (r__2 =
747 r_imag(&bcoeff), dabs(r__2));
748 r__3 = scale, r__4 = 1.f / (safmin * dmax(r__5,r__6));
749 scale = dmin(r__3,r__4);
750 if (lsa) {
751 acoeff = ascale * (scale * sbeta);
752 } else {
753 acoeff = scale * acoeff;
754 }
755 if (lsb) {
756 q__2.r = scale * salpha.r, q__2.i = scale * salpha.i;
757 q__1.r = bscale * q__2.r, q__1.i = bscale * q__2.i;
758 bcoeff.r = q__1.r, bcoeff.i = q__1.i;
759 } else {
760 q__1.r = scale * bcoeff.r, q__1.i = scale * bcoeff.i;
761 bcoeff.r = q__1.r, bcoeff.i = q__1.i;
762 }
763 }
764
765 acoefa = dabs(acoeff);
766 bcoefa = (r__1 = bcoeff.r, dabs(r__1)) + (r__2 = r_imag(&
767 bcoeff), dabs(r__2));
768 xmax = 1.f;
769 i__1 = *n;
770 for (jr = 1; jr <= i__1; ++jr) {
771 i__2 = jr;
772 work[i__2].r = 0.f, work[i__2].i = 0.f;
773 /* L160: */
774 }
775 i__1 = je;
776 work[i__1].r = 1.f, work[i__1].i = 0.f;
777 /* Computing MAX */
778 r__1 = ulp * acoefa * anorm, r__2 = ulp * bcoefa * bnorm,
779 r__1 = max(r__1,r__2);
780 dmin__ = dmax(r__1,safmin);
781
782 /* Triangular solve of (a A - b B) x = 0 (columnwise) */
783
784 /* WORK(1:j-1) contains sums w, */
785 /* WORK(j+1:JE) contains x */
786
787 i__1 = je - 1;
788 for (jr = 1; jr <= i__1; ++jr) {
789 i__2 = jr;
790 i__3 = jr + je * a_dim1;
791 q__2.r = acoeff * a[i__3].r, q__2.i = acoeff * a[i__3].i;
792 i__4 = jr + je * b_dim1;
793 q__3.r = bcoeff.r * b[i__4].r - bcoeff.i * b[i__4].i,
794 q__3.i = bcoeff.r * b[i__4].i + bcoeff.i * b[i__4]
795 .r;
796 q__1.r = q__2.r - q__3.r, q__1.i = q__2.i - q__3.i;
797 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
798 /* L170: */
799 }
800 i__1 = je;
801 work[i__1].r = 1.f, work[i__1].i = 0.f;
802
803 for (j = je - 1; j >= 1; --j) {
804
805 /* Form x(j) := - w(j) / d */
806 /* with scaling and perturbation of the denominator */
807
808 i__1 = j + j * a_dim1;
809 q__2.r = acoeff * a[i__1].r, q__2.i = acoeff * a[i__1].i;
810 i__2 = j + j * b_dim1;
811 q__3.r = bcoeff.r * b[i__2].r - bcoeff.i * b[i__2].i,
812 q__3.i = bcoeff.r * b[i__2].i + bcoeff.i * b[i__2]
813 .r;
814 q__1.r = q__2.r - q__3.r, q__1.i = q__2.i - q__3.i;
815 d__.r = q__1.r, d__.i = q__1.i;
816 if ((r__1 = d__.r, dabs(r__1)) + (r__2 = r_imag(&d__),
817 dabs(r__2)) <= dmin__) {
818 q__1.r = dmin__, q__1.i = 0.f;
819 d__.r = q__1.r, d__.i = q__1.i;
820 }
821
822 if ((r__1 = d__.r, dabs(r__1)) + (r__2 = r_imag(&d__),
823 dabs(r__2)) < 1.f) {
824 i__1 = j;
825 if ((r__1 = work[i__1].r, dabs(r__1)) + (r__2 =
826 r_imag(&work[j]), dabs(r__2)) >= bignum * ((
827 r__3 = d__.r, dabs(r__3)) + (r__4 = r_imag(&
828 d__), dabs(r__4)))) {
829 i__1 = j;
830 temp = 1.f / ((r__1 = work[i__1].r, dabs(r__1)) +
831 (r__2 = r_imag(&work[j]), dabs(r__2)));
832 i__1 = je;
833 for (jr = 1; jr <= i__1; ++jr) {
834 i__2 = jr;
835 i__3 = jr;
836 q__1.r = temp * work[i__3].r, q__1.i = temp *
837 work[i__3].i;
838 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
839 /* L180: */
840 }
841 }
842 }
843
844 i__1 = j;
845 i__2 = j;
846 q__2.r = -work[i__2].r, q__2.i = -work[i__2].i;
847 cladiv_(&q__1, &q__2, &d__);
848 work[i__1].r = q__1.r, work[i__1].i = q__1.i;
849
850 if (j > 1) {
851
852 /* w = w + x(j)*(a A(*,j) - b B(*,j) ) with scaling */
853
854 i__1 = j;
855 if ((r__1 = work[i__1].r, dabs(r__1)) + (r__2 =
856 r_imag(&work[j]), dabs(r__2)) > 1.f) {
857 i__1 = j;
858 temp = 1.f / ((r__1 = work[i__1].r, dabs(r__1)) +
859 (r__2 = r_imag(&work[j]), dabs(r__2)));
860 if (acoefa * rwork[j] + bcoefa * rwork[*n + j] >=
861 bignum * temp) {
862 i__1 = je;
863 for (jr = 1; jr <= i__1; ++jr) {
864 i__2 = jr;
865 i__3 = jr;
866 q__1.r = temp * work[i__3].r, q__1.i =
867 temp * work[i__3].i;
868 work[i__2].r = q__1.r, work[i__2].i =
869 q__1.i;
870 /* L190: */
871 }
872 }
873 }
874
875 i__1 = j;
876 q__1.r = acoeff * work[i__1].r, q__1.i = acoeff *
877 work[i__1].i;
878 ca.r = q__1.r, ca.i = q__1.i;
879 i__1 = j;
880 q__1.r = bcoeff.r * work[i__1].r - bcoeff.i * work[
881 i__1].i, q__1.i = bcoeff.r * work[i__1].i +
882 bcoeff.i * work[i__1].r;
883 cb.r = q__1.r, cb.i = q__1.i;
884 i__1 = j - 1;
885 for (jr = 1; jr <= i__1; ++jr) {
886 i__2 = jr;
887 i__3 = jr;
888 i__4 = jr + j * a_dim1;
889 q__3.r = ca.r * a[i__4].r - ca.i * a[i__4].i,
890 q__3.i = ca.r * a[i__4].i + ca.i * a[i__4]
891 .r;
892 q__2.r = work[i__3].r + q__3.r, q__2.i = work[
893 i__3].i + q__3.i;
894 i__5 = jr + j * b_dim1;
895 q__4.r = cb.r * b[i__5].r - cb.i * b[i__5].i,
896 q__4.i = cb.r * b[i__5].i + cb.i * b[i__5]
897 .r;
898 q__1.r = q__2.r - q__4.r, q__1.i = q__2.i -
899 q__4.i;
900 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
901 /* L200: */
902 }
903 }
904 /* L210: */
905 }
906
907 /* Back transform eigenvector if HOWMNY='B'. */
908
909 if (ilback) {
910 cgemv_("N", n, &je, &c_b2, &vr[vr_offset], ldvr, &work[1],
911 &c__1, &c_b1, &work[*n + 1], &c__1, (ftnlen)1);
912 isrc = 2;
913 iend = *n;
914 } else {
915 isrc = 1;
916 iend = je;
917 }
918
919 /* Copy and scale eigenvector into column of VR */
920
921 xmax = 0.f;
922 i__1 = iend;
923 for (jr = 1; jr <= i__1; ++jr) {
924 /* Computing MAX */
925 i__2 = (isrc - 1) * *n + jr;
926 r__3 = xmax, r__4 = (r__1 = work[i__2].r, dabs(r__1)) + (
927 r__2 = r_imag(&work[(isrc - 1) * *n + jr]), dabs(
928 r__2));
929 xmax = dmax(r__3,r__4);
930 /* L220: */
931 }
932
933 if (xmax > safmin) {
934 temp = 1.f / xmax;
935 i__1 = iend;
936 for (jr = 1; jr <= i__1; ++jr) {
937 i__2 = jr + ieig * vr_dim1;
938 i__3 = (isrc - 1) * *n + jr;
939 q__1.r = temp * work[i__3].r, q__1.i = temp * work[
940 i__3].i;
941 vr[i__2].r = q__1.r, vr[i__2].i = q__1.i;
942 /* L230: */
943 }
944 } else {
945 iend = 0;
946 }
947
948 i__1 = *n;
949 for (jr = iend + 1; jr <= i__1; ++jr) {
950 i__2 = jr + ieig * vr_dim1;
951 vr[i__2].r = 0.f, vr[i__2].i = 0.f;
952 /* L240: */
953 }
954
955 }
956 L250:
957 ;
958 }
959 }
960
961 return 0;
962
963 /* End of CTGEVC */
964
965 } /* ctgevc_ */
966
967