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