Return to ztrevc.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1: *> \brief \b ZTREVC
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZTREVC + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrevc.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrevc.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrevc.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
   22: *                          LDVR, MM, M, WORK, RWORK, INFO )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          HOWMNY, SIDE
   26: *       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       LOGICAL            SELECT( * )
   30: *       DOUBLE PRECISION   RWORK( * )
   31: *       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
   32: *      $                   WORK( * )
   33: *       ..
   34: *  
   35: *
   36: *> \par Purpose:
   37: *  =============
   38: *>
   39: *> \verbatim
   40: *>
   41: *> ZTREVC computes some or all of the right and/or left eigenvectors of
   42: *> a complex upper triangular matrix T.
   43: *> Matrices of this type are produced by the Schur factorization of
   44: *> a complex general matrix:  A = Q*T*Q**H, as computed by ZHSEQR.
   45: *> 
   46: *> The right eigenvector x and the left eigenvector y of T corresponding
   47: *> to an eigenvalue w are defined by:
   48: *> 
   49: *>              T*x = w*x,     (y**H)*T = w*(y**H)
   50: *> 
   51: *> where y**H denotes the conjugate transpose of the vector y.
   52: *> The eigenvalues are not input to this routine, but are read directly
   53: *> from the diagonal of T.
   54: *> 
   55: *> This routine returns the matrices X and/or Y of right and left
   56: *> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
   57: *> input matrix.  If Q is the unitary factor that reduces a matrix A to
   58: *> Schur form T, then Q*X and Q*Y are the matrices of right and left
   59: *> eigenvectors of A.
   60: *> \endverbatim
   61: *
   62: *  Arguments:
   63: *  ==========
   64: *
   65: *> \param[in] SIDE
   66: *> \verbatim
   67: *>          SIDE is CHARACTER*1
   68: *>          = 'R':  compute right eigenvectors only;
   69: *>          = 'L':  compute left eigenvectors only;
   70: *>          = 'B':  compute both right and left eigenvectors.
   71: *> \endverbatim
   72: *>
   73: *> \param[in] HOWMNY
   74: *> \verbatim
   75: *>          HOWMNY is CHARACTER*1
   76: *>          = 'A':  compute all right and/or left eigenvectors;
   77: *>          = 'B':  compute all right and/or left eigenvectors,
   78: *>                  backtransformed using the matrices supplied in
   79: *>                  VR and/or VL;
   80: *>          = 'S':  compute selected right and/or left eigenvectors,
   81: *>                  as indicated by the logical array SELECT.
   82: *> \endverbatim
   83: *>
   84: *> \param[in] SELECT
   85: *> \verbatim
   86: *>          SELECT is LOGICAL array, dimension (N)
   87: *>          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
   88: *>          computed.
   89: *>          The eigenvector corresponding to the j-th eigenvalue is
   90: *>          computed if SELECT(j) = .TRUE..
   91: *>          Not referenced if HOWMNY = 'A' or 'B'.
   92: *> \endverbatim
   93: *>
   94: *> \param[in] N
   95: *> \verbatim
   96: *>          N is INTEGER
   97: *>          The order of the matrix T. N >= 0.
   98: *> \endverbatim
   99: *>
  100: *> \param[in,out] T
  101: *> \verbatim
  102: *>          T is COMPLEX*16 array, dimension (LDT,N)
  103: *>          The upper triangular matrix T.  T is modified, but restored
  104: *>          on exit.
  105: *> \endverbatim
  106: *>
  107: *> \param[in] LDT
  108: *> \verbatim
  109: *>          LDT is INTEGER
  110: *>          The leading dimension of the array T. LDT >= max(1,N).
  111: *> \endverbatim
  112: *>
  113: *> \param[in,out] VL
  114: *> \verbatim
  115: *>          VL is COMPLEX*16 array, dimension (LDVL,MM)
  116: *>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
  117: *>          contain an N-by-N matrix Q (usually the unitary matrix Q of
  118: *>          Schur vectors returned by ZHSEQR).
  119: *>          On exit, if SIDE = 'L' or 'B', VL contains:
  120: *>          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
  121: *>          if HOWMNY = 'B', the matrix Q*Y;
  122: *>          if HOWMNY = 'S', the left eigenvectors of T specified by
  123: *>                           SELECT, stored consecutively in the columns
  124: *>                           of VL, in the same order as their
  125: *>                           eigenvalues.
  126: *>          Not referenced if SIDE = 'R'.
  127: *> \endverbatim
  128: *>
  129: *> \param[in] LDVL
  130: *> \verbatim
  131: *>          LDVL is INTEGER
  132: *>          The leading dimension of the array VL.  LDVL >= 1, and if
  133: *>          SIDE = 'L' or 'B', LDVL >= N.
  134: *> \endverbatim
  135: *>
  136: *> \param[in,out] VR
  137: *> \verbatim
  138: *>          VR is COMPLEX*16 array, dimension (LDVR,MM)
  139: *>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
  140: *>          contain an N-by-N matrix Q (usually the unitary matrix Q of
  141: *>          Schur vectors returned by ZHSEQR).
  142: *>          On exit, if SIDE = 'R' or 'B', VR contains:
  143: *>          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
  144: *>          if HOWMNY = 'B', the matrix Q*X;
  145: *>          if HOWMNY = 'S', the right eigenvectors of T specified by
  146: *>                           SELECT, stored consecutively in the columns
  147: *>                           of VR, in the same order as their
  148: *>                           eigenvalues.
  149: *>          Not referenced if SIDE = 'L'.
  150: *> \endverbatim
  151: *>
  152: *> \param[in] LDVR
  153: *> \verbatim
  154: *>          LDVR is INTEGER
  155: *>          The leading dimension of the array VR.  LDVR >= 1, and if
  156: *>          SIDE = 'R' or 'B'; LDVR >= N.
  157: *> \endverbatim
  158: *>
  159: *> \param[in] MM
  160: *> \verbatim
  161: *>          MM is INTEGER
  162: *>          The number of columns in the arrays VL and/or VR. MM >= M.
  163: *> \endverbatim
  164: *>
  165: *> \param[out] M
  166: *> \verbatim
  167: *>          M is INTEGER
  168: *>          The number of columns in the arrays VL and/or VR actually
  169: *>          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
  170: *>          is set to N.  Each selected eigenvector occupies one
  171: *>          column.
  172: *> \endverbatim
  173: *>
  174: *> \param[out] WORK
  175: *> \verbatim
  176: *>          WORK is COMPLEX*16 array, dimension (2*N)
  177: *> \endverbatim
  178: *>
  179: *> \param[out] RWORK
  180: *> \verbatim
  181: *>          RWORK is DOUBLE PRECISION array, dimension (N)
  182: *> \endverbatim
  183: *>
  184: *> \param[out] INFO
  185: *> \verbatim
  186: *>          INFO is INTEGER
  187: *>          = 0:  successful exit
  188: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  189: *> \endverbatim
  190: *
  191: *  Authors:
  192: *  ========
  193: *
  194: *> \author Univ. of Tennessee 
  195: *> \author Univ. of California Berkeley 
  196: *> \author Univ. of Colorado Denver 
  197: *> \author NAG Ltd. 
  198: *
  199: *> \date November 2011
  200: *
  201: *> \ingroup complex16OTHERcomputational
  202: *
  203: *> \par Further Details:
  204: *  =====================
  205: *>
  206: *> \verbatim
  207: *>
  208: *>  The algorithm used in this program is basically backward (forward)
  209: *>  substitution, with scaling to make the the code robust against
  210: *>  possible overflow.
  211: *>
  212: *>  Each eigenvector is normalized so that the element of largest
  213: *>  magnitude has magnitude 1; here the magnitude of a complex number
  214: *>  (x,y) is taken to be |x| + |y|.
  215: *> \endverbatim
  216: *>
  217: *  =====================================================================
  218:       SUBROUTINE ZTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
  219:      $                   LDVR, MM, M, WORK, RWORK, INFO )
  220: *
  221: *  -- LAPACK computational routine (version 3.4.0) --
  222: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  223: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  224: *     November 2011
  225: *
  226: *     .. Scalar Arguments ..
  227:       CHARACTER          HOWMNY, SIDE
  228:       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
  229: *     ..
  230: *     .. Array Arguments ..
  231:       LOGICAL            SELECT( * )
  232:       DOUBLE PRECISION   RWORK( * )
  233:       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
  234:      $                   WORK( * )
  235: *     ..
  236: *
  237: *  =====================================================================
  238: *
  239: *     .. Parameters ..
  240:       DOUBLE PRECISION   ZERO, ONE
  241:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  242:       COMPLEX*16         CMZERO, CMONE
  243:       PARAMETER          ( CMZERO = ( 0.0D+0, 0.0D+0 ),
  244:      $                   CMONE = ( 1.0D+0, 0.0D+0 ) )
  245: *     ..
  246: *     .. Local Scalars ..
  247:       LOGICAL            ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
  248:       INTEGER            I, II, IS, J, K, KI
  249:       DOUBLE PRECISION   OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
  250:       COMPLEX*16         CDUM
  251: *     ..
  252: *     .. External Functions ..
  253:       LOGICAL            LSAME
  254:       INTEGER            IZAMAX
  255:       DOUBLE PRECISION   DLAMCH, DZASUM
  256:       EXTERNAL           LSAME, IZAMAX, DLAMCH, DZASUM
  257: *     ..
  258: *     .. External Subroutines ..
  259:       EXTERNAL           XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLATRS
  260: *     ..
  261: *     .. Intrinsic Functions ..
  262:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX
  263: *     ..
  264: *     .. Statement Functions ..
  265:       DOUBLE PRECISION   CABS1
  266: *     ..
  267: *     .. Statement Function definitions ..
  268:       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
  269: *     ..
  270: *     .. Executable Statements ..
  271: *
  272: *     Decode and test the input parameters
  273: *
  274:       BOTHV = LSAME( SIDE, 'B' )
  275:       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
  276:       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
  277: *
  278:       ALLV = LSAME( HOWMNY, 'A' )
  279:       OVER = LSAME( HOWMNY, 'B' )
  280:       SOMEV = LSAME( HOWMNY, 'S' )
  281: *
  282: *     Set M to the number of columns required to store the selected
  283: *     eigenvectors.
  284: *
  285:       IF( SOMEV ) THEN
  286:          M = 0
  287:          DO 10 J = 1, N
  288:             IF( SELECT( J ) )
  289:      $         M = M + 1
  290:    10    CONTINUE
  291:       ELSE
  292:          M = N
  293:       END IF
  294: *
  295:       INFO = 0
  296:       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
  297:          INFO = -1
  298:       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
  299:          INFO = -2
  300:       ELSE IF( N.LT.0 ) THEN
  301:          INFO = -4
  302:       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
  303:          INFO = -6
  304:       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
  305:          INFO = -8
  306:       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
  307:          INFO = -10
  308:       ELSE IF( MM.LT.M ) THEN
  309:          INFO = -11
  310:       END IF
  311:       IF( INFO.NE.0 ) THEN
  312:          CALL XERBLA( 'ZTREVC', -INFO )
  313:          RETURN
  314:       END IF
  315: *
  316: *     Quick return if possible.
  317: *
  318:       IF( N.EQ.0 )
  319:      $   RETURN
  320: *
  321: *     Set the constants to control overflow.
  322: *
  323:       UNFL = DLAMCH( 'Safe minimum' )
  324:       OVFL = ONE / UNFL
  325:       CALL DLABAD( UNFL, OVFL )
  326:       ULP = DLAMCH( 'Precision' )
  327:       SMLNUM = UNFL*( N / ULP )
  328: *
  329: *     Store the diagonal elements of T in working array WORK.
  330: *
  331:       DO 20 I = 1, N
  332:          WORK( I+N ) = T( I, I )
  333:    20 CONTINUE
  334: *
  335: *     Compute 1-norm of each column of strictly upper triangular
  336: *     part of T to control overflow in triangular solver.
  337: *
  338:       RWORK( 1 ) = ZERO
  339:       DO 30 J = 2, N
  340:          RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 )
  341:    30 CONTINUE
  342: *
  343:       IF( RIGHTV ) THEN
  344: *
  345: *        Compute right eigenvectors.
  346: *
  347:          IS = M
  348:          DO 80 KI = N, 1, -1
  349: *
  350:             IF( SOMEV ) THEN
  351:                IF( .NOT.SELECT( KI ) )
  352:      $            GO TO 80
  353:             END IF
  354:             SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
  355: *
  356:             WORK( 1 ) = CMONE
  357: *
  358: *           Form right-hand side.
  359: *
  360:             DO 40 K = 1, KI - 1
  361:                WORK( K ) = -T( K, KI )
  362:    40       CONTINUE
  363: *
  364: *           Solve the triangular system:
  365: *              (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK.
  366: *
  367:             DO 50 K = 1, KI - 1
  368:                T( K, K ) = T( K, K ) - T( KI, KI )
  369:                IF( CABS1( T( K, K ) ).LT.SMIN )
  370:      $            T( K, K ) = SMIN
  371:    50       CONTINUE
  372: *
  373:             IF( KI.GT.1 ) THEN
  374:                CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
  375:      $                      KI-1, T, LDT, WORK( 1 ), SCALE, RWORK,
  376:      $                      INFO )
  377:                WORK( KI ) = SCALE
  378:             END IF
  379: *
  380: *           Copy the vector x or Q*x to VR and normalize.
  381: *
  382:             IF( .NOT.OVER ) THEN
  383:                CALL ZCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 )
  384: *
  385:                II = IZAMAX( KI, VR( 1, IS ), 1 )
  386:                REMAX = ONE / CABS1( VR( II, IS ) )
  387:                CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 )
  388: *
  389:                DO 60 K = KI + 1, N
  390:                   VR( K, IS ) = CMZERO
  391:    60          CONTINUE
  392:             ELSE
  393:                IF( KI.GT.1 )
  394:      $            CALL ZGEMV( 'N', N, KI-1, CMONE, VR, LDVR, WORK( 1 ),
  395:      $                        1, DCMPLX( SCALE ), VR( 1, KI ), 1 )
  396: *
  397:                II = IZAMAX( N, VR( 1, KI ), 1 )
  398:                REMAX = ONE / CABS1( VR( II, KI ) )
  399:                CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 )
  400:             END IF
  401: *
  402: *           Set back the original diagonal elements of T.
  403: *
  404:             DO 70 K = 1, KI - 1
  405:                T( K, K ) = WORK( K+N )
  406:    70       CONTINUE
  407: *
  408:             IS = IS - 1
  409:    80    CONTINUE
  410:       END IF
  411: *
  412:       IF( LEFTV ) THEN
  413: *
  414: *        Compute left eigenvectors.
  415: *
  416:          IS = 1
  417:          DO 130 KI = 1, N
  418: *
  419:             IF( SOMEV ) THEN
  420:                IF( .NOT.SELECT( KI ) )
  421:      $            GO TO 130
  422:             END IF
  423:             SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
  424: *
  425:             WORK( N ) = CMONE
  426: *
  427: *           Form right-hand side.
  428: *
  429:             DO 90 K = KI + 1, N
  430:                WORK( K ) = -DCONJG( T( KI, K ) )
  431:    90       CONTINUE
  432: *
  433: *           Solve the triangular system:
  434: *              (T(KI+1:N,KI+1:N) - T(KI,KI))**H * X = SCALE*WORK.
  435: *
  436:             DO 100 K = KI + 1, N
  437:                T( K, K ) = T( K, K ) - T( KI, KI )
  438:                IF( CABS1( T( K, K ) ).LT.SMIN )
  439:      $            T( K, K ) = SMIN
  440:   100       CONTINUE
  441: *
  442:             IF( KI.LT.N ) THEN
  443:                CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
  444:      $                      'Y', N-KI, T( KI+1, KI+1 ), LDT,
  445:      $                      WORK( KI+1 ), SCALE, RWORK, INFO )
  446:                WORK( KI ) = SCALE
  447:             END IF
  448: *
  449: *           Copy the vector x or Q*x to VL and normalize.
  450: *
  451:             IF( .NOT.OVER ) THEN
  452:                CALL ZCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 )
  453: *
  454:                II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
  455:                REMAX = ONE / CABS1( VL( II, IS ) )
  456:                CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
  457: *
  458:                DO 110 K = 1, KI - 1
  459:                   VL( K, IS ) = CMZERO
  460:   110          CONTINUE
  461:             ELSE
  462:                IF( KI.LT.N )
  463:      $            CALL ZGEMV( 'N', N, N-KI, CMONE, VL( 1, KI+1 ), LDVL,
  464:      $                        WORK( KI+1 ), 1, DCMPLX( SCALE ),
  465:      $                        VL( 1, KI ), 1 )
  466: *
  467:                II = IZAMAX( N, VL( 1, KI ), 1 )
  468:                REMAX = ONE / CABS1( VL( II, KI ) )
  469:                CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 )
  470:             END IF
  471: *
  472: *           Set back the original diagonal elements of T.
  473: *
  474:             DO 120 K = KI + 1, N
  475:                T( K, K ) = WORK( K+N )
  476:   120       CONTINUE
  477: *
  478:             IS = IS + 1
  479:   130    CONTINUE
  480:       END IF
  481: *
  482:       RETURN
  483: *
  484: *     End of ZTREVC
  485: *
  486:       END