SRC/dggev.f

   1 *> \brief <b> DGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGGEV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggev.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggev.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggev.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
  22 *                         BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBVL, JOBVR
  26 *       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
  30 *      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
  31 *      $                   VR( LDVR, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> DGGEV computes for a pair of N-by-N real nonsymmetric matrices (A,B)
  41 *> the generalized eigenvalues, and optionally, the left and/or right
  42 *> generalized eigenvectors.
  43 *>
  44 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
  45 *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
  46 *> singular. It is usually represented as the pair (alpha,beta), as
  47 *> there is a reasonable interpretation for beta=0, and even for both
  48 *> being zero.
  49 *>
  50 *> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
  51 *> of (A,B) satisfies
  52 *>
  53 *>                  A * v(j) = lambda(j) * B * v(j).
  54 *>
  55 *> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
  56 *> of (A,B) satisfies
  57 *>
  58 *>                  u(j)**H * A  = lambda(j) * u(j)**H * B .
  59 *>
  60 *> where u(j)**H is the conjugate-transpose of u(j).
  61 *>
  62 *> \endverbatim
  63 *
  64 *  Arguments:
  65 *  ==========
  66 *
  67 *> \param[in] JOBVL
  68 *> \verbatim
  69 *>          JOBVL is CHARACTER*1
  70 *>          = 'N':  do not compute the left generalized eigenvectors;
  71 *>          = 'V':  compute the left generalized eigenvectors.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] JOBVR
  75 *> \verbatim
  76 *>          JOBVR is CHARACTER*1
  77 *>          = 'N':  do not compute the right generalized eigenvectors;
  78 *>          = 'V':  compute the right generalized eigenvectors.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] N
  82 *> \verbatim
  83 *>          N is INTEGER
  84 *>          The order of the matrices A, B, VL, and VR.  N >= 0.
  85 *> \endverbatim
  86 *>
  87 *> \param[in,out] A
  88 *> \verbatim
  89 *>          A is DOUBLE PRECISION array, dimension (LDA, N)
  90 *>          On entry, the matrix A in the pair (A,B).
  91 *>          On exit, A has been overwritten.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] LDA
  95 *> \verbatim
  96 *>          LDA is INTEGER
  97 *>          The leading dimension of A.  LDA >= max(1,N).
  98 *> \endverbatim
  99 *>
 100 *> \param[in,out] B
 101 *> \verbatim
 102 *>          B is DOUBLE PRECISION array, dimension (LDB, N)
 103 *>          On entry, the matrix B in the pair (A,B).
 104 *>          On exit, B has been overwritten.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] LDB
 108 *> \verbatim
 109 *>          LDB is INTEGER
 110 *>          The leading dimension of B.  LDB >= max(1,N).
 111 *> \endverbatim
 112 *>
 113 *> \param[out] ALPHAR
 114 *> \verbatim
 115 *>          ALPHAR is DOUBLE PRECISION array, dimension (N)
 116 *> \endverbatim
 117 *>
 118 *> \param[out] ALPHAI
 119 *> \verbatim
 120 *>          ALPHAI is DOUBLE PRECISION array, dimension (N)
 121 *> \endverbatim
 122 *>
 123 *> \param[out] BETA
 124 *> \verbatim
 125 *>          BETA is DOUBLE PRECISION array, dimension (N)
 126 *>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
 127 *>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
 128 *>          the j-th eigenvalue is real; if positive, then the j-th and
 129 *>          (j+1)-st eigenvalues are a complex conjugate pair, with
 130 *>          ALPHAI(j+1) negative.
 131 *>
 132 *>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
 133 *>          may easily over- or underflow, and BETA(j) may even be zero.
 134 *>          Thus, the user should avoid naively computing the ratio
 135 *>          alpha/beta.  However, ALPHAR and ALPHAI will be always less
 136 *>          than and usually comparable with norm(A) in magnitude, and
 137 *>          BETA always less than and usually comparable with norm(B).
 138 *> \endverbatim
 139 *>
 140 *> \param[out] VL
 141 *> \verbatim
 142 *>          VL is DOUBLE PRECISION array, dimension (LDVL,N)
 143 *>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
 144 *>          after another in the columns of VL, in the same order as
 145 *>          their eigenvalues. If the j-th eigenvalue is real, then
 146 *>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
 147 *>          (j+1)-th eigenvalues form a complex conjugate pair, then
 148 *>          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
 149 *>          Each eigenvector is scaled so the largest component has
 150 *>          abs(real part)+abs(imag. part)=1.
 151 *>          Not referenced if JOBVL = 'N'.
 152 *> \endverbatim
 153 *>
 154 *> \param[in] LDVL
 155 *> \verbatim
 156 *>          LDVL is INTEGER
 157 *>          The leading dimension of the matrix VL. LDVL >= 1, and
 158 *>          if JOBVL = 'V', LDVL >= N.
 159 *> \endverbatim
 160 *>
 161 *> \param[out] VR
 162 *> \verbatim
 163 *>          VR is DOUBLE PRECISION array, dimension (LDVR,N)
 164 *>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
 165 *>          after another in the columns of VR, in the same order as
 166 *>          their eigenvalues. If the j-th eigenvalue is real, then
 167 *>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
 168 *>          (j+1)-th eigenvalues form a complex conjugate pair, then
 169 *>          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
 170 *>          Each eigenvector is scaled so the largest component has
 171 *>          abs(real part)+abs(imag. part)=1.
 172 *>          Not referenced if JOBVR = 'N'.
 173 *> \endverbatim
 174 *>
 175 *> \param[in] LDVR
 176 *> \verbatim
 177 *>          LDVR is INTEGER
 178 *>          The leading dimension of the matrix VR. LDVR >= 1, and
 179 *>          if JOBVR = 'V', LDVR >= N.
 180 *> \endverbatim
 181 *>
 182 *> \param[out] WORK
 183 *> \verbatim
 184 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 185 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 186 *> \endverbatim
 187 *>
 188 *> \param[in] LWORK
 189 *> \verbatim
 190 *>          LWORK is INTEGER
 191 *>          The dimension of the array WORK.  LWORK >= max(1,8*N).
 192 *>          For good performance, LWORK must generally be larger.
 193 *>
 194 *>          If LWORK = -1, then a workspace query is assumed; the routine
 195 *>          only calculates the optimal size of the WORK array, returns
 196 *>          this value as the first entry of the WORK array, and no error
 197 *>          message related to LWORK is issued by XERBLA.
 198 *> \endverbatim
 199 *>
 200 *> \param[out] INFO
 201 *> \verbatim
 202 *>          INFO is INTEGER
 203 *>          = 0:  successful exit
 204 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 205 *>          = 1,...,N:
 206 *>                The QZ iteration failed.  No eigenvectors have been
 207 *>                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
 208 *>                should be correct for j=INFO+1,...,N.
 209 *>          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
 210 *>                =N+2: error return from DTGEVC.
 211 *> \endverbatim
 212 *
 213 *  Authors:
 214 *  ========
 215 *
 216 *> \author Univ. of Tennessee
 217 *> \author Univ. of California Berkeley
 218 *> \author Univ. of Colorado Denver
 219 *> \author NAG Ltd.
 220 *
 221 *> \date April 2012
 222 *
 223 *> \ingroup doubleGEeigen
 224 *
 225 *  =====================================================================
 226       SUBROUTINE DGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
 227      $                  BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
 228 *
 229 *  -- LAPACK driver routine (version 3.4.1) --
 230 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 231 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 232 *     April 2012
 233 *
 234 *     .. Scalar Arguments ..
 235       CHARACTER          JOBVL, JOBVR
 236       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
 237 *     ..
 238 *     .. Array Arguments ..
 239       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
 240      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
 241      $                   VR( LDVR, * ), WORK( * )
 242 *     ..
 243 *
 244 *  =====================================================================
 245 *
 246 *     .. Parameters ..
 247       DOUBLE PRECISION   ZERO, ONE
 248       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 249 *     ..
 250 *     .. Local Scalars ..
 251       LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
 252       CHARACTER          CHTEMP
 253       INTEGER            ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
 254      $                   IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, MAXWRK,
 255      $                   MINWRK
 256       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
 257      $                   SMLNUM, TEMP
 258 *     ..
 259 *     .. Local Arrays ..
 260       LOGICAL            LDUMMA( 1 )
 261 *     ..
 262 *     .. External Subroutines ..
 263       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
 264      $                   DLACPY,DLASCL, DLASET, DORGQR, DORMQR, DTGEVC,
 265      $                   XERBLA
 266 *     ..
 267 *     .. External Functions ..
 268       LOGICAL            LSAME
 269       INTEGER            ILAENV
 270       DOUBLE PRECISION   DLAMCH, DLANGE
 271       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
 272 *     ..
 273 *     .. Intrinsic Functions ..
 274       INTRINSIC          ABS, MAX, SQRT
 275 *     ..
 276 *     .. Executable Statements ..
 277 *
 278 *     Decode the input arguments
 279 *
 280       IF( LSAME( JOBVL, 'N' ) ) THEN
 281          IJOBVL = 1
 282          ILVL = .FALSE.
 283       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
 284          IJOBVL = 2
 285          ILVL = .TRUE.
 286       ELSE
 287          IJOBVL = -1
 288          ILVL = .FALSE.
 289       END IF
 290 *
 291       IF( LSAME( JOBVR, 'N' ) ) THEN
 292          IJOBVR = 1
 293          ILVR = .FALSE.
 294       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
 295          IJOBVR = 2
 296          ILVR = .TRUE.
 297       ELSE
 298          IJOBVR = -1
 299          ILVR = .FALSE.
 300       END IF
 301       ILV = ILVL .OR. ILVR
 302 *
 303 *     Test the input arguments
 304 *
 305       INFO = 0
 306       LQUERY = ( LWORK.EQ.-1 )
 307       IF( IJOBVL.LE.0 ) THEN
 308          INFO = -1
 309       ELSE IF( IJOBVR.LE.0 ) THEN
 310          INFO = -2
 311       ELSE IF( N.LT.0 ) THEN
 312          INFO = -3
 313       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 314          INFO = -5
 315       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 316          INFO = -7
 317       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
 318          INFO = -12
 319       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
 320          INFO = -14
 321       END IF
 322 *
 323 *     Compute workspace
 324 *      (Note: Comments in the code beginning "Workspace:" describe the
 325 *       minimal amount of workspace needed at that point in the code,
 326 *       as well as the preferred amount for good performance.
 327 *       NB refers to the optimal block size for the immediately
 328 *       following subroutine, as returned by ILAENV. The workspace is
 329 *       computed assuming ILO = 1 and IHI = N, the worst case.)
 330 *
 331       IF( INFO.EQ.0 ) THEN
 332          MINWRK = MAX( 1, 8*N )
 333          MAXWRK = MAX( 1, N*( 7 +
 334      $                 ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 ) ) )
 335          MAXWRK = MAX( MAXWRK, N*( 7 +
 336      $                 ILAENV( 1, 'DORMQR', ' ', N, 1, N, 0 ) ) )
 337          IF( ILVL ) THEN
 338             MAXWRK = MAX( MAXWRK, N*( 7 +
 339      $                 ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) ) )
 340          END IF
 341          WORK( 1 ) = MAXWRK
 342 *
 343          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
 344      $      INFO = -16
 345       END IF
 346 *
 347       IF( INFO.NE.0 ) THEN
 348          CALL XERBLA( 'DGGEV ', -INFO )
 349          RETURN
 350       ELSE IF( LQUERY ) THEN
 351          RETURN
 352       END IF
 353 *
 354 *     Quick return if possible
 355 *
 356       IF( N.EQ.0 )
 357      $   RETURN
 358 *
 359 *     Get machine constants
 360 *
 361       EPS = DLAMCH( 'P' )
 362       SMLNUM = DLAMCH( 'S' )
 363       BIGNUM = ONE / SMLNUM
 364       CALL DLABAD( SMLNUM, BIGNUM )
 365       SMLNUM = SQRT( SMLNUM ) / EPS
 366       BIGNUM = ONE / SMLNUM
 367 *
 368 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 369 *
 370       ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
 371       ILASCL = .FALSE.
 372       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 373          ANRMTO = SMLNUM
 374          ILASCL = .TRUE.
 375       ELSE IF( ANRM.GT.BIGNUM ) THEN
 376          ANRMTO = BIGNUM
 377          ILASCL = .TRUE.
 378       END IF
 379       IF( ILASCL )
 380      $   CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
 381 *
 382 *     Scale B if max element outside range [SMLNUM,BIGNUM]
 383 *
 384       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
 385       ILBSCL = .FALSE.
 386       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 387          BNRMTO = SMLNUM
 388          ILBSCL = .TRUE.
 389       ELSE IF( BNRM.GT.BIGNUM ) THEN
 390          BNRMTO = BIGNUM
 391          ILBSCL = .TRUE.
 392       END IF
 393       IF( ILBSCL )
 394      $   CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
 395 *
 396 *     Permute the matrices A, B to isolate eigenvalues if possible
 397 *     (Workspace: need 6*N)
 398 *
 399       ILEFT = 1
 400       IRIGHT = N + 1
 401       IWRK = IRIGHT + N
 402       CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
 403      $             WORK( IRIGHT ), WORK( IWRK ), IERR )
 404 *
 405 *     Reduce B to triangular form (QR decomposition of B)
 406 *     (Workspace: need N, prefer N*NB)
 407 *
 408       IROWS = IHI + 1 - ILO
 409       IF( ILV ) THEN
 410          ICOLS = N + 1 - ILO
 411       ELSE
 412          ICOLS = IROWS
 413       END IF
 414       ITAU = IWRK
 415       IWRK = ITAU + IROWS
 416       CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
 417      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
 418 *
 419 *     Apply the orthogonal transformation to matrix A
 420 *     (Workspace: need N, prefer N*NB)
 421 *
 422       CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
 423      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
 424      $             LWORK+1-IWRK, IERR )
 425 *
 426 *     Initialize VL
 427 *     (Workspace: need N, prefer N*NB)
 428 *
 429       IF( ILVL ) THEN
 430          CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
 431          IF( IROWS.GT.1 ) THEN
 432             CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
 433      $                   VL( ILO+1, ILO ), LDVL )
 434          END IF
 435          CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
 436      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
 437       END IF
 438 *
 439 *     Initialize VR
 440 *
 441       IF( ILVR )
 442      $   CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
 443 *
 444 *     Reduce to generalized Hessenberg form
 445 *     (Workspace: none needed)
 446 *
 447       IF( ILV ) THEN
 448 *
 449 *        Eigenvectors requested -- work on whole matrix.
 450 *
 451          CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
 452      $                LDVL, VR, LDVR, IERR )
 453       ELSE
 454          CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
 455      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
 456       END IF
 457 *
 458 *     Perform QZ algorithm (Compute eigenvalues, and optionally, the
 459 *     Schur forms and Schur vectors)
 460 *     (Workspace: need N)
 461 *
 462       IWRK = ITAU
 463       IF( ILV ) THEN
 464          CHTEMP = 'S'
 465       ELSE
 466          CHTEMP = 'E'
 467       END IF
 468       CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
 469      $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
 470      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
 471       IF( IERR.NE.0 ) THEN
 472          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
 473             INFO = IERR
 474          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
 475             INFO = IERR - N
 476          ELSE
 477             INFO = N + 1
 478          END IF
 479          GO TO 110
 480       END IF
 481 *
 482 *     Compute Eigenvectors
 483 *     (Workspace: need 6*N)
 484 *
 485       IF( ILV ) THEN
 486          IF( ILVL ) THEN
 487             IF( ILVR ) THEN
 488                CHTEMP = 'B'
 489             ELSE
 490                CHTEMP = 'L'
 491             END IF
 492          ELSE
 493             CHTEMP = 'R'
 494          END IF
 495          CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
 496      $                VR, LDVR, N, IN, WORK( IWRK ), IERR )
 497          IF( IERR.NE.0 ) THEN
 498             INFO = N + 2
 499             GO TO 110
 500          END IF
 501 *
 502 *        Undo balancing on VL and VR and normalization
 503 *        (Workspace: none needed)
 504 *
 505          IF( ILVL ) THEN
 506             CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
 507      $                   WORK( IRIGHT ), N, VL, LDVL, IERR )
 508             DO 50 JC = 1, N
 509                IF( ALPHAI( JC ).LT.ZERO )
 510      $            GO TO 50
 511                TEMP = ZERO
 512                IF( ALPHAI( JC ).EQ.ZERO ) THEN
 513                   DO 10 JR = 1, N
 514                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
 515    10             CONTINUE
 516                ELSE
 517                   DO 20 JR = 1, N
 518                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
 519      $                      ABS( VL( JR, JC+1 ) ) )
 520    20             CONTINUE
 521                END IF
 522                IF( TEMP.LT.SMLNUM )
 523      $            GO TO 50
 524                TEMP = ONE / TEMP
 525                IF( ALPHAI( JC ).EQ.ZERO ) THEN
 526                   DO 30 JR = 1, N
 527                      VL( JR, JC ) = VL( JR, JC )*TEMP
 528    30             CONTINUE
 529                ELSE
 530                   DO 40 JR = 1, N
 531                      VL( JR, JC ) = VL( JR, JC )*TEMP
 532                      VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
 533    40             CONTINUE
 534                END IF
 535    50       CONTINUE
 536          END IF
 537          IF( ILVR ) THEN
 538             CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
 539      $                   WORK( IRIGHT ), N, VR, LDVR, IERR )
 540             DO 100 JC = 1, N
 541                IF( ALPHAI( JC ).LT.ZERO )
 542      $            GO TO 100
 543                TEMP = ZERO
 544                IF( ALPHAI( JC ).EQ.ZERO ) THEN
 545                   DO 60 JR = 1, N
 546                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
 547    60             CONTINUE
 548                ELSE
 549                   DO 70 JR = 1, N
 550                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
 551      $                      ABS( VR( JR, JC+1 ) ) )
 552    70             CONTINUE
 553                END IF
 554                IF( TEMP.LT.SMLNUM )
 555      $            GO TO 100
 556                TEMP = ONE / TEMP
 557                IF( ALPHAI( JC ).EQ.ZERO ) THEN
 558                   DO 80 JR = 1, N
 559                      VR( JR, JC ) = VR( JR, JC )*TEMP
 560    80             CONTINUE
 561                ELSE
 562                   DO 90 JR = 1, N
 563                      VR( JR, JC ) = VR( JR, JC )*TEMP
 564                      VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
 565    90             CONTINUE
 566                END IF
 567   100       CONTINUE
 568          END IF
 569 *
 570 *        End of eigenvector calculation
 571 *
 572       END IF
 573 *
 574 *     Undo scaling if necessary
 575 *
 576   110 CONTINUE
 577 *
 578       IF( ILASCL ) THEN
 579          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
 580          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
 581       END IF
 582 *
 583       IF( ILBSCL ) THEN
 584          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
 585       END IF
 586 *
 587       WORK( 1 ) = MAXWRK
 588       RETURN
 589 *
 590 *     End of DGGEV
 591 *
 592       END