LAPACK: SRC/cgegv.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
00002      $                  VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
00003 *
00004 *  -- LAPACK driver routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          JOBVL, JOBVR
00011       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               RWORK( * )
00015       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
00016      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
00017      $                   WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  This routine is deprecated and has been replaced by routine CGGEV.
00024 *
00025 *  CGEGV computes the eigenvalues and, optionally, the left and/or right
00026 *  eigenvectors of a complex matrix pair (A,B).
00027 *  Given two square matrices A and B,
00028 *  the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
00029 *  eigenvalues lambda and corresponding (non-zero) eigenvectors x such
00030 *  that
00031 *     A*x = lambda*B*x.
00032 *
00033 *  An alternate form is to find the eigenvalues mu and corresponding
00034 *  eigenvectors y such that
00035 *     mu*A*y = B*y.
00036 *
00037 *  These two forms are equivalent with mu = 1/lambda and x = y if
00038 *  neither lambda nor mu is zero.  In order to deal with the case that
00039 *  lambda or mu is zero or small, two values alpha and beta are returned
00040 *  for each eigenvalue, such that lambda = alpha/beta and
00041 *  mu = beta/alpha.
00042 *  
00043 *  The vectors x and y in the above equations are right eigenvectors of
00044 *  the matrix pair (A,B).  Vectors u and v satisfying
00045 *     u**H*A = lambda*u**H*B  or  mu*v**H*A = v**H*B
00046 *  are left eigenvectors of (A,B).
00047 *
00048 *  Note: this routine performs "full balancing" on A and B -- see
00049 *  "Further Details", below.
00050 *
00051 *  Arguments
00052 *  =========
00053 *
00054 *  JOBVL   (input) CHARACTER*1
00055 *          = 'N':  do not compute the left generalized eigenvectors;
00056 *          = 'V':  compute the left generalized eigenvectors (returned
00057 *                  in VL).
00058 *
00059 *  JOBVR   (input) CHARACTER*1
00060 *          = 'N':  do not compute the right generalized eigenvectors;
00061 *          = 'V':  compute the right generalized eigenvectors (returned
00062 *                  in VR).
00063 *
00064 *  N       (input) INTEGER
00065 *          The order of the matrices A, B, VL, and VR.  N >= 0.
00066 *
00067 *  A       (input/output) COMPLEX array, dimension (LDA, N)
00068 *          On entry, the matrix A.
00069 *          If JOBVL = 'V' or JOBVR = 'V', then on exit A
00070 *          contains the Schur form of A from the generalized Schur
00071 *          factorization of the pair (A,B) after balancing.  If no
00072 *          eigenvectors were computed, then only the diagonal elements
00073 *          of the Schur form will be correct.  See CGGHRD and CHGEQZ
00074 *          for details.
00075 *
00076 *  LDA     (input) INTEGER
00077 *          The leading dimension of A.  LDA >= max(1,N).
00078 *
00079 *  B       (input/output) COMPLEX array, dimension (LDB, N)
00080 *          On entry, the matrix B.
00081 *          If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
00082 *          upper triangular matrix obtained from B in the generalized
00083 *          Schur factorization of the pair (A,B) after balancing.
00084 *          If no eigenvectors were computed, then only the diagonal
00085 *          elements of B will be correct.  See CGGHRD and CHGEQZ for
00086 *          details.
00087 *
00088 *  LDB     (input) INTEGER
00089 *          The leading dimension of B.  LDB >= max(1,N).
00090 *
00091 *  ALPHA   (output) COMPLEX array, dimension (N)
00092 *          The complex scalars alpha that define the eigenvalues of
00093 *          GNEP.
00094 *
00095 *  BETA    (output) COMPLEX array, dimension (N)
00096 *          The complex scalars beta that define the eigenvalues of GNEP.
00097 *          
00098 *          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
00099 *          represent the j-th eigenvalue of the matrix pair (A,B), in
00100 *          one of the forms lambda = alpha/beta or mu = beta/alpha.
00101 *          Since either lambda or mu may overflow, they should not,
00102 *          in general, be computed.
00103 
00104 *
00105 *  VL      (output) COMPLEX array, dimension (LDVL,N)
00106 *          If JOBVL = 'V', the left eigenvectors u(j) are stored
00107 *          in the columns of VL, in the same order as their eigenvalues.
00108 *          Each eigenvector is scaled so that its largest component has
00109 *          abs(real part) + abs(imag. part) = 1, except for eigenvectors
00110 *          corresponding to an eigenvalue with alpha = beta = 0, which
00111 *          are set to zero.
00112 *          Not referenced if JOBVL = 'N'.
00113 *
00114 *  LDVL    (input) INTEGER
00115 *          The leading dimension of the matrix VL. LDVL >= 1, and
00116 *          if JOBVL = 'V', LDVL >= N.
00117 *
00118 *  VR      (output) COMPLEX array, dimension (LDVR,N)
00119 *          If JOBVR = 'V', the right eigenvectors x(j) are stored
00120 *          in the columns of VR, in the same order as their eigenvalues.
00121 *          Each eigenvector is scaled so that its largest component has
00122 *          abs(real part) + abs(imag. part) = 1, except for eigenvectors
00123 *          corresponding to an eigenvalue with alpha = beta = 0, which
00124 *          are set to zero.
00125 *          Not referenced if JOBVR = 'N'.
00126 *
00127 *  LDVR    (input) INTEGER
00128 *          The leading dimension of the matrix VR. LDVR >= 1, and
00129 *          if JOBVR = 'V', LDVR >= N.
00130 *
00131 *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
00132 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00133 *
00134 *  LWORK   (input) INTEGER
00135 *          The dimension of the array WORK.  LWORK >= max(1,2*N).
00136 *          For good performance, LWORK must generally be larger.
00137 *          To compute the optimal value of LWORK, call ILAENV to get
00138 *          blocksizes (for CGEQRF, CUNMQR, and CUNGQR.)  Then compute:
00139 *          NB  -- MAX of the blocksizes for CGEQRF, CUNMQR, and CUNGQR;
00140 *          The optimal LWORK is  MAX( 2*N, N*(NB+1) ).
00141 *
00142 *          If LWORK = -1, then a workspace query is assumed; the routine
00143 *          only calculates the optimal size of the WORK array, returns
00144 *          this value as the first entry of the WORK array, and no error
00145 *          message related to LWORK is issued by XERBLA.
00146 *
00147 *  RWORK   (workspace/output) REAL array, dimension (8*N)
00148 *
00149 *  INFO    (output) INTEGER
00150 *          = 0:  successful exit
00151 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00152 *          =1,...,N:
00153 *                The QZ iteration failed.  No eigenvectors have been
00154 *                calculated, but ALPHA(j) and BETA(j) should be
00155 *                correct for j=INFO+1,...,N.
00156 *          > N:  errors that usually indicate LAPACK problems:
00157 *                =N+1: error return from CGGBAL
00158 *                =N+2: error return from CGEQRF
00159 *                =N+3: error return from CUNMQR
00160 *                =N+4: error return from CUNGQR
00161 *                =N+5: error return from CGGHRD
00162 *                =N+6: error return from CHGEQZ (other than failed
00163 *                                               iteration)
00164 *                =N+7: error return from CTGEVC
00165 *                =N+8: error return from CGGBAK (computing VL)
00166 *                =N+9: error return from CGGBAK (computing VR)
00167 *                =N+10: error return from CLASCL (various calls)
00168 *
00169 *  Further Details
00170 *  ===============
00171 *
00172 *  Balancing
00173 *  ---------
00174 *
00175 *  This driver calls CGGBAL to both permute and scale rows and columns
00176 *  of A and B.  The permutations PL and PR are chosen so that PL*A*PR
00177 *  and PL*B*R will be upper triangular except for the diagonal blocks
00178 *  A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
00179 *  possible.  The diagonal scaling matrices DL and DR are chosen so
00180 *  that the pair  DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
00181 *  one (except for the elements that start out zero.)
00182 *
00183 *  After the eigenvalues and eigenvectors of the balanced matrices
00184 *  have been computed, CGGBAK transforms the eigenvectors back to what
00185 *  they would have been (in perfect arithmetic) if they had not been
00186 *  balanced.
00187 *
00188 *  Contents of A and B on Exit
00189 *  -------- -- - --- - -- ----
00190 *
00191 *  If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
00192 *  both), then on exit the arrays A and B will contain the complex Schur
00193 *  form[*] of the "balanced" versions of A and B.  If no eigenvectors
00194 *  are computed, then only the diagonal blocks will be correct.
00195 *
00196 *  [*] In other words, upper triangular form.
00197 *
00198 *  =====================================================================
00199 *
00200 *     .. Parameters ..
00201       REAL               ZERO, ONE
00202       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00203       COMPLEX            CZERO, CONE
00204       PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ),
00205      $                   CONE = ( 1.0E0, 0.0E0 ) )
00206 *     ..
00207 *     .. Local Scalars ..
00208       LOGICAL            ILIMIT, ILV, ILVL, ILVR, LQUERY
00209       CHARACTER          CHTEMP
00210       INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
00211      $                   IN, IRIGHT, IROWS, IRWORK, ITAU, IWORK, JC, JR,
00212      $                   LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3
00213       REAL               ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
00214      $                   BNRM1, BNRM2, EPS, SAFMAX, SAFMIN, SALFAI,
00215      $                   SALFAR, SBETA, SCALE, TEMP
00216       COMPLEX            X
00217 *     ..
00218 *     .. Local Arrays ..
00219       LOGICAL            LDUMMA( 1 )
00220 *     ..
00221 *     .. External Subroutines ..
00222       EXTERNAL           CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY,
00223      $                   CLASCL, CLASET, CTGEVC, CUNGQR, CUNMQR, XERBLA
00224 *     ..
00225 *     .. External Functions ..
00226       LOGICAL            LSAME
00227       INTEGER            ILAENV
00228       REAL               CLANGE, SLAMCH
00229       EXTERNAL           ILAENV, LSAME, CLANGE, SLAMCH
00230 *     ..
00231 *     .. Intrinsic Functions ..
00232       INTRINSIC          ABS, AIMAG, CMPLX, INT, MAX, REAL
00233 *     ..
00234 *     .. Statement Functions ..
00235       REAL               ABS1
00236 *     ..
00237 *     .. Statement Function definitions ..
00238       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
00239 *     ..
00240 *     .. Executable Statements ..
00241 *
00242 *     Decode the input arguments
00243 *
00244       IF( LSAME( JOBVL, 'N' ) ) THEN
00245          IJOBVL = 1
00246          ILVL = .FALSE.
00247       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
00248          IJOBVL = 2
00249          ILVL = .TRUE.
00250       ELSE
00251          IJOBVL = -1
00252          ILVL = .FALSE.
00253       END IF
00254 *
00255       IF( LSAME( JOBVR, 'N' ) ) THEN
00256          IJOBVR = 1
00257          ILVR = .FALSE.
00258       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
00259          IJOBVR = 2
00260          ILVR = .TRUE.
00261       ELSE
00262          IJOBVR = -1
00263          ILVR = .FALSE.
00264       END IF
00265       ILV = ILVL .OR. ILVR
00266 *
00267 *     Test the input arguments
00268 *
00269       LWKMIN = MAX( 2*N, 1 )
00270       LWKOPT = LWKMIN
00271       WORK( 1 ) = LWKOPT
00272       LQUERY = ( LWORK.EQ.-1 )
00273       INFO = 0
00274       IF( IJOBVL.LE.0 ) THEN
00275          INFO = -1
00276       ELSE IF( IJOBVR.LE.0 ) THEN
00277          INFO = -2
00278       ELSE IF( N.LT.0 ) THEN
00279          INFO = -3
00280       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00281          INFO = -5
00282       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00283          INFO = -7
00284       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
00285          INFO = -11
00286       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
00287          INFO = -13
00288       ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
00289          INFO = -15
00290       END IF
00291 *
00292       IF( INFO.EQ.0 ) THEN
00293          NB1 = ILAENV( 1, 'CGEQRF', ' ', N, N, -1, -1 )
00294          NB2 = ILAENV( 1, 'CUNMQR', ' ', N, N, N, -1 )
00295          NB3 = ILAENV( 1, 'CUNGQR', ' ', N, N, N, -1 )
00296          NB = MAX( NB1, NB2, NB3 )
00297          LOPT = MAX( 2*N, N*(NB+1) )
00298          WORK( 1 ) = LOPT
00299       END IF
00300 *
00301       IF( INFO.NE.0 ) THEN
00302          CALL XERBLA( 'CGEGV ', -INFO )
00303          RETURN
00304       ELSE IF( LQUERY ) THEN
00305          RETURN
00306       END IF
00307 *
00308 *     Quick return if possible
00309 *
00310       IF( N.EQ.0 )
00311      $   RETURN
00312 *
00313 *     Get machine constants
00314 *
00315       EPS = SLAMCH( 'E' )*SLAMCH( 'B' )
00316       SAFMIN = SLAMCH( 'S' )
00317       SAFMIN = SAFMIN + SAFMIN
00318       SAFMAX = ONE / SAFMIN
00319 *
00320 *     Scale A
00321 *
00322       ANRM = CLANGE( 'M', N, N, A, LDA, RWORK )
00323       ANRM1 = ANRM
00324       ANRM2 = ONE
00325       IF( ANRM.LT.ONE ) THEN
00326          IF( SAFMAX*ANRM.LT.ONE ) THEN
00327             ANRM1 = SAFMIN
00328             ANRM2 = SAFMAX*ANRM
00329          END IF
00330       END IF
00331 *
00332       IF( ANRM.GT.ZERO ) THEN
00333          CALL CLASCL( 'G', -1, -1, ANRM, ONE, N, N, A, LDA, IINFO )
00334          IF( IINFO.NE.0 ) THEN
00335             INFO = N + 10
00336             RETURN
00337          END IF
00338       END IF
00339 *
00340 *     Scale B
00341 *
00342       BNRM = CLANGE( 'M', N, N, B, LDB, RWORK )
00343       BNRM1 = BNRM
00344       BNRM2 = ONE
00345       IF( BNRM.LT.ONE ) THEN
00346          IF( SAFMAX*BNRM.LT.ONE ) THEN
00347             BNRM1 = SAFMIN
00348             BNRM2 = SAFMAX*BNRM
00349          END IF
00350       END IF
00351 *
00352       IF( BNRM.GT.ZERO ) THEN
00353          CALL CLASCL( 'G', -1, -1, BNRM, ONE, N, N, B, LDB, IINFO )
00354          IF( IINFO.NE.0 ) THEN
00355             INFO = N + 10
00356             RETURN
00357          END IF
00358       END IF
00359 *
00360 *     Permute the matrix to make it more nearly triangular
00361 *     Also "balance" the matrix.
00362 *
00363       ILEFT = 1
00364       IRIGHT = N + 1
00365       IRWORK = IRIGHT + N
00366       CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
00367      $             RWORK( IRIGHT ), RWORK( IRWORK ), IINFO )
00368       IF( IINFO.NE.0 ) THEN
00369          INFO = N + 1
00370          GO TO 80
00371       END IF
00372 *
00373 *     Reduce B to triangular form, and initialize VL and/or VR
00374 *
00375       IROWS = IHI + 1 - ILO
00376       IF( ILV ) THEN
00377          ICOLS = N + 1 - ILO
00378       ELSE
00379          ICOLS = IROWS
00380       END IF
00381       ITAU = 1
00382       IWORK = ITAU + IROWS
00383       CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00384      $             WORK( IWORK ), LWORK+1-IWORK, IINFO )
00385       IF( IINFO.GE.0 )
00386      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00387       IF( IINFO.NE.0 ) THEN
00388          INFO = N + 2
00389          GO TO 80
00390       END IF
00391 *
00392       CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00393      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
00394      $             LWORK+1-IWORK, IINFO )
00395       IF( IINFO.GE.0 )
00396      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00397       IF( IINFO.NE.0 ) THEN
00398          INFO = N + 3
00399          GO TO 80
00400       END IF
00401 *
00402       IF( ILVL ) THEN
00403          CALL CLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
00404          CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00405      $                VL( ILO+1, ILO ), LDVL )
00406          CALL CUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
00407      $                WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
00408      $                IINFO )
00409          IF( IINFO.GE.0 )
00410      $      LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00411          IF( IINFO.NE.0 ) THEN
00412             INFO = N + 4
00413             GO TO 80
00414          END IF
00415       END IF
00416 *
00417       IF( ILVR )
00418      $   CALL CLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
00419 *
00420 *     Reduce to generalized Hessenberg form
00421 *
00422       IF( ILV ) THEN
00423 *
00424 *        Eigenvectors requested -- work on whole matrix.
00425 *
00426          CALL CGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
00427      $                LDVL, VR, LDVR, IINFO )
00428       ELSE
00429          CALL CGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
00430      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IINFO )
00431       END IF
00432       IF( IINFO.NE.0 ) THEN
00433          INFO = N + 5
00434          GO TO 80
00435       END IF
00436 *
00437 *     Perform QZ algorithm
00438 *
00439       IWORK = ITAU
00440       IF( ILV ) THEN
00441          CHTEMP = 'S'
00442       ELSE
00443          CHTEMP = 'E'
00444       END IF
00445       CALL CHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
00446      $             ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWORK ),
00447      $             LWORK+1-IWORK, RWORK( IRWORK ), IINFO )
00448       IF( IINFO.GE.0 )
00449      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00450       IF( IINFO.NE.0 ) THEN
00451          IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
00452             INFO = IINFO
00453          ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
00454             INFO = IINFO - N
00455          ELSE
00456             INFO = N + 6
00457          END IF
00458          GO TO 80
00459       END IF
00460 *
00461       IF( ILV ) THEN
00462 *
00463 *        Compute Eigenvectors
00464 *
00465          IF( ILVL ) THEN
00466             IF( ILVR ) THEN
00467                CHTEMP = 'B'
00468             ELSE
00469                CHTEMP = 'L'
00470             END IF
00471          ELSE
00472             CHTEMP = 'R'
00473          END IF
00474 *
00475          CALL CTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
00476      $                VR, LDVR, N, IN, WORK( IWORK ), RWORK( IRWORK ),
00477      $                IINFO )
00478          IF( IINFO.NE.0 ) THEN
00479             INFO = N + 7
00480             GO TO 80
00481          END IF
00482 *
00483 *        Undo balancing on VL and VR, rescale
00484 *
00485          IF( ILVL ) THEN
00486             CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
00487      $                   RWORK( IRIGHT ), N, VL, LDVL, IINFO )
00488             IF( IINFO.NE.0 ) THEN
00489                INFO = N + 8
00490                GO TO 80
00491             END IF
00492             DO 30 JC = 1, N
00493                TEMP = ZERO
00494                DO 10 JR = 1, N
00495                   TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
00496    10          CONTINUE
00497                IF( TEMP.LT.SAFMIN )
00498      $            GO TO 30
00499                TEMP = ONE / TEMP
00500                DO 20 JR = 1, N
00501                   VL( JR, JC ) = VL( JR, JC )*TEMP
00502    20          CONTINUE
00503    30       CONTINUE
00504          END IF
00505          IF( ILVR ) THEN
00506             CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
00507      $                   RWORK( IRIGHT ), N, VR, LDVR, IINFO )
00508             IF( IINFO.NE.0 ) THEN
00509                INFO = N + 9
00510                GO TO 80
00511             END IF
00512             DO 60 JC = 1, N
00513                TEMP = ZERO
00514                DO 40 JR = 1, N
00515                   TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
00516    40          CONTINUE
00517                IF( TEMP.LT.SAFMIN )
00518      $            GO TO 60
00519                TEMP = ONE / TEMP
00520                DO 50 JR = 1, N
00521                   VR( JR, JC ) = VR( JR, JC )*TEMP
00522    50          CONTINUE
00523    60       CONTINUE
00524          END IF
00525 *
00526 *        End of eigenvector calculation
00527 *
00528       END IF
00529 *
00530 *     Undo scaling in alpha, beta
00531 *
00532 *     Note: this does not give the alpha and beta for the unscaled
00533 *     problem.
00534 *
00535 *     Un-scaling is limited to avoid underflow in alpha and beta
00536 *     if they are significant.
00537 *
00538       DO 70 JC = 1, N
00539          ABSAR = ABS( REAL( ALPHA( JC ) ) )
00540          ABSAI = ABS( AIMAG( ALPHA( JC ) ) )
00541          ABSB = ABS( REAL( BETA( JC ) ) )
00542          SALFAR = ANRM*REAL( ALPHA( JC ) )
00543          SALFAI = ANRM*AIMAG( ALPHA( JC ) )
00544          SBETA = BNRM*REAL( BETA( JC ) )
00545          ILIMIT = .FALSE.
00546          SCALE = ONE
00547 *
00548 *        Check for significant underflow in imaginary part of ALPHA
00549 *
00550          IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
00551      $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
00552             ILIMIT = .TRUE.
00553             SCALE = ( SAFMIN / ANRM1 ) / MAX( SAFMIN, ANRM2*ABSAI )
00554          END IF
00555 *
00556 *        Check for significant underflow in real part of ALPHA
00557 *
00558          IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
00559      $       MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
00560             ILIMIT = .TRUE.
00561             SCALE = MAX( SCALE, ( SAFMIN / ANRM1 ) /
00562      $              MAX( SAFMIN, ANRM2*ABSAR ) )
00563          END IF
00564 *
00565 *        Check for significant underflow in BETA
00566 *
00567          IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE.
00568      $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN
00569             ILIMIT = .TRUE.
00570             SCALE = MAX( SCALE, ( SAFMIN / BNRM1 ) /
00571      $              MAX( SAFMIN, BNRM2*ABSB ) )
00572          END IF
00573 *
00574 *        Check for possible overflow when limiting scaling
00575 *
00576          IF( ILIMIT ) THEN
00577             TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ),
00578      $             ABS( SBETA ) )
00579             IF( TEMP.GT.ONE )
00580      $         SCALE = SCALE / TEMP
00581             IF( SCALE.LT.ONE )
00582      $         ILIMIT = .FALSE.
00583          END IF
00584 *
00585 *        Recompute un-scaled ALPHA, BETA if necessary.
00586 *
00587          IF( ILIMIT ) THEN
00588             SALFAR = ( SCALE*REAL( ALPHA( JC ) ) )*ANRM
00589             SALFAI = ( SCALE*AIMAG( ALPHA( JC ) ) )*ANRM
00590             SBETA = ( SCALE*BETA( JC ) )*BNRM
00591          END IF
00592          ALPHA( JC ) = CMPLX( SALFAR, SALFAI )
00593          BETA( JC ) = SBETA
00594    70 CONTINUE
00595 *
00596    80 CONTINUE
00597       WORK( 1 ) = LWKOPT
00598 *
00599       RETURN
00600 *
00601 *     End of CGEGV
00602 *
00603       END