1 /* ./src_f77/cggesx.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 static integer c__0 = 0;
14 static integer c_n1 = -1;
15
cggesx_(char * jobvsl,char * jobvsr,char * sort,L_fp selctg,char * sense,integer * n,complex * a,integer * lda,complex * b,integer * ldb,integer * sdim,complex * alpha,complex * beta,complex * vsl,integer * ldvsl,complex * vsr,integer * ldvsr,real * rconde,real * rcondv,complex * work,integer * lwork,real * rwork,integer * iwork,integer * liwork,logical * bwork,integer * info,ftnlen jobvsl_len,ftnlen jobvsr_len,ftnlen sort_len,ftnlen sense_len)16 /* Subroutine */ int cggesx_(char *jobvsl, char *jobvsr, char *sort, L_fp
17 selctg, char *sense, integer *n, complex *a, integer *lda, complex *b,
18 integer *ldb, integer *sdim, complex *alpha, complex *beta, complex *
19 vsl, integer *ldvsl, complex *vsr, integer *ldvsr, real *rconde, real
20 *rcondv, complex *work, integer *lwork, real *rwork, integer *iwork,
21 integer *liwork, logical *bwork, integer *info, ftnlen jobvsl_len,
22 ftnlen jobvsr_len, ftnlen sort_len, ftnlen sense_len)
23 {
24 /* System generated locals */
25 integer a_dim1, a_offset, b_dim1, b_offset, vsl_dim1, vsl_offset,
26 vsr_dim1, vsr_offset, i__1, i__2;
27
28 /* Builtin functions */
29 double sqrt(doublereal);
30
31 /* Local variables */
32 static integer i__;
33 static real pl, pr, dif[2];
34 static integer ihi, ilo;
35 static real eps;
36 static integer ijob;
37 static real anrm, bnrm;
38 static integer ierr, itau, iwrk;
39 extern logical lsame_(char *, char *, ftnlen, ftnlen);
40 static integer ileft, icols;
41 static logical cursl, ilvsl, ilvsr;
42 static integer irwrk, irows;
43 extern /* Subroutine */ int cggbak_(char *, char *, integer *, integer *,
44 integer *, real *, real *, integer *, complex *, integer *,
45 integer *, ftnlen, ftnlen), cggbal_(char *, integer *, complex *,
46 integer *, complex *, integer *, integer *, integer *, real *,
47 real *, real *, integer *, ftnlen), slabad_(real *, real *);
48 extern doublereal clange_(char *, integer *, integer *, complex *,
49 integer *, real *, ftnlen);
50 extern /* Subroutine */ int cgghrd_(char *, char *, integer *, integer *,
51 integer *, complex *, integer *, complex *, integer *, complex *,
52 integer *, complex *, integer *, integer *, ftnlen, ftnlen),
53 clascl_(char *, integer *, integer *, real *, real *, integer *,
54 integer *, complex *, integer *, integer *, ftnlen);
55 static logical ilascl, ilbscl;
56 extern /* Subroutine */ int cgeqrf_(integer *, integer *, complex *,
57 integer *, complex *, complex *, integer *, integer *), clacpy_(
58 char *, integer *, integer *, complex *, integer *, complex *,
59 integer *, ftnlen), claset_(char *, integer *, integer *, complex
60 *, complex *, complex *, integer *, ftnlen), xerbla_(char *,
61 integer *, ftnlen);
62 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
63 integer *, integer *, ftnlen, ftnlen);
64 extern doublereal slamch_(char *, ftnlen);
65 static real bignum;
66 extern /* Subroutine */ int chgeqz_(char *, char *, char *, integer *,
67 integer *, integer *, complex *, integer *, complex *, integer *,
68 complex *, complex *, complex *, integer *, complex *, integer *,
69 complex *, integer *, real *, integer *, ftnlen, ftnlen, ftnlen),
70 ctgsen_(integer *, logical *, logical *, logical *, integer *,
71 complex *, integer *, complex *, integer *, complex *, complex *,
72 complex *, integer *, complex *, integer *, integer *, real *,
73 real *, real *, complex *, integer *, integer *, integer *,
74 integer *);
75 static integer ijobvl, iright, ijobvr;
76 static logical wantsb;
77 static integer liwmin;
78 static logical wantse, lastsl;
79 static real anrmto, bnrmto;
80 extern /* Subroutine */ int cungqr_(integer *, integer *, integer *,
81 complex *, integer *, complex *, complex *, integer *, integer *);
82 static integer minwrk, maxwrk;
83 static logical wantsn;
84 static real smlnum;
85 extern /* Subroutine */ int cunmqr_(char *, char *, integer *, integer *,
86 integer *, complex *, integer *, complex *, complex *, integer *,
87 complex *, integer *, integer *, ftnlen, ftnlen);
88 static logical wantst, wantsv;
89
90
91 /* -- LAPACK driver routine (version 3.0) -- */
92 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
93 /* Courant Institute, Argonne National Lab, and Rice University */
94 /* June 30, 1999 */
95
96 /* .. Scalar Arguments .. */
97 /* .. */
98 /* .. Array Arguments .. */
99 /* .. */
100 /* .. Function Arguments .. */
101 /* .. */
102
103 /* Purpose */
104 /* ======= */
105
106 /* CGGESX computes for a pair of N-by-N complex nonsymmetric matrices */
107 /* (A,B), the generalized eigenvalues, the complex Schur form (S,T), */
108 /* and, optionally, the left and/or right matrices of Schur vectors (VSL */
109 /* and VSR). This gives the generalized Schur factorization */
110
111 /* (A,B) = ( (VSL) S (VSR)**H, (VSL) T (VSR)**H ) */
112
113 /* where (VSR)**H is the conjugate-transpose of VSR. */
114
115 /* Optionally, it also orders the eigenvalues so that a selected cluster */
116 /* of eigenvalues appears in the leading diagonal blocks of the upper */
117 /* triangular matrix S and the upper triangular matrix T; computes */
118 /* a reciprocal condition number for the average of the selected */
119 /* eigenvalues (RCONDE); and computes a reciprocal condition number for */
120 /* the right and left deflating subspaces corresponding to the selected */
121 /* eigenvalues (RCONDV). The leading columns of VSL and VSR then form */
122 /* an orthonormal basis for the corresponding left and right eigenspaces */
123 /* (deflating subspaces). */
124
125 /* A generalized eigenvalue for a pair of matrices (A,B) is a scalar w */
126 /* or a ratio alpha/beta = w, such that A - w*B is singular. It is */
127 /* usually represented as the pair (alpha,beta), as there is a */
128 /* reasonable interpretation for beta=0 or for both being zero. */
129
130 /* A pair of matrices (S,T) is in generalized complex Schur form if T is */
131 /* upper triangular with non-negative diagonal and S is upper */
132 /* triangular. */
133
134 /* Arguments */
135 /* ========= */
136
137 /* JOBVSL (input) CHARACTER*1 */
138 /* = 'N': do not compute the left Schur vectors; */
139 /* = 'V': compute the left Schur vectors. */
140
141 /* JOBVSR (input) CHARACTER*1 */
142 /* = 'N': do not compute the right Schur vectors; */
143 /* = 'V': compute the right Schur vectors. */
144
145 /* SORT (input) CHARACTER*1 */
146 /* Specifies whether or not to order the eigenvalues on the */
147 /* diagonal of the generalized Schur form. */
148 /* = 'N': Eigenvalues are not ordered; */
149 /* = 'S': Eigenvalues are ordered (see SELCTG). */
150
151 /* SELCTG (input) LOGICAL FUNCTION of two COMPLEX arguments */
152 /* SELCTG must be declared EXTERNAL in the calling subroutine. */
153 /* If SORT = 'N', SELCTG is not referenced. */
154 /* If SORT = 'S', SELCTG is used to select eigenvalues to sort */
155 /* to the top left of the Schur form. */
156 /* Note that a selected complex eigenvalue may no longer satisfy */
157 /* SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since */
158 /* ordering may change the value of complex eigenvalues */
159 /* (especially if the eigenvalue is ill-conditioned), in this */
160 /* case INFO is set to N+3 see INFO below). */
161
162 /* SENSE (input) CHARACTER */
163 /* Determines which reciprocal condition numbers are computed. */
164 /* = 'N' : None are computed; */
165 /* = 'E' : Computed for average of selected eigenvalues only; */
166 /* = 'V' : Computed for selected deflating subspaces only; */
167 /* = 'B' : Computed for both. */
168 /* If SENSE = 'E', 'V', or 'B', SORT must equal 'S'. */
169
170 /* N (input) INTEGER */
171 /* The order of the matrices A, B, VSL, and VSR. N >= 0. */
172
173 /* A (input/output) COMPLEX array, dimension (LDA, N) */
174 /* On entry, the first of the pair of matrices. */
175 /* On exit, A has been overwritten by its generalized Schur */
176 /* form S. */
177
178 /* LDA (input) INTEGER */
179 /* The leading dimension of A. LDA >= max(1,N). */
180
181 /* B (input/output) COMPLEX array, dimension (LDB, N) */
182 /* On entry, the second of the pair of matrices. */
183 /* On exit, B has been overwritten by its generalized Schur */
184 /* form T. */
185
186 /* LDB (input) INTEGER */
187 /* The leading dimension of B. LDB >= max(1,N). */
188
189 /* SDIM (output) INTEGER */
190 /* If SORT = 'N', SDIM = 0. */
191 /* If SORT = 'S', SDIM = number of eigenvalues (after sorting) */
192 /* for which SELCTG is true. */
193
194 /* ALPHA (output) COMPLEX array, dimension (N) */
195 /* BETA (output) COMPLEX array, dimension (N) */
196 /* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the */
197 /* generalized eigenvalues. ALPHA(j) and BETA(j),j=1,...,N are */
198 /* the diagonals of the complex Schur form (S,T). BETA(j) will */
199 /* be non-negative real. */
200
201 /* Note: the quotients ALPHA(j)/BETA(j) may easily over- or */
202 /* underflow, and BETA(j) may even be zero. Thus, the user */
203 /* should avoid naively computing the ratio alpha/beta. */
204 /* However, ALPHA will be always less than and usually */
205 /* comparable with norm(A) in magnitude, and BETA always less */
206 /* than and usually comparable with norm(B). */
207
208 /* VSL (output) COMPLEX array, dimension (LDVSL,N) */
209 /* If JOBVSL = 'V', VSL will contain the left Schur vectors. */
210 /* Not referenced if JOBVSL = 'N'. */
211
212 /* LDVSL (input) INTEGER */
213 /* The leading dimension of the matrix VSL. LDVSL >=1, and */
214 /* if JOBVSL = 'V', LDVSL >= N. */
215
216 /* VSR (output) COMPLEX array, dimension (LDVSR,N) */
217 /* If JOBVSR = 'V', VSR will contain the right Schur vectors. */
218 /* Not referenced if JOBVSR = 'N'. */
219
220 /* LDVSR (input) INTEGER */
221 /* The leading dimension of the matrix VSR. LDVSR >= 1, and */
222 /* if JOBVSR = 'V', LDVSR >= N. */
223
224 /* RCONDE (output) REAL array, dimension ( 2 ) */
225 /* If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the */
226 /* reciprocal condition numbers for the average of the selected */
227 /* eigenvalues. */
228 /* Not referenced if SENSE = 'N' or 'V'. */
229
230 /* RCONDV (output) REAL array, dimension ( 2 ) */
231 /* If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the */
232 /* reciprocal condition number for the selected deflating */
233 /* subspaces. */
234 /* Not referenced if SENSE = 'N' or 'E'. */
235
236 /* WORK (workspace/output) COMPLEX array, dimension (LWORK) */
237 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
238
239 /* LWORK (input) INTEGER */
240 /* The dimension of the array WORK. LWORK >= 2*N. */
241 /* If SENSE = 'E', 'V', or 'B', */
242 /* LWORK >= MAX(2*N, 2*SDIM*(N-SDIM)). */
243
244 /* RWORK (workspace) REAL array, dimension ( 8*N ) */
245 /* Real workspace. */
246
247 /* IWORK (workspace/output) INTEGER array, dimension (LIWORK) */
248 /* Not referenced if SENSE = 'N'. */
249 /* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */
250
251 /* LIWORK (input) INTEGER */
252 /* The dimension of the array WORK. LIWORK >= N+2. */
253
254 /* BWORK (workspace) LOGICAL array, dimension (N) */
255 /* Not referenced if SORT = 'N'. */
256
257 /* INFO (output) INTEGER */
258 /* = 0: successful exit */
259 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
260 /* = 1,...,N: */
261 /* The QZ iteration failed. (A,B) are not in Schur */
262 /* form, but ALPHA(j) and BETA(j) should be correct for */
263 /* j=INFO+1,...,N. */
264 /* > N: =N+1: other than QZ iteration failed in CHGEQZ */
265 /* =N+2: after reordering, roundoff changed values of */
266 /* some complex eigenvalues so that leading */
267 /* eigenvalues in the Generalized Schur form no */
268 /* longer satisfy SELCTG=.TRUE. This could also */
269 /* be caused due to scaling. */
270 /* =N+3: reordering failed in CTGSEN. */
271
272 /* ===================================================================== */
273
274 /* .. Parameters .. */
275 /* .. */
276 /* .. Local Scalars .. */
277 /* .. */
278 /* .. Local Arrays .. */
279 /* .. */
280 /* .. External Subroutines .. */
281 /* .. */
282 /* .. External Functions .. */
283 /* .. */
284 /* .. Intrinsic Functions .. */
285 /* .. */
286 /* .. Executable Statements .. */
287
288 /* Decode the input arguments */
289
290 /* Parameter adjustments */
291 a_dim1 = *lda;
292 a_offset = 1 + a_dim1;
293 a -= a_offset;
294 b_dim1 = *ldb;
295 b_offset = 1 + b_dim1;
296 b -= b_offset;
297 --alpha;
298 --beta;
299 vsl_dim1 = *ldvsl;
300 vsl_offset = 1 + vsl_dim1;
301 vsl -= vsl_offset;
302 vsr_dim1 = *ldvsr;
303 vsr_offset = 1 + vsr_dim1;
304 vsr -= vsr_offset;
305 --rconde;
306 --rcondv;
307 --work;
308 --rwork;
309 --iwork;
310 --bwork;
311
312 /* Function Body */
313 if (lsame_(jobvsl, "N", (ftnlen)1, (ftnlen)1)) {
314 ijobvl = 1;
315 ilvsl = FALSE_;
316 } else if (lsame_(jobvsl, "V", (ftnlen)1, (ftnlen)1)) {
317 ijobvl = 2;
318 ilvsl = TRUE_;
319 } else {
320 ijobvl = -1;
321 ilvsl = FALSE_;
322 }
323
324 if (lsame_(jobvsr, "N", (ftnlen)1, (ftnlen)1)) {
325 ijobvr = 1;
326 ilvsr = FALSE_;
327 } else if (lsame_(jobvsr, "V", (ftnlen)1, (ftnlen)1)) {
328 ijobvr = 2;
329 ilvsr = TRUE_;
330 } else {
331 ijobvr = -1;
332 ilvsr = FALSE_;
333 }
334
335 wantst = lsame_(sort, "S", (ftnlen)1, (ftnlen)1);
336 wantsn = lsame_(sense, "N", (ftnlen)1, (ftnlen)1);
337 wantse = lsame_(sense, "E", (ftnlen)1, (ftnlen)1);
338 wantsv = lsame_(sense, "V", (ftnlen)1, (ftnlen)1);
339 wantsb = lsame_(sense, "B", (ftnlen)1, (ftnlen)1);
340 if (wantsn) {
341 ijob = 0;
342 iwork[1] = 1;
343 } else if (wantse) {
344 ijob = 1;
345 } else if (wantsv) {
346 ijob = 2;
347 } else if (wantsb) {
348 ijob = 4;
349 }
350
351 /* Test the input arguments */
352
353 *info = 0;
354 if (ijobvl <= 0) {
355 *info = -1;
356 } else if (ijobvr <= 0) {
357 *info = -2;
358 } else if (! wantst && ! lsame_(sort, "N", (ftnlen)1, (ftnlen)1)) {
359 *info = -3;
360 } else if (! (wantsn || wantse || wantsv || wantsb) || ! wantst && !
361 wantsn) {
362 *info = -5;
363 } else if (*n < 0) {
364 *info = -6;
365 } else if (*lda < max(1,*n)) {
366 *info = -8;
367 } else if (*ldb < max(1,*n)) {
368 *info = -10;
369 } else if (*ldvsl < 1 || ilvsl && *ldvsl < *n) {
370 *info = -15;
371 } else if (*ldvsr < 1 || ilvsr && *ldvsr < *n) {
372 *info = -17;
373 }
374
375 /* Compute workspace */
376 /* (Note: Comments in the code beginning "Workspace:" describe the */
377 /* minimal amount of workspace needed at that point in the code, */
378 /* as well as the preferred amount for good performance. */
379 /* NB refers to the optimal block size for the immediately */
380 /* following subroutine, as returned by ILAENV.) */
381
382 minwrk = 1;
383 if (*info == 0 && *lwork >= 1) {
384 /* Computing MAX */
385 i__1 = 1, i__2 = *n << 1;
386 minwrk = max(i__1,i__2);
387 maxwrk = *n + *n * ilaenv_(&c__1, "CGEQRF", " ", n, &c__1, n, &c__0, (
388 ftnlen)6, (ftnlen)1);
389 if (ilvsl) {
390 /* Computing MAX */
391 i__1 = maxwrk, i__2 = *n + *n * ilaenv_(&c__1, "CUNGQR", " ", n, &
392 c__1, n, &c_n1, (ftnlen)6, (ftnlen)1);
393 maxwrk = max(i__1,i__2);
394 }
395 work[1].r = (real) maxwrk, work[1].i = 0.f;
396 }
397 if (! wantsn) {
398 liwmin = *n + 2;
399 } else {
400 liwmin = 1;
401 }
402 iwork[1] = liwmin;
403
404 if (*info == 0 && *lwork < minwrk) {
405 *info = -21;
406 } else if (*info == 0 && ijob >= 1) {
407 if (*liwork < liwmin) {
408 *info = -24;
409 }
410 }
411
412 if (*info != 0) {
413 i__1 = -(*info);
414 xerbla_("CGGESX", &i__1, (ftnlen)6);
415 return 0;
416 }
417
418 /* Quick return if possible */
419
420 if (*n == 0) {
421 *sdim = 0;
422 return 0;
423 }
424
425 /* Get machine constants */
426
427 eps = slamch_("P", (ftnlen)1);
428 smlnum = slamch_("S", (ftnlen)1);
429 bignum = 1.f / smlnum;
430 slabad_(&smlnum, &bignum);
431 smlnum = sqrt(smlnum) / eps;
432 bignum = 1.f / smlnum;
433
434 /* Scale A if max element outside range [SMLNUM,BIGNUM] */
435
436 anrm = clange_("M", n, n, &a[a_offset], lda, &rwork[1], (ftnlen)1);
437 ilascl = FALSE_;
438 if (anrm > 0.f && anrm < smlnum) {
439 anrmto = smlnum;
440 ilascl = TRUE_;
441 } else if (anrm > bignum) {
442 anrmto = bignum;
443 ilascl = TRUE_;
444 }
445 if (ilascl) {
446 clascl_("G", &c__0, &c__0, &anrm, &anrmto, n, n, &a[a_offset], lda, &
447 ierr, (ftnlen)1);
448 }
449
450 /* Scale B if max element outside range [SMLNUM,BIGNUM] */
451
452 bnrm = clange_("M", n, n, &b[b_offset], ldb, &rwork[1], (ftnlen)1);
453 ilbscl = FALSE_;
454 if (bnrm > 0.f && bnrm < smlnum) {
455 bnrmto = smlnum;
456 ilbscl = TRUE_;
457 } else if (bnrm > bignum) {
458 bnrmto = bignum;
459 ilbscl = TRUE_;
460 }
461 if (ilbscl) {
462 clascl_("G", &c__0, &c__0, &bnrm, &bnrmto, n, n, &b[b_offset], ldb, &
463 ierr, (ftnlen)1);
464 }
465
466 /* Permute the matrix to make it more nearly triangular */
467 /* (Real Workspace: need 6*N) */
468
469 ileft = 1;
470 iright = *n + 1;
471 irwrk = iright + *n;
472 cggbal_("P", n, &a[a_offset], lda, &b[b_offset], ldb, &ilo, &ihi, &rwork[
473 ileft], &rwork[iright], &rwork[irwrk], &ierr, (ftnlen)1);
474
475 /* Reduce B to triangular form (QR decomposition of B) */
476 /* (Complex Workspace: need N, prefer N*NB) */
477
478 irows = ihi + 1 - ilo;
479 icols = *n + 1 - ilo;
480 itau = 1;
481 iwrk = itau + irows;
482 i__1 = *lwork + 1 - iwrk;
483 cgeqrf_(&irows, &icols, &b[ilo + ilo * b_dim1], ldb, &work[itau], &work[
484 iwrk], &i__1, &ierr);
485
486 /* Apply the unitary transformation to matrix A */
487 /* (Complex Workspace: need N, prefer N*NB) */
488
489 i__1 = *lwork + 1 - iwrk;
490 cunmqr_("L", "C", &irows, &icols, &irows, &b[ilo + ilo * b_dim1], ldb, &
491 work[itau], &a[ilo + ilo * a_dim1], lda, &work[iwrk], &i__1, &
492 ierr, (ftnlen)1, (ftnlen)1);
493
494 /* Initialize VSL */
495 /* (Complex Workspace: need N, prefer N*NB) */
496
497 if (ilvsl) {
498 claset_("Full", n, n, &c_b1, &c_b2, &vsl[vsl_offset], ldvsl, (ftnlen)
499 4);
500 i__1 = irows - 1;
501 i__2 = irows - 1;
502 clacpy_("L", &i__1, &i__2, &b[ilo + 1 + ilo * b_dim1], ldb, &vsl[ilo
503 + 1 + ilo * vsl_dim1], ldvsl, (ftnlen)1);
504 i__1 = *lwork + 1 - iwrk;
505 cungqr_(&irows, &irows, &irows, &vsl[ilo + ilo * vsl_dim1], ldvsl, &
506 work[itau], &work[iwrk], &i__1, &ierr);
507 }
508
509 /* Initialize VSR */
510
511 if (ilvsr) {
512 claset_("Full", n, n, &c_b1, &c_b2, &vsr[vsr_offset], ldvsr, (ftnlen)
513 4);
514 }
515
516 /* Reduce to generalized Hessenberg form */
517 /* (Workspace: none needed) */
518
519 cgghrd_(jobvsl, jobvsr, n, &ilo, &ihi, &a[a_offset], lda, &b[b_offset],
520 ldb, &vsl[vsl_offset], ldvsl, &vsr[vsr_offset], ldvsr, &ierr, (
521 ftnlen)1, (ftnlen)1);
522
523 *sdim = 0;
524
525 /* Perform QZ algorithm, computing Schur vectors if desired */
526 /* (Complex Workspace: need N) */
527 /* (Real Workspace: need N) */
528
529 iwrk = itau;
530 i__1 = *lwork + 1 - iwrk;
531 chgeqz_("S", jobvsl, jobvsr, n, &ilo, &ihi, &a[a_offset], lda, &b[
532 b_offset], ldb, &alpha[1], &beta[1], &vsl[vsl_offset], ldvsl, &
533 vsr[vsr_offset], ldvsr, &work[iwrk], &i__1, &rwork[irwrk], &ierr,
534 (ftnlen)1, (ftnlen)1, (ftnlen)1);
535 if (ierr != 0) {
536 if (ierr > 0 && ierr <= *n) {
537 *info = ierr;
538 } else if (ierr > *n && ierr <= *n << 1) {
539 *info = ierr - *n;
540 } else {
541 *info = *n + 1;
542 }
543 goto L40;
544 }
545
546 /* Sort eigenvalues ALPHA/BETA and compute the reciprocal of */
547 /* condition number(s) */
548
549 if (wantst) {
550
551 /* Undo scaling on eigenvalues before SELCTGing */
552
553 if (ilascl) {
554 clascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alpha[1], n,
555 &ierr, (ftnlen)1);
556 }
557 if (ilbscl) {
558 clascl_("G", &c__0, &c__0, &bnrmto, &bnrm, n, &c__1, &beta[1], n,
559 &ierr, (ftnlen)1);
560 }
561
562 /* Select eigenvalues */
563
564 i__1 = *n;
565 for (i__ = 1; i__ <= i__1; ++i__) {
566 bwork[i__] = (*selctg)(&alpha[i__], &beta[i__]);
567 /* L10: */
568 }
569
570 /* Reorder eigenvalues, transform Generalized Schur vectors, and */
571 /* compute reciprocal condition numbers */
572 /* (Complex Workspace: If IJOB >= 1, need MAX(1, 2*SDIM*(N-SDIM)) */
573 /* otherwise, need 1 ) */
574
575 i__1 = *lwork - iwrk + 1;
576 ctgsen_(&ijob, &ilvsl, &ilvsr, &bwork[1], n, &a[a_offset], lda, &b[
577 b_offset], ldb, &alpha[1], &beta[1], &vsl[vsl_offset], ldvsl,
578 &vsr[vsr_offset], ldvsr, sdim, &pl, &pr, dif, &work[iwrk], &
579 i__1, &iwork[1], liwork, &ierr);
580
581 if (ijob >= 1) {
582 /* Computing MAX */
583 i__1 = maxwrk, i__2 = (*sdim << 1) * (*n - *sdim);
584 maxwrk = max(i__1,i__2);
585 }
586 if (ierr == -21) {
587
588 /* not enough complex workspace */
589
590 *info = -21;
591 } else {
592 rconde[1] = pl;
593 rconde[2] = pl;
594 rcondv[1] = dif[0];
595 rcondv[2] = dif[1];
596 if (ierr == 1) {
597 *info = *n + 3;
598 }
599 }
600
601 }
602
603 /* Apply permutation to VSL and VSR */
604 /* (Workspace: none needed) */
605
606 if (ilvsl) {
607 cggbak_("P", "L", n, &ilo, &ihi, &rwork[ileft], &rwork[iright], n, &
608 vsl[vsl_offset], ldvsl, &ierr, (ftnlen)1, (ftnlen)1);
609 }
610
611 if (ilvsr) {
612 cggbak_("P", "R", n, &ilo, &ihi, &rwork[ileft], &rwork[iright], n, &
613 vsr[vsr_offset], ldvsr, &ierr, (ftnlen)1, (ftnlen)1);
614 }
615
616 /* Undo scaling */
617
618 if (ilascl) {
619 clascl_("U", &c__0, &c__0, &anrmto, &anrm, n, n, &a[a_offset], lda, &
620 ierr, (ftnlen)1);
621 clascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alpha[1], n, &
622 ierr, (ftnlen)1);
623 }
624
625 if (ilbscl) {
626 clascl_("U", &c__0, &c__0, &bnrmto, &bnrm, n, n, &b[b_offset], ldb, &
627 ierr, (ftnlen)1);
628 clascl_("G", &c__0, &c__0, &bnrmto, &bnrm, n, &c__1, &beta[1], n, &
629 ierr, (ftnlen)1);
630 }
631
632 /* L20: */
633
634 if (wantst) {
635
636 /* Check if reordering is correct */
637
638 lastsl = TRUE_;
639 *sdim = 0;
640 i__1 = *n;
641 for (i__ = 1; i__ <= i__1; ++i__) {
642 cursl = (*selctg)(&alpha[i__], &beta[i__]);
643 if (cursl) {
644 ++(*sdim);
645 }
646 if (cursl && ! lastsl) {
647 *info = *n + 2;
648 }
649 lastsl = cursl;
650 /* L30: */
651 }
652
653 }
654
655 L40:
656
657 work[1].r = (real) maxwrk, work[1].i = 0.f;
658 iwork[1] = liwmin;
659
660 return 0;
661
662 /* End of CGGESX */
663
664 } /* cggesx_ */
665
666