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