SRC/dgeevx.f

   1 *> \brief <b> DGEEVX 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 DGEEVX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeevx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeevx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeevx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
  22 *                          VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
  23 *                          RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          BALANC, JOBVL, JOBVR, SENSE
  27 *       INTEGER            IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
  28 *       DOUBLE PRECISION   ABNRM
  29 *       ..
  30 *       .. Array Arguments ..
  31 *       INTEGER            IWORK( * )
  32 *       DOUBLE PRECISION   A( LDA, * ), RCONDE( * ), RCONDV( * ),
  33 *      $                   SCALE( * ), VL( LDVL, * ), VR( LDVR, * ),
  34 *      $                   WI( * ), WORK( * ), WR( * )
  35 *       ..
  36 *
  37 *
  38 *> \par Purpose:
  39 *  =============
  40 *>
  41 *> \verbatim
  42 *>
  43 *> DGEEVX computes for an N-by-N real nonsymmetric matrix A, the
  44 *> eigenvalues and, optionally, the left and/or right eigenvectors.
  45 *>
  46 *> Optionally also, it computes a balancing transformation to improve
  47 *> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
  48 *> SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
  49 *> (RCONDE), and reciprocal condition numbers for the right
  50 *> eigenvectors (RCONDV).
  51 *>
  52 *> The right eigenvector v(j) of A satisfies
  53 *>                  A * v(j) = lambda(j) * v(j)
  54 *> where lambda(j) is its eigenvalue.
  55 *> The left eigenvector u(j) of A satisfies
  56 *>               u(j)**H * A = lambda(j) * u(j)**H
  57 *> where u(j)**H denotes the conjugate-transpose of u(j).
  58 *>
  59 *> The computed eigenvectors are normalized to have Euclidean norm
  60 *> equal to 1 and largest component real.
  61 *>
  62 *> Balancing a matrix means permuting the rows and columns to make it
  63 *> more nearly upper triangular, and applying a diagonal similarity
  64 *> transformation D * A * D**(-1), where D is a diagonal matrix, to
  65 *> make its rows and columns closer in norm and the condition numbers
  66 *> of its eigenvalues and eigenvectors smaller.  The computed
  67 *> reciprocal condition numbers correspond to the balanced matrix.
  68 *> Permuting rows and columns will not change the condition numbers
  69 *> (in exact arithmetic) but diagonal scaling will.  For further
  70 *> explanation of balancing, see section 4.10.2 of the LAPACK
  71 *> Users' Guide.
  72 *> \endverbatim
  73 *
  74 *  Arguments:
  75 *  ==========
  76 *
  77 *> \param[in] BALANC
  78 *> \verbatim
  79 *>          BALANC is CHARACTER*1
  80 *>          Indicates how the input matrix should be diagonally scaled
  81 *>          and/or permuted to improve the conditioning of its
  82 *>          eigenvalues.
  83 *>          = 'N': Do not diagonally scale or permute;
  84 *>          = 'P': Perform permutations to make the matrix more nearly
  85 *>                 upper triangular. Do not diagonally scale;
  86 *>          = 'S': Diagonally scale the matrix, i.e. replace A by
  87 *>                 D*A*D**(-1), where D is a diagonal matrix chosen
  88 *>                 to make the rows and columns of A more equal in
  89 *>                 norm. Do not permute;
  90 *>          = 'B': Both diagonally scale and permute A.
  91 *>
  92 *>          Computed reciprocal condition numbers will be for the matrix
  93 *>          after balancing and/or permuting. Permuting does not change
  94 *>          condition numbers (in exact arithmetic), but balancing does.
  95 *> \endverbatim
  96 *>
  97 *> \param[in] JOBVL
  98 *> \verbatim
  99 *>          JOBVL is CHARACTER*1
 100 *>          = 'N': left eigenvectors of A are not computed;
 101 *>          = 'V': left eigenvectors of A are computed.
 102 *>          If SENSE = 'E' or 'B', JOBVL must = 'V'.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] JOBVR
 106 *> \verbatim
 107 *>          JOBVR is CHARACTER*1
 108 *>          = 'N': right eigenvectors of A are not computed;
 109 *>          = 'V': right eigenvectors of A are computed.
 110 *>          If SENSE = 'E' or 'B', JOBVR must = 'V'.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] SENSE
 114 *> \verbatim
 115 *>          SENSE is CHARACTER*1
 116 *>          Determines which reciprocal condition numbers are computed.
 117 *>          = 'N': None are computed;
 118 *>          = 'E': Computed for eigenvalues only;
 119 *>          = 'V': Computed for right eigenvectors only;
 120 *>          = 'B': Computed for eigenvalues and right eigenvectors.
 121 *>
 122 *>          If SENSE = 'E' or 'B', both left and right eigenvectors
 123 *>          must also be computed (JOBVL = 'V' and JOBVR = 'V').
 124 *> \endverbatim
 125 *>
 126 *> \param[in] N
 127 *> \verbatim
 128 *>          N is INTEGER
 129 *>          The order of the matrix A. N >= 0.
 130 *> \endverbatim
 131 *>
 132 *> \param[in,out] A
 133 *> \verbatim
 134 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
 135 *>          On entry, the N-by-N matrix A.
 136 *>          On exit, A has been overwritten.  If JOBVL = 'V' or
 137 *>          JOBVR = 'V', A contains the real Schur form of the balanced
 138 *>          version of the input matrix A.
 139 *> \endverbatim
 140 *>
 141 *> \param[in] LDA
 142 *> \verbatim
 143 *>          LDA is INTEGER
 144 *>          The leading dimension of the array A.  LDA >= max(1,N).
 145 *> \endverbatim
 146 *>
 147 *> \param[out] WR
 148 *> \verbatim
 149 *>          WR is DOUBLE PRECISION array, dimension (N)
 150 *> \endverbatim
 151 *>
 152 *> \param[out] WI
 153 *> \verbatim
 154 *>          WI is DOUBLE PRECISION array, dimension (N)
 155 *>          WR and WI contain the real and imaginary parts,
 156 *>          respectively, of the computed eigenvalues.  Complex
 157 *>          conjugate pairs of eigenvalues will appear consecutively
 158 *>          with the eigenvalue having the positive imaginary part
 159 *>          first.
 160 *> \endverbatim
 161 *>
 162 *> \param[out] VL
 163 *> \verbatim
 164 *>          VL is DOUBLE PRECISION array, dimension (LDVL,N)
 165 *>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
 166 *>          after another in the columns of VL, in the same order
 167 *>          as their eigenvalues.
 168 *>          If JOBVL = 'N', VL is not referenced.
 169 *>          If the j-th eigenvalue is real, then u(j) = VL(:,j),
 170 *>          the j-th column of VL.
 171 *>          If the j-th and (j+1)-st eigenvalues form a complex
 172 *>          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
 173 *>          u(j+1) = VL(:,j) - i*VL(:,j+1).
 174 *> \endverbatim
 175 *>
 176 *> \param[in] LDVL
 177 *> \verbatim
 178 *>          LDVL is INTEGER
 179 *>          The leading dimension of the array VL.  LDVL >= 1; if
 180 *>          JOBVL = 'V', LDVL >= N.
 181 *> \endverbatim
 182 *>
 183 *> \param[out] VR
 184 *> \verbatim
 185 *>          VR is DOUBLE PRECISION array, dimension (LDVR,N)
 186 *>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
 187 *>          after another in the columns of VR, in the same order
 188 *>          as their eigenvalues.
 189 *>          If JOBVR = 'N', VR is not referenced.
 190 *>          If the j-th eigenvalue is real, then v(j) = VR(:,j),
 191 *>          the j-th column of VR.
 192 *>          If the j-th and (j+1)-st eigenvalues form a complex
 193 *>          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
 194 *>          v(j+1) = VR(:,j) - i*VR(:,j+1).
 195 *> \endverbatim
 196 *>
 197 *> \param[in] LDVR
 198 *> \verbatim
 199 *>          LDVR is INTEGER
 200 *>          The leading dimension of the array VR.  LDVR >= 1, and if
 201 *>          JOBVR = 'V', LDVR >= N.
 202 *> \endverbatim
 203 *>
 204 *> \param[out] ILO
 205 *> \verbatim
 206 *>          ILO is INTEGER
 207 *> \endverbatim
 208 *>
 209 *> \param[out] IHI
 210 *> \verbatim
 211 *>          IHI is INTEGER
 212 *>          ILO and IHI are integer values determined when A was
 213 *>          balanced.  The balanced A(i,j) = 0 if I > J and
 214 *>          J = 1,...,ILO-1 or I = IHI+1,...,N.
 215 *> \endverbatim
 216 *>
 217 *> \param[out] SCALE
 218 *> \verbatim
 219 *>          SCALE is DOUBLE PRECISION array, dimension (N)
 220 *>          Details of the permutations and scaling factors applied
 221 *>          when balancing A.  If P(j) is the index of the row and column
 222 *>          interchanged with row and column j, and D(j) is the scaling
 223 *>          factor applied to row and column j, then
 224 *>          SCALE(J) = P(J),    for J = 1,...,ILO-1
 225 *>                   = D(J),    for J = ILO,...,IHI
 226 *>                   = P(J)     for J = IHI+1,...,N.
 227 *>          The order in which the interchanges are made is N to IHI+1,
 228 *>          then 1 to ILO-1.
 229 *> \endverbatim
 230 *>
 231 *> \param[out] ABNRM
 232 *> \verbatim
 233 *>          ABNRM is DOUBLE PRECISION
 234 *>          The one-norm of the balanced matrix (the maximum
 235 *>          of the sum of absolute values of elements of any column).
 236 *> \endverbatim
 237 *>
 238 *> \param[out] RCONDE
 239 *> \verbatim
 240 *>          RCONDE is DOUBLE PRECISION array, dimension (N)
 241 *>          RCONDE(j) is the reciprocal condition number of the j-th
 242 *>          eigenvalue.
 243 *> \endverbatim
 244 *>
 245 *> \param[out] RCONDV
 246 *> \verbatim
 247 *>          RCONDV is DOUBLE PRECISION array, dimension (N)
 248 *>          RCONDV(j) is the reciprocal condition number of the j-th
 249 *>          right eigenvector.
 250 *> \endverbatim
 251 *>
 252 *> \param[out] WORK
 253 *> \verbatim
 254 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 255 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 256 *> \endverbatim
 257 *>
 258 *> \param[in] LWORK
 259 *> \verbatim
 260 *>          LWORK is INTEGER
 261 *>          The dimension of the array WORK.   If SENSE = 'N' or 'E',
 262 *>          LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V',
 263 *>          LWORK >= 3*N.  If SENSE = 'V' or 'B', LWORK >= N*(N+6).
 264 *>          For good performance, LWORK must generally be larger.
 265 *>
 266 *>          If LWORK = -1, then a workspace query is assumed; the routine
 267 *>          only calculates the optimal size of the WORK array, returns
 268 *>          this value as the first entry of the WORK array, and no error
 269 *>          message related to LWORK is issued by XERBLA.
 270 *> \endverbatim
 271 *>
 272 *> \param[out] IWORK
 273 *> \verbatim
 274 *>          IWORK is INTEGER array, dimension (2*N-2)
 275 *>          If SENSE = 'N' or 'E', not referenced.
 276 *> \endverbatim
 277 *>
 278 *> \param[out] INFO
 279 *> \verbatim
 280 *>          INFO is INTEGER
 281 *>          = 0:  successful exit
 282 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 283 *>          > 0:  if INFO = i, the QR algorithm failed to compute all the
 284 *>                eigenvalues, and no eigenvectors or condition numbers
 285 *>                have been computed; elements 1:ILO-1 and i+1:N of WR
 286 *>                and WI contain eigenvalues which have converged.
 287 *> \endverbatim
 288 *
 289 *  Authors:
 290 *  ========
 291 *
 292 *> \author Univ. of Tennessee
 293 *> \author Univ. of California Berkeley
 294 *> \author Univ. of Colorado Denver
 295 *> \author NAG Ltd.
 296 *
 297 *> \date June 2016
 298 *
 299 *  @precisions fortran d -> s
 300 *
 301 *> \ingroup doubleGEeigen
 302 *
 303 *  =====================================================================
 304       SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
 305      $                   VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
 306      $                   RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
 307       implicit none
 308 *
 309 *  -- LAPACK driver routine (version 3.6.1) --
 310 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 311 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 312 *     June 2016
 313 *
 314 *     .. Scalar Arguments ..
 315       CHARACTER          BALANC, JOBVL, JOBVR, SENSE
 316       INTEGER            IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
 317       DOUBLE PRECISION   ABNRM
 318 *     ..
 319 *     .. Array Arguments ..
 320       INTEGER            IWORK( * )
 321       DOUBLE PRECISION   A( LDA, * ), RCONDE( * ), RCONDV( * ),
 322      $                   SCALE( * ), VL( LDVL, * ), VR( LDVR, * ),
 323      $                   WI( * ), WORK( * ), WR( * )
 324 *     ..
 325 *
 326 *  =====================================================================
 327 *
 328 *     .. Parameters ..
 329       DOUBLE PRECISION   ZERO, ONE
 330       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 331 *     ..
 332 *     .. Local Scalars ..
 333       LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
 334      $                   WNTSNN, WNTSNV
 335       CHARACTER          JOB, SIDE
 336       INTEGER            HSWORK, I, ICOND, IERR, ITAU, IWRK, K,
 337      $                   LWORK_TREVC, MAXWRK, MINWRK, NOUT
 338       DOUBLE PRECISION   ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
 339      $                   SN
 340 *     ..
 341 *     .. Local Arrays ..
 342       LOGICAL            SELECT( 1 )
 343       DOUBLE PRECISION   DUM( 1 )
 344 *     ..
 345 *     .. External Subroutines ..
 346       EXTERNAL           DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
 347      $                   DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC3,
 348      $                   DTRSNA, XERBLA
 349 *     ..
 350 *     .. External Functions ..
 351       LOGICAL            LSAME
 352       INTEGER            IDAMAX, ILAENV
 353       DOUBLE PRECISION   DLAMCH, DLANGE, DLAPY2, DNRM2
 354       EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
 355      $                   DNRM2
 356 *     ..
 357 *     .. Intrinsic Functions ..
 358       INTRINSIC          MAX, SQRT
 359 *     ..
 360 *     .. Executable Statements ..
 361 *
 362 *     Test the input arguments
 363 *
 364       INFO = 0
 365       LQUERY = ( LWORK.EQ.-1 )
 366       WANTVL = LSAME( JOBVL, 'V' )
 367       WANTVR = LSAME( JOBVR, 'V' )
 368       WNTSNN = LSAME( SENSE, 'N' )
 369       WNTSNE = LSAME( SENSE, 'E' )
 370       WNTSNV = LSAME( SENSE, 'V' )
 371       WNTSNB = LSAME( SENSE, 'B' )
 372       IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'S' )
 373      $      .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) )
 374      $     THEN
 375          INFO = -1
 376       ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
 377          INFO = -2
 378       ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
 379          INFO = -3
 380       ELSE IF( .NOT.( WNTSNN .OR. WNTSNE .OR. WNTSNB .OR. WNTSNV ) .OR.
 381      $         ( ( WNTSNE .OR. WNTSNB ) .AND. .NOT.( WANTVL .AND.
 382      $         WANTVR ) ) ) THEN
 383          INFO = -4
 384       ELSE IF( N.LT.0 ) THEN
 385          INFO = -5
 386       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 387          INFO = -7
 388       ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
 389          INFO = -11
 390       ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
 391          INFO = -13
 392       END IF
 393 *
 394 *     Compute workspace
 395 *      (Note: Comments in the code beginning "Workspace:" describe the
 396 *       minimal amount of workspace needed at that point in the code,
 397 *       as well as the preferred amount for good performance.
 398 *       NB refers to the optimal block size for the immediately
 399 *       following subroutine, as returned by ILAENV.
 400 *       HSWORK refers to the workspace preferred by DHSEQR, as
 401 *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
 402 *       the worst case.)
 403 *
 404       IF( INFO.EQ.0 ) THEN
 405          IF( N.EQ.0 ) THEN
 406             MINWRK = 1
 407             MAXWRK = 1
 408          ELSE
 409             MAXWRK = N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
 410 *
 411             IF( WANTVL ) THEN
 412                CALL DTREVC3( 'L', 'B', SELECT, N, A, LDA,
 413      $                       VL, LDVL, VR, LDVR,
 414      $                       N, NOUT, WORK, -1, IERR )
 415                LWORK_TREVC = INT( WORK(1) )
 416                MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
 417                CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
 418      $                WORK, -1, INFO )
 419             ELSE IF( WANTVR ) THEN
 420                CALL DTREVC3( 'R', 'B', SELECT, N, A, LDA,
 421      $                       VL, LDVL, VR, LDVR,
 422      $                       N, NOUT, WORK, -1, IERR )
 423                LWORK_TREVC = INT( WORK(1) )
 424                MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
 425                CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
 426      $                WORK, -1, INFO )
 427             ELSE
 428                IF( WNTSNN ) THEN
 429                   CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR,
 430      $                LDVR, WORK, -1, INFO )
 431                ELSE
 432                   CALL DHSEQR( 'S', 'N', N, 1, N, A, LDA, WR, WI, VR,
 433      $                LDVR, WORK, -1, INFO )
 434                END IF
 435             END IF
 436             HSWORK = INT( WORK(1) )
 437 *
 438             IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
 439                MINWRK = 2*N
 440                IF( .NOT.WNTSNN )
 441      $            MINWRK = MAX( MINWRK, N*N+6*N )
 442                MAXWRK = MAX( MAXWRK, HSWORK )
 443                IF( .NOT.WNTSNN )
 444      $            MAXWRK = MAX( MAXWRK, N*N + 6*N )
 445             ELSE
 446                MINWRK = 3*N
 447                IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) )
 448      $            MINWRK = MAX( MINWRK, N*N + 6*N )
 449                MAXWRK = MAX( MAXWRK, HSWORK )
 450                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'DORGHR',
 451      $                       ' ', N, 1, N, -1 ) )
 452                IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) )
 453      $            MAXWRK = MAX( MAXWRK, N*N + 6*N )
 454                MAXWRK = MAX( MAXWRK, 3*N )
 455             END IF
 456             MAXWRK = MAX( MAXWRK, MINWRK )
 457          END IF
 458          WORK( 1 ) = MAXWRK
 459 *
 460          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
 461             INFO = -21
 462          END IF
 463       END IF
 464 *
 465       IF( INFO.NE.0 ) THEN
 466          CALL XERBLA( 'DGEEVX', -INFO )
 467          RETURN
 468       ELSE IF( LQUERY ) THEN
 469          RETURN
 470       END IF
 471 *
 472 *     Quick return if possible
 473 *
 474       IF( N.EQ.0 )
 475      $   RETURN
 476 *
 477 *     Get machine constants
 478 *
 479       EPS = DLAMCH( 'P' )
 480       SMLNUM = DLAMCH( 'S' )
 481       BIGNUM = ONE / SMLNUM
 482       CALL DLABAD( SMLNUM, BIGNUM )
 483       SMLNUM = SQRT( SMLNUM ) / EPS
 484       BIGNUM = ONE / SMLNUM
 485 *
 486 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 487 *
 488       ICOND = 0
 489       ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
 490       SCALEA = .FALSE.
 491       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 492          SCALEA = .TRUE.
 493          CSCALE = SMLNUM
 494       ELSE IF( ANRM.GT.BIGNUM ) THEN
 495          SCALEA = .TRUE.
 496          CSCALE = BIGNUM
 497       END IF
 498       IF( SCALEA )
 499      $   CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
 500 *
 501 *     Balance the matrix and compute ABNRM
 502 *
 503       CALL DGEBAL( BALANC, N, A, LDA, ILO, IHI, SCALE, IERR )
 504       ABNRM = DLANGE( '1', N, N, A, LDA, DUM )
 505       IF( SCALEA ) THEN
 506          DUM( 1 ) = ABNRM
 507          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
 508          ABNRM = DUM( 1 )
 509       END IF
 510 *
 511 *     Reduce to upper Hessenberg form
 512 *     (Workspace: need 2*N, prefer N+N*NB)
 513 *
 514       ITAU = 1
 515       IWRK = ITAU + N
 516       CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
 517      $             LWORK-IWRK+1, IERR )
 518 *
 519       IF( WANTVL ) THEN
 520 *
 521 *        Want left eigenvectors
 522 *        Copy Householder vectors to VL
 523 *
 524          SIDE = 'L'
 525          CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
 526 *
 527 *        Generate orthogonal matrix in VL
 528 *        (Workspace: need 2*N-1, prefer N+(N-1)*NB)
 529 *
 530          CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
 531      $                LWORK-IWRK+1, IERR )
 532 *
 533 *        Perform QR iteration, accumulating Schur vectors in VL
 534 *        (Workspace: need 1, prefer HSWORK (see comments) )
 535 *
 536          IWRK = ITAU
 537          CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
 538      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
 539 *
 540          IF( WANTVR ) THEN
 541 *
 542 *           Want left and right eigenvectors
 543 *           Copy Schur vectors to VR
 544 *
 545             SIDE = 'B'
 546             CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
 547          END IF
 548 *
 549       ELSE IF( WANTVR ) THEN
 550 *
 551 *        Want right eigenvectors
 552 *        Copy Householder vectors to VR
 553 *
 554          SIDE = 'R'
 555          CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
 556 *
 557 *        Generate orthogonal matrix in VR
 558 *        (Workspace: need 2*N-1, prefer N+(N-1)*NB)
 559 *
 560          CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
 561      $                LWORK-IWRK+1, IERR )
 562 *
 563 *        Perform QR iteration, accumulating Schur vectors in VR
 564 *        (Workspace: need 1, prefer HSWORK (see comments) )
 565 *
 566          IWRK = ITAU
 567          CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
 568      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
 569 *
 570       ELSE
 571 *
 572 *        Compute eigenvalues only
 573 *        If condition numbers desired, compute Schur form
 574 *
 575          IF( WNTSNN ) THEN
 576             JOB = 'E'
 577          ELSE
 578             JOB = 'S'
 579          END IF
 580 *
 581 *        (Workspace: need 1, prefer HSWORK (see comments) )
 582 *
 583          IWRK = ITAU
 584          CALL DHSEQR( JOB, 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
 585      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
 586       END IF
 587 *
 588 *     If INFO .NE. 0 from DHSEQR, then quit
 589 *
 590       IF( INFO.NE.0 )
 591      $   GO TO 50
 592 *
 593       IF( WANTVL .OR. WANTVR ) THEN
 594 *
 595 *        Compute left and/or right eigenvectors
 596 *        (Workspace: need 3*N, prefer N + 2*N*NB)
 597 *
 598          CALL DTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
 599      $                 N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR )
 600       END IF
 601 *
 602 *     Compute condition numbers if desired
 603 *     (Workspace: need N*N+6*N unless SENSE = 'E')
 604 *
 605       IF( .NOT.WNTSNN ) THEN
 606          CALL DTRSNA( SENSE, 'A', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
 607      $                RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, IWORK,
 608      $                ICOND )
 609       END IF
 610 *
 611       IF( WANTVL ) THEN
 612 *
 613 *        Undo balancing of left eigenvectors
 614 *
 615          CALL DGEBAK( BALANC, 'L', N, ILO, IHI, SCALE, N, VL, LDVL,
 616      $                IERR )
 617 *
 618 *        Normalize left eigenvectors and make largest component real
 619 *
 620          DO 20 I = 1, N
 621             IF( WI( I ).EQ.ZERO ) THEN
 622                SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
 623                CALL DSCAL( N, SCL, VL( 1, I ), 1 )
 624             ELSE IF( WI( I ).GT.ZERO ) THEN
 625                SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
 626      $               DNRM2( N, VL( 1, I+1 ), 1 ) )
 627                CALL DSCAL( N, SCL, VL( 1, I ), 1 )
 628                CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
 629                DO 10 K = 1, N
 630                   WORK( K ) = VL( K, I )**2 + VL( K, I+1 )**2
 631    10          CONTINUE
 632                K = IDAMAX( N, WORK, 1 )
 633                CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
 634                CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
 635                VL( K, I+1 ) = ZERO
 636             END IF
 637    20    CONTINUE
 638       END IF
 639 *
 640       IF( WANTVR ) THEN
 641 *
 642 *        Undo balancing of right eigenvectors
 643 *
 644          CALL DGEBAK( BALANC, 'R', N, ILO, IHI, SCALE, N, VR, LDVR,
 645      $                IERR )
 646 *
 647 *        Normalize right eigenvectors and make largest component real
 648 *
 649          DO 40 I = 1, N
 650             IF( WI( I ).EQ.ZERO ) THEN
 651                SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
 652                CALL DSCAL( N, SCL, VR( 1, I ), 1 )
 653             ELSE IF( WI( I ).GT.ZERO ) THEN
 654                SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
 655      $               DNRM2( N, VR( 1, I+1 ), 1 ) )
 656                CALL DSCAL( N, SCL, VR( 1, I ), 1 )
 657                CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
 658                DO 30 K = 1, N
 659                   WORK( K ) = VR( K, I )**2 + VR( K, I+1 )**2
 660    30          CONTINUE
 661                K = IDAMAX( N, WORK, 1 )
 662                CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
 663                CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
 664                VR( K, I+1 ) = ZERO
 665             END IF
 666    40    CONTINUE
 667       END IF
 668 *
 669 *     Undo scaling if necessary
 670 *
 671    50 CONTINUE
 672       IF( SCALEA ) THEN
 673          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
 674      $                MAX( N-INFO, 1 ), IERR )
 675          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
 676      $                MAX( N-INFO, 1 ), IERR )
 677          IF( INFO.EQ.0 ) THEN
 678             IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 )
 679      $         CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, RCONDV, N,
 680      $                      IERR )
 681          ELSE
 682             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
 683      $                   IERR )
 684             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
 685      $                   IERR )
 686          END IF
 687       END IF
 688 *
 689       WORK( 1 ) = MAXWRK
 690       RETURN
 691 *
 692 *     End of DGEEVX
 693 *
 694       END