001:       SUBROUTINE DGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
002:      $                   B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
003:      $                   VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK,
004:      $                   LIWORK, BWORK, INFO )
005: *
006: *  -- LAPACK driver routine (version 3.2.1)                           --
007: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
008: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
009: *  -- April 2009                                                      --
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   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
020:      $                   B( LDB, * ), BETA( * ), RCONDE( 2 ),
021:      $                   RCONDV( 2 ), VSL( LDVSL, * ), VSR( LDVSR, * ),
022:      $                   WORK( * )
023: *     ..
024: *     .. Function Arguments ..
025:       LOGICAL            SELCTG
026:       EXTERNAL           SELCTG
027: *     ..
028: *
029: *  Purpose
030: *  =======
031: *
032: *  DGGESX computes for a pair of N-by-N real nonsymmetric matrices
033: *  (A,B), the generalized eigenvalues, the real Schur form (S,T), and,
034: *  optionally, the left and/or right matrices of Schur vectors (VSL and
035: *  VSR).  This gives the generalized Schur factorization
036: *
037: *       (A,B) = ( (VSL) S (VSR)**T, (VSL) T (VSR)**T )
038: *
039: *  Optionally, it also orders the eigenvalues so that a selected cluster
040: *  of eigenvalues appears in the leading diagonal blocks of the upper
041: *  quasi-triangular matrix S and the upper triangular matrix T; computes
042: *  a reciprocal condition number for the average of the selected
043: *  eigenvalues (RCONDE); and computes a reciprocal condition number for
044: *  the right and left deflating subspaces corresponding to the selected
045: *  eigenvalues (RCONDV). The leading columns of VSL and VSR then form
046: *  an orthonormal basis for the corresponding left and right eigenspaces
047: *  (deflating subspaces).
048: *
049: *  A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
050: *  or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
051: *  usually represented as the pair (alpha,beta), as there is a
052: *  reasonable interpretation for beta=0 or for both being zero.
053: *
054: *  A pair of matrices (S,T) is in generalized real Schur form if T is
055: *  upper triangular with non-negative diagonal and S is block upper
056: *  triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond
057: *  to real generalized eigenvalues, while 2-by-2 blocks of S will be
058: *  "standardized" by making the corresponding elements of T have the
059: *  form:
060: *          [  a  0  ]
061: *          [  0  b  ]
062: *
063: *  and the pair of corresponding 2-by-2 blocks in S and T will have a
064: *  complex conjugate pair of generalized eigenvalues.
065: *
066: *
067: *  Arguments
068: *  =========
069: *
070: *  JOBVSL  (input) CHARACTER*1
071: *          = 'N':  do not compute the left Schur vectors;
072: *          = 'V':  compute the left Schur vectors.
073: *
074: *  JOBVSR  (input) CHARACTER*1
075: *          = 'N':  do not compute the right Schur vectors;
076: *          = 'V':  compute the right Schur vectors.
077: *
078: *  SORT    (input) CHARACTER*1
079: *          Specifies whether or not to order the eigenvalues on the
080: *          diagonal of the generalized Schur form.
081: *          = 'N':  Eigenvalues are not ordered;
082: *          = 'S':  Eigenvalues are ordered (see SELCTG).
083: *
084: *  SELCTG  (external procedure) LOGICAL FUNCTION of three DOUBLE PRECISION arguments
085: *          SELCTG must be declared EXTERNAL in the calling subroutine.
086: *          If SORT = 'N', SELCTG is not referenced.
087: *          If SORT = 'S', SELCTG is used to select eigenvalues to sort
088: *          to the top left of the Schur form.
089: *          An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
090: *          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
091: *          one of a complex conjugate pair of eigenvalues is selected,
092: *          then both complex eigenvalues are selected.
093: *          Note that a selected complex eigenvalue may no longer satisfy
094: *          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) = .TRUE. after ordering,
095: *          since ordering may change the value of complex eigenvalues
096: *          (especially if the eigenvalue is ill-conditioned), in this
097: *          case INFO is set to N+3.
098: *
099: *  SENSE   (input) CHARACTER*1
100: *          Determines which reciprocal condition numbers are computed.
101: *          = 'N' : None are computed;
102: *          = 'E' : Computed for average of selected eigenvalues only;
103: *          = 'V' : Computed for selected deflating subspaces only;
104: *          = 'B' : Computed for both.
105: *          If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
106: *
107: *  N       (input) INTEGER
108: *          The order of the matrices A, B, VSL, and VSR.  N >= 0.
109: *
110: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
111: *          On entry, the first of the pair of matrices.
112: *          On exit, A has been overwritten by its generalized Schur
113: *          form S.
114: *
115: *  LDA     (input) INTEGER
116: *          The leading dimension of A.  LDA >= max(1,N).
117: *
118: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)
119: *          On entry, the second of the pair of matrices.
120: *          On exit, B has been overwritten by its generalized Schur
121: *          form T.
122: *
123: *  LDB     (input) INTEGER
124: *          The leading dimension of B.  LDB >= max(1,N).
125: *
126: *  SDIM    (output) INTEGER
127: *          If SORT = 'N', SDIM = 0.
128: *          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
129: *          for which SELCTG is true.  (Complex conjugate pairs for which
130: *          SELCTG is true for either eigenvalue count as 2.)
131: *
132: *  ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
133: *  ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
134: *  BETA    (output) DOUBLE PRECISION array, dimension (N)
135: *          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
136: *          be the generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i
137: *          and BETA(j),j=1,...,N  are the diagonals of the complex Schur
138: *          form (S,T) that would result if the 2-by-2 diagonal blocks of
139: *          the real Schur form of (A,B) were further reduced to
140: *          triangular form using 2-by-2 complex unitary transformations.
141: *          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
142: *          positive, then the j-th and (j+1)-st eigenvalues are a
143: *          complex conjugate pair, with ALPHAI(j+1) negative.
144: *
145: *          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
146: *          may easily over- or underflow, and BETA(j) may even be zero.
147: *          Thus, the user should avoid naively computing the ratio.
148: *          However, ALPHAR and ALPHAI will be always less than and
149: *          usually comparable with norm(A) in magnitude, and BETA always
150: *          less than and usually comparable with norm(B).
151: *
152: *  VSL     (output) DOUBLE PRECISION array, dimension (LDVSL,N)
153: *          If JOBVSL = 'V', VSL will contain the left Schur vectors.
154: *          Not referenced if JOBVSL = 'N'.
155: *
156: *  LDVSL   (input) INTEGER
157: *          The leading dimension of the matrix VSL. LDVSL >=1, and
158: *          if JOBVSL = 'V', LDVSL >= N.
159: *
160: *  VSR     (output) DOUBLE PRECISION array, dimension (LDVSR,N)
161: *          If JOBVSR = 'V', VSR will contain the right Schur vectors.
162: *          Not referenced if JOBVSR = 'N'.
163: *
164: *  LDVSR   (input) INTEGER
165: *          The leading dimension of the matrix VSR. LDVSR >= 1, and
166: *          if JOBVSR = 'V', LDVSR >= N.
167: *
168: *  RCONDE  (output) DOUBLE PRECISION array, dimension ( 2 )
169: *          If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
170: *          reciprocal condition numbers for the average of the selected
171: *          eigenvalues.
172: *          Not referenced if SENSE = 'N' or 'V'.
173: *
174: *  RCONDV  (output) DOUBLE PRECISION array, dimension ( 2 )
175: *          If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
176: *          reciprocal condition numbers for the selected deflating
177: *          subspaces.
178: *          Not referenced if SENSE = 'N' or 'E'.
179: *
180: *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
181: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
182: *
183: *  LWORK   (input) INTEGER
184: *          The dimension of the array WORK.
185: *          If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
186: *          LWORK >= max( 8*N, 6*N+16, 2*SDIM*(N-SDIM) ), else
187: *          LWORK >= max( 8*N, 6*N+16 ).
188: *          Note that 2*SDIM*(N-SDIM) <= N*N/2.
189: *          Note also that an error is only returned if
190: *          LWORK < max( 8*N, 6*N+16), but if SENSE = 'E' or 'V' or 'B'
191: *          this may not be large enough.
192: *
193: *          If LWORK = -1, then a workspace query is assumed; the routine
194: *          only calculates the bound on the optimal size of the WORK
195: *          array and the minimum size of the IWORK array, returns these
196: *          values as the first entries of the WORK and IWORK arrays, and
197: *          no error message related to LWORK or LIWORK is issued by
198: *          XERBLA.
199: *
200: *  IWORK   (workspace) INTEGER array, dimension (MAX(1,LIWORK))
201: *          On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
202: *
203: *  LIWORK  (input) INTEGER
204: *          The dimension of the array IWORK.
205: *          If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
206: *          LIWORK >= N+6.
207: *
208: *          If LIWORK = -1, then a workspace query is assumed; the
209: *          routine only calculates the bound on the optimal size of the
210: *          WORK array and the minimum size of the IWORK array, returns
211: *          these values as the first entries of the WORK and IWORK
212: *          arrays, and no error message related to LWORK or LIWORK is
213: *          issued by XERBLA.
214: *
215: *  BWORK   (workspace) LOGICAL array, dimension (N)
216: *          Not referenced if SORT = 'N'.
217: *
218: *  INFO    (output) INTEGER
219: *          = 0:  successful exit
220: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
221: *          = 1,...,N:
222: *                The QZ iteration failed.  (A,B) are not in Schur
223: *                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
224: *                be correct for j=INFO+1,...,N.
225: *          > N:  =N+1: other than QZ iteration failed in DHGEQZ
226: *                =N+2: after reordering, roundoff changed values of
227: *                      some complex eigenvalues so that leading
228: *                      eigenvalues in the Generalized Schur form no
229: *                      longer satisfy SELCTG=.TRUE.  This could also
230: *                      be caused due to scaling.
231: *                =N+3: reordering failed in DTGSEN.
232: *
233: *  Further Details
234: *  ===============
235: *
236: *  An approximate (asymptotic) bound on the average absolute error of
237: *  the selected eigenvalues is
238: *
239: *       EPS * norm((A, B)) / RCONDE( 1 ).
240: *
241: *  An approximate (asymptotic) bound on the maximum angular error in
242: *  the computed deflating subspaces is
243: *
244: *       EPS * norm((A, B)) / RCONDV( 2 ).
245: *
246: *  See LAPACK User's Guide, section 4.11 for more information.
247: *
248: *  =====================================================================
249: *
250: *     .. Parameters ..
251:       DOUBLE PRECISION   ZERO, ONE
252:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
253: *     ..
254: *     .. Local Scalars ..
255:       LOGICAL            CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
256:      $                   LQUERY, LST2SL, WANTSB, WANTSE, WANTSN, WANTST,
257:      $                   WANTSV
258:       INTEGER            I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
259:      $                   ILEFT, ILO, IP, IRIGHT, IROWS, ITAU, IWRK,
260:      $                   LIWMIN, LWRK, MAXWRK, MINWRK
261:       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
262:      $                   PR, SAFMAX, SAFMIN, SMLNUM
263: *     ..
264: *     .. Local Arrays ..
265:       DOUBLE PRECISION   DIF( 2 )
266: *     ..
267: *     .. External Subroutines ..
268:       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
269:      $                   DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN,
270:      $                   XERBLA
271: *     ..
272: *     .. External Functions ..
273:       LOGICAL            LSAME
274:       INTEGER            ILAENV
275:       DOUBLE PRECISION   DLAMCH, DLANGE
276:       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
277: *     ..
278: *     .. Intrinsic Functions ..
279:       INTRINSIC          ABS, MAX, SQRT
280: *     ..
281: *     .. Executable Statements ..
282: *
283: *     Decode the input arguments
284: *
285:       IF( LSAME( JOBVSL, 'N' ) ) THEN
286:          IJOBVL = 1
287:          ILVSL = .FALSE.
288:       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
289:          IJOBVL = 2
290:          ILVSL = .TRUE.
291:       ELSE
292:          IJOBVL = -1
293:          ILVSL = .FALSE.
294:       END IF
295: *
296:       IF( LSAME( JOBVSR, 'N' ) ) THEN
297:          IJOBVR = 1
298:          ILVSR = .FALSE.
299:       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
300:          IJOBVR = 2
301:          ILVSR = .TRUE.
302:       ELSE
303:          IJOBVR = -1
304:          ILVSR = .FALSE.
305:       END IF
306: *
307:       WANTST = LSAME( SORT, 'S' )
308:       WANTSN = LSAME( SENSE, 'N' )
309:       WANTSE = LSAME( SENSE, 'E' )
310:       WANTSV = LSAME( SENSE, 'V' )
311:       WANTSB = LSAME( SENSE, 'B' )
312:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
313:       IF( WANTSN ) THEN
314:          IJOB = 0
315:       ELSE IF( WANTSE ) THEN
316:          IJOB = 1
317:       ELSE IF( WANTSV ) THEN
318:          IJOB = 2
319:       ELSE IF( WANTSB ) THEN
320:          IJOB = 4
321:       END IF
322: *
323: *     Test the input arguments
324: *
325:       INFO = 0
326:       IF( IJOBVL.LE.0 ) THEN
327:          INFO = -1
328:       ELSE IF( IJOBVR.LE.0 ) THEN
329:          INFO = -2
330:       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
331:          INFO = -3
332:       ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
333:      $         ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
334:          INFO = -5
335:       ELSE IF( N.LT.0 ) THEN
336:          INFO = -6
337:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
338:          INFO = -8
339:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
340:          INFO = -10
341:       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
342:          INFO = -16
343:       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
344:          INFO = -18
345:       END IF
346: *
347: *     Compute workspace
348: *      (Note: Comments in the code beginning "Workspace:" describe the
349: *       minimal amount of workspace needed at that point in the code,
350: *       as well as the preferred amount for good performance.
351: *       NB refers to the optimal block size for the immediately
352: *       following subroutine, as returned by ILAENV.)
353: *
354:       IF( INFO.EQ.0 ) THEN
355:          IF( N.GT.0) THEN
356:             MINWRK = MAX( 8*N, 6*N + 16 )
357:             MAXWRK = MINWRK - N +
358:      $               N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 )
359:             MAXWRK = MAX( MAXWRK, MINWRK - N +
360:      $               N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, -1 ) )
361:             IF( ILVSL ) THEN
362:                MAXWRK = MAX( MAXWRK, MINWRK - N +
363:      $                  N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) )
364:             END IF
365:             LWRK = MAXWRK
366:             IF( IJOB.GE.1 )
367:      $         LWRK = MAX( LWRK, N*N/2 )
368:          ELSE
369:             MINWRK = 1
370:             MAXWRK = 1
371:             LWRK   = 1
372:          END IF
373:          WORK( 1 ) = LWRK
374:          IF( WANTSN .OR. N.EQ.0 ) THEN
375:             LIWMIN = 1
376:          ELSE
377:             LIWMIN = N + 6
378:          END IF
379:          IWORK( 1 ) = LIWMIN
380: *
381:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
382:             INFO = -22
383:          ELSE IF( LIWORK.LT.LIWMIN  .AND. .NOT.LQUERY ) THEN
384:             INFO = -24
385:          END IF
386:       END IF
387: *
388:       IF( INFO.NE.0 ) THEN
389:          CALL XERBLA( 'DGGESX', -INFO )
390:          RETURN
391:       ELSE IF (LQUERY) THEN
392:          RETURN
393:       END IF
394: *
395: *     Quick return if possible
396: *
397:       IF( N.EQ.0 ) THEN
398:          SDIM = 0
399:          RETURN
400:       END IF
401: *
402: *     Get machine constants
403: *
404:       EPS = DLAMCH( 'P' )
405:       SAFMIN = DLAMCH( 'S' )
406:       SAFMAX = ONE / SAFMIN
407:       CALL DLABAD( SAFMIN, SAFMAX )
408:       SMLNUM = SQRT( SAFMIN ) / EPS
409:       BIGNUM = ONE / SMLNUM
410: *
411: *     Scale A if max element outside range [SMLNUM,BIGNUM]
412: *
413:       ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
414:       ILASCL = .FALSE.
415:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
416:          ANRMTO = SMLNUM
417:          ILASCL = .TRUE.
418:       ELSE IF( ANRM.GT.BIGNUM ) THEN
419:          ANRMTO = BIGNUM
420:          ILASCL = .TRUE.
421:       END IF
422:       IF( ILASCL )
423:      $   CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
424: *
425: *     Scale B if max element outside range [SMLNUM,BIGNUM]
426: *
427:       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
428:       ILBSCL = .FALSE.
429:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
430:          BNRMTO = SMLNUM
431:          ILBSCL = .TRUE.
432:       ELSE IF( BNRM.GT.BIGNUM ) THEN
433:          BNRMTO = BIGNUM
434:          ILBSCL = .TRUE.
435:       END IF
436:       IF( ILBSCL )
437:      $   CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
438: *
439: *     Permute the matrix to make it more nearly triangular
440: *     (Workspace: need 6*N + 2*N for permutation parameters)
441: *
442:       ILEFT = 1
443:       IRIGHT = N + 1
444:       IWRK = IRIGHT + N
445:       CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
446:      $             WORK( IRIGHT ), WORK( IWRK ), IERR )
447: *
448: *     Reduce B to triangular form (QR decomposition of B)
449: *     (Workspace: need N, prefer N*NB)
450: *
451:       IROWS = IHI + 1 - ILO
452:       ICOLS = N + 1 - ILO
453:       ITAU = IWRK
454:       IWRK = ITAU + IROWS
455:       CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
456:      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
457: *
458: *     Apply the orthogonal transformation to matrix A
459: *     (Workspace: need N, prefer N*NB)
460: *
461:       CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
462:      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
463:      $             LWORK+1-IWRK, IERR )
464: *
465: *     Initialize VSL
466: *     (Workspace: need N, prefer N*NB)
467: *
468:       IF( ILVSL ) THEN
469:          CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
470:          IF( IROWS.GT.1 ) THEN
471:             CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
472:      $                   VSL( ILO+1, ILO ), LDVSL )
473:          END IF
474:          CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
475:      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
476:       END IF
477: *
478: *     Initialize VSR
479: *
480:       IF( ILVSR )
481:      $   CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
482: *
483: *     Reduce to generalized Hessenberg form
484: *     (Workspace: none needed)
485: *
486:       CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
487:      $             LDVSL, VSR, LDVSR, IERR )
488: *
489:       SDIM = 0
490: *
491: *     Perform QZ algorithm, computing Schur vectors if desired
492: *     (Workspace: need N)
493: *
494:       IWRK = ITAU
495:       CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
496:      $             ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
497:      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
498:       IF( IERR.NE.0 ) THEN
499:          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
500:             INFO = IERR
501:          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
502:             INFO = IERR - N
503:          ELSE
504:             INFO = N + 1
505:          END IF
506:          GO TO 60
507:       END IF
508: *
509: *     Sort eigenvalues ALPHA/BETA and compute the reciprocal of
510: *     condition number(s)
511: *     (Workspace: If IJOB >= 1, need MAX( 8*(N+1), 2*SDIM*(N-SDIM) )
512: *                 otherwise, need 8*(N+1) )
513: *
514:       IF( WANTST ) THEN
515: *
516: *        Undo scaling on eigenvalues before SELCTGing
517: *
518:          IF( ILASCL ) THEN
519:             CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
520:      $                   IERR )
521:             CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
522:      $                   IERR )
523:          END IF
524:          IF( ILBSCL )
525:      $      CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
526: *
527: *        Select eigenvalues
528: *
529:          DO 10 I = 1, N
530:             BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
531:    10    CONTINUE
532: *
533: *        Reorder eigenvalues, transform Generalized Schur vectors, and
534: *        compute reciprocal condition numbers
535: *
536:          CALL DTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
537:      $                ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
538:      $                SDIM, PL, PR, DIF, WORK( IWRK ), LWORK-IWRK+1,
539:      $                IWORK, LIWORK, IERR )
540: *
541:          IF( IJOB.GE.1 )
542:      $      MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
543:          IF( IERR.EQ.-22 ) THEN
544: *
545: *            not enough real workspace
546: *
547:             INFO = -22
548:          ELSE
549:             IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN
550:                RCONDE( 1 ) = PL
551:                RCONDE( 2 ) = PR
552:             END IF
553:             IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
554:                RCONDV( 1 ) = DIF( 1 )
555:                RCONDV( 2 ) = DIF( 2 )
556:             END IF
557:             IF( IERR.EQ.1 )
558:      $         INFO = N + 3
559:          END IF
560: *
561:       END IF
562: *
563: *     Apply permutation to VSL and VSR
564: *     (Workspace: none needed)
565: *
566:       IF( ILVSL )
567:      $   CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
568:      $                WORK( IRIGHT ), N, VSL, LDVSL, IERR )
569: *
570:       IF( ILVSR )
571:      $   CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
572:      $                WORK( IRIGHT ), N, VSR, LDVSR, IERR )
573: *
574: *     Check if unscaling would cause over/underflow, if so, rescale
575: *     (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
576: *     B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
577: *
578:       IF( ILASCL ) THEN
579:          DO 20 I = 1, N
580:             IF( ALPHAI( I ).NE.ZERO ) THEN
581:                IF( ( ALPHAR( I ) / SAFMAX ).GT.( ANRMTO / ANRM ) .OR.
582:      $             ( SAFMIN / ALPHAR( I ) ).GT.( ANRM / ANRMTO ) ) THEN
583:                   WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) )
584:                   BETA( I ) = BETA( I )*WORK( 1 )
585:                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
586:                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
587:                ELSE IF( ( ALPHAI( I ) / SAFMAX ).GT.
588:      $                  ( ANRMTO / ANRM ) .OR.
589:      $                  ( SAFMIN / ALPHAI( I ) ).GT.( ANRM / ANRMTO ) )
590:      $                   THEN
591:                   WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) )
592:                   BETA( I ) = BETA( I )*WORK( 1 )
593:                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
594:                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
595:                END IF
596:             END IF
597:    20    CONTINUE
598:       END IF
599: *
600:       IF( ILBSCL ) THEN
601:          DO 30 I = 1, N
602:             IF( ALPHAI( I ).NE.ZERO ) THEN
603:                IF( ( BETA( I ) / SAFMAX ).GT.( BNRMTO / BNRM ) .OR.
604:      $             ( SAFMIN / BETA( I ) ).GT.( BNRM / BNRMTO ) ) THEN
605:                   WORK( 1 ) = ABS( B( I, I ) / BETA( I ) )
606:                   BETA( I ) = BETA( I )*WORK( 1 )
607:                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
608:                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
609:                END IF
610:             END IF
611:    30    CONTINUE
612:       END IF
613: *
614: *     Undo scaling
615: *
616:       IF( ILASCL ) THEN
617:          CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
618:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
619:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
620:       END IF
621: *
622:       IF( ILBSCL ) THEN
623:          CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
624:          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
625:       END IF
626: *
627:       IF( WANTST ) THEN
628: *
629: *        Check if reordering is correct
630: *
631:          LASTSL = .TRUE.
632:          LST2SL = .TRUE.
633:          SDIM = 0
634:          IP = 0
635:          DO 50 I = 1, N
636:             CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
637:             IF( ALPHAI( I ).EQ.ZERO ) THEN
638:                IF( CURSL )
639:      $            SDIM = SDIM + 1
640:                IP = 0
641:                IF( CURSL .AND. .NOT.LASTSL )
642:      $            INFO = N + 2
643:             ELSE
644:                IF( IP.EQ.1 ) THEN
645: *
646: *                 Last eigenvalue of conjugate pair
647: *
648:                   CURSL = CURSL .OR. LASTSL
649:                   LASTSL = CURSL
650:                   IF( CURSL )
651:      $               SDIM = SDIM + 2
652:                   IP = -1
653:                   IF( CURSL .AND. .NOT.LST2SL )
654:      $               INFO = N + 2
655:                ELSE
656: *
657: *                 First eigenvalue of conjugate pair
658: *
659:                   IP = 1
660:                END IF
661:             END IF
662:             LST2SL = LASTSL
663:             LASTSL = CURSL
664:    50    CONTINUE
665: *
666:       END IF
667: *
668:    60 CONTINUE
669: *
670:       WORK( 1 ) = MAXWRK
671:       IWORK( 1 ) = LIWMIN
672: *
673:       RETURN
674: *
675: *     End of DGGESX
676: *
677:       END
678: