001:       SUBROUTINE ZGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
002:      $                  VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          JOBVL, JOBVR
010:       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   RWORK( * )
014:       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
015:      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
016:      $                   WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  This routine is deprecated and has been replaced by routine ZGGEV.
023: *
024: *  ZGEGV computes the eigenvalues and, optionally, the left and/or right
025: *  eigenvectors of a complex matrix pair (A,B).
026: *  Given two square matrices A and B,
027: *  the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
028: *  eigenvalues lambda and corresponding (non-zero) eigenvectors x such
029: *  that
030: *     A*x = lambda*B*x.
031: *
032: *  An alternate form is to find the eigenvalues mu and corresponding
033: *  eigenvectors y such that
034: *     mu*A*y = B*y.
035: *
036: *  These two forms are equivalent with mu = 1/lambda and x = y if
037: *  neither lambda nor mu is zero.  In order to deal with the case that
038: *  lambda or mu is zero or small, two values alpha and beta are returned
039: *  for each eigenvalue, such that lambda = alpha/beta and
040: *  mu = beta/alpha.
041: *
042: *  The vectors x and y in the above equations are right eigenvectors of
043: *  the matrix pair (A,B).  Vectors u and v satisfying
044: *     u**H*A = lambda*u**H*B  or  mu*v**H*A = v**H*B
045: *  are left eigenvectors of (A,B).
046: *
047: *  Note: this routine performs "full balancing" on A and B -- see
048: *  "Further Details", below.
049: *
050: *  Arguments
051: *  =========
052: *
053: *  JOBVL   (input) CHARACTER*1
054: *          = 'N':  do not compute the left generalized eigenvectors;
055: *          = 'V':  compute the left generalized eigenvectors (returned
056: *                  in VL).
057: *
058: *  JOBVR   (input) CHARACTER*1
059: *          = 'N':  do not compute the right generalized eigenvectors;
060: *          = 'V':  compute the right generalized eigenvectors (returned
061: *                  in VR).
062: *
063: *  N       (input) INTEGER
064: *          The order of the matrices A, B, VL, and VR.  N >= 0.
065: *
066: *  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
067: *          On entry, the matrix A.
068: *          If JOBVL = 'V' or JOBVR = 'V', then on exit A
069: *          contains the Schur form of A from the generalized Schur
070: *          factorization of the pair (A,B) after balancing.  If no
071: *          eigenvectors were computed, then only the diagonal elements
072: *          of the Schur form will be correct.  See ZGGHRD and ZHGEQZ
073: *          for details.
074: *
075: *  LDA     (input) INTEGER
076: *          The leading dimension of A.  LDA >= max(1,N).
077: *
078: *  B       (input/output) COMPLEX*16 array, dimension (LDB, N)
079: *          On entry, the matrix B.
080: *          If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
081: *          upper triangular matrix obtained from B in the generalized
082: *          Schur factorization of the pair (A,B) after balancing.
083: *          If no eigenvectors were computed, then only the diagonal
084: *          elements of B will be correct.  See ZGGHRD and ZHGEQZ for
085: *          details.
086: *
087: *  LDB     (input) INTEGER
088: *          The leading dimension of B.  LDB >= max(1,N).
089: *
090: *  ALPHA   (output) COMPLEX*16 array, dimension (N)
091: *          The complex scalars alpha that define the eigenvalues of
092: *          GNEP.
093: *
094: *  BETA    (output) COMPLEX*16 array, dimension (N)
095: *          The complex scalars beta that define the eigenvalues of GNEP.
096: *          
097: *          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
098: *          represent the j-th eigenvalue of the matrix pair (A,B), in
099: *          one of the forms lambda = alpha/beta or mu = beta/alpha.
100: *          Since either lambda or mu may overflow, they should not,
101: *          in general, be computed.
102: *
103: *  VL      (output) COMPLEX*16 array, dimension (LDVL,N)
104: *          If JOBVL = 'V', the left eigenvectors u(j) are stored
105: *          in the columns of VL, in the same order as their eigenvalues.
106: *          Each eigenvector is scaled so that its largest component has
107: *          abs(real part) + abs(imag. part) = 1, except for eigenvectors
108: *          corresponding to an eigenvalue with alpha = beta = 0, which
109: *          are set to zero.
110: *          Not referenced if JOBVL = 'N'.
111: *
112: *  LDVL    (input) INTEGER
113: *          The leading dimension of the matrix VL. LDVL >= 1, and
114: *          if JOBVL = 'V', LDVL >= N.
115: *
116: *  VR      (output) COMPLEX*16 array, dimension (LDVR,N)
117: *          If JOBVR = 'V', the right eigenvectors x(j) are stored
118: *          in the columns of VR, in the same order as their eigenvalues.
119: *          Each eigenvector is scaled so that its largest component has
120: *          abs(real part) + abs(imag. part) = 1, except for eigenvectors
121: *          corresponding to an eigenvalue with alpha = beta = 0, which
122: *          are set to zero.
123: *          Not referenced if JOBVR = 'N'.
124: *
125: *  LDVR    (input) INTEGER
126: *          The leading dimension of the matrix VR. LDVR >= 1, and
127: *          if JOBVR = 'V', LDVR >= N.
128: *
129: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
130: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
131: *
132: *  LWORK   (input) INTEGER
133: *          The dimension of the array WORK.  LWORK >= max(1,2*N).
134: *          For good performance, LWORK must generally be larger.
135: *          To compute the optimal value of LWORK, call ILAENV to get
136: *          blocksizes (for ZGEQRF, ZUNMQR, and ZUNGQR.)  Then compute:
137: *          NB  -- MAX of the blocksizes for ZGEQRF, ZUNMQR, and ZUNGQR;
138: *          The optimal LWORK is  MAX( 2*N, N*(NB+1) ).
139: *
140: *          If LWORK = -1, then a workspace query is assumed; the routine
141: *          only calculates the optimal size of the WORK array, returns
142: *          this value as the first entry of the WORK array, and no error
143: *          message related to LWORK is issued by XERBLA.
144: *
145: *  RWORK   (workspace/output) DOUBLE PRECISION array, dimension (8*N)
146: *
147: *  INFO    (output) INTEGER
148: *          = 0:  successful exit
149: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
150: *          =1,...,N:
151: *                The QZ iteration failed.  No eigenvectors have been
152: *                calculated, but ALPHA(j) and BETA(j) should be
153: *                correct for j=INFO+1,...,N.
154: *          > N:  errors that usually indicate LAPACK problems:
155: *                =N+1: error return from ZGGBAL
156: *                =N+2: error return from ZGEQRF
157: *                =N+3: error return from ZUNMQR
158: *                =N+4: error return from ZUNGQR
159: *                =N+5: error return from ZGGHRD
160: *                =N+6: error return from ZHGEQZ (other than failed
161: *                                               iteration)
162: *                =N+7: error return from ZTGEVC
163: *                =N+8: error return from ZGGBAK (computing VL)
164: *                =N+9: error return from ZGGBAK (computing VR)
165: *                =N+10: error return from ZLASCL (various calls)
166: *
167: *  Further Details
168: *  ===============
169: *
170: *  Balancing
171: *  ---------
172: *
173: *  This driver calls ZGGBAL to both permute and scale rows and columns
174: *  of A and B.  The permutations PL and PR are chosen so that PL*A*PR
175: *  and PL*B*R will be upper triangular except for the diagonal blocks
176: *  A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
177: *  possible.  The diagonal scaling matrices DL and DR are chosen so
178: *  that the pair  DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
179: *  one (except for the elements that start out zero.)
180: *
181: *  After the eigenvalues and eigenvectors of the balanced matrices
182: *  have been computed, ZGGBAK transforms the eigenvectors back to what
183: *  they would have been (in perfect arithmetic) if they had not been
184: *  balanced.
185: *
186: *  Contents of A and B on Exit
187: *  -------- -- - --- - -- ----
188: *
189: *  If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
190: *  both), then on exit the arrays A and B will contain the complex Schur
191: *  form[*] of the "balanced" versions of A and B.  If no eigenvectors
192: *  are computed, then only the diagonal blocks will be correct.
193: *
194: *  [*] In other words, upper triangular form.
195: *
196: *  =====================================================================
197: *
198: *     .. Parameters ..
199:       DOUBLE PRECISION   ZERO, ONE
200:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
201:       COMPLEX*16         CZERO, CONE
202:       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
203:      $                   CONE = ( 1.0D0, 0.0D0 ) )
204: *     ..
205: *     .. Local Scalars ..
206:       LOGICAL            ILIMIT, ILV, ILVL, ILVR, LQUERY
207:       CHARACTER          CHTEMP
208:       INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
209:      $                   IN, IRIGHT, IROWS, IRWORK, ITAU, IWORK, JC, JR,
210:      $                   LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3
211:       DOUBLE PRECISION   ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
212:      $                   BNRM1, BNRM2, EPS, SAFMAX, SAFMIN, SALFAI,
213:      $                   SALFAR, SBETA, SCALE, TEMP
214:       COMPLEX*16         X
215: *     ..
216: *     .. Local Arrays ..
217:       LOGICAL            LDUMMA( 1 )
218: *     ..
219: *     .. External Subroutines ..
220:       EXTERNAL           XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD, ZHGEQZ,
221:      $                   ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR, ZUNMQR
222: *     ..
223: *     .. External Functions ..
224:       LOGICAL            LSAME
225:       INTEGER            ILAENV
226:       DOUBLE PRECISION   DLAMCH, ZLANGE
227:       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
228: *     ..
229: *     .. Intrinsic Functions ..
230:       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, INT, MAX
231: *     ..
232: *     .. Statement Functions ..
233:       DOUBLE PRECISION   ABS1
234: *     ..
235: *     .. Statement Function definitions ..
236:       ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
237: *     ..
238: *     .. Executable Statements ..
239: *
240: *     Decode the input arguments
241: *
242:       IF( LSAME( JOBVL, 'N' ) ) THEN
243:          IJOBVL = 1
244:          ILVL = .FALSE.
245:       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
246:          IJOBVL = 2
247:          ILVL = .TRUE.
248:       ELSE
249:          IJOBVL = -1
250:          ILVL = .FALSE.
251:       END IF
252: *
253:       IF( LSAME( JOBVR, 'N' ) ) THEN
254:          IJOBVR = 1
255:          ILVR = .FALSE.
256:       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
257:          IJOBVR = 2
258:          ILVR = .TRUE.
259:       ELSE
260:          IJOBVR = -1
261:          ILVR = .FALSE.
262:       END IF
263:       ILV = ILVL .OR. ILVR
264: *
265: *     Test the input arguments
266: *
267:       LWKMIN = MAX( 2*N, 1 )
268:       LWKOPT = LWKMIN
269:       WORK( 1 ) = LWKOPT
270:       LQUERY = ( LWORK.EQ.-1 )
271:       INFO = 0
272:       IF( IJOBVL.LE.0 ) THEN
273:          INFO = -1
274:       ELSE IF( IJOBVR.LE.0 ) THEN
275:          INFO = -2
276:       ELSE IF( N.LT.0 ) THEN
277:          INFO = -3
278:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
279:          INFO = -5
280:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
281:          INFO = -7
282:       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
283:          INFO = -11
284:       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
285:          INFO = -13
286:       ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
287:          INFO = -15
288:       END IF
289: *
290:       IF( INFO.EQ.0 ) THEN
291:          NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, N, -1, -1 )
292:          NB2 = ILAENV( 1, 'ZUNMQR', ' ', N, N, N, -1 )
293:          NB3 = ILAENV( 1, 'ZUNGQR', ' ', N, N, N, -1 )
294:          NB = MAX( NB1, NB2, NB3 )
295:          LOPT = MAX( 2*N, N*( NB+1 ) )
296:          WORK( 1 ) = LOPT
297:       END IF
298: *
299:       IF( INFO.NE.0 ) THEN
300:          CALL XERBLA( 'ZGEGV ', -INFO )
301:          RETURN
302:       ELSE IF( LQUERY ) THEN
303:          RETURN
304:       END IF
305: *
306: *     Quick return if possible
307: *
308:       IF( N.EQ.0 )
309:      $   RETURN
310: *
311: *     Get machine constants
312: *
313:       EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
314:       SAFMIN = DLAMCH( 'S' )
315:       SAFMIN = SAFMIN + SAFMIN
316:       SAFMAX = ONE / SAFMIN
317: *
318: *     Scale A
319: *
320:       ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
321:       ANRM1 = ANRM
322:       ANRM2 = ONE
323:       IF( ANRM.LT.ONE ) THEN
324:          IF( SAFMAX*ANRM.LT.ONE ) THEN
325:             ANRM1 = SAFMIN
326:             ANRM2 = SAFMAX*ANRM
327:          END IF
328:       END IF
329: *
330:       IF( ANRM.GT.ZERO ) THEN
331:          CALL ZLASCL( 'G', -1, -1, ANRM, ONE, N, N, A, LDA, IINFO )
332:          IF( IINFO.NE.0 ) THEN
333:             INFO = N + 10
334:             RETURN
335:          END IF
336:       END IF
337: *
338: *     Scale B
339: *
340:       BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
341:       BNRM1 = BNRM
342:       BNRM2 = ONE
343:       IF( BNRM.LT.ONE ) THEN
344:          IF( SAFMAX*BNRM.LT.ONE ) THEN
345:             BNRM1 = SAFMIN
346:             BNRM2 = SAFMAX*BNRM
347:          END IF
348:       END IF
349: *
350:       IF( BNRM.GT.ZERO ) THEN
351:          CALL ZLASCL( 'G', -1, -1, BNRM, ONE, N, N, B, LDB, IINFO )
352:          IF( IINFO.NE.0 ) THEN
353:             INFO = N + 10
354:             RETURN
355:          END IF
356:       END IF
357: *
358: *     Permute the matrix to make it more nearly triangular
359: *     Also "balance" the matrix.
360: *
361:       ILEFT = 1
362:       IRIGHT = N + 1
363:       IRWORK = IRIGHT + N
364:       CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
365:      $             RWORK( IRIGHT ), RWORK( IRWORK ), IINFO )
366:       IF( IINFO.NE.0 ) THEN
367:          INFO = N + 1
368:          GO TO 80
369:       END IF
370: *
371: *     Reduce B to triangular form, and initialize VL and/or VR
372: *
373:       IROWS = IHI + 1 - ILO
374:       IF( ILV ) THEN
375:          ICOLS = N + 1 - ILO
376:       ELSE
377:          ICOLS = IROWS
378:       END IF
379:       ITAU = 1
380:       IWORK = ITAU + IROWS
381:       CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
382:      $             WORK( IWORK ), LWORK+1-IWORK, IINFO )
383:       IF( IINFO.GE.0 )
384:      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
385:       IF( IINFO.NE.0 ) THEN
386:          INFO = N + 2
387:          GO TO 80
388:       END IF
389: *
390:       CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
391:      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
392:      $             LWORK+1-IWORK, IINFO )
393:       IF( IINFO.GE.0 )
394:      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
395:       IF( IINFO.NE.0 ) THEN
396:          INFO = N + 3
397:          GO TO 80
398:       END IF
399: *
400:       IF( ILVL ) THEN
401:          CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
402:          CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
403:      $                VL( ILO+1, ILO ), LDVL )
404:          CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
405:      $                WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
406:      $                IINFO )
407:          IF( IINFO.GE.0 )
408:      $      LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
409:          IF( IINFO.NE.0 ) THEN
410:             INFO = N + 4
411:             GO TO 80
412:          END IF
413:       END IF
414: *
415:       IF( ILVR )
416:      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
417: *
418: *     Reduce to generalized Hessenberg form
419: *
420:       IF( ILV ) THEN
421: *
422: *        Eigenvectors requested -- work on whole matrix.
423: *
424:          CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
425:      $                LDVL, VR, LDVR, IINFO )
426:       ELSE
427:          CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
428:      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IINFO )
429:       END IF
430:       IF( IINFO.NE.0 ) THEN
431:          INFO = N + 5
432:          GO TO 80
433:       END IF
434: *
435: *     Perform QZ algorithm
436: *
437:       IWORK = ITAU
438:       IF( ILV ) THEN
439:          CHTEMP = 'S'
440:       ELSE
441:          CHTEMP = 'E'
442:       END IF
443:       CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
444:      $             ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWORK ),
445:      $             LWORK+1-IWORK, RWORK( IRWORK ), IINFO )
446:       IF( IINFO.GE.0 )
447:      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
448:       IF( IINFO.NE.0 ) THEN
449:          IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
450:             INFO = IINFO
451:          ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
452:             INFO = IINFO - N
453:          ELSE
454:             INFO = N + 6
455:          END IF
456:          GO TO 80
457:       END IF
458: *
459:       IF( ILV ) THEN
460: *
461: *        Compute Eigenvectors
462: *
463:          IF( ILVL ) THEN
464:             IF( ILVR ) THEN
465:                CHTEMP = 'B'
466:             ELSE
467:                CHTEMP = 'L'
468:             END IF
469:          ELSE
470:             CHTEMP = 'R'
471:          END IF
472: *
473:          CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
474:      $                VR, LDVR, N, IN, WORK( IWORK ), RWORK( IRWORK ),
475:      $                IINFO )
476:          IF( IINFO.NE.0 ) THEN
477:             INFO = N + 7
478:             GO TO 80
479:          END IF
480: *
481: *        Undo balancing on VL and VR, rescale
482: *
483:          IF( ILVL ) THEN
484:             CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
485:      $                   RWORK( IRIGHT ), N, VL, LDVL, IINFO )
486:             IF( IINFO.NE.0 ) THEN
487:                INFO = N + 8
488:                GO TO 80
489:             END IF
490:             DO 30 JC = 1, N
491:                TEMP = ZERO
492:                DO 10 JR = 1, N
493:                   TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
494:    10          CONTINUE
495:                IF( TEMP.LT.SAFMIN )
496:      $            GO TO 30
497:                TEMP = ONE / TEMP
498:                DO 20 JR = 1, N
499:                   VL( JR, JC ) = VL( JR, JC )*TEMP
500:    20          CONTINUE
501:    30       CONTINUE
502:          END IF
503:          IF( ILVR ) THEN
504:             CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
505:      $                   RWORK( IRIGHT ), N, VR, LDVR, IINFO )
506:             IF( IINFO.NE.0 ) THEN
507:                INFO = N + 9
508:                GO TO 80
509:             END IF
510:             DO 60 JC = 1, N
511:                TEMP = ZERO
512:                DO 40 JR = 1, N
513:                   TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
514:    40          CONTINUE
515:                IF( TEMP.LT.SAFMIN )
516:      $            GO TO 60
517:                TEMP = ONE / TEMP
518:                DO 50 JR = 1, N
519:                   VR( JR, JC ) = VR( JR, JC )*TEMP
520:    50          CONTINUE
521:    60       CONTINUE
522:          END IF
523: *
524: *        End of eigenvector calculation
525: *
526:       END IF
527: *
528: *     Undo scaling in alpha, beta
529: *
530: *     Note: this does not give the alpha and beta for the unscaled
531: *     problem.
532: *
533: *     Un-scaling is limited to avoid underflow in alpha and beta
534: *     if they are significant.
535: *
536:       DO 70 JC = 1, N
537:          ABSAR = ABS( DBLE( ALPHA( JC ) ) )
538:          ABSAI = ABS( DIMAG( ALPHA( JC ) ) )
539:          ABSB = ABS( DBLE( BETA( JC ) ) )
540:          SALFAR = ANRM*DBLE( ALPHA( JC ) )
541:          SALFAI = ANRM*DIMAG( ALPHA( JC ) )
542:          SBETA = BNRM*DBLE( BETA( JC ) )
543:          ILIMIT = .FALSE.
544:          SCALE = ONE
545: *
546: *        Check for significant underflow in imaginary part of ALPHA
547: *
548:          IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
549:      $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
550:             ILIMIT = .TRUE.
551:             SCALE = ( SAFMIN / ANRM1 ) / MAX( SAFMIN, ANRM2*ABSAI )
552:          END IF
553: *
554: *        Check for significant underflow in real part of ALPHA
555: *
556:          IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
557:      $       MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
558:             ILIMIT = .TRUE.
559:             SCALE = MAX( SCALE, ( SAFMIN / ANRM1 ) /
560:      $              MAX( SAFMIN, ANRM2*ABSAR ) )
561:          END IF
562: *
563: *        Check for significant underflow in BETA
564: *
565:          IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE.
566:      $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN
567:             ILIMIT = .TRUE.
568:             SCALE = MAX( SCALE, ( SAFMIN / BNRM1 ) /
569:      $              MAX( SAFMIN, BNRM2*ABSB ) )
570:          END IF
571: *
572: *        Check for possible overflow when limiting scaling
573: *
574:          IF( ILIMIT ) THEN
575:             TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ),
576:      $             ABS( SBETA ) )
577:             IF( TEMP.GT.ONE )
578:      $         SCALE = SCALE / TEMP
579:             IF( SCALE.LT.ONE )
580:      $         ILIMIT = .FALSE.
581:          END IF
582: *
583: *        Recompute un-scaled ALPHA, BETA if necessary.
584: *
585:          IF( ILIMIT ) THEN
586:             SALFAR = ( SCALE*DBLE( ALPHA( JC ) ) )*ANRM
587:             SALFAI = ( SCALE*DIMAG( ALPHA( JC ) ) )*ANRM
588:             SBETA = ( SCALE*BETA( JC ) )*BNRM
589:          END IF
590:          ALPHA( JC ) = DCMPLX( SALFAR, SALFAI )
591:          BETA( JC ) = SBETA
592:    70 CONTINUE
593: *
594:    80 CONTINUE
595:       WORK( 1 ) = LWKOPT
596: *
597:       RETURN
598: *
599: *     End of ZGEGV
600: *
601:       END
602: