SRC/zggevx.f

   1 *> \brief <b> ZGGEVX 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 ZGGEVX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggevx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggevx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggevx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
  22 *                          ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI,
  23 *                          LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV,
  24 *                          WORK, LWORK, RWORK, IWORK, BWORK, INFO )
  25 *
  26 *       .. Scalar Arguments ..
  27 *       CHARACTER          BALANC, JOBVL, JOBVR, SENSE
  28 *       INTEGER            IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
  29 *       DOUBLE PRECISION   ABNRM, BBNRM
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       LOGICAL            BWORK( * )
  33 *       INTEGER            IWORK( * )
  34 *       DOUBLE PRECISION   LSCALE( * ), RCONDE( * ), RCONDV( * ),
  35 *      $                   RSCALE( * ), RWORK( * )
  36 *       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
  37 *      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
  38 *      $                   WORK( * )
  39 *       ..
  40 *
  41 *
  42 *> \par Purpose:
  43 *  =============
  44 *>
  45 *> \verbatim
  46 *>
  47 *> ZGGEVX computes for a pair of N-by-N complex nonsymmetric matrices
  48 *> (A,B) the generalized eigenvalues, and optionally, the left and/or
  49 *> right generalized eigenvectors.
  50 *>
  51 *> Optionally, it also computes a balancing transformation to improve
  52 *> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
  53 *> LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
  54 *> the eigenvalues (RCONDE), and reciprocal condition numbers for the
  55 *> right eigenvectors (RCONDV).
  56 *>
  57 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
  58 *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
  59 *> singular. It is usually represented as the pair (alpha,beta), as
  60 *> there is a reasonable interpretation for beta=0, and even for both
  61 *> being zero.
  62 *>
  63 *> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
  64 *> of (A,B) satisfies
  65 *>                  A * v(j) = lambda(j) * B * v(j) .
  66 *> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
  67 *> of (A,B) satisfies
  68 *>                  u(j)**H * A  = lambda(j) * u(j)**H * B.
  69 *> where u(j)**H is the conjugate-transpose of u(j).
  70 *>
  71 *> \endverbatim
  72 *
  73 *  Arguments:
  74 *  ==========
  75 *
  76 *> \param[in] BALANC
  77 *> \verbatim
  78 *>          BALANC is CHARACTER*1
  79 *>          Specifies the balance option to be performed:
  80 *>          = 'N':  do not diagonally scale or permute;
  81 *>          = 'P':  permute only;
  82 *>          = 'S':  scale only;
  83 *>          = 'B':  both permute and scale.
  84 *>          Computed reciprocal condition numbers will be for the
  85 *>          matrices after permuting and/or balancing. Permuting does
  86 *>          not change condition numbers (in exact arithmetic), but
  87 *>          balancing does.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] JOBVL
  91 *> \verbatim
  92 *>          JOBVL is CHARACTER*1
  93 *>          = 'N':  do not compute the left generalized eigenvectors;
  94 *>          = 'V':  compute the left generalized eigenvectors.
  95 *> \endverbatim
  96 *>
  97 *> \param[in] JOBVR
  98 *> \verbatim
  99 *>          JOBVR is CHARACTER*1
 100 *>          = 'N':  do not compute the right generalized eigenvectors;
 101 *>          = 'V':  compute the right generalized eigenvectors.
 102 *> \endverbatim
 103 *>
 104 *> \param[in] SENSE
 105 *> \verbatim
 106 *>          SENSE is CHARACTER*1
 107 *>          Determines which reciprocal condition numbers are computed.
 108 *>          = 'N': none are computed;
 109 *>          = 'E': computed for eigenvalues only;
 110 *>          = 'V': computed for eigenvectors only;
 111 *>          = 'B': computed for eigenvalues and eigenvectors.
 112 *> \endverbatim
 113 *>
 114 *> \param[in] N
 115 *> \verbatim
 116 *>          N is INTEGER
 117 *>          The order of the matrices A, B, VL, and VR.  N >= 0.
 118 *> \endverbatim
 119 *>
 120 *> \param[in,out] A
 121 *> \verbatim
 122 *>          A is COMPLEX*16 array, dimension (LDA, N)
 123 *>          On entry, the matrix A in the pair (A,B).
 124 *>          On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
 125 *>          or both, then A contains the first part of the complex Schur
 126 *>          form of the "balanced" versions of the input A and B.
 127 *> \endverbatim
 128 *>
 129 *> \param[in] LDA
 130 *> \verbatim
 131 *>          LDA is INTEGER
 132 *>          The leading dimension of A.  LDA >= max(1,N).
 133 *> \endverbatim
 134 *>
 135 *> \param[in,out] B
 136 *> \verbatim
 137 *>          B is COMPLEX*16 array, dimension (LDB, N)
 138 *>          On entry, the matrix B in the pair (A,B).
 139 *>          On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
 140 *>          or both, then B contains the second part of the complex
 141 *>          Schur form of the "balanced" versions of the input A and B.
 142 *> \endverbatim
 143 *>
 144 *> \param[in] LDB
 145 *> \verbatim
 146 *>          LDB is INTEGER
 147 *>          The leading dimension of B.  LDB >= max(1,N).
 148 *> \endverbatim
 149 *>
 150 *> \param[out] ALPHA
 151 *> \verbatim
 152 *>          ALPHA is COMPLEX*16 array, dimension (N)
 153 *> \endverbatim
 154 *>
 155 *> \param[out] BETA
 156 *> \verbatim
 157 *>          BETA is COMPLEX*16 array, dimension (N)
 158 *>          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the generalized
 159 *>          eigenvalues.
 160 *>
 161 *>          Note: the quotient ALPHA(j)/BETA(j) ) may easily over- or
 162 *>          underflow, and BETA(j) may even be zero.  Thus, the user
 163 *>          should avoid naively computing the ratio ALPHA/BETA.
 164 *>          However, ALPHA will be always less than and usually
 165 *>          comparable with norm(A) in magnitude, and BETA always less
 166 *>          than and usually comparable with norm(B).
 167 *> \endverbatim
 168 *>
 169 *> \param[out] VL
 170 *> \verbatim
 171 *>          VL is COMPLEX*16 array, dimension (LDVL,N)
 172 *>          If JOBVL = 'V', the left generalized eigenvectors u(j) are
 173 *>          stored one after another in the columns of VL, in the same
 174 *>          order as their eigenvalues.
 175 *>          Each eigenvector will be scaled so the largest component
 176 *>          will have abs(real part) + abs(imag. part) = 1.
 177 *>          Not referenced if JOBVL = 'N'.
 178 *> \endverbatim
 179 *>
 180 *> \param[in] LDVL
 181 *> \verbatim
 182 *>          LDVL is INTEGER
 183 *>          The leading dimension of the matrix VL. LDVL >= 1, and
 184 *>          if JOBVL = 'V', LDVL >= N.
 185 *> \endverbatim
 186 *>
 187 *> \param[out] VR
 188 *> \verbatim
 189 *>          VR is COMPLEX*16 array, dimension (LDVR,N)
 190 *>          If JOBVR = 'V', the right generalized eigenvectors v(j) are
 191 *>          stored one after another in the columns of VR, in the same
 192 *>          order as their eigenvalues.
 193 *>          Each eigenvector will be scaled so the largest component
 194 *>          will have abs(real part) + abs(imag. part) = 1.
 195 *>          Not referenced if JOBVR = 'N'.
 196 *> \endverbatim
 197 *>
 198 *> \param[in] LDVR
 199 *> \verbatim
 200 *>          LDVR is INTEGER
 201 *>          The leading dimension of the matrix VR. LDVR >= 1, and
 202 *>          if JOBVR = 'V', LDVR >= N.
 203 *> \endverbatim
 204 *>
 205 *> \param[out] ILO
 206 *> \verbatim
 207 *>          ILO is INTEGER
 208 *> \endverbatim
 209 *>
 210 *> \param[out] IHI
 211 *> \verbatim
 212 *>          IHI is INTEGER
 213 *>          ILO and IHI are integer values such that on exit
 214 *>          A(i,j) = 0 and B(i,j) = 0 if i > j and
 215 *>          j = 1,...,ILO-1 or i = IHI+1,...,N.
 216 *>          If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
 217 *> \endverbatim
 218 *>
 219 *> \param[out] LSCALE
 220 *> \verbatim
 221 *>          LSCALE is DOUBLE PRECISION array, dimension (N)
 222 *>          Details of the permutations and scaling factors applied
 223 *>          to the left side of A and B.  If PL(j) is the index of the
 224 *>          row interchanged with row j, and DL(j) is the scaling
 225 *>          factor applied to row j, then
 226 *>            LSCALE(j) = PL(j)  for j = 1,...,ILO-1
 227 *>                      = DL(j)  for j = ILO,...,IHI
 228 *>                      = PL(j)  for j = IHI+1,...,N.
 229 *>          The order in which the interchanges are made is N to IHI+1,
 230 *>          then 1 to ILO-1.
 231 *> \endverbatim
 232 *>
 233 *> \param[out] RSCALE
 234 *> \verbatim
 235 *>          RSCALE is DOUBLE PRECISION array, dimension (N)
 236 *>          Details of the permutations and scaling factors applied
 237 *>          to the right side of A and B.  If PR(j) is the index of the
 238 *>          column interchanged with column j, and DR(j) is the scaling
 239 *>          factor applied to column j, then
 240 *>            RSCALE(j) = PR(j)  for j = 1,...,ILO-1
 241 *>                      = DR(j)  for j = ILO,...,IHI
 242 *>                      = PR(j)  for j = IHI+1,...,N
 243 *>          The order in which the interchanges are made is N to IHI+1,
 244 *>          then 1 to ILO-1.
 245 *> \endverbatim
 246 *>
 247 *> \param[out] ABNRM
 248 *> \verbatim
 249 *>          ABNRM is DOUBLE PRECISION
 250 *>          The one-norm of the balanced matrix A.
 251 *> \endverbatim
 252 *>
 253 *> \param[out] BBNRM
 254 *> \verbatim
 255 *>          BBNRM is DOUBLE PRECISION
 256 *>          The one-norm of the balanced matrix B.
 257 *> \endverbatim
 258 *>
 259 *> \param[out] RCONDE
 260 *> \verbatim
 261 *>          RCONDE is DOUBLE PRECISION array, dimension (N)
 262 *>          If SENSE = 'E' or 'B', the reciprocal condition numbers of
 263 *>          the eigenvalues, stored in consecutive elements of the array.
 264 *>          If SENSE = 'N' or 'V', RCONDE is not referenced.
 265 *> \endverbatim
 266 *>
 267 *> \param[out] RCONDV
 268 *> \verbatim
 269 *>          RCONDV is DOUBLE PRECISION array, dimension (N)
 270 *>          If JOB = 'V' or 'B', the estimated reciprocal condition
 271 *>          numbers of the eigenvectors, stored in consecutive elements
 272 *>          of the array. If the eigenvalues cannot be reordered to
 273 *>          compute RCONDV(j), RCONDV(j) is set to 0; this can only occur
 274 *>          when the true value would be very small anyway.
 275 *>          If SENSE = 'N' or 'E', RCONDV is not referenced.
 276 *> \endverbatim
 277 *>
 278 *> \param[out] WORK
 279 *> \verbatim
 280 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 281 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 282 *> \endverbatim
 283 *>
 284 *> \param[in] LWORK
 285 *> \verbatim
 286 *>          LWORK is INTEGER
 287 *>          The dimension of the array WORK. LWORK >= max(1,2*N).
 288 *>          If SENSE = 'E', LWORK >= max(1,4*N).
 289 *>          If SENSE = 'V' or 'B', LWORK >= max(1,2*N*N+2*N).
 290 *>
 291 *>          If LWORK = -1, then a workspace query is assumed; the routine
 292 *>          only calculates the optimal size of the WORK array, returns
 293 *>          this value as the first entry of the WORK array, and no error
 294 *>          message related to LWORK is issued by XERBLA.
 295 *> \endverbatim
 296 *>
 297 *> \param[out] RWORK
 298 *> \verbatim
 299 *>          RWORK is DOUBLE PRECISION array, dimension (lrwork)
 300 *>          lrwork must be at least max(1,6*N) if BALANC = 'S' or 'B',
 301 *>          and at least max(1,2*N) otherwise.
 302 *>          Real workspace.
 303 *> \endverbatim
 304 *>
 305 *> \param[out] IWORK
 306 *> \verbatim
 307 *>          IWORK is INTEGER array, dimension (N+2)
 308 *>          If SENSE = 'E', IWORK is not referenced.
 309 *> \endverbatim
 310 *>
 311 *> \param[out] BWORK
 312 *> \verbatim
 313 *>          BWORK is LOGICAL array, dimension (N)
 314 *>          If SENSE = 'N', BWORK is not referenced.
 315 *> \endverbatim
 316 *>
 317 *> \param[out] INFO
 318 *> \verbatim
 319 *>          INFO is INTEGER
 320 *>          = 0:  successful exit
 321 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 322 *>          = 1,...,N:
 323 *>                The QZ iteration failed.  No eigenvectors have been
 324 *>                calculated, but ALPHA(j) and BETA(j) should be correct
 325 *>                for j=INFO+1,...,N.
 326 *>          > N:  =N+1: other than QZ iteration failed in ZHGEQZ.
 327 *>                =N+2: error return from ZTGEVC.
 328 *> \endverbatim
 329 *
 330 *  Authors:
 331 *  ========
 332 *
 333 *> \author Univ. of Tennessee
 334 *> \author Univ. of California Berkeley
 335 *> \author Univ. of Colorado Denver
 336 *> \author NAG Ltd.
 337 *
 338 *> \date April 2012
 339 *
 340 *> \ingroup complex16GEeigen
 341 *
 342 *> \par Further Details:
 343 *  =====================
 344 *>
 345 *> \verbatim
 346 *>
 347 *>  Balancing a matrix pair (A,B) includes, first, permuting rows and
 348 *>  columns to isolate eigenvalues, second, applying diagonal similarity
 349 *>  transformation to the rows and columns to make the rows and columns
 350 *>  as close in norm as possible. The computed reciprocal condition
 351 *>  numbers correspond to the balanced matrix. Permuting rows and columns
 352 *>  will not change the condition numbers (in exact arithmetic) but
 353 *>  diagonal scaling will.  For further explanation of balancing, see
 354 *>  section 4.11.1.2 of LAPACK Users' Guide.
 355 *>
 356 *>  An approximate error bound on the chordal distance between the i-th
 357 *>  computed generalized eigenvalue w and the corresponding exact
 358 *>  eigenvalue lambda is
 359 *>
 360 *>       chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
 361 *>
 362 *>  An approximate error bound for the angle between the i-th computed
 363 *>  eigenvector VL(i) or VR(i) is given by
 364 *>
 365 *>       EPS * norm(ABNRM, BBNRM) / DIF(i).
 366 *>
 367 *>  For further explanation of the reciprocal condition numbers RCONDE
 368 *>  and RCONDV, see section 4.11 of LAPACK User's Guide.
 369 *> \endverbatim
 370 *>
 371 *  =====================================================================
 372       SUBROUTINE ZGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
 373      $                   ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI,
 374      $                   LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV,
 375      $                   WORK, LWORK, RWORK, IWORK, BWORK, INFO )
 376 *
 377 *  -- LAPACK driver routine (version 3.4.1) --
 378 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 379 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 380 *     April 2012
 381 *
 382 *     .. Scalar Arguments ..
 383       CHARACTER          BALANC, JOBVL, JOBVR, SENSE
 384       INTEGER            IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
 385       DOUBLE PRECISION   ABNRM, BBNRM
 386 *     ..
 387 *     .. Array Arguments ..
 388       LOGICAL            BWORK( * )
 389       INTEGER            IWORK( * )
 390       DOUBLE PRECISION   LSCALE( * ), RCONDE( * ), RCONDV( * ),
 391      $                   RSCALE( * ), RWORK( * )
 392       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
 393      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
 394      $                   WORK( * )
 395 *     ..
 396 *
 397 *  =====================================================================
 398 *
 399 *     .. Parameters ..
 400       DOUBLE PRECISION   ZERO, ONE
 401       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 402       COMPLEX*16         CZERO, CONE
 403       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
 404      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
 405 *     ..
 406 *     .. Local Scalars ..
 407       LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
 408      $                   WANTSB, WANTSE, WANTSN, WANTSV
 409       CHARACTER          CHTEMP
 410       INTEGER            I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
 411      $                   ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK, MINWRK
 412       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
 413      $                   SMLNUM, TEMP
 414       COMPLEX*16         X
 415 *     ..
 416 *     .. Local Arrays ..
 417       LOGICAL            LDUMMA( 1 )
 418 *     ..
 419 *     .. External Subroutines ..
 420       EXTERNAL           DLABAD, DLASCL, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL,
 421      $                   ZGGHRD, ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC,
 422      $                   ZTGSNA, ZUNGQR, ZUNMQR
 423 *     ..
 424 *     .. External Functions ..
 425       LOGICAL            LSAME
 426       INTEGER            ILAENV
 427       DOUBLE PRECISION   DLAMCH, ZLANGE
 428       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
 429 *     ..
 430 *     .. Intrinsic Functions ..
 431       INTRINSIC          ABS, DBLE, DIMAG, MAX, SQRT
 432 *     ..
 433 *     .. Statement Functions ..
 434       DOUBLE PRECISION   ABS1
 435 *     ..
 436 *     .. Statement Function definitions ..
 437       ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
 438 *     ..
 439 *     .. Executable Statements ..
 440 *
 441 *     Decode the input arguments
 442 *
 443       IF( LSAME( JOBVL, 'N' ) ) THEN
 444          IJOBVL = 1
 445          ILVL = .FALSE.
 446       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
 447          IJOBVL = 2
 448          ILVL = .TRUE.
 449       ELSE
 450          IJOBVL = -1
 451          ILVL = .FALSE.
 452       END IF
 453 *
 454       IF( LSAME( JOBVR, 'N' ) ) THEN
 455          IJOBVR = 1
 456          ILVR = .FALSE.
 457       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
 458          IJOBVR = 2
 459          ILVR = .TRUE.
 460       ELSE
 461          IJOBVR = -1
 462          ILVR = .FALSE.
 463       END IF
 464       ILV = ILVL .OR. ILVR
 465 *
 466       NOSCL  = LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'P' )
 467       WANTSN = LSAME( SENSE, 'N' )
 468       WANTSE = LSAME( SENSE, 'E' )
 469       WANTSV = LSAME( SENSE, 'V' )
 470       WANTSB = LSAME( SENSE, 'B' )
 471 *
 472 *     Test the input arguments
 473 *
 474       INFO = 0
 475       LQUERY = ( LWORK.EQ.-1 )
 476       IF( .NOT.( NOSCL .OR. LSAME( BALANC,'S' ) .OR.
 477      $    LSAME( BALANC, 'B' ) ) ) THEN
 478          INFO = -1
 479       ELSE IF( IJOBVL.LE.0 ) THEN
 480          INFO = -2
 481       ELSE IF( IJOBVR.LE.0 ) THEN
 482          INFO = -3
 483       ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
 484      $          THEN
 485          INFO = -4
 486       ELSE IF( N.LT.0 ) THEN
 487          INFO = -5
 488       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 489          INFO = -7
 490       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 491          INFO = -9
 492       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
 493          INFO = -13
 494       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
 495          INFO = -15
 496       END IF
 497 *
 498 *     Compute workspace
 499 *      (Note: Comments in the code beginning "Workspace:" describe the
 500 *       minimal amount of workspace needed at that point in the code,
 501 *       as well as the preferred amount for good performance.
 502 *       NB refers to the optimal block size for the immediately
 503 *       following subroutine, as returned by ILAENV. The workspace is
 504 *       computed assuming ILO = 1 and IHI = N, the worst case.)
 505 *
 506       IF( INFO.EQ.0 ) THEN
 507          IF( N.EQ.0 ) THEN
 508             MINWRK = 1
 509             MAXWRK = 1
 510          ELSE
 511             MINWRK = 2*N
 512             IF( WANTSE ) THEN
 513                MINWRK = 4*N
 514             ELSE IF( WANTSV .OR. WANTSB ) THEN
 515                MINWRK = 2*N*( N + 1)
 516             END IF
 517             MAXWRK = MINWRK
 518             MAXWRK = MAX( MAXWRK,
 519      $                    N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
 520             MAXWRK = MAX( MAXWRK,
 521      $                    N + N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) )
 522             IF( ILVL ) THEN
 523                MAXWRK = MAX( MAXWRK, N +
 524      $                       N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, 0 ) )
 525             END IF
 526          END IF
 527          WORK( 1 ) = MAXWRK
 528 *
 529          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
 530             INFO = -25
 531          END IF
 532       END IF
 533 *
 534       IF( INFO.NE.0 ) THEN
 535          CALL XERBLA( 'ZGGEVX', -INFO )
 536          RETURN
 537       ELSE IF( LQUERY ) THEN
 538          RETURN
 539       END IF
 540 *
 541 *     Quick return if possible
 542 *
 543       IF( N.EQ.0 )
 544      $   RETURN
 545 *
 546 *     Get machine constants
 547 *
 548       EPS = DLAMCH( 'P' )
 549       SMLNUM = DLAMCH( 'S' )
 550       BIGNUM = ONE / SMLNUM
 551       CALL DLABAD( SMLNUM, BIGNUM )
 552       SMLNUM = SQRT( SMLNUM ) / EPS
 553       BIGNUM = ONE / SMLNUM
 554 *
 555 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 556 *
 557       ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
 558       ILASCL = .FALSE.
 559       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 560          ANRMTO = SMLNUM
 561          ILASCL = .TRUE.
 562       ELSE IF( ANRM.GT.BIGNUM ) THEN
 563          ANRMTO = BIGNUM
 564          ILASCL = .TRUE.
 565       END IF
 566       IF( ILASCL )
 567      $   CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
 568 *
 569 *     Scale B if max element outside range [SMLNUM,BIGNUM]
 570 *
 571       BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
 572       ILBSCL = .FALSE.
 573       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 574          BNRMTO = SMLNUM
 575          ILBSCL = .TRUE.
 576       ELSE IF( BNRM.GT.BIGNUM ) THEN
 577          BNRMTO = BIGNUM
 578          ILBSCL = .TRUE.
 579       END IF
 580       IF( ILBSCL )
 581      $   CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
 582 *
 583 *     Permute and/or balance the matrix pair (A,B)
 584 *     (Real Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
 585 *
 586       CALL ZGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
 587      $             RWORK, IERR )
 588 *
 589 *     Compute ABNRM and BBNRM
 590 *
 591       ABNRM = ZLANGE( '1', N, N, A, LDA, RWORK( 1 ) )
 592       IF( ILASCL ) THEN
 593          RWORK( 1 ) = ABNRM
 594          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, RWORK( 1 ), 1,
 595      $                IERR )
 596          ABNRM = RWORK( 1 )
 597       END IF
 598 *
 599       BBNRM = ZLANGE( '1', N, N, B, LDB, RWORK( 1 ) )
 600       IF( ILBSCL ) THEN
 601          RWORK( 1 ) = BBNRM
 602          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, RWORK( 1 ), 1,
 603      $                IERR )
 604          BBNRM = RWORK( 1 )
 605       END IF
 606 *
 607 *     Reduce B to triangular form (QR decomposition of B)
 608 *     (Complex Workspace: need N, prefer N*NB )
 609 *
 610       IROWS = IHI + 1 - ILO
 611       IF( ILV .OR. .NOT.WANTSN ) THEN
 612          ICOLS = N + 1 - ILO
 613       ELSE
 614          ICOLS = IROWS
 615       END IF
 616       ITAU = 1
 617       IWRK = ITAU + IROWS
 618       CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
 619      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
 620 *
 621 *     Apply the unitary transformation to A
 622 *     (Complex Workspace: need N, prefer N*NB)
 623 *
 624       CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
 625      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
 626      $             LWORK+1-IWRK, IERR )
 627 *
 628 *     Initialize VL and/or VR
 629 *     (Workspace: need N, prefer N*NB)
 630 *
 631       IF( ILVL ) THEN
 632          CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
 633          IF( IROWS.GT.1 ) THEN
 634             CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
 635      $                   VL( ILO+1, ILO ), LDVL )
 636          END IF
 637          CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
 638      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
 639       END IF
 640 *
 641       IF( ILVR )
 642      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
 643 *
 644 *     Reduce to generalized Hessenberg form
 645 *     (Workspace: none needed)
 646 *
 647       IF( ILV .OR. .NOT.WANTSN ) THEN
 648 *
 649 *        Eigenvectors requested -- work on whole matrix.
 650 *
 651          CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
 652      $                LDVL, VR, LDVR, IERR )
 653       ELSE
 654          CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
 655      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
 656       END IF
 657 *
 658 *     Perform QZ algorithm (Compute eigenvalues, and optionally, the
 659 *     Schur forms and Schur vectors)
 660 *     (Complex Workspace: need N)
 661 *     (Real Workspace: need N)
 662 *
 663       IWRK = ITAU
 664       IF( ILV .OR. .NOT.WANTSN ) THEN
 665          CHTEMP = 'S'
 666       ELSE
 667          CHTEMP = 'E'
 668       END IF
 669 *
 670       CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
 671      $             ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
 672      $             LWORK+1-IWRK, RWORK, IERR )
 673       IF( IERR.NE.0 ) THEN
 674          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
 675             INFO = IERR
 676          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
 677             INFO = IERR - N
 678          ELSE
 679             INFO = N + 1
 680          END IF
 681          GO TO 90
 682       END IF
 683 *
 684 *     Compute Eigenvectors and estimate condition numbers if desired
 685 *     ZTGEVC: (Complex Workspace: need 2*N )
 686 *             (Real Workspace:    need 2*N )
 687 *     ZTGSNA: (Complex Workspace: need 2*N*N if SENSE='V' or 'B')
 688 *             (Integer Workspace: need N+2 )
 689 *
 690       IF( ILV .OR. .NOT.WANTSN ) THEN
 691          IF( ILV ) THEN
 692             IF( ILVL ) THEN
 693                IF( ILVR ) THEN
 694                   CHTEMP = 'B'
 695                ELSE
 696                   CHTEMP = 'L'
 697                END IF
 698             ELSE
 699                CHTEMP = 'R'
 700             END IF
 701 *
 702             CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
 703      $                   LDVL, VR, LDVR, N, IN, WORK( IWRK ), RWORK,
 704      $                   IERR )
 705             IF( IERR.NE.0 ) THEN
 706                INFO = N + 2
 707                GO TO 90
 708             END IF
 709          END IF
 710 *
 711          IF( .NOT.WANTSN ) THEN
 712 *
 713 *           compute eigenvectors (DTGEVC) and estimate condition
 714 *           numbers (DTGSNA). Note that the definition of the condition
 715 *           number is not invariant under transformation (u,v) to
 716 *           (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
 717 *           Schur form (S,T), Q and Z are orthogonal matrices. In order
 718 *           to avoid using extra 2*N*N workspace, we have to
 719 *           re-calculate eigenvectors and estimate the condition numbers
 720 *           one at a time.
 721 *
 722             DO 20 I = 1, N
 723 *
 724                DO 10 J = 1, N
 725                   BWORK( J ) = .FALSE.
 726    10          CONTINUE
 727                BWORK( I ) = .TRUE.
 728 *
 729                IWRK = N + 1
 730                IWRK1 = IWRK + N
 731 *
 732                IF( WANTSE .OR. WANTSB ) THEN
 733                   CALL ZTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
 734      $                         WORK( 1 ), N, WORK( IWRK ), N, 1, M,
 735      $                         WORK( IWRK1 ), RWORK, IERR )
 736                   IF( IERR.NE.0 ) THEN
 737                      INFO = N + 2
 738                      GO TO 90
 739                   END IF
 740                END IF
 741 *
 742                CALL ZTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
 743      $                      WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
 744      $                      RCONDV( I ), 1, M, WORK( IWRK1 ),
 745      $                      LWORK-IWRK1+1, IWORK, IERR )
 746 *
 747    20       CONTINUE
 748          END IF
 749       END IF
 750 *
 751 *     Undo balancing on VL and VR and normalization
 752 *     (Workspace: none needed)
 753 *
 754       IF( ILVL ) THEN
 755          CALL ZGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
 756      $                LDVL, IERR )
 757 *
 758          DO 50 JC = 1, N
 759             TEMP = ZERO
 760             DO 30 JR = 1, N
 761                TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
 762    30       CONTINUE
 763             IF( TEMP.LT.SMLNUM )
 764      $         GO TO 50
 765             TEMP = ONE / TEMP
 766             DO 40 JR = 1, N
 767                VL( JR, JC ) = VL( JR, JC )*TEMP
 768    40       CONTINUE
 769    50    CONTINUE
 770       END IF
 771 *
 772       IF( ILVR ) THEN
 773          CALL ZGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
 774      $                LDVR, IERR )
 775          DO 80 JC = 1, N
 776             TEMP = ZERO
 777             DO 60 JR = 1, N
 778                TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
 779    60       CONTINUE
 780             IF( TEMP.LT.SMLNUM )
 781      $         GO TO 80
 782             TEMP = ONE / TEMP
 783             DO 70 JR = 1, N
 784                VR( JR, JC ) = VR( JR, JC )*TEMP
 785    70       CONTINUE
 786    80    CONTINUE
 787       END IF
 788 *
 789 *     Undo scaling if necessary
 790 *
 791    90 CONTINUE
 792 *
 793       IF( ILASCL )
 794      $   CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
 795 *
 796       IF( ILBSCL )
 797      $   CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
 798 *
 799       WORK( 1 ) = MAXWRK
 800       RETURN
 801 *
 802 *     End of ZGGEVX
 803 *
 804       END