Annotation of rpl/lapack/lapack/zggevx.f, revision 1.11

1.8       bertrand    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
1.10      bertrand  299: *>          RWORK is DOUBLE PRECISION array, dimension (lrwork)
1.8       bertrand  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: *
1.10      bertrand  338: *> \date April 2012
1.8       bertrand  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: *  =====================================================================
1.1       bertrand  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: *
1.10      bertrand  377: *  -- LAPACK driver routine (version 3.4.1) --
1.1       bertrand  378: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    379: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.10      bertrand  380: *     April 2012
1.1       bertrand  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: *
1.8       bertrand  397: *  =====================================================================
1.1       bertrand  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: *
1.10      bertrand  791:    90 CONTINUE
                    792: *
1.1       bertrand  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

CVSweb interface <joel.bertrand@systella.fr>