Annotation of rpl/lapack/lapack/zgeevx.f, revision 1.17

1.8       bertrand    1: *> \brief <b> ZGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.15      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.8       bertrand    7: *
                      8: *> \htmlonly
1.15      bertrand    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">
1.8       bertrand   15: *> [TXT]</a>
1.15      bertrand   16: *> \endhtmlonly
1.8       bertrand   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 )
1.15      bertrand   24: *
1.8       bertrand   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: *       ..
1.15      bertrand   36: *
1.8       bertrand   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: *
1.15      bertrand  274: *> \author Univ. of Tennessee
                    275: *> \author Univ. of California Berkeley
                    276: *> \author Univ. of Colorado Denver
                    277: *> \author NAG Ltd.
1.8       bertrand  278: *
1.13      bertrand  279: *> \date June 2016
                    280: *
                    281: *  @precisions fortran z -> c
1.8       bertrand  282: *
                    283: *> \ingroup complex16GEeigen
                    284: *
                    285: *  =====================================================================
1.1       bertrand  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 )
1.13      bertrand  289:       implicit none
1.1       bertrand  290: *
1.15      bertrand  291: *  -- LAPACK driver routine (version 3.7.0) --
1.1       bertrand  292: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    293: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.13      bertrand  294: *     June 2016
1.1       bertrand  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
1.13      bertrand  318:       INTEGER            HSWORK, I, ICOND, IERR, ITAU, IWRK, K,
                    319:      $                   LWORK_TREVC, MAXWRK, MINWRK, NOUT
1.1       bertrand  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,
1.13      bertrand  329:      $                   ZGEHRD, ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC3,
1.1       bertrand  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 ..
1.13      bertrand  339:       INTRINSIC          DBLE, DCMPLX, CONJG, AIMAG, MAX, SQRT
1.1       bertrand  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
1.13      bertrand  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 )
1.1       bertrand  398:                CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
                    399:      $                WORK, -1, INFO )
                    400:             ELSE IF( WANTVR ) THEN
1.13      bertrand  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 )
1.1       bertrand  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
1.13      bertrand  417:             HSWORK = INT( WORK(1) )
1.1       bertrand  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: *
1.13      bertrand  575: *     If INFO .NE. 0 from ZHSEQR, then quit
1.1       bertrand  576: *
1.13      bertrand  577:       IF( INFO.NE.0 )
1.1       bertrand  578:      $   GO TO 50
                    579: *
                    580:       IF( WANTVL .OR. WANTVR ) THEN
                    581: *
                    582: *        Compute left and/or right eigenvectors
1.13      bertrand  583: *        (CWorkspace: need 2*N, prefer N + 2*N*NB)
1.1       bertrand  584: *        (RWorkspace: need N)
                    585: *
1.13      bertrand  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 )
1.1       bertrand  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 +
1.13      bertrand  615:      $                      AIMAG( VL( K, I ) )**2
1.1       bertrand  616:    10       CONTINUE
                    617:             K = IDAMAX( N, RWORK, 1 )
1.13      bertrand  618:             TMP = CONJG( VL( K, I ) ) / SQRT( RWORK( K ) )
1.1       bertrand  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 +
1.13      bertrand  638:      $                      AIMAG( VR( K, I ) )**2
1.1       bertrand  639:    30       CONTINUE
                    640:             K = IDAMAX( N, RWORK, 1 )
1.13      bertrand  641:             TMP = CONJG( VR( K, I ) ) / SQRT( RWORK( K ) )
1.1       bertrand  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

CVSweb interface <joel.bertrand@systella.fr>