1 /* ../netlib/ctgevc.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
3 #include "FLA_f2c.h" /* Table of constant values */
4 static complex c_b1 =
5 {
6 0.f,0.f
7 }
8 ;
9 static complex c_b2 =
10 {
11 1.f,0.f
12 }
13 ;
14 static integer c__1 = 1;
15 /* > \brief \b CTGEVC */
16 /* =========== DOCUMENTATION =========== */
17 /* Online html documentation available at */
18 /* http://www.netlib.org/lapack/explore-html/ */
19 /* > \htmlonly */
20 /* > Download CTGEVC + dependencies */
21 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgevc. f"> */
22 /* > [TGZ]</a> */
23 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgevc. f"> */
24 /* > [ZIP]</a> */
25 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgevc. f"> */
26 /* > [TXT]</a> */
27 /* > \endhtmlonly */
28 /* Definition: */
29 /* =========== */
30 /* SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, */
31 /* LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO ) */
32 /* .. Scalar Arguments .. */
33 /* CHARACTER HOWMNY, SIDE */
34 /* INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N */
35 /* .. */
36 /* .. Array Arguments .. */
37 /* LOGICAL SELECT( * ) */
38 /* REAL RWORK( * ) */
39 /* COMPLEX P( LDP, * ), S( LDS, * ), VL( LDVL, * ), */
40 /* $ VR( LDVR, * ), WORK( * ) */
41 /* .. */
42 /* > \par Purpose: */
43 /* ============= */
44 /* > */
45 /* > \verbatim */
46 /* > */
47 /* > CTGEVC computes some or all of the right and/or left eigenvectors of */
48 /* > a pair of complex matrices (S,P), where S and P are upper triangular. */
49 /* > Matrix pairs of this type are produced by the generalized Schur */
50 /* > factorization of a complex matrix pair (A,B): */
51 /* > */
52 /* > A = Q*S*Z**H, B = Q*P*Z**H */
53 /* > */
54 /* > as computed by CGGHRD + CHGEQZ. */
55 /* > */
56 /* > The right eigenvector x and the left eigenvector y of (S,P) */
57 /* > corresponding to an eigenvalue w are defined by: */
58 /* > */
59 /* > S*x = w*P*x, (y**H)*S = w*(y**H)*P, */
60 /* > */
61 /* > where y**H denotes the conjugate tranpose of y. */
62 /* > The eigenvalues are not input to this routine, but are computed */
63 /* > directly from the diagonal elements of S and P. */
64 /* > */
65 /* > This routine returns the matrices X and/or Y of right and left */
66 /* > eigenvectors of (S,P), or the products Z*X and/or Q*Y, */
67 /* > where Z and Q are input matrices. */
68 /* > If Q and Z are the unitary factors from the generalized Schur */
69 /* > factorization of a matrix pair (A,B), then Z*X and Q*Y */
70 /* > are the matrices of right and left eigenvectors of (A,B). */
71 /* > \endverbatim */
72 /* Arguments: */
73 /* ========== */
74 /* > \param[in] SIDE */
75 /* > \verbatim */
76 /* > SIDE is CHARACTER*1 */
77 /* > = 'R': compute right eigenvectors only;
78 */
79 /* > = 'L': compute left eigenvectors only;
80 */
81 /* > = 'B': compute both right and left eigenvectors. */
82 /* > \endverbatim */
83 /* > */
84 /* > \param[in] HOWMNY */
85 /* > \verbatim */
86 /* > HOWMNY is CHARACTER*1 */
87 /* > = 'A': compute all right and/or left eigenvectors;
88 */
89 /* > = 'B': compute all right and/or left eigenvectors, */
90 /* > backtransformed by the matrices in VR and/or VL;
91 */
92 /* > = 'S': compute selected right and/or left eigenvectors, */
93 /* > specified by the logical array SELECT. */
94 /* > \endverbatim */
95 /* > */
96 /* > \param[in] SELECT */
97 /* > \verbatim */
98 /* > SELECT is LOGICAL array, dimension (N) */
99 /* > If HOWMNY='S', SELECT specifies the eigenvectors to be */
100 /* > computed. The eigenvector corresponding to the j-th */
101 /* > eigenvalue is computed if SELECT(j) = .TRUE.. */
102 /* > Not referenced if HOWMNY = 'A' or 'B'. */
103 /* > \endverbatim */
104 /* > */
105 /* > \param[in] N */
106 /* > \verbatim */
107 /* > N is INTEGER */
108 /* > The order of the matrices S and P. N >= 0. */
109 /* > \endverbatim */
110 /* > */
111 /* > \param[in] S */
112 /* > \verbatim */
113 /* > S is COMPLEX array, dimension (LDS,N) */
114 /* > The upper triangular matrix S from a generalized Schur */
115 /* > factorization, as computed by CHGEQZ. */
116 /* > \endverbatim */
117 /* > */
118 /* > \param[in] LDS */
119 /* > \verbatim */
120 /* > LDS is INTEGER */
121 /* > The leading dimension of array S. LDS >= max(1,N). */
122 /* > \endverbatim */
123 /* > */
124 /* > \param[in] P */
125 /* > \verbatim */
126 /* > P is COMPLEX array, dimension (LDP,N) */
127 /* > The upper triangular matrix P from a generalized Schur */
128 /* > factorization, as computed by CHGEQZ. P must have real */
129 /* > diagonal elements. */
130 /* > \endverbatim */
131 /* > */
132 /* > \param[in] LDP */
133 /* > \verbatim */
134 /* > LDP is INTEGER */
135 /* > The leading dimension of array P. LDP >= max(1,N). */
136 /* > \endverbatim */
137 /* > */
138 /* > \param[in,out] VL */
139 /* > \verbatim */
140 /* > VL is COMPLEX array, dimension (LDVL,MM) */
141 /* > On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must */
142 /* > contain an N-by-N matrix Q (usually the unitary matrix Q */
143 /* > of left Schur vectors returned by CHGEQZ). */
144 /* > On exit, if SIDE = 'L' or 'B', VL contains: */
145 /* > if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
146 */
147 /* > if HOWMNY = 'B', the matrix Q*Y;
148 */
149 /* > if HOWMNY = 'S', the left eigenvectors of (S,P) specified by */
150 /* > SELECT, stored consecutively in the columns of */
151 /* > VL, in the same order as their eigenvalues. */
152 /* > Not referenced if SIDE = 'R'. */
153 /* > \endverbatim */
154 /* > */
155 /* > \param[in] LDVL */
156 /* > \verbatim */
157 /* > LDVL is INTEGER */
158 /* > The leading dimension of array VL. LDVL >= 1, and if */
159 /* > SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N. */
160 /* > \endverbatim */
161 /* > */
162 /* > \param[in,out] VR */
163 /* > \verbatim */
164 /* > VR is COMPLEX array, dimension (LDVR,MM) */
165 /* > On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must */
166 /* > contain an N-by-N matrix Q (usually the unitary matrix Z */
167 /* > of right Schur vectors returned by CHGEQZ). */
168 /* > On exit, if SIDE = 'R' or 'B', VR contains: */
169 /* > if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
170 */
171 /* > if HOWMNY = 'B', the matrix Z*X;
172 */
173 /* > if HOWMNY = 'S', the right eigenvectors of (S,P) specified by */
174 /* > SELECT, stored consecutively in the columns of */
175 /* > VR, in the same order as their eigenvalues. */
176 /* > Not referenced if SIDE = 'L'. */
177 /* > \endverbatim */
178 /* > */
179 /* > \param[in] LDVR */
180 /* > \verbatim */
181 /* > LDVR is INTEGER */
182 /* > The leading dimension of the array VR. LDVR >= 1, and if */
183 /* > SIDE = 'R' or 'B', LDVR >= N. */
184 /* > \endverbatim */
185 /* > */
186 /* > \param[in] MM */
187 /* > \verbatim */
188 /* > MM is INTEGER */
189 /* > The number of columns in the arrays VL and/or VR. MM >= M. */
190 /* > \endverbatim */
191 /* > */
192 /* > \param[out] M */
193 /* > \verbatim */
194 /* > M is INTEGER */
195 /* > The number of columns in the arrays VL and/or VR actually */
196 /* > used to store the eigenvectors. If HOWMNY = 'A' or 'B', M */
197 /* > is set to N. Each selected eigenvector occupies one column. */
198 /* > \endverbatim */
199 /* > */
200 /* > \param[out] WORK */
201 /* > \verbatim */
202 /* > WORK is COMPLEX array, dimension (2*N) */
203 /* > \endverbatim */
204 /* > */
205 /* > \param[out] RWORK */
206 /* > \verbatim */
207 /* > RWORK is REAL array, dimension (2*N) */
208 /* > \endverbatim */
209 /* > */
210 /* > \param[out] INFO */
211 /* > \verbatim */
212 /* > INFO is INTEGER */
213 /* > = 0: successful exit. */
214 /* > < 0: if INFO = -i, the i-th argument had an illegal value. */
215 /* > \endverbatim */
216 /* Authors: */
217 /* ======== */
218 /* > \author Univ. of Tennessee */
219 /* > \author Univ. of California Berkeley */
220 /* > \author Univ. of Colorado Denver */
221 /* > \author NAG Ltd. */
222 /* > \date November 2011 */
223 /* > \ingroup complexGEcomputational */
224 /* ===================================================================== */
225 /* Subroutine */
ctgevc_(char * side,char * howmny,logical * select,integer * n,complex * s,integer * lds,complex * p,integer * ldp,complex * vl,integer * ldvl,complex * vr,integer * ldvr,integer * mm,integer * m,complex * work,real * rwork,integer * info)226 int ctgevc_(char *side, char *howmny, logical *select, integer *n, complex *s, integer *lds, complex *p, integer *ldp, complex *vl, integer *ldvl, complex *vr, integer *ldvr, integer *mm, integer *m, complex *work, real *rwork, integer *info)
227 {
228 /* System generated locals */
229 integer p_dim1, p_offset, s_dim1, s_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, i__2, i__3, i__4, i__5;
230 real r__1, r__2, r__3, r__4, r__5, r__6;
231 complex q__1, q__2, q__3, q__4;
232 /* Builtin functions */
233 double r_imag(complex *);
234 void r_cnjg(complex *, complex *);
235 /* Local variables */
236 complex d__;
237 integer i__, j;
238 complex ca, cb;
239 integer je, im, jr;
240 real big;
241 logical lsa, lsb;
242 real ulp;
243 complex sum;
244 integer ibeg, ieig, iend;
245 real dmin__;
246 integer isrc;
247 real temp;
248 complex suma, sumb;
249 real xmax, scale;
250 logical ilall;
251 integer iside;
252 real sbeta;
253 extern logical lsame_(char *, char *);
254 extern /* Subroutine */
255 int cgemv_(char *, integer *, integer *, complex * , complex *, integer *, complex *, integer *, complex *, complex * , integer *);
256 real small;
257 logical compl;
258 real anorm, bnorm;
259 logical compr, ilbbad;
260 real acoefa, bcoefa, acoeff;
261 complex bcoeff;
262 logical ilback;
263 extern /* Subroutine */
264 int slabad_(real *, real *);
265 real ascale, bscale;
266 extern /* Complex */
267 VOID cladiv_(complex *, complex *, complex *);
268 extern real slamch_(char *);
269 complex salpha;
270 real safmin;
271 extern /* Subroutine */
272 int xerbla_(char *, integer *);
273 real bignum;
274 logical ilcomp;
275 integer ihwmny;
276 /* -- LAPACK computational routine (version 3.4.0) -- */
277 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
278 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
279 /* November 2011 */
280 /* .. Scalar Arguments .. */
281 /* .. */
282 /* .. Array Arguments .. */
283 /* .. */
284 /* ===================================================================== */
285 /* .. Parameters .. */
286 /* .. */
287 /* .. Local Scalars .. */
288 /* .. */
289 /* .. External Functions .. */
290 /* .. */
291 /* .. External Subroutines .. */
292 /* .. */
293 /* .. Intrinsic Functions .. */
294 /* .. */
295 /* .. Statement Functions .. */
296 /* .. */
297 /* .. Statement Function definitions .. */
298 /* .. */
299 /* .. Executable Statements .. */
300 /* Decode and Test the input parameters */
301 /* Parameter adjustments */
302 --select;
303 s_dim1 = *lds;
304 s_offset = 1 + s_dim1;
305 s -= s_offset;
306 p_dim1 = *ldp;
307 p_offset = 1 + p_dim1;
308 p -= p_offset;
309 vl_dim1 = *ldvl;
310 vl_offset = 1 + vl_dim1;
311 vl -= vl_offset;
312 vr_dim1 = *ldvr;
313 vr_offset = 1 + vr_dim1;
314 vr -= vr_offset;
315 --work;
316 --rwork;
317 /* Function Body */
318 if (lsame_(howmny, "A"))
319 {
320 ihwmny = 1;
321 ilall = TRUE_;
322 ilback = FALSE_;
323 }
324 else if (lsame_(howmny, "S"))
325 {
326 ihwmny = 2;
327 ilall = FALSE_;
328 ilback = FALSE_;
329 }
330 else if (lsame_(howmny, "B"))
331 {
332 ihwmny = 3;
333 ilall = TRUE_;
334 ilback = TRUE_;
335 }
336 else
337 {
338 ihwmny = -1;
339 }
340 if (lsame_(side, "R"))
341 {
342 iside = 1;
343 compl = FALSE_;
344 compr = TRUE_;
345 }
346 else if (lsame_(side, "L"))
347 {
348 iside = 2;
349 compl = TRUE_;
350 compr = FALSE_;
351 }
352 else if (lsame_(side, "B"))
353 {
354 iside = 3;
355 compl = TRUE_;
356 compr = TRUE_;
357 }
358 else
359 {
360 iside = -1;
361 }
362 *info = 0;
363 if (iside < 0)
364 {
365 *info = -1;
366 }
367 else if (ihwmny < 0)
368 {
369 *info = -2;
370 }
371 else if (*n < 0)
372 {
373 *info = -4;
374 }
375 else if (*lds < max(1,*n))
376 {
377 *info = -6;
378 }
379 else if (*ldp < max(1,*n))
380 {
381 *info = -8;
382 }
383 if (*info != 0)
384 {
385 i__1 = -(*info);
386 xerbla_("CTGEVC", &i__1);
387 return 0;
388 }
389 /* Count the number of eigenvectors */
390 if (! ilall)
391 {
392 im = 0;
393 i__1 = *n;
394 for (j = 1;
395 j <= i__1;
396 ++j)
397 {
398 if (select[j])
399 {
400 ++im;
401 }
402 /* L10: */
403 }
404 }
405 else
406 {
407 im = *n;
408 }
409 /* Check diagonal of B */
410 ilbbad = FALSE_;
411 i__1 = *n;
412 for (j = 1;
413 j <= i__1;
414 ++j)
415 {
416 if (r_imag(&p[j + j * p_dim1]) != 0.f)
417 {
418 ilbbad = TRUE_;
419 }
420 /* L20: */
421 }
422 if (ilbbad)
423 {
424 *info = -7;
425 }
426 else if (compl && *ldvl < *n || *ldvl < 1)
427 {
428 *info = -10;
429 }
430 else if (compr && *ldvr < *n || *ldvr < 1)
431 {
432 *info = -12;
433 }
434 else if (*mm < im)
435 {
436 *info = -13;
437 }
438 if (*info != 0)
439 {
440 i__1 = -(*info);
441 xerbla_("CTGEVC", &i__1);
442 return 0;
443 }
444 /* Quick return if possible */
445 *m = im;
446 if (*n == 0)
447 {
448 return 0;
449 }
450 /* Machine Constants */
451 safmin = slamch_("Safe minimum");
452 big = 1.f / safmin;
453 slabad_(&safmin, &big);
454 ulp = slamch_("Epsilon") * slamch_("Base");
455 small = safmin * *n / ulp;
456 big = 1.f / small;
457 bignum = 1.f / (safmin * *n);
458 /* Compute the 1-norm of each column of the strictly upper triangular */
459 /* part of A and B to check for possible overflow in the triangular */
460 /* solver. */
461 i__1 = s_dim1 + 1;
462 anorm = (r__1 = s[i__1].r, f2c_abs(r__1)) + (r__2 = r_imag(&s[s_dim1 + 1]), f2c_abs(r__2));
463 i__1 = p_dim1 + 1;
464 bnorm = (r__1 = p[i__1].r, f2c_abs(r__1)) + (r__2 = r_imag(&p[p_dim1 + 1]), f2c_abs(r__2));
465 rwork[1] = 0.f;
466 rwork[*n + 1] = 0.f;
467 i__1 = *n;
468 for (j = 2;
469 j <= i__1;
470 ++j)
471 {
472 rwork[j] = 0.f;
473 rwork[*n + j] = 0.f;
474 i__2 = j - 1;
475 for (i__ = 1;
476 i__ <= i__2;
477 ++i__)
478 {
479 i__3 = i__ + j * s_dim1;
480 rwork[j] += (r__1 = s[i__3].r, f2c_abs(r__1)) + (r__2 = r_imag(&s[i__ + j * s_dim1]), f2c_abs(r__2));
481 i__3 = i__ + j * p_dim1;
482 rwork[*n + j] += (r__1 = p[i__3].r, f2c_abs(r__1)) + (r__2 = r_imag(& p[i__ + j * p_dim1]), f2c_abs(r__2));
483 /* L30: */
484 }
485 /* Computing MAX */
486 i__2 = j + j * s_dim1;
487 r__3 = anorm;
488 r__4 = rwork[j] + ((r__1 = s[i__2].r, f2c_abs(r__1)) + ( r__2 = r_imag(&s[j + j * s_dim1]), f2c_abs(r__2))); // , expr subst
489 anorm = max(r__3,r__4);
490 /* Computing MAX */
491 i__2 = j + j * p_dim1;
492 r__3 = bnorm;
493 r__4 = rwork[*n + j] + ((r__1 = p[i__2].r, f2c_abs(r__1)) + (r__2 = r_imag(&p[j + j * p_dim1]), f2c_abs(r__2))); // , expr subst
494 bnorm = max(r__3,r__4);
495 /* L40: */
496 }
497 ascale = 1.f / max(anorm,safmin);
498 bscale = 1.f / max(bnorm,safmin);
499 /* Left eigenvectors */
500 if (compl)
501 {
502 ieig = 0;
503 /* Main loop over eigenvalues */
504 i__1 = *n;
505 for (je = 1;
506 je <= i__1;
507 ++je)
508 {
509 if (ilall)
510 {
511 ilcomp = TRUE_;
512 }
513 else
514 {
515 ilcomp = select[je];
516 }
517 if (ilcomp)
518 {
519 ++ieig;
520 i__2 = je + je * s_dim1;
521 i__3 = je + je * p_dim1;
522 if ((r__2 = s[i__2].r, f2c_abs(r__2)) + (r__3 = r_imag(&s[je + je * s_dim1]), f2c_abs(r__3)) <= safmin && (r__1 = p[i__3].r, f2c_abs(r__1)) <= safmin)
523 {
524 /* Singular matrix pencil -- return unit eigenvector */
525 i__2 = *n;
526 for (jr = 1;
527 jr <= i__2;
528 ++jr)
529 {
530 i__3 = jr + ieig * vl_dim1;
531 vl[i__3].r = 0.f;
532 vl[i__3].i = 0.f; // , expr subst
533 /* L50: */
534 }
535 i__2 = ieig + ieig * vl_dim1;
536 vl[i__2].r = 1.f;
537 vl[i__2].i = 0.f; // , expr subst
538 goto L140;
539 }
540 /* Non-singular eigenvalue: */
541 /* Compute coefficients a and b in */
542 /* H */
543 /* y ( a A - b B ) = 0 */
544 /* Computing MAX */
545 i__2 = je + je * s_dim1;
546 i__3 = je + je * p_dim1;
547 r__4 = ((r__2 = s[i__2].r, f2c_abs(r__2)) + (r__3 = r_imag(&s[je + je * s_dim1]), f2c_abs(r__3))) * ascale;
548 r__5 = (r__1 = p[i__3].r, f2c_abs(r__1)) * bscale;
549 r__4 = max(r__4,r__5); // ; expr subst
550 temp = 1.f / max(r__4,safmin);
551 i__2 = je + je * s_dim1;
552 q__2.r = temp * s[i__2].r;
553 q__2.i = temp * s[i__2].i; // , expr subst
554 q__1.r = ascale * q__2.r;
555 q__1.i = ascale * q__2.i; // , expr subst
556 salpha.r = q__1.r;
557 salpha.i = q__1.i; // , expr subst
558 i__2 = je + je * p_dim1;
559 sbeta = temp * p[i__2].r * bscale;
560 acoeff = sbeta * ascale;
561 q__1.r = bscale * salpha.r;
562 q__1.i = bscale * salpha.i; // , expr subst
563 bcoeff.r = q__1.r;
564 bcoeff.i = q__1.i; // , expr subst
565 /* Scale to avoid underflow */
566 lsa = f2c_abs(sbeta) >= safmin && f2c_abs(acoeff) < small;
567 lsb = (r__1 = salpha.r, f2c_abs(r__1)) + (r__2 = r_imag(&salpha), f2c_abs(r__2)) >= safmin && (r__3 = bcoeff.r, f2c_abs(r__3)) + (r__4 = r_imag(&bcoeff), f2c_abs(r__4)) < small;
568 scale = 1.f;
569 if (lsa)
570 {
571 scale = small / f2c_abs(sbeta) * min(anorm,big);
572 }
573 if (lsb)
574 {
575 /* Computing MAX */
576 r__3 = scale;
577 r__4 = small / ((r__1 = salpha.r, f2c_abs(r__1)) + (r__2 = r_imag(&salpha), f2c_abs(r__2))) * min( bnorm,big); // , expr subst
578 scale = max(r__3,r__4);
579 }
580 if (lsa || lsb)
581 {
582 /* Computing MIN */
583 /* Computing MAX */
584 r__5 = 1.f, r__6 = f2c_abs(acoeff);
585 r__5 = max(r__5,r__6);
586 r__6 = (r__1 = bcoeff.r, f2c_abs(r__1)) + (r__2 = r_imag(&bcoeff), f2c_abs(r__2)); // ; expr subst
587 r__3 = scale;
588 r__4 = 1.f / (safmin * max(r__5,r__6)); // , expr subst
589 scale = min(r__3,r__4);
590 if (lsa)
591 {
592 acoeff = ascale * (scale * sbeta);
593 }
594 else
595 {
596 acoeff = scale * acoeff;
597 }
598 if (lsb)
599 {
600 q__2.r = scale * salpha.r;
601 q__2.i = scale * salpha.i; // , expr subst
602 q__1.r = bscale * q__2.r;
603 q__1.i = bscale * q__2.i; // , expr subst
604 bcoeff.r = q__1.r;
605 bcoeff.i = q__1.i; // , expr subst
606 }
607 else
608 {
609 q__1.r = scale * bcoeff.r;
610 q__1.i = scale * bcoeff.i; // , expr subst
611 bcoeff.r = q__1.r;
612 bcoeff.i = q__1.i; // , expr subst
613 }
614 }
615 acoefa = f2c_abs(acoeff);
616 bcoefa = (r__1 = bcoeff.r, f2c_abs(r__1)) + (r__2 = r_imag(& bcoeff), f2c_abs(r__2));
617 xmax = 1.f;
618 i__2 = *n;
619 for (jr = 1;
620 jr <= i__2;
621 ++jr)
622 {
623 i__3 = jr;
624 work[i__3].r = 0.f;
625 work[i__3].i = 0.f; // , expr subst
626 /* L60: */
627 }
628 i__2 = je;
629 work[i__2].r = 1.f;
630 work[i__2].i = 0.f; // , expr subst
631 /* Computing MAX */
632 r__1 = ulp * acoefa * anorm;
633 r__2 = ulp * bcoefa * bnorm;
634 r__1 = max(r__1,r__2); // ; expr subst
635 dmin__ = max(r__1,safmin);
636 /* H */
637 /* Triangular solve of (a A - b B) y = 0 */
638 /* H */
639 /* (rowwise in (a A - b B) , or columnwise in a A - b B) */
640 i__2 = *n;
641 for (j = je + 1;
642 j <= i__2;
643 ++j)
644 {
645 /* Compute */
646 /* j-1 */
647 /* SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k) */
648 /* k=je */
649 /* (Scale if necessary) */
650 temp = 1.f / xmax;
651 if (acoefa * rwork[j] + bcoefa * rwork[*n + j] > bignum * temp)
652 {
653 i__3 = j - 1;
654 for (jr = je;
655 jr <= i__3;
656 ++jr)
657 {
658 i__4 = jr;
659 i__5 = jr;
660 q__1.r = temp * work[i__5].r;
661 q__1.i = temp * work[i__5].i; // , expr subst
662 work[i__4].r = q__1.r;
663 work[i__4].i = q__1.i; // , expr subst
664 /* L70: */
665 }
666 xmax = 1.f;
667 }
668 suma.r = 0.f;
669 suma.i = 0.f; // , expr subst
670 sumb.r = 0.f;
671 sumb.i = 0.f; // , expr subst
672 i__3 = j - 1;
673 for (jr = je;
674 jr <= i__3;
675 ++jr)
676 {
677 r_cnjg(&q__3, &s[jr + j * s_dim1]);
678 i__4 = jr;
679 q__2.r = q__3.r * work[i__4].r - q__3.i * work[i__4] .i;
680 q__2.i = q__3.r * work[i__4].i + q__3.i * work[i__4].r; // , expr subst
681 q__1.r = suma.r + q__2.r;
682 q__1.i = suma.i + q__2.i; // , expr subst
683 suma.r = q__1.r;
684 suma.i = q__1.i; // , expr subst
685 r_cnjg(&q__3, &p[jr + j * p_dim1]);
686 i__4 = jr;
687 q__2.r = q__3.r * work[i__4].r - q__3.i * work[i__4] .i;
688 q__2.i = q__3.r * work[i__4].i + q__3.i * work[i__4].r; // , expr subst
689 q__1.r = sumb.r + q__2.r;
690 q__1.i = sumb.i + q__2.i; // , expr subst
691 sumb.r = q__1.r;
692 sumb.i = q__1.i; // , expr subst
693 /* L80: */
694 }
695 q__2.r = acoeff * suma.r;
696 q__2.i = acoeff * suma.i; // , expr subst
697 r_cnjg(&q__4, &bcoeff);
698 q__3.r = q__4.r * sumb.r - q__4.i * sumb.i;
699 q__3.i = q__4.r * sumb.i + q__4.i * sumb.r; // , expr subst
700 q__1.r = q__2.r - q__3.r;
701 q__1.i = q__2.i - q__3.i; // , expr subst
702 sum.r = q__1.r;
703 sum.i = q__1.i; // , expr subst
704 /* Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) ) */
705 /* with scaling and perturbation of the denominator */
706 i__3 = j + j * s_dim1;
707 q__3.r = acoeff * s[i__3].r;
708 q__3.i = acoeff * s[i__3].i; // , expr subst
709 i__4 = j + j * p_dim1;
710 q__4.r = bcoeff.r * p[i__4].r - bcoeff.i * p[i__4].i;
711 q__4.i = bcoeff.r * p[i__4].i + bcoeff.i * p[i__4] .r; // , expr subst
712 q__2.r = q__3.r - q__4.r;
713 q__2.i = q__3.i - q__4.i; // , expr subst
714 r_cnjg(&q__1, &q__2);
715 d__.r = q__1.r;
716 d__.i = q__1.i; // , expr subst
717 if ((r__1 = d__.r, f2c_abs(r__1)) + (r__2 = r_imag(&d__), f2c_abs( r__2)) <= dmin__)
718 {
719 q__1.r = dmin__;
720 q__1.i = 0.f; // , expr subst
721 d__.r = q__1.r;
722 d__.i = q__1.i; // , expr subst
723 }
724 if ((r__1 = d__.r, f2c_abs(r__1)) + (r__2 = r_imag(&d__), f2c_abs( r__2)) < 1.f)
725 {
726 if ((r__1 = sum.r, f2c_abs(r__1)) + (r__2 = r_imag(&sum), f2c_abs(r__2)) >= bignum * ((r__3 = d__.r, f2c_abs( r__3)) + (r__4 = r_imag(&d__), f2c_abs(r__4))))
727 {
728 temp = 1.f / ((r__1 = sum.r, f2c_abs(r__1)) + (r__2 = r_imag(&sum), f2c_abs(r__2)));
729 i__3 = j - 1;
730 for (jr = je;
731 jr <= i__3;
732 ++jr)
733 {
734 i__4 = jr;
735 i__5 = jr;
736 q__1.r = temp * work[i__5].r;
737 q__1.i = temp * work[i__5].i; // , expr subst
738 work[i__4].r = q__1.r;
739 work[i__4].i = q__1.i; // , expr subst
740 /* L90: */
741 }
742 xmax = temp * xmax;
743 q__1.r = temp * sum.r;
744 q__1.i = temp * sum.i; // , expr subst
745 sum.r = q__1.r;
746 sum.i = q__1.i; // , expr subst
747 }
748 }
749 i__3 = j;
750 q__2.r = -sum.r;
751 q__2.i = -sum.i; // , expr subst
752 cladiv_(&q__1, &q__2, &d__);
753 work[i__3].r = q__1.r;
754 work[i__3].i = q__1.i; // , expr subst
755 /* Computing MAX */
756 i__3 = j;
757 r__3 = xmax;
758 r__4 = (r__1 = work[i__3].r, f2c_abs(r__1)) + ( r__2 = r_imag(&work[j]), f2c_abs(r__2)); // , expr subst
759 xmax = max(r__3,r__4);
760 /* L100: */
761 }
762 /* Back transform eigenvector if HOWMNY='B'. */
763 if (ilback)
764 {
765 i__2 = *n + 1 - je;
766 cgemv_("N", n, &i__2, &c_b2, &vl[je * vl_dim1 + 1], ldvl, &work[je], &c__1, &c_b1, &work[*n + 1], &c__1);
767 isrc = 2;
768 ibeg = 1;
769 }
770 else
771 {
772 isrc = 1;
773 ibeg = je;
774 }
775 /* Copy and scale eigenvector into column of VL */
776 xmax = 0.f;
777 i__2 = *n;
778 for (jr = ibeg;
779 jr <= i__2;
780 ++jr)
781 {
782 /* Computing MAX */
783 i__3 = (isrc - 1) * *n + jr;
784 r__3 = xmax;
785 r__4 = (r__1 = work[i__3].r, f2c_abs(r__1)) + ( r__2 = r_imag(&work[(isrc - 1) * *n + jr]), f2c_abs( r__2)); // , expr subst
786 xmax = max(r__3,r__4);
787 /* L110: */
788 }
789 if (xmax > safmin)
790 {
791 temp = 1.f / xmax;
792 i__2 = *n;
793 for (jr = ibeg;
794 jr <= i__2;
795 ++jr)
796 {
797 i__3 = jr + ieig * vl_dim1;
798 i__4 = (isrc - 1) * *n + jr;
799 q__1.r = temp * work[i__4].r;
800 q__1.i = temp * work[ i__4].i; // , expr subst
801 vl[i__3].r = q__1.r;
802 vl[i__3].i = q__1.i; // , expr subst
803 /* L120: */
804 }
805 }
806 else
807 {
808 ibeg = *n + 1;
809 }
810 i__2 = ibeg - 1;
811 for (jr = 1;
812 jr <= i__2;
813 ++jr)
814 {
815 i__3 = jr + ieig * vl_dim1;
816 vl[i__3].r = 0.f;
817 vl[i__3].i = 0.f; // , expr subst
818 /* L130: */
819 }
820 }
821 L140:
822 ;
823 }
824 }
825 /* Right eigenvectors */
826 if (compr)
827 {
828 ieig = im + 1;
829 /* Main loop over eigenvalues */
830 for (je = *n;
831 je >= 1;
832 --je)
833 {
834 if (ilall)
835 {
836 ilcomp = TRUE_;
837 }
838 else
839 {
840 ilcomp = select[je];
841 }
842 if (ilcomp)
843 {
844 --ieig;
845 i__1 = je + je * s_dim1;
846 i__2 = je + je * p_dim1;
847 if ((r__2 = s[i__1].r, f2c_abs(r__2)) + (r__3 = r_imag(&s[je + je * s_dim1]), f2c_abs(r__3)) <= safmin && (r__1 = p[i__2].r, f2c_abs(r__1)) <= safmin)
848 {
849 /* Singular matrix pencil -- return unit eigenvector */
850 i__1 = *n;
851 for (jr = 1;
852 jr <= i__1;
853 ++jr)
854 {
855 i__2 = jr + ieig * vr_dim1;
856 vr[i__2].r = 0.f;
857 vr[i__2].i = 0.f; // , expr subst
858 /* L150: */
859 }
860 i__1 = ieig + ieig * vr_dim1;
861 vr[i__1].r = 1.f;
862 vr[i__1].i = 0.f; // , expr subst
863 goto L250;
864 }
865 /* Non-singular eigenvalue: */
866 /* Compute coefficients a and b in */
867 /* ( a A - b B ) x = 0 */
868 /* Computing MAX */
869 i__1 = je + je * s_dim1;
870 i__2 = je + je * p_dim1;
871 r__4 = ((r__2 = s[i__1].r, f2c_abs(r__2)) + (r__3 = r_imag(&s[je + je * s_dim1]), f2c_abs(r__3))) * ascale;
872 r__5 = (r__1 = p[i__2].r, f2c_abs(r__1)) * bscale;
873 r__4 = max(r__4,r__5); // ; expr subst
874 temp = 1.f / max(r__4,safmin);
875 i__1 = je + je * s_dim1;
876 q__2.r = temp * s[i__1].r;
877 q__2.i = temp * s[i__1].i; // , expr subst
878 q__1.r = ascale * q__2.r;
879 q__1.i = ascale * q__2.i; // , expr subst
880 salpha.r = q__1.r;
881 salpha.i = q__1.i; // , expr subst
882 i__1 = je + je * p_dim1;
883 sbeta = temp * p[i__1].r * bscale;
884 acoeff = sbeta * ascale;
885 q__1.r = bscale * salpha.r;
886 q__1.i = bscale * salpha.i; // , expr subst
887 bcoeff.r = q__1.r;
888 bcoeff.i = q__1.i; // , expr subst
889 /* Scale to avoid underflow */
890 lsa = f2c_abs(sbeta) >= safmin && f2c_abs(acoeff) < small;
891 lsb = (r__1 = salpha.r, f2c_abs(r__1)) + (r__2 = r_imag(&salpha), f2c_abs(r__2)) >= safmin && (r__3 = bcoeff.r, f2c_abs(r__3)) + (r__4 = r_imag(&bcoeff), f2c_abs(r__4)) < small;
892 scale = 1.f;
893 if (lsa)
894 {
895 scale = small / f2c_abs(sbeta) * min(anorm,big);
896 }
897 if (lsb)
898 {
899 /* Computing MAX */
900 r__3 = scale;
901 r__4 = small / ((r__1 = salpha.r, f2c_abs(r__1)) + (r__2 = r_imag(&salpha), f2c_abs(r__2))) * min( bnorm,big); // , expr subst
902 scale = max(r__3,r__4);
903 }
904 if (lsa || lsb)
905 {
906 /* Computing MIN */
907 /* Computing MAX */
908 r__5 = 1.f, r__6 = f2c_abs(acoeff);
909 r__5 = max(r__5,r__6);
910 r__6 = (r__1 = bcoeff.r, f2c_abs(r__1)) + (r__2 = r_imag(&bcoeff), f2c_abs(r__2)); // ; expr subst
911 r__3 = scale;
912 r__4 = 1.f / (safmin * max(r__5,r__6)); // , expr subst
913 scale = min(r__3,r__4);
914 if (lsa)
915 {
916 acoeff = ascale * (scale * sbeta);
917 }
918 else
919 {
920 acoeff = scale * acoeff;
921 }
922 if (lsb)
923 {
924 q__2.r = scale * salpha.r;
925 q__2.i = scale * salpha.i; // , expr subst
926 q__1.r = bscale * q__2.r;
927 q__1.i = bscale * q__2.i; // , expr subst
928 bcoeff.r = q__1.r;
929 bcoeff.i = q__1.i; // , expr subst
930 }
931 else
932 {
933 q__1.r = scale * bcoeff.r;
934 q__1.i = scale * bcoeff.i; // , expr subst
935 bcoeff.r = q__1.r;
936 bcoeff.i = q__1.i; // , expr subst
937 }
938 }
939 acoefa = f2c_abs(acoeff);
940 bcoefa = (r__1 = bcoeff.r, f2c_abs(r__1)) + (r__2 = r_imag(& bcoeff), f2c_abs(r__2));
941 xmax = 1.f;
942 i__1 = *n;
943 for (jr = 1;
944 jr <= i__1;
945 ++jr)
946 {
947 i__2 = jr;
948 work[i__2].r = 0.f;
949 work[i__2].i = 0.f; // , expr subst
950 /* L160: */
951 }
952 i__1 = je;
953 work[i__1].r = 1.f;
954 work[i__1].i = 0.f; // , expr subst
955 /* Computing MAX */
956 r__1 = ulp * acoefa * anorm;
957 r__2 = ulp * bcoefa * bnorm;
958 r__1 = max(r__1,r__2); // ; expr subst
959 dmin__ = max(r__1,safmin);
960 /* Triangular solve of (a A - b B) x = 0 (columnwise) */
961 /* WORK(1:j-1) contains sums w, */
962 /* WORK(j+1:JE) contains x */
963 i__1 = je - 1;
964 for (jr = 1;
965 jr <= i__1;
966 ++jr)
967 {
968 i__2 = jr;
969 i__3 = jr + je * s_dim1;
970 q__2.r = acoeff * s[i__3].r;
971 q__2.i = acoeff * s[i__3].i; // , expr subst
972 i__4 = jr + je * p_dim1;
973 q__3.r = bcoeff.r * p[i__4].r - bcoeff.i * p[i__4].i;
974 q__3.i = bcoeff.r * p[i__4].i + bcoeff.i * p[i__4] .r; // , expr subst
975 q__1.r = q__2.r - q__3.r;
976 q__1.i = q__2.i - q__3.i; // , expr subst
977 work[i__2].r = q__1.r;
978 work[i__2].i = q__1.i; // , expr subst
979 /* L170: */
980 }
981 i__1 = je;
982 work[i__1].r = 1.f;
983 work[i__1].i = 0.f; // , expr subst
984 for (j = je - 1;
985 j >= 1;
986 --j)
987 {
988 /* Form x(j) := - w(j) / d */
989 /* with scaling and perturbation of the denominator */
990 i__1 = j + j * s_dim1;
991 q__2.r = acoeff * s[i__1].r;
992 q__2.i = acoeff * s[i__1].i; // , expr subst
993 i__2 = j + j * p_dim1;
994 q__3.r = bcoeff.r * p[i__2].r - bcoeff.i * p[i__2].i;
995 q__3.i = bcoeff.r * p[i__2].i + bcoeff.i * p[i__2] .r; // , expr subst
996 q__1.r = q__2.r - q__3.r;
997 q__1.i = q__2.i - q__3.i; // , expr subst
998 d__.r = q__1.r;
999 d__.i = q__1.i; // , expr subst
1000 if ((r__1 = d__.r, f2c_abs(r__1)) + (r__2 = r_imag(&d__), f2c_abs( r__2)) <= dmin__)
1001 {
1002 q__1.r = dmin__;
1003 q__1.i = 0.f; // , expr subst
1004 d__.r = q__1.r;
1005 d__.i = q__1.i; // , expr subst
1006 }
1007 if ((r__1 = d__.r, f2c_abs(r__1)) + (r__2 = r_imag(&d__), f2c_abs( r__2)) < 1.f)
1008 {
1009 i__1 = j;
1010 if ((r__1 = work[i__1].r, f2c_abs(r__1)) + (r__2 = r_imag( &work[j]), f2c_abs(r__2)) >= bignum * ((r__3 = d__.r, f2c_abs(r__3)) + (r__4 = r_imag(&d__), f2c_abs( r__4))))
1011 {
1012 i__1 = j;
1013 temp = 1.f / ((r__1 = work[i__1].r, f2c_abs(r__1)) + ( r__2 = r_imag(&work[j]), f2c_abs(r__2)));
1014 i__1 = je;
1015 for (jr = 1;
1016 jr <= i__1;
1017 ++jr)
1018 {
1019 i__2 = jr;
1020 i__3 = jr;
1021 q__1.r = temp * work[i__3].r;
1022 q__1.i = temp * work[i__3].i; // , expr subst
1023 work[i__2].r = q__1.r;
1024 work[i__2].i = q__1.i; // , expr subst
1025 /* L180: */
1026 }
1027 }
1028 }
1029 i__1 = j;
1030 i__2 = j;
1031 q__2.r = -work[i__2].r;
1032 q__2.i = -work[i__2].i; // , expr subst
1033 cladiv_(&q__1, &q__2, &d__);
1034 work[i__1].r = q__1.r;
1035 work[i__1].i = q__1.i; // , expr subst
1036 if (j > 1)
1037 {
1038 /* w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling */
1039 i__1 = j;
1040 if ((r__1 = work[i__1].r, f2c_abs(r__1)) + (r__2 = r_imag( &work[j]), f2c_abs(r__2)) > 1.f)
1041 {
1042 i__1 = j;
1043 temp = 1.f / ((r__1 = work[i__1].r, f2c_abs(r__1)) + ( r__2 = r_imag(&work[j]), f2c_abs(r__2)));
1044 if (acoefa * rwork[j] + bcoefa * rwork[*n + j] >= bignum * temp)
1045 {
1046 i__1 = je;
1047 for (jr = 1;
1048 jr <= i__1;
1049 ++jr)
1050 {
1051 i__2 = jr;
1052 i__3 = jr;
1053 q__1.r = temp * work[i__3].r;
1054 q__1.i = temp * work[i__3].i; // , expr subst
1055 work[i__2].r = q__1.r;
1056 work[i__2].i = q__1.i; // , expr subst
1057 /* L190: */
1058 }
1059 }
1060 }
1061 i__1 = j;
1062 q__1.r = acoeff * work[i__1].r;
1063 q__1.i = acoeff * work[i__1].i; // , expr subst
1064 ca.r = q__1.r;
1065 ca.i = q__1.i; // , expr subst
1066 i__1 = j;
1067 q__1.r = bcoeff.r * work[i__1].r - bcoeff.i * work[ i__1].i;
1068 q__1.i = bcoeff.r * work[i__1].i + bcoeff.i * work[i__1].r; // , expr subst
1069 cb.r = q__1.r;
1070 cb.i = q__1.i; // , expr subst
1071 i__1 = j - 1;
1072 for (jr = 1;
1073 jr <= i__1;
1074 ++jr)
1075 {
1076 i__2 = jr;
1077 i__3 = jr;
1078 i__4 = jr + j * s_dim1;
1079 q__3.r = ca.r * s[i__4].r - ca.i * s[i__4].i;
1080 q__3.i = ca.r * s[i__4].i + ca.i * s[i__4] .r; // , expr subst
1081 q__2.r = work[i__3].r + q__3.r;
1082 q__2.i = work[ i__3].i + q__3.i; // , expr subst
1083 i__5 = jr + j * p_dim1;
1084 q__4.r = cb.r * p[i__5].r - cb.i * p[i__5].i;
1085 q__4.i = cb.r * p[i__5].i + cb.i * p[i__5] .r; // , expr subst
1086 q__1.r = q__2.r - q__4.r;
1087 q__1.i = q__2.i - q__4.i; // , expr subst
1088 work[i__2].r = q__1.r;
1089 work[i__2].i = q__1.i; // , expr subst
1090 /* L200: */
1091 }
1092 }
1093 /* L210: */
1094 }
1095 /* Back transform eigenvector if HOWMNY='B'. */
1096 if (ilback)
1097 {
1098 cgemv_("N", n, &je, &c_b2, &vr[vr_offset], ldvr, &work[1], &c__1, &c_b1, &work[*n + 1], &c__1);
1099 isrc = 2;
1100 iend = *n;
1101 }
1102 else
1103 {
1104 isrc = 1;
1105 iend = je;
1106 }
1107 /* Copy and scale eigenvector into column of VR */
1108 xmax = 0.f;
1109 i__1 = iend;
1110 for (jr = 1;
1111 jr <= i__1;
1112 ++jr)
1113 {
1114 /* Computing MAX */
1115 i__2 = (isrc - 1) * *n + jr;
1116 r__3 = xmax;
1117 r__4 = (r__1 = work[i__2].r, f2c_abs(r__1)) + ( r__2 = r_imag(&work[(isrc - 1) * *n + jr]), f2c_abs( r__2)); // , expr subst
1118 xmax = max(r__3,r__4);
1119 /* L220: */
1120 }
1121 if (xmax > safmin)
1122 {
1123 temp = 1.f / xmax;
1124 i__1 = iend;
1125 for (jr = 1;
1126 jr <= i__1;
1127 ++jr)
1128 {
1129 i__2 = jr + ieig * vr_dim1;
1130 i__3 = (isrc - 1) * *n + jr;
1131 q__1.r = temp * work[i__3].r;
1132 q__1.i = temp * work[ i__3].i; // , expr subst
1133 vr[i__2].r = q__1.r;
1134 vr[i__2].i = q__1.i; // , expr subst
1135 /* L230: */
1136 }
1137 }
1138 else
1139 {
1140 iend = 0;
1141 }
1142 i__1 = *n;
1143 for (jr = iend + 1;
1144 jr <= i__1;
1145 ++jr)
1146 {
1147 i__2 = jr + ieig * vr_dim1;
1148 vr[i__2].r = 0.f;
1149 vr[i__2].i = 0.f; // , expr subst
1150 /* L240: */
1151 }
1152 }
1153 L250:
1154 ;
1155 }
1156 }
1157 return 0;
1158 /* End of CTGEVC */
1159 }
1160 /* ctgevc_ */
1161