LAPACK: SRC/zgegv.f Source File

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