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