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