1 SUBROUTINE ZGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
2 $ B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR,
3 $ LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK,
4 $ IWORK, LIWORK, BWORK, INFO )
5 *
6 * -- LAPACK driver routine (version 3.2) --
7 * -- LAPACK is a software package provided by Univ. of Tennessee, --
8 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9 * November 2006
10 *
11 * .. Scalar Arguments ..
12 CHARACTER JOBVSL, JOBVSR, SENSE, SORT
13 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
14 $ SDIM
15 * ..
16 * .. Array Arguments ..
17 LOGICAL BWORK( * )
18 INTEGER IWORK( * )
19 DOUBLE PRECISION RCONDE( 2 ), RCONDV( 2 ), RWORK( * )
20 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
21 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
22 $ WORK( * )
23 * ..
24 * .. Function Arguments ..
25 LOGICAL SELCTG
26 EXTERNAL SELCTG
27 * ..
28 *
29 * Purpose
30 * =======
31 *
32 * ZGGESX computes for a pair of N-by-N complex nonsymmetric matrices
33 * (A,B), the generalized eigenvalues, the complex Schur form (S,T),
34 * and, optionally, the left and/or right matrices of Schur vectors (VSL
35 * and VSR). This gives the generalized Schur factorization
36 *
37 * (A,B) = ( (VSL) S (VSR)**H, (VSL) T (VSR)**H )
38 *
39 * where (VSR)**H is the conjugate-transpose of VSR.
40 *
41 * Optionally, it also orders the eigenvalues so that a selected cluster
42 * of eigenvalues appears in the leading diagonal blocks of the upper
43 * triangular matrix S and the upper triangular matrix T; computes
44 * a reciprocal condition number for the average of the selected
45 * eigenvalues (RCONDE); and computes a reciprocal condition number for
46 * the right and left deflating subspaces corresponding to the selected
47 * eigenvalues (RCONDV). The leading columns of VSL and VSR then form
48 * an orthonormal basis for the corresponding left and right eigenspaces
49 * (deflating subspaces).
50 *
51 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
52 * or a ratio alpha/beta = w, such that A - w*B is singular. It is
53 * usually represented as the pair (alpha,beta), as there is a
54 * reasonable interpretation for beta=0 or for both being zero.
55 *
56 * A pair of matrices (S,T) is in generalized complex Schur form if T is
57 * upper triangular with non-negative diagonal and S is upper
58 * triangular.
59 *
60 * Arguments
61 * =========
62 *
63 * JOBVSL (input) CHARACTER*1
64 * = 'N': do not compute the left Schur vectors;
65 * = 'V': compute the left Schur vectors.
66 *
67 * JOBVSR (input) CHARACTER*1
68 * = 'N': do not compute the right Schur vectors;
69 * = 'V': compute the right Schur vectors.
70 *
71 * SORT (input) CHARACTER*1
72 * Specifies whether or not to order the eigenvalues on the
73 * diagonal of the generalized Schur form.
74 * = 'N': Eigenvalues are not ordered;
75 * = 'S': Eigenvalues are ordered (see SELCTG).
76 *
77 * SELCTG (external procedure) LOGICAL FUNCTION of two COMPLEX*16 arguments
78 * SELCTG must be declared EXTERNAL in the calling subroutine.
79 * If SORT = 'N', SELCTG is not referenced.
80 * If SORT = 'S', SELCTG is used to select eigenvalues to sort
81 * to the top left of the Schur form.
82 * Note that a selected complex eigenvalue may no longer satisfy
83 * SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
84 * ordering may change the value of complex eigenvalues
85 * (especially if the eigenvalue is ill-conditioned), in this
86 * case INFO is set to N+3 see INFO below).
87 *
88 * SENSE (input) CHARACTER*1
89 * Determines which reciprocal condition numbers are computed.
90 * = 'N' : None are computed;
91 * = 'E' : Computed for average of selected eigenvalues only;
92 * = 'V' : Computed for selected deflating subspaces only;
93 * = 'B' : Computed for both.
94 * If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
95 *
96 * N (input) INTEGER
97 * The order of the matrices A, B, VSL, and VSR. N >= 0.
98 *
99 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
100 * On entry, the first of the pair of matrices.
101 * On exit, A has been overwritten by its generalized Schur
102 * form S.
103 *
104 * LDA (input) INTEGER
105 * The leading dimension of A. LDA >= max(1,N).
106 *
107 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
108 * On entry, the second of the pair of matrices.
109 * On exit, B has been overwritten by its generalized Schur
110 * form T.
111 *
112 * LDB (input) INTEGER
113 * The leading dimension of B. LDB >= max(1,N).
114 *
115 * SDIM (output) INTEGER
116 * If SORT = 'N', SDIM = 0.
117 * If SORT = 'S', SDIM = number of eigenvalues (after sorting)
118 * for which SELCTG is true.
119 *
120 * ALPHA (output) COMPLEX*16 array, dimension (N)
121 * BETA (output) COMPLEX*16 array, dimension (N)
122 * On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
123 * generalized eigenvalues. ALPHA(j) and BETA(j),j=1,...,N are
124 * the diagonals of the complex Schur form (S,T). BETA(j) will
125 * be non-negative real.
126 *
127 * Note: the quotients ALPHA(j)/BETA(j) may easily over- or
128 * underflow, and BETA(j) may even be zero. Thus, the user
129 * should avoid naively computing the ratio alpha/beta.
130 * However, ALPHA will be always less than and usually
131 * comparable with norm(A) in magnitude, and BETA always less
132 * than and usually comparable with norm(B).
133 *
134 * VSL (output) COMPLEX*16 array, dimension (LDVSL,N)
135 * If JOBVSL = 'V', VSL will contain the left Schur vectors.
136 * Not referenced if JOBVSL = 'N'.
137 *
138 * LDVSL (input) INTEGER
139 * The leading dimension of the matrix VSL. LDVSL >=1, and
140 * if JOBVSL = 'V', LDVSL >= N.
141 *
142 * VSR (output) COMPLEX*16 array, dimension (LDVSR,N)
143 * If JOBVSR = 'V', VSR will contain the right Schur vectors.
144 * Not referenced if JOBVSR = 'N'.
145 *
146 * LDVSR (input) INTEGER
147 * The leading dimension of the matrix VSR. LDVSR >= 1, and
148 * if JOBVSR = 'V', LDVSR >= N.
149 *
150 * RCONDE (output) DOUBLE PRECISION array, dimension ( 2 )
151 * If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
152 * reciprocal condition numbers for the average of the selected
153 * eigenvalues.
154 * Not referenced if SENSE = 'N' or 'V'.
155 *
156 * RCONDV (output) DOUBLE PRECISION array, dimension ( 2 )
157 * If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
158 * reciprocal condition number for the selected deflating
159 * subspaces.
160 * Not referenced if SENSE = 'N' or 'E'.
161 *
162 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
163 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
164 *
165 * LWORK (input) INTEGER
166 * The dimension of the array WORK.
167 * If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
168 * LWORK >= MAX(1,2*N,2*SDIM*(N-SDIM)), else
169 * LWORK >= MAX(1,2*N). Note that 2*SDIM*(N-SDIM) <= N*N/2.
170 * Note also that an error is only returned if
171 * LWORK < MAX(1,2*N), but if SENSE = 'E' or 'V' or 'B' this may
172 * not be large enough.
173 *
174 * If LWORK = -1, then a workspace query is assumed; the routine
175 * only calculates the bound on the optimal size of the WORK
176 * array and the minimum size of the IWORK array, returns these
177 * values as the first entries of the WORK and IWORK arrays, and
178 * no error message related to LWORK or LIWORK is issued by
179 * XERBLA.
180 *
181 * RWORK (workspace) DOUBLE PRECISION array, dimension ( 8*N )
182 * Real workspace.
183 *
184 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
185 * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
186 *
187 * LIWORK (input) INTEGER
188 * The dimension of the array IWORK.
189 * If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
190 * LIWORK >= N+2.
191 *
192 * If LIWORK = -1, then a workspace query is assumed; the
193 * routine only calculates the bound on the optimal size of the
194 * WORK array and the minimum size of the IWORK array, returns
195 * these values as the first entries of the WORK and IWORK
196 * arrays, and no error message related to LWORK or LIWORK is
197 * issued by XERBLA.
198 *
199 * BWORK (workspace) LOGICAL array, dimension (N)
200 * Not referenced if SORT = 'N'.
201 *
202 * INFO (output) INTEGER
203 * = 0: successful exit
204 * < 0: if INFO = -i, the i-th argument had an illegal value.
205 * = 1,...,N:
206 * The QZ iteration failed. (A,B) are not in Schur
207 * form, but ALPHA(j) and BETA(j) should be correct for
208 * j=INFO+1,...,N.
209 * > N: =N+1: other than QZ iteration failed in ZHGEQZ
210 * =N+2: after reordering, roundoff changed values of
211 * some complex eigenvalues so that leading
212 * eigenvalues in the Generalized Schur form no
213 * longer satisfy SELCTG=.TRUE. This could also
214 * be caused due to scaling.
215 * =N+3: reordering failed in ZTGSEN.
216 *
217 * =====================================================================
218 *
219 * .. Parameters ..
220 DOUBLE PRECISION ZERO, ONE
221 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
222 COMPLEX*16 CZERO, CONE
223 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
224 $ CONE = ( 1.0D+0, 0.0D+0 ) )
225 * ..
226 * .. Local Scalars ..
227 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
228 $ LQUERY, WANTSB, WANTSE, WANTSN, WANTST, WANTSV
229 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
230 $ ILEFT, ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK,
231 $ LIWMIN, LWRK, MAXWRK, MINWRK
232 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
233 $ PR, SMLNUM
234 * ..
235 * .. Local Arrays ..
236 DOUBLE PRECISION DIF( 2 )
237 * ..
238 * .. External Subroutines ..
239 EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
240 $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR,
241 $ ZUNMQR
242 * ..
243 * .. External Functions ..
244 LOGICAL LSAME
245 INTEGER ILAENV
246 DOUBLE PRECISION DLAMCH, ZLANGE
247 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
248 * ..
249 * .. Intrinsic Functions ..
250 INTRINSIC MAX, SQRT
251 * ..
252 * .. Executable Statements ..
253 *
254 * Decode the input arguments
255 *
256 IF( LSAME( JOBVSL, 'N' ) ) THEN
257 IJOBVL = 1
258 ILVSL = .FALSE.
259 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
260 IJOBVL = 2
261 ILVSL = .TRUE.
262 ELSE
263 IJOBVL = -1
264 ILVSL = .FALSE.
265 END IF
266 *
267 IF( LSAME( JOBVSR, 'N' ) ) THEN
268 IJOBVR = 1
269 ILVSR = .FALSE.
270 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
271 IJOBVR = 2
272 ILVSR = .TRUE.
273 ELSE
274 IJOBVR = -1
275 ILVSR = .FALSE.
276 END IF
277 *
278 WANTST = LSAME( SORT, 'S' )
279 WANTSN = LSAME( SENSE, 'N' )
280 WANTSE = LSAME( SENSE, 'E' )
281 WANTSV = LSAME( SENSE, 'V' )
282 WANTSB = LSAME( SENSE, 'B' )
283 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
284 IF( WANTSN ) THEN
285 IJOB = 0
286 ELSE IF( WANTSE ) THEN
287 IJOB = 1
288 ELSE IF( WANTSV ) THEN
289 IJOB = 2
290 ELSE IF( WANTSB ) THEN
291 IJOB = 4
292 END IF
293 *
294 * Test the input arguments
295 *
296 INFO = 0
297 IF( IJOBVL.LE.0 ) THEN
298 INFO = -1
299 ELSE IF( IJOBVR.LE.0 ) THEN
300 INFO = -2
301 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
302 INFO = -3
303 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
304 $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
305 INFO = -5
306 ELSE IF( N.LT.0 ) THEN
307 INFO = -6
308 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
309 INFO = -8
310 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
311 INFO = -10
312 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
313 INFO = -15
314 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
315 INFO = -17
316 END IF
317 *
318 * Compute workspace
319 * (Note: Comments in the code beginning "Workspace:" describe the
320 * minimal amount of workspace needed at that point in the code,
321 * as well as the preferred amount for good performance.
322 * NB refers to the optimal block size for the immediately
323 * following subroutine, as returned by ILAENV.)
324 *
325 IF( INFO.EQ.0 ) THEN
326 IF( N.GT.0) THEN
327 MINWRK = 2*N
328 MAXWRK = N*(1 + ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
329 MAXWRK = MAX( MAXWRK, N*( 1 +
330 $ ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, -1 ) ) )
331 IF( ILVSL ) THEN
332 MAXWRK = MAX( MAXWRK, N*( 1 +
333 $ ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) ) )
334 END IF
335 LWRK = MAXWRK
336 IF( IJOB.GE.1 )
337 $ LWRK = MAX( LWRK, N*N/2 )
338 ELSE
339 MINWRK = 1
340 MAXWRK = 1
341 LWRK = 1
342 END IF
343 WORK( 1 ) = LWRK
344 IF( WANTSN .OR. N.EQ.0 ) THEN
345 LIWMIN = 1
346 ELSE
347 LIWMIN = N + 2
348 END IF
349 IWORK( 1 ) = LIWMIN
350 *
351 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
352 INFO = -21
353 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY) THEN
354 INFO = -24
355 END IF
356 END IF
357 *
358 IF( INFO.NE.0 ) THEN
359 CALL XERBLA( 'ZGGESX', -INFO )
360 RETURN
361 ELSE IF (LQUERY) THEN
362 RETURN
363 END IF
364 *
365 * Quick return if possible
366 *
367 IF( N.EQ.0 ) THEN
368 SDIM = 0
369 RETURN
370 END IF
371 *
372 * Get machine constants
373 *
374 EPS = DLAMCH( 'P' )
375 SMLNUM = DLAMCH( 'S' )
376 BIGNUM = ONE / SMLNUM
377 CALL DLABAD( SMLNUM, BIGNUM )
378 SMLNUM = SQRT( SMLNUM ) / EPS
379 BIGNUM = ONE / SMLNUM
380 *
381 * Scale A if max element outside range [SMLNUM,BIGNUM]
382 *
383 ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
384 ILASCL = .FALSE.
385 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
386 ANRMTO = SMLNUM
387 ILASCL = .TRUE.
388 ELSE IF( ANRM.GT.BIGNUM ) THEN
389 ANRMTO = BIGNUM
390 ILASCL = .TRUE.
391 END IF
392 IF( ILASCL )
393 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
394 *
395 * Scale B if max element outside range [SMLNUM,BIGNUM]
396 *
397 BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
398 ILBSCL = .FALSE.
399 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
400 BNRMTO = SMLNUM
401 ILBSCL = .TRUE.
402 ELSE IF( BNRM.GT.BIGNUM ) THEN
403 BNRMTO = BIGNUM
404 ILBSCL = .TRUE.
405 END IF
406 IF( ILBSCL )
407 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
408 *
409 * Permute the matrix to make it more nearly triangular
410 * (Real Workspace: need 6*N)
411 *
412 ILEFT = 1
413 IRIGHT = N + 1
414 IRWRK = IRIGHT + N
415 CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
416 $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
417 *
418 * Reduce B to triangular form (QR decomposition of B)
419 * (Complex Workspace: need N, prefer N*NB)
420 *
421 IROWS = IHI + 1 - ILO
422 ICOLS = N + 1 - ILO
423 ITAU = 1
424 IWRK = ITAU + IROWS
425 CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
426 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
427 *
428 * Apply the unitary transformation to matrix A
429 * (Complex Workspace: need N, prefer N*NB)
430 *
431 CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
432 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
433 $ LWORK+1-IWRK, IERR )
434 *
435 * Initialize VSL
436 * (Complex Workspace: need N, prefer N*NB)
437 *
438 IF( ILVSL ) THEN
439 CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
440 IF( IROWS.GT.1 ) THEN
441 CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
442 $ VSL( ILO+1, ILO ), LDVSL )
443 END IF
444 CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
445 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
446 END IF
447 *
448 * Initialize VSR
449 *
450 IF( ILVSR )
451 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
452 *
453 * Reduce to generalized Hessenberg form
454 * (Workspace: none needed)
455 *
456 CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
457 $ LDVSL, VSR, LDVSR, IERR )
458 *
459 SDIM = 0
460 *
461 * Perform QZ algorithm, computing Schur vectors if desired
462 * (Complex Workspace: need N)
463 * (Real Workspace: need N)
464 *
465 IWRK = ITAU
466 CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
467 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
468 $ LWORK+1-IWRK, RWORK( IRWRK ), IERR )
469 IF( IERR.NE.0 ) THEN
470 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
471 INFO = IERR
472 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
473 INFO = IERR - N
474 ELSE
475 INFO = N + 1
476 END IF
477 GO TO 40
478 END IF
479 *
480 * Sort eigenvalues ALPHA/BETA and compute the reciprocal of
481 * condition number(s)
482 *
483 IF( WANTST ) THEN
484 *
485 * Undo scaling on eigenvalues before SELCTGing
486 *
487 IF( ILASCL )
488 $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
489 IF( ILBSCL )
490 $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
491 *
492 * Select eigenvalues
493 *
494 DO 10 I = 1, N
495 BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
496 10 CONTINUE
497 *
498 * Reorder eigenvalues, transform Generalized Schur vectors, and
499 * compute reciprocal condition numbers
500 * (Complex Workspace: If IJOB >= 1, need MAX(1, 2*SDIM*(N-SDIM))
501 * otherwise, need 1 )
502 *
503 CALL ZTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
504 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PL, PR,
505 $ DIF, WORK( IWRK ), LWORK-IWRK+1, IWORK, LIWORK,
506 $ IERR )
507 *
508 IF( IJOB.GE.1 )
509 $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
510 IF( IERR.EQ.-21 ) THEN
511 *
512 * not enough complex workspace
513 *
514 INFO = -21
515 ELSE
516 IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN
517 RCONDE( 1 ) = PL
518 RCONDE( 2 ) = PR
519 END IF
520 IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
521 RCONDV( 1 ) = DIF( 1 )
522 RCONDV( 2 ) = DIF( 2 )
523 END IF
524 IF( IERR.EQ.1 )
525 $ INFO = N + 3
526 END IF
527 *
528 END IF
529 *
530 * Apply permutation to VSL and VSR
531 * (Workspace: none needed)
532 *
533 IF( ILVSL )
534 $ CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
535 $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
536 *
537 IF( ILVSR )
538 $ CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
539 $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
540 *
541 * Undo scaling
542 *
543 IF( ILASCL ) THEN
544 CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
545 CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
546 END IF
547 *
548 IF( ILBSCL ) THEN
549 CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
550 CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
551 END IF
552 *
553 IF( WANTST ) THEN
554 *
555 * Check if reordering is correct
556 *
557 LASTSL = .TRUE.
558 SDIM = 0
559 DO 30 I = 1, N
560 CURSL = SELCTG( ALPHA( I ), BETA( I ) )
561 IF( CURSL )
562 $ SDIM = SDIM + 1
563 IF( CURSL .AND. .NOT.LASTSL )
564 $ INFO = N + 2
565 LASTSL = CURSL
566 30 CONTINUE
567 *
568 END IF
569 *
570 40 CONTINUE
571 *
572 WORK( 1 ) = MAXWRK
573 IWORK( 1 ) = LIWMIN
574 *
575 RETURN
576 *
577 * End of ZGGESX
578 *
579 END
2 $ B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR,
3 $ LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK,
4 $ IWORK, LIWORK, BWORK, INFO )
5 *
6 * -- LAPACK driver routine (version 3.2) --
7 * -- LAPACK is a software package provided by Univ. of Tennessee, --
8 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9 * November 2006
10 *
11 * .. Scalar Arguments ..
12 CHARACTER JOBVSL, JOBVSR, SENSE, SORT
13 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
14 $ SDIM
15 * ..
16 * .. Array Arguments ..
17 LOGICAL BWORK( * )
18 INTEGER IWORK( * )
19 DOUBLE PRECISION RCONDE( 2 ), RCONDV( 2 ), RWORK( * )
20 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
21 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
22 $ WORK( * )
23 * ..
24 * .. Function Arguments ..
25 LOGICAL SELCTG
26 EXTERNAL SELCTG
27 * ..
28 *
29 * Purpose
30 * =======
31 *
32 * ZGGESX computes for a pair of N-by-N complex nonsymmetric matrices
33 * (A,B), the generalized eigenvalues, the complex Schur form (S,T),
34 * and, optionally, the left and/or right matrices of Schur vectors (VSL
35 * and VSR). This gives the generalized Schur factorization
36 *
37 * (A,B) = ( (VSL) S (VSR)**H, (VSL) T (VSR)**H )
38 *
39 * where (VSR)**H is the conjugate-transpose of VSR.
40 *
41 * Optionally, it also orders the eigenvalues so that a selected cluster
42 * of eigenvalues appears in the leading diagonal blocks of the upper
43 * triangular matrix S and the upper triangular matrix T; computes
44 * a reciprocal condition number for the average of the selected
45 * eigenvalues (RCONDE); and computes a reciprocal condition number for
46 * the right and left deflating subspaces corresponding to the selected
47 * eigenvalues (RCONDV). The leading columns of VSL and VSR then form
48 * an orthonormal basis for the corresponding left and right eigenspaces
49 * (deflating subspaces).
50 *
51 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
52 * or a ratio alpha/beta = w, such that A - w*B is singular. It is
53 * usually represented as the pair (alpha,beta), as there is a
54 * reasonable interpretation for beta=0 or for both being zero.
55 *
56 * A pair of matrices (S,T) is in generalized complex Schur form if T is
57 * upper triangular with non-negative diagonal and S is upper
58 * triangular.
59 *
60 * Arguments
61 * =========
62 *
63 * JOBVSL (input) CHARACTER*1
64 * = 'N': do not compute the left Schur vectors;
65 * = 'V': compute the left Schur vectors.
66 *
67 * JOBVSR (input) CHARACTER*1
68 * = 'N': do not compute the right Schur vectors;
69 * = 'V': compute the right Schur vectors.
70 *
71 * SORT (input) CHARACTER*1
72 * Specifies whether or not to order the eigenvalues on the
73 * diagonal of the generalized Schur form.
74 * = 'N': Eigenvalues are not ordered;
75 * = 'S': Eigenvalues are ordered (see SELCTG).
76 *
77 * SELCTG (external procedure) LOGICAL FUNCTION of two COMPLEX*16 arguments
78 * SELCTG must be declared EXTERNAL in the calling subroutine.
79 * If SORT = 'N', SELCTG is not referenced.
80 * If SORT = 'S', SELCTG is used to select eigenvalues to sort
81 * to the top left of the Schur form.
82 * Note that a selected complex eigenvalue may no longer satisfy
83 * SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
84 * ordering may change the value of complex eigenvalues
85 * (especially if the eigenvalue is ill-conditioned), in this
86 * case INFO is set to N+3 see INFO below).
87 *
88 * SENSE (input) CHARACTER*1
89 * Determines which reciprocal condition numbers are computed.
90 * = 'N' : None are computed;
91 * = 'E' : Computed for average of selected eigenvalues only;
92 * = 'V' : Computed for selected deflating subspaces only;
93 * = 'B' : Computed for both.
94 * If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
95 *
96 * N (input) INTEGER
97 * The order of the matrices A, B, VSL, and VSR. N >= 0.
98 *
99 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
100 * On entry, the first of the pair of matrices.
101 * On exit, A has been overwritten by its generalized Schur
102 * form S.
103 *
104 * LDA (input) INTEGER
105 * The leading dimension of A. LDA >= max(1,N).
106 *
107 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
108 * On entry, the second of the pair of matrices.
109 * On exit, B has been overwritten by its generalized Schur
110 * form T.
111 *
112 * LDB (input) INTEGER
113 * The leading dimension of B. LDB >= max(1,N).
114 *
115 * SDIM (output) INTEGER
116 * If SORT = 'N', SDIM = 0.
117 * If SORT = 'S', SDIM = number of eigenvalues (after sorting)
118 * for which SELCTG is true.
119 *
120 * ALPHA (output) COMPLEX*16 array, dimension (N)
121 * BETA (output) COMPLEX*16 array, dimension (N)
122 * On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
123 * generalized eigenvalues. ALPHA(j) and BETA(j),j=1,...,N are
124 * the diagonals of the complex Schur form (S,T). BETA(j) will
125 * be non-negative real.
126 *
127 * Note: the quotients ALPHA(j)/BETA(j) may easily over- or
128 * underflow, and BETA(j) may even be zero. Thus, the user
129 * should avoid naively computing the ratio alpha/beta.
130 * However, ALPHA will be always less than and usually
131 * comparable with norm(A) in magnitude, and BETA always less
132 * than and usually comparable with norm(B).
133 *
134 * VSL (output) COMPLEX*16 array, dimension (LDVSL,N)
135 * If JOBVSL = 'V', VSL will contain the left Schur vectors.
136 * Not referenced if JOBVSL = 'N'.
137 *
138 * LDVSL (input) INTEGER
139 * The leading dimension of the matrix VSL. LDVSL >=1, and
140 * if JOBVSL = 'V', LDVSL >= N.
141 *
142 * VSR (output) COMPLEX*16 array, dimension (LDVSR,N)
143 * If JOBVSR = 'V', VSR will contain the right Schur vectors.
144 * Not referenced if JOBVSR = 'N'.
145 *
146 * LDVSR (input) INTEGER
147 * The leading dimension of the matrix VSR. LDVSR >= 1, and
148 * if JOBVSR = 'V', LDVSR >= N.
149 *
150 * RCONDE (output) DOUBLE PRECISION array, dimension ( 2 )
151 * If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
152 * reciprocal condition numbers for the average of the selected
153 * eigenvalues.
154 * Not referenced if SENSE = 'N' or 'V'.
155 *
156 * RCONDV (output) DOUBLE PRECISION array, dimension ( 2 )
157 * If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
158 * reciprocal condition number for the selected deflating
159 * subspaces.
160 * Not referenced if SENSE = 'N' or 'E'.
161 *
162 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
163 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
164 *
165 * LWORK (input) INTEGER
166 * The dimension of the array WORK.
167 * If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
168 * LWORK >= MAX(1,2*N,2*SDIM*(N-SDIM)), else
169 * LWORK >= MAX(1,2*N). Note that 2*SDIM*(N-SDIM) <= N*N/2.
170 * Note also that an error is only returned if
171 * LWORK < MAX(1,2*N), but if SENSE = 'E' or 'V' or 'B' this may
172 * not be large enough.
173 *
174 * If LWORK = -1, then a workspace query is assumed; the routine
175 * only calculates the bound on the optimal size of the WORK
176 * array and the minimum size of the IWORK array, returns these
177 * values as the first entries of the WORK and IWORK arrays, and
178 * no error message related to LWORK or LIWORK is issued by
179 * XERBLA.
180 *
181 * RWORK (workspace) DOUBLE PRECISION array, dimension ( 8*N )
182 * Real workspace.
183 *
184 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
185 * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
186 *
187 * LIWORK (input) INTEGER
188 * The dimension of the array IWORK.
189 * If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
190 * LIWORK >= N+2.
191 *
192 * If LIWORK = -1, then a workspace query is assumed; the
193 * routine only calculates the bound on the optimal size of the
194 * WORK array and the minimum size of the IWORK array, returns
195 * these values as the first entries of the WORK and IWORK
196 * arrays, and no error message related to LWORK or LIWORK is
197 * issued by XERBLA.
198 *
199 * BWORK (workspace) LOGICAL array, dimension (N)
200 * Not referenced if SORT = 'N'.
201 *
202 * INFO (output) INTEGER
203 * = 0: successful exit
204 * < 0: if INFO = -i, the i-th argument had an illegal value.
205 * = 1,...,N:
206 * The QZ iteration failed. (A,B) are not in Schur
207 * form, but ALPHA(j) and BETA(j) should be correct for
208 * j=INFO+1,...,N.
209 * > N: =N+1: other than QZ iteration failed in ZHGEQZ
210 * =N+2: after reordering, roundoff changed values of
211 * some complex eigenvalues so that leading
212 * eigenvalues in the Generalized Schur form no
213 * longer satisfy SELCTG=.TRUE. This could also
214 * be caused due to scaling.
215 * =N+3: reordering failed in ZTGSEN.
216 *
217 * =====================================================================
218 *
219 * .. Parameters ..
220 DOUBLE PRECISION ZERO, ONE
221 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
222 COMPLEX*16 CZERO, CONE
223 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
224 $ CONE = ( 1.0D+0, 0.0D+0 ) )
225 * ..
226 * .. Local Scalars ..
227 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
228 $ LQUERY, WANTSB, WANTSE, WANTSN, WANTST, WANTSV
229 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
230 $ ILEFT, ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK,
231 $ LIWMIN, LWRK, MAXWRK, MINWRK
232 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
233 $ PR, SMLNUM
234 * ..
235 * .. Local Arrays ..
236 DOUBLE PRECISION DIF( 2 )
237 * ..
238 * .. External Subroutines ..
239 EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
240 $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR,
241 $ ZUNMQR
242 * ..
243 * .. External Functions ..
244 LOGICAL LSAME
245 INTEGER ILAENV
246 DOUBLE PRECISION DLAMCH, ZLANGE
247 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
248 * ..
249 * .. Intrinsic Functions ..
250 INTRINSIC MAX, SQRT
251 * ..
252 * .. Executable Statements ..
253 *
254 * Decode the input arguments
255 *
256 IF( LSAME( JOBVSL, 'N' ) ) THEN
257 IJOBVL = 1
258 ILVSL = .FALSE.
259 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
260 IJOBVL = 2
261 ILVSL = .TRUE.
262 ELSE
263 IJOBVL = -1
264 ILVSL = .FALSE.
265 END IF
266 *
267 IF( LSAME( JOBVSR, 'N' ) ) THEN
268 IJOBVR = 1
269 ILVSR = .FALSE.
270 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
271 IJOBVR = 2
272 ILVSR = .TRUE.
273 ELSE
274 IJOBVR = -1
275 ILVSR = .FALSE.
276 END IF
277 *
278 WANTST = LSAME( SORT, 'S' )
279 WANTSN = LSAME( SENSE, 'N' )
280 WANTSE = LSAME( SENSE, 'E' )
281 WANTSV = LSAME( SENSE, 'V' )
282 WANTSB = LSAME( SENSE, 'B' )
283 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
284 IF( WANTSN ) THEN
285 IJOB = 0
286 ELSE IF( WANTSE ) THEN
287 IJOB = 1
288 ELSE IF( WANTSV ) THEN
289 IJOB = 2
290 ELSE IF( WANTSB ) THEN
291 IJOB = 4
292 END IF
293 *
294 * Test the input arguments
295 *
296 INFO = 0
297 IF( IJOBVL.LE.0 ) THEN
298 INFO = -1
299 ELSE IF( IJOBVR.LE.0 ) THEN
300 INFO = -2
301 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
302 INFO = -3
303 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
304 $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
305 INFO = -5
306 ELSE IF( N.LT.0 ) THEN
307 INFO = -6
308 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
309 INFO = -8
310 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
311 INFO = -10
312 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
313 INFO = -15
314 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
315 INFO = -17
316 END IF
317 *
318 * Compute workspace
319 * (Note: Comments in the code beginning "Workspace:" describe the
320 * minimal amount of workspace needed at that point in the code,
321 * as well as the preferred amount for good performance.
322 * NB refers to the optimal block size for the immediately
323 * following subroutine, as returned by ILAENV.)
324 *
325 IF( INFO.EQ.0 ) THEN
326 IF( N.GT.0) THEN
327 MINWRK = 2*N
328 MAXWRK = N*(1 + ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
329 MAXWRK = MAX( MAXWRK, N*( 1 +
330 $ ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, -1 ) ) )
331 IF( ILVSL ) THEN
332 MAXWRK = MAX( MAXWRK, N*( 1 +
333 $ ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) ) )
334 END IF
335 LWRK = MAXWRK
336 IF( IJOB.GE.1 )
337 $ LWRK = MAX( LWRK, N*N/2 )
338 ELSE
339 MINWRK = 1
340 MAXWRK = 1
341 LWRK = 1
342 END IF
343 WORK( 1 ) = LWRK
344 IF( WANTSN .OR. N.EQ.0 ) THEN
345 LIWMIN = 1
346 ELSE
347 LIWMIN = N + 2
348 END IF
349 IWORK( 1 ) = LIWMIN
350 *
351 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
352 INFO = -21
353 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY) THEN
354 INFO = -24
355 END IF
356 END IF
357 *
358 IF( INFO.NE.0 ) THEN
359 CALL XERBLA( 'ZGGESX', -INFO )
360 RETURN
361 ELSE IF (LQUERY) THEN
362 RETURN
363 END IF
364 *
365 * Quick return if possible
366 *
367 IF( N.EQ.0 ) THEN
368 SDIM = 0
369 RETURN
370 END IF
371 *
372 * Get machine constants
373 *
374 EPS = DLAMCH( 'P' )
375 SMLNUM = DLAMCH( 'S' )
376 BIGNUM = ONE / SMLNUM
377 CALL DLABAD( SMLNUM, BIGNUM )
378 SMLNUM = SQRT( SMLNUM ) / EPS
379 BIGNUM = ONE / SMLNUM
380 *
381 * Scale A if max element outside range [SMLNUM,BIGNUM]
382 *
383 ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
384 ILASCL = .FALSE.
385 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
386 ANRMTO = SMLNUM
387 ILASCL = .TRUE.
388 ELSE IF( ANRM.GT.BIGNUM ) THEN
389 ANRMTO = BIGNUM
390 ILASCL = .TRUE.
391 END IF
392 IF( ILASCL )
393 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
394 *
395 * Scale B if max element outside range [SMLNUM,BIGNUM]
396 *
397 BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
398 ILBSCL = .FALSE.
399 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
400 BNRMTO = SMLNUM
401 ILBSCL = .TRUE.
402 ELSE IF( BNRM.GT.BIGNUM ) THEN
403 BNRMTO = BIGNUM
404 ILBSCL = .TRUE.
405 END IF
406 IF( ILBSCL )
407 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
408 *
409 * Permute the matrix to make it more nearly triangular
410 * (Real Workspace: need 6*N)
411 *
412 ILEFT = 1
413 IRIGHT = N + 1
414 IRWRK = IRIGHT + N
415 CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
416 $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
417 *
418 * Reduce B to triangular form (QR decomposition of B)
419 * (Complex Workspace: need N, prefer N*NB)
420 *
421 IROWS = IHI + 1 - ILO
422 ICOLS = N + 1 - ILO
423 ITAU = 1
424 IWRK = ITAU + IROWS
425 CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
426 $ WORK( IWRK ), LWORK+1-IWRK, IERR )
427 *
428 * Apply the unitary transformation to matrix A
429 * (Complex Workspace: need N, prefer N*NB)
430 *
431 CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
432 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
433 $ LWORK+1-IWRK, IERR )
434 *
435 * Initialize VSL
436 * (Complex Workspace: need N, prefer N*NB)
437 *
438 IF( ILVSL ) THEN
439 CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
440 IF( IROWS.GT.1 ) THEN
441 CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
442 $ VSL( ILO+1, ILO ), LDVSL )
443 END IF
444 CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
445 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
446 END IF
447 *
448 * Initialize VSR
449 *
450 IF( ILVSR )
451 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
452 *
453 * Reduce to generalized Hessenberg form
454 * (Workspace: none needed)
455 *
456 CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
457 $ LDVSL, VSR, LDVSR, IERR )
458 *
459 SDIM = 0
460 *
461 * Perform QZ algorithm, computing Schur vectors if desired
462 * (Complex Workspace: need N)
463 * (Real Workspace: need N)
464 *
465 IWRK = ITAU
466 CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
467 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
468 $ LWORK+1-IWRK, RWORK( IRWRK ), IERR )
469 IF( IERR.NE.0 ) THEN
470 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
471 INFO = IERR
472 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
473 INFO = IERR - N
474 ELSE
475 INFO = N + 1
476 END IF
477 GO TO 40
478 END IF
479 *
480 * Sort eigenvalues ALPHA/BETA and compute the reciprocal of
481 * condition number(s)
482 *
483 IF( WANTST ) THEN
484 *
485 * Undo scaling on eigenvalues before SELCTGing
486 *
487 IF( ILASCL )
488 $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
489 IF( ILBSCL )
490 $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
491 *
492 * Select eigenvalues
493 *
494 DO 10 I = 1, N
495 BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
496 10 CONTINUE
497 *
498 * Reorder eigenvalues, transform Generalized Schur vectors, and
499 * compute reciprocal condition numbers
500 * (Complex Workspace: If IJOB >= 1, need MAX(1, 2*SDIM*(N-SDIM))
501 * otherwise, need 1 )
502 *
503 CALL ZTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
504 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PL, PR,
505 $ DIF, WORK( IWRK ), LWORK-IWRK+1, IWORK, LIWORK,
506 $ IERR )
507 *
508 IF( IJOB.GE.1 )
509 $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
510 IF( IERR.EQ.-21 ) THEN
511 *
512 * not enough complex workspace
513 *
514 INFO = -21
515 ELSE
516 IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN
517 RCONDE( 1 ) = PL
518 RCONDE( 2 ) = PR
519 END IF
520 IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
521 RCONDV( 1 ) = DIF( 1 )
522 RCONDV( 2 ) = DIF( 2 )
523 END IF
524 IF( IERR.EQ.1 )
525 $ INFO = N + 3
526 END IF
527 *
528 END IF
529 *
530 * Apply permutation to VSL and VSR
531 * (Workspace: none needed)
532 *
533 IF( ILVSL )
534 $ CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
535 $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
536 *
537 IF( ILVSR )
538 $ CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
539 $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
540 *
541 * Undo scaling
542 *
543 IF( ILASCL ) THEN
544 CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
545 CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
546 END IF
547 *
548 IF( ILBSCL ) THEN
549 CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
550 CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
551 END IF
552 *
553 IF( WANTST ) THEN
554 *
555 * Check if reordering is correct
556 *
557 LASTSL = .TRUE.
558 SDIM = 0
559 DO 30 I = 1, N
560 CURSL = SELCTG( ALPHA( I ), BETA( I ) )
561 IF( CURSL )
562 $ SDIM = SDIM + 1
563 IF( CURSL .AND. .NOT.LASTSL )
564 $ INFO = N + 2
565 LASTSL = CURSL
566 30 CONTINUE
567 *
568 END IF
569 *
570 40 CONTINUE
571 *
572 WORK( 1 ) = MAXWRK
573 IWORK( 1 ) = LIWMIN
574 *
575 RETURN
576 *
577 * End of ZGGESX
578 *
579 END