SRC/dggev3.f

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