001:       SUBROUTINE ZGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
002:      $                   B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR,
003:      $                   LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK,
004:      $                   IWORK, LIWORK, BWORK, INFO )
005: *
006: *  -- LAPACK driver routine (version 3.2) --
007: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
008: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
009: *     November 2006
010: *
011: *     .. Scalar Arguments ..
012:       CHARACTER          JOBVSL, JOBVSR, SENSE, SORT
013:       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
014:      $                   SDIM
015: *     ..
016: *     .. Array Arguments ..
017:       LOGICAL            BWORK( * )
018:       INTEGER            IWORK( * )
019:       DOUBLE PRECISION   RCONDE( 2 ), RCONDV( 2 ), RWORK( * )
020:       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
021:      $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
022:      $                   WORK( * )
023: *     ..
024: *     .. Function Arguments ..
025:       LOGICAL            SELCTG
026:       EXTERNAL           SELCTG
027: *     ..
028: *
029: *  Purpose
030: *  =======
031: *
032: *  ZGGESX computes for a pair of N-by-N complex nonsymmetric matrices
033: *  (A,B), the generalized eigenvalues, the complex Schur form (S,T),
034: *  and, optionally, the left and/or right matrices of Schur vectors (VSL
035: *  and VSR).  This gives the generalized Schur factorization
036: *
037: *       (A,B) = ( (VSL) S (VSR)**H, (VSL) T (VSR)**H )
038: *
039: *  where (VSR)**H is the conjugate-transpose of VSR.
040: *
041: *  Optionally, it also orders the eigenvalues so that a selected cluster
042: *  of eigenvalues appears in the leading diagonal blocks of the upper
043: *  triangular matrix S and the upper triangular matrix T; computes
044: *  a reciprocal condition number for the average of the selected
045: *  eigenvalues (RCONDE); and computes a reciprocal condition number for
046: *  the right and left deflating subspaces corresponding to the selected
047: *  eigenvalues (RCONDV). The leading columns of VSL and VSR then form
048: *  an orthonormal basis for the corresponding left and right eigenspaces
049: *  (deflating subspaces).
050: *
051: *  A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
052: *  or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
053: *  usually represented as the pair (alpha,beta), as there is a
054: *  reasonable interpretation for beta=0 or for both being zero.
055: *
056: *  A pair of matrices (S,T) is in generalized complex Schur form if T is
057: *  upper triangular with non-negative diagonal and S is upper
058: *  triangular.
059: *
060: *  Arguments
061: *  =========
062: *
063: *  JOBVSL  (input) CHARACTER*1
064: *          = 'N':  do not compute the left Schur vectors;
065: *          = 'V':  compute the left Schur vectors.
066: *
067: *  JOBVSR  (input) CHARACTER*1
068: *          = 'N':  do not compute the right Schur vectors;
069: *          = 'V':  compute the right Schur vectors.
070: *
071: *  SORT    (input) CHARACTER*1
072: *          Specifies whether or not to order the eigenvalues on the
073: *          diagonal of the generalized Schur form.
074: *          = 'N':  Eigenvalues are not ordered;
075: *          = 'S':  Eigenvalues are ordered (see SELCTG).
076: *
077: *  SELCTG  (external procedure) LOGICAL FUNCTION of two COMPLEX*16 arguments
078: *          SELCTG must be declared EXTERNAL in the calling subroutine.
079: *          If SORT = 'N', SELCTG is not referenced.
080: *          If SORT = 'S', SELCTG is used to select eigenvalues to sort
081: *          to the top left of the Schur form.
082: *          Note that a selected complex eigenvalue may no longer satisfy
083: *          SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
084: *          ordering may change the value of complex eigenvalues
085: *          (especially if the eigenvalue is ill-conditioned), in this
086: *          case INFO is set to N+3 see INFO below).
087: *
088: *  SENSE   (input) CHARACTER*1
089: *          Determines which reciprocal condition numbers are computed.
090: *          = 'N' : None are computed;
091: *          = 'E' : Computed for average of selected eigenvalues only;
092: *          = 'V' : Computed for selected deflating subspaces only;
093: *          = 'B' : Computed for both.
094: *          If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
095: *
096: *  N       (input) INTEGER
097: *          The order of the matrices A, B, VSL, and VSR.  N >= 0.
098: *
099: *  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
580: