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