SRC/zgeevx.f

   1 *> \brief <b> ZGEEVX 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 ZGEEVX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeevx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeevx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeevx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL,
  22 *                          LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE,
  23 *                          RCONDV, WORK, LWORK, RWORK, 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 *       DOUBLE PRECISION   RCONDE( * ), RCONDV( * ), RWORK( * ),
  32 *      $                   SCALE( * )
  33 *       COMPLEX*16         A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
  34 *      $                   W( * ), WORK( * )
  35 *       ..
  36 *
  37 *
  38 *> \par Purpose:
  39 *  =============
  40 *>
  41 *> \verbatim
  42 *>
  43 *> ZGEEVX computes for an N-by-N complex 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, ie. 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 COMPLEX*16 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 Schur form of the balanced
 138 *>          version of the 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] W
 148 *> \verbatim
 149 *>          W is COMPLEX*16 array, dimension (N)
 150 *>          W contains the computed eigenvalues.
 151 *> \endverbatim
 152 *>
 153 *> \param[out] VL
 154 *> \verbatim
 155 *>          VL is COMPLEX*16 array, dimension (LDVL,N)
 156 *>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
 157 *>          after another in the columns of VL, in the same order
 158 *>          as their eigenvalues.
 159 *>          If JOBVL = 'N', VL is not referenced.
 160 *>          u(j) = VL(:,j), the j-th column of VL.
 161 *> \endverbatim
 162 *>
 163 *> \param[in] LDVL
 164 *> \verbatim
 165 *>          LDVL is INTEGER
 166 *>          The leading dimension of the array VL.  LDVL >= 1; if
 167 *>          JOBVL = 'V', LDVL >= N.
 168 *> \endverbatim
 169 *>
 170 *> \param[out] VR
 171 *> \verbatim
 172 *>          VR is COMPLEX*16 array, dimension (LDVR,N)
 173 *>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
 174 *>          after another in the columns of VR, in the same order
 175 *>          as their eigenvalues.
 176 *>          If JOBVR = 'N', VR is not referenced.
 177 *>          v(j) = VR(:,j), the j-th column of VR.
 178 *> \endverbatim
 179 *>
 180 *> \param[in] LDVR
 181 *> \verbatim
 182 *>          LDVR is INTEGER
 183 *>          The leading dimension of the array VR.  LDVR >= 1; if
 184 *>          JOBVR = 'V', LDVR >= N.
 185 *> \endverbatim
 186 *>
 187 *> \param[out] ILO
 188 *> \verbatim
 189 *>          ILO is INTEGER
 190 *> \endverbatim
 191 *>
 192 *> \param[out] IHI
 193 *> \verbatim
 194 *>          IHI is INTEGER
 195 *>          ILO and IHI are integer values determined when A was
 196 *>          balanced.  The balanced A(i,j) = 0 if I > J and
 197 *>          J = 1,...,ILO-1 or I = IHI+1,...,N.
 198 *> \endverbatim
 199 *>
 200 *> \param[out] SCALE
 201 *> \verbatim
 202 *>          SCALE is DOUBLE PRECISION array, dimension (N)
 203 *>          Details of the permutations and scaling factors applied
 204 *>          when balancing A.  If P(j) is the index of the row and column
 205 *>          interchanged with row and column j, and D(j) is the scaling
 206 *>          factor applied to row and column j, then
 207 *>          SCALE(J) = P(J),    for J = 1,...,ILO-1
 208 *>                   = D(J),    for J = ILO,...,IHI
 209 *>                   = P(J)     for J = IHI+1,...,N.
 210 *>          The order in which the interchanges are made is N to IHI+1,
 211 *>          then 1 to ILO-1.
 212 *> \endverbatim
 213 *>
 214 *> \param[out] ABNRM
 215 *> \verbatim
 216 *>          ABNRM is DOUBLE PRECISION
 217 *>          The one-norm of the balanced matrix (the maximum
 218 *>          of the sum of absolute values of elements of any column).
 219 *> \endverbatim
 220 *>
 221 *> \param[out] RCONDE
 222 *> \verbatim
 223 *>          RCONDE is DOUBLE PRECISION array, dimension (N)
 224 *>          RCONDE(j) is the reciprocal condition number of the j-th
 225 *>          eigenvalue.
 226 *> \endverbatim
 227 *>
 228 *> \param[out] RCONDV
 229 *> \verbatim
 230 *>          RCONDV is DOUBLE PRECISION array, dimension (N)
 231 *>          RCONDV(j) is the reciprocal condition number of the j-th
 232 *>          right eigenvector.
 233 *> \endverbatim
 234 *>
 235 *> \param[out] WORK
 236 *> \verbatim
 237 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 238 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 239 *> \endverbatim
 240 *>
 241 *> \param[in] LWORK
 242 *> \verbatim
 243 *>          LWORK is INTEGER
 244 *>          The dimension of the array WORK.  If SENSE = 'N' or 'E',
 245 *>          LWORK >= max(1,2*N), and if SENSE = 'V' or 'B',
 246 *>          LWORK >= N*N+2*N.
 247 *>          For good performance, LWORK must generally be larger.
 248 *>
 249 *>          If LWORK = -1, then a workspace query is assumed; the routine
 250 *>          only calculates the optimal size of the WORK array, returns
 251 *>          this value as the first entry of the WORK array, and no error
 252 *>          message related to LWORK is issued by XERBLA.
 253 *> \endverbatim
 254 *>
 255 *> \param[out] RWORK
 256 *> \verbatim
 257 *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
 258 *> \endverbatim
 259 *>
 260 *> \param[out] INFO
 261 *> \verbatim
 262 *>          INFO is INTEGER
 263 *>          = 0:  successful exit
 264 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 265 *>          > 0:  if INFO = i, the QR algorithm failed to compute all the
 266 *>                eigenvalues, and no eigenvectors or condition numbers
 267 *>                have been computed; elements 1:ILO-1 and i+1:N of W
 268 *>                contain eigenvalues which have converged.
 269 *> \endverbatim
 270 *
 271 *  Authors:
 272 *  ========
 273 *
 274 *> \author Univ. of Tennessee
 275 *> \author Univ. of California Berkeley
 276 *> \author Univ. of Colorado Denver
 277 *> \author NAG Ltd.
 278 *
 279 *> \date June 2016
 280 *
 281 *  @precisions fortran z -> c
 282 *
 283 *> \ingroup complex16GEeigen
 284 *
 285 *  =====================================================================
 286       SUBROUTINE ZGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL,
 287      $                   LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE,
 288      $                   RCONDV, WORK, LWORK, RWORK, INFO )
 289       implicit none
 290 *
 291 *  -- LAPACK driver routine (version 3.6.1) --
 292 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 293 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 294 *     June 2016
 295 *
 296 *     .. Scalar Arguments ..
 297       CHARACTER          BALANC, JOBVL, JOBVR, SENSE
 298       INTEGER            IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
 299       DOUBLE PRECISION   ABNRM
 300 *     ..
 301 *     .. Array Arguments ..
 302       DOUBLE PRECISION   RCONDE( * ), RCONDV( * ), RWORK( * ),
 303      $                   SCALE( * )
 304       COMPLEX*16         A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
 305      $                   W( * ), WORK( * )
 306 *     ..
 307 *
 308 *  =====================================================================
 309 *
 310 *     .. Parameters ..
 311       DOUBLE PRECISION   ZERO, ONE
 312       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 313 *     ..
 314 *     .. Local Scalars ..
 315       LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
 316      $                   WNTSNN, WNTSNV
 317       CHARACTER          JOB, SIDE
 318       INTEGER            HSWORK, I, ICOND, IERR, ITAU, IWRK, K,
 319      $                   LWORK_TREVC, MAXWRK, MINWRK, NOUT
 320       DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
 321       COMPLEX*16         TMP
 322 *     ..
 323 *     .. Local Arrays ..
 324       LOGICAL            SELECT( 1 )
 325       DOUBLE PRECISION   DUM( 1 )
 326 *     ..
 327 *     .. External Subroutines ..
 328       EXTERNAL           DLABAD, DLASCL, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL,
 329      $                   ZGEHRD, ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC3,
 330      $                   ZTRSNA, ZUNGHR
 331 *     ..
 332 *     .. External Functions ..
 333       LOGICAL            LSAME
 334       INTEGER            IDAMAX, ILAENV
 335       DOUBLE PRECISION   DLAMCH, DZNRM2, ZLANGE
 336       EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DZNRM2, ZLANGE
 337 *     ..
 338 *     .. Intrinsic Functions ..
 339       INTRINSIC          DBLE, DCMPLX, CONJG, AIMAG, MAX, SQRT
 340 *     ..
 341 *     .. Executable Statements ..
 342 *
 343 *     Test the input arguments
 344 *
 345       INFO = 0
 346       LQUERY = ( LWORK.EQ.-1 )
 347       WANTVL = LSAME( JOBVL, 'V' )
 348       WANTVR = LSAME( JOBVR, 'V' )
 349       WNTSNN = LSAME( SENSE, 'N' )
 350       WNTSNE = LSAME( SENSE, 'E' )
 351       WNTSNV = LSAME( SENSE, 'V' )
 352       WNTSNB = LSAME( SENSE, 'B' )
 353       IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'S' ) .OR.
 354      $    LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) ) THEN
 355          INFO = -1
 356       ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
 357          INFO = -2
 358       ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
 359          INFO = -3
 360       ELSE IF( .NOT.( WNTSNN .OR. WNTSNE .OR. WNTSNB .OR. WNTSNV ) .OR.
 361      $         ( ( WNTSNE .OR. WNTSNB ) .AND. .NOT.( WANTVL .AND.
 362      $         WANTVR ) ) ) THEN
 363          INFO = -4
 364       ELSE IF( N.LT.0 ) THEN
 365          INFO = -5
 366       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 367          INFO = -7
 368       ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
 369          INFO = -10
 370       ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
 371          INFO = -12
 372       END IF
 373 *
 374 *     Compute workspace
 375 *      (Note: Comments in the code beginning "Workspace:" describe the
 376 *       minimal amount of workspace needed at that point in the code,
 377 *       as well as the preferred amount for good performance.
 378 *       CWorkspace refers to complex workspace, and RWorkspace to real
 379 *       workspace. NB refers to the optimal block size for the
 380 *       immediately following subroutine, as returned by ILAENV.
 381 *       HSWORK refers to the workspace preferred by ZHSEQR, as
 382 *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
 383 *       the worst case.)
 384 *
 385       IF( INFO.EQ.0 ) THEN
 386          IF( N.EQ.0 ) THEN
 387             MINWRK = 1
 388             MAXWRK = 1
 389          ELSE
 390             MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
 391 *
 392             IF( WANTVL ) THEN
 393                CALL ZTREVC3( 'L', 'B', SELECT, N, A, LDA,
 394      $                       VL, LDVL, VR, LDVR,
 395      $                       N, NOUT, WORK, -1, RWORK, -1, IERR )
 396                LWORK_TREVC = INT( WORK(1) )
 397                MAXWRK = MAX( MAXWRK, LWORK_TREVC )
 398                CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
 399      $                WORK, -1, INFO )
 400             ELSE IF( WANTVR ) THEN
 401                CALL ZTREVC3( 'R', 'B', SELECT, N, A, LDA,
 402      $                       VL, LDVL, VR, LDVR,
 403      $                       N, NOUT, WORK, -1, RWORK, -1, IERR )
 404                LWORK_TREVC = INT( WORK(1) )
 405                MAXWRK = MAX( MAXWRK, LWORK_TREVC )
 406                CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
 407      $                WORK, -1, INFO )
 408             ELSE
 409                IF( WNTSNN ) THEN
 410                   CALL ZHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
 411      $                WORK, -1, INFO )
 412                ELSE
 413                   CALL ZHSEQR( 'S', 'N', N, 1, N, A, LDA, W, VR, LDVR,
 414      $                WORK, -1, INFO )
 415                END IF
 416             END IF
 417             HSWORK = INT( WORK(1) )
 418 *
 419             IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
 420                MINWRK = 2*N
 421                IF( .NOT.( WNTSNN .OR. WNTSNE ) )
 422      $            MINWRK = MAX( MINWRK, N*N + 2*N )
 423                MAXWRK = MAX( MAXWRK, HSWORK )
 424                IF( .NOT.( WNTSNN .OR. WNTSNE ) )
 425      $            MAXWRK = MAX( MAXWRK, N*N + 2*N )
 426             ELSE
 427                MINWRK = 2*N
 428                IF( .NOT.( WNTSNN .OR. WNTSNE ) )
 429      $            MINWRK = MAX( MINWRK, N*N + 2*N )
 430                MAXWRK = MAX( MAXWRK, HSWORK )
 431                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
 432      $                       ' ', N, 1, N, -1 ) )
 433                IF( .NOT.( WNTSNN .OR. WNTSNE ) )
 434      $            MAXWRK = MAX( MAXWRK, N*N + 2*N )
 435                MAXWRK = MAX( MAXWRK, 2*N )
 436             END IF
 437             MAXWRK = MAX( MAXWRK, MINWRK )
 438          END IF
 439          WORK( 1 ) = MAXWRK
 440 *
 441          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
 442             INFO = -20
 443          END IF
 444       END IF
 445 *
 446       IF( INFO.NE.0 ) THEN
 447          CALL XERBLA( 'ZGEEVX', -INFO )
 448          RETURN
 449       ELSE IF( LQUERY ) THEN
 450          RETURN
 451       END IF
 452 *
 453 *     Quick return if possible
 454 *
 455       IF( N.EQ.0 )
 456      $   RETURN
 457 *
 458 *     Get machine constants
 459 *
 460       EPS = DLAMCH( 'P' )
 461       SMLNUM = DLAMCH( 'S' )
 462       BIGNUM = ONE / SMLNUM
 463       CALL DLABAD( SMLNUM, BIGNUM )
 464       SMLNUM = SQRT( SMLNUM ) / EPS
 465       BIGNUM = ONE / SMLNUM
 466 *
 467 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 468 *
 469       ICOND = 0
 470       ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
 471       SCALEA = .FALSE.
 472       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 473          SCALEA = .TRUE.
 474          CSCALE = SMLNUM
 475       ELSE IF( ANRM.GT.BIGNUM ) THEN
 476          SCALEA = .TRUE.
 477          CSCALE = BIGNUM
 478       END IF
 479       IF( SCALEA )
 480      $   CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
 481 *
 482 *     Balance the matrix and compute ABNRM
 483 *
 484       CALL ZGEBAL( BALANC, N, A, LDA, ILO, IHI, SCALE, IERR )
 485       ABNRM = ZLANGE( '1', N, N, A, LDA, DUM )
 486       IF( SCALEA ) THEN
 487          DUM( 1 ) = ABNRM
 488          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
 489          ABNRM = DUM( 1 )
 490       END IF
 491 *
 492 *     Reduce to upper Hessenberg form
 493 *     (CWorkspace: need 2*N, prefer N+N*NB)
 494 *     (RWorkspace: none)
 495 *
 496       ITAU = 1
 497       IWRK = ITAU + N
 498       CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
 499      $             LWORK-IWRK+1, IERR )
 500 *
 501       IF( WANTVL ) THEN
 502 *
 503 *        Want left eigenvectors
 504 *        Copy Householder vectors to VL
 505 *
 506          SIDE = 'L'
 507          CALL ZLACPY( 'L', N, N, A, LDA, VL, LDVL )
 508 *
 509 *        Generate unitary matrix in VL
 510 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
 511 *        (RWorkspace: none)
 512 *
 513          CALL ZUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
 514      $                LWORK-IWRK+1, IERR )
 515 *
 516 *        Perform QR iteration, accumulating Schur vectors in VL
 517 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
 518 *        (RWorkspace: none)
 519 *
 520          IWRK = ITAU
 521          CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
 522      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
 523 *
 524          IF( WANTVR ) THEN
 525 *
 526 *           Want left and right eigenvectors
 527 *           Copy Schur vectors to VR
 528 *
 529             SIDE = 'B'
 530             CALL ZLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
 531          END IF
 532 *
 533       ELSE IF( WANTVR ) THEN
 534 *
 535 *        Want right eigenvectors
 536 *        Copy Householder vectors to VR
 537 *
 538          SIDE = 'R'
 539          CALL ZLACPY( 'L', N, N, A, LDA, VR, LDVR )
 540 *
 541 *        Generate unitary matrix in VR
 542 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
 543 *        (RWorkspace: none)
 544 *
 545          CALL ZUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
 546      $                LWORK-IWRK+1, IERR )
 547 *
 548 *        Perform QR iteration, accumulating Schur vectors in VR
 549 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
 550 *        (RWorkspace: none)
 551 *
 552          IWRK = ITAU
 553          CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
 554      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
 555 *
 556       ELSE
 557 *
 558 *        Compute eigenvalues only
 559 *        If condition numbers desired, compute Schur form
 560 *
 561          IF( WNTSNN ) THEN
 562             JOB = 'E'
 563          ELSE
 564             JOB = 'S'
 565          END IF
 566 *
 567 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
 568 *        (RWorkspace: none)
 569 *
 570          IWRK = ITAU
 571          CALL ZHSEQR( JOB, 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
 572      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
 573       END IF
 574 *
 575 *     If INFO .NE. 0 from ZHSEQR, then quit
 576 *
 577       IF( INFO.NE.0 )
 578      $   GO TO 50
 579 *
 580       IF( WANTVL .OR. WANTVR ) THEN
 581 *
 582 *        Compute left and/or right eigenvectors
 583 *        (CWorkspace: need 2*N, prefer N + 2*N*NB)
 584 *        (RWorkspace: need N)
 585 *
 586          CALL ZTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
 587      $                 N, NOUT, WORK( IWRK ), LWORK-IWRK+1,
 588      $                 RWORK, N, IERR )
 589       END IF
 590 *
 591 *     Compute condition numbers if desired
 592 *     (CWorkspace: need N*N+2*N unless SENSE = 'E')
 593 *     (RWorkspace: need 2*N unless SENSE = 'E')
 594 *
 595       IF( .NOT.WNTSNN ) THEN
 596          CALL ZTRSNA( SENSE, 'A', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
 597      $                RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, RWORK,
 598      $                ICOND )
 599       END IF
 600 *
 601       IF( WANTVL ) THEN
 602 *
 603 *        Undo balancing of left eigenvectors
 604 *
 605          CALL ZGEBAK( BALANC, 'L', N, ILO, IHI, SCALE, N, VL, LDVL,
 606      $                IERR )
 607 *
 608 *        Normalize left eigenvectors and make largest component real
 609 *
 610          DO 20 I = 1, N
 611             SCL = ONE / DZNRM2( N, VL( 1, I ), 1 )
 612             CALL ZDSCAL( N, SCL, VL( 1, I ), 1 )
 613             DO 10 K = 1, N
 614                RWORK( K ) = DBLE( VL( K, I ) )**2 +
 615      $                      AIMAG( VL( K, I ) )**2
 616    10       CONTINUE
 617             K = IDAMAX( N, RWORK, 1 )
 618             TMP = CONJG( VL( K, I ) ) / SQRT( RWORK( K ) )
 619             CALL ZSCAL( N, TMP, VL( 1, I ), 1 )
 620             VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO )
 621    20    CONTINUE
 622       END IF
 623 *
 624       IF( WANTVR ) THEN
 625 *
 626 *        Undo balancing of right eigenvectors
 627 *
 628          CALL ZGEBAK( BALANC, 'R', N, ILO, IHI, SCALE, N, VR, LDVR,
 629      $                IERR )
 630 *
 631 *        Normalize right eigenvectors and make largest component real
 632 *
 633          DO 40 I = 1, N
 634             SCL = ONE / DZNRM2( N, VR( 1, I ), 1 )
 635             CALL ZDSCAL( N, SCL, VR( 1, I ), 1 )
 636             DO 30 K = 1, N
 637                RWORK( K ) = DBLE( VR( K, I ) )**2 +
 638      $                      AIMAG( VR( K, I ) )**2
 639    30       CONTINUE
 640             K = IDAMAX( N, RWORK, 1 )
 641             TMP = CONJG( VR( K, I ) ) / SQRT( RWORK( K ) )
 642             CALL ZSCAL( N, TMP, VR( 1, I ), 1 )
 643             VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO )
 644    40    CONTINUE
 645       END IF
 646 *
 647 *     Undo scaling if necessary
 648 *
 649    50 CONTINUE
 650       IF( SCALEA ) THEN
 651          CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
 652      $                MAX( N-INFO, 1 ), IERR )
 653          IF( INFO.EQ.0 ) THEN
 654             IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 )
 655      $         CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, RCONDV, N,
 656      $                      IERR )
 657          ELSE
 658             CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
 659          END IF
 660       END IF
 661 *
 662       WORK( 1 ) = MAXWRK
 663       RETURN
 664 *
 665 *     End of ZGEEVX
 666 *
 667       END