Return to dgeevx.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

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