Return to dhsein.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *> \brief \b DHSEIN
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DHSEIN + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhsein.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhsein.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhsein.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI,
   22: *                          VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL,
   23: *                          IFAILR, INFO )
   24: * 
   25: *       .. Scalar Arguments ..
   26: *       CHARACTER          EIGSRC, INITV, SIDE
   27: *       INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
   28: *       ..
   29: *       .. Array Arguments ..
   30: *       LOGICAL            SELECT( * )
   31: *       INTEGER            IFAILL( * ), IFAILR( * )
   32: *       DOUBLE PRECISION   H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
   33: *      $                   WI( * ), WORK( * ), WR( * )
   34: *       ..
   35: *  
   36: *
   37: *> \par Purpose:
   38: *  =============
   39: *>
   40: *> \verbatim
   41: *>
   42: *> DHSEIN uses inverse iteration to find specified right and/or left
   43: *> eigenvectors of a real upper Hessenberg matrix H.
   44: *>
   45: *> The right eigenvector x and the left eigenvector y of the matrix H
   46: *> corresponding to an eigenvalue w are defined by:
   47: *>
   48: *>              H * x = w * x,     y**h * H = w * y**h
   49: *>
   50: *> where y**h denotes the conjugate transpose of the vector y.
   51: *> \endverbatim
   52: *
   53: *  Arguments:
   54: *  ==========
   55: *
   56: *> \param[in] SIDE
   57: *> \verbatim
   58: *>          SIDE is CHARACTER*1
   59: *>          = 'R': compute right eigenvectors only;
   60: *>          = 'L': compute left eigenvectors only;
   61: *>          = 'B': compute both right and left eigenvectors.
   62: *> \endverbatim
   63: *>
   64: *> \param[in] EIGSRC
   65: *> \verbatim
   66: *>          EIGSRC is CHARACTER*1
   67: *>          Specifies the source of eigenvalues supplied in (WR,WI):
   68: *>          = 'Q': the eigenvalues were found using DHSEQR; thus, if
   69: *>                 H has zero subdiagonal elements, and so is
   70: *>                 block-triangular, then the j-th eigenvalue can be
   71: *>                 assumed to be an eigenvalue of the block containing
   72: *>                 the j-th row/column.  This property allows DHSEIN to
   73: *>                 perform inverse iteration on just one diagonal block.
   74: *>          = 'N': no assumptions are made on the correspondence
   75: *>                 between eigenvalues and diagonal blocks.  In this
   76: *>                 case, DHSEIN must always perform inverse iteration
   77: *>                 using the whole matrix H.
   78: *> \endverbatim
   79: *>
   80: *> \param[in] INITV
   81: *> \verbatim
   82: *>          INITV is CHARACTER*1
   83: *>          = 'N': no initial vectors are supplied;
   84: *>          = 'U': user-supplied initial vectors are stored in the arrays
   85: *>                 VL and/or VR.
   86: *> \endverbatim
   87: *>
   88: *> \param[in,out] SELECT
   89: *> \verbatim
   90: *>          SELECT is LOGICAL array, dimension (N)
   91: *>          Specifies the eigenvectors to be computed. To select the
   92: *>          real eigenvector corresponding to a real eigenvalue WR(j),
   93: *>          SELECT(j) must be set to .TRUE.. To select the complex
   94: *>          eigenvector corresponding to a complex eigenvalue
   95: *>          (WR(j),WI(j)), with complex conjugate (WR(j+1),WI(j+1)),
   96: *>          either SELECT(j) or SELECT(j+1) or both must be set to
   97: *>          .TRUE.; then on exit SELECT(j) is .TRUE. and SELECT(j+1) is
   98: *>          .FALSE..
   99: *> \endverbatim
  100: *>
  101: *> \param[in] N
  102: *> \verbatim
  103: *>          N is INTEGER
  104: *>          The order of the matrix H.  N >= 0.
  105: *> \endverbatim
  106: *>
  107: *> \param[in] H
  108: *> \verbatim
  109: *>          H is DOUBLE PRECISION array, dimension (LDH,N)
  110: *>          The upper Hessenberg matrix H.
  111: *>          If a NaN is detected in H, the routine will return with INFO=-6.
  112: *> \endverbatim
  113: *>
  114: *> \param[in] LDH
  115: *> \verbatim
  116: *>          LDH is INTEGER
  117: *>          The leading dimension of the array H.  LDH >= max(1,N).
  118: *> \endverbatim
  119: *>
  120: *> \param[in,out] WR
  121: *> \verbatim
  122: *>          WR is DOUBLE PRECISION array, dimension (N)
  123: *> \endverbatim
  124: *>
  125: *> \param[in] WI
  126: *> \verbatim
  127: *>          WI is DOUBLE PRECISION array, dimension (N)
  128: *>
  129: *>          On entry, the real and imaginary parts of the eigenvalues of
  130: *>          H; a complex conjugate pair of eigenvalues must be stored in
  131: *>          consecutive elements of WR and WI.
  132: *>          On exit, WR may have been altered since close eigenvalues
  133: *>          are perturbed slightly in searching for independent
  134: *>          eigenvectors.
  135: *> \endverbatim
  136: *>
  137: *> \param[in,out] VL
  138: *> \verbatim
  139: *>          VL is DOUBLE PRECISION array, dimension (LDVL,MM)
  140: *>          On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
  141: *>          contain starting vectors for the inverse iteration for the
  142: *>          left eigenvectors; the starting vector for each eigenvector
  143: *>          must be in the same column(s) in which the eigenvector will
  144: *>          be stored.
  145: *>          On exit, if SIDE = 'L' or 'B', the left eigenvectors
  146: *>          specified by SELECT will be stored consecutively in the
  147: *>          columns of VL, in the same order as their eigenvalues. A
  148: *>          complex eigenvector corresponding to a complex eigenvalue is
  149: *>          stored in two consecutive columns, the first holding the real
  150: *>          part and the second the imaginary part.
  151: *>          If SIDE = 'R', VL is not referenced.
  152: *> \endverbatim
  153: *>
  154: *> \param[in] LDVL
  155: *> \verbatim
  156: *>          LDVL is INTEGER
  157: *>          The leading dimension of the array VL.
  158: *>          LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
  159: *> \endverbatim
  160: *>
  161: *> \param[in,out] VR
  162: *> \verbatim
  163: *>          VR is DOUBLE PRECISION array, dimension (LDVR,MM)
  164: *>          On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
  165: *>          contain starting vectors for the inverse iteration for the
  166: *>          right eigenvectors; the starting vector for each eigenvector
  167: *>          must be in the same column(s) in which the eigenvector will
  168: *>          be stored.
  169: *>          On exit, if SIDE = 'R' or 'B', the right eigenvectors
  170: *>          specified by SELECT will be stored consecutively in the
  171: *>          columns of VR, in the same order as their eigenvalues. A
  172: *>          complex eigenvector corresponding to a complex eigenvalue is
  173: *>          stored in two consecutive columns, the first holding the real
  174: *>          part and the second the imaginary part.
  175: *>          If SIDE = 'L', VR is not referenced.
  176: *> \endverbatim
  177: *>
  178: *> \param[in] LDVR
  179: *> \verbatim
  180: *>          LDVR is INTEGER
  181: *>          The leading dimension of the array VR.
  182: *>          LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
  183: *> \endverbatim
  184: *>
  185: *> \param[in] MM
  186: *> \verbatim
  187: *>          MM is INTEGER
  188: *>          The number of columns in the arrays VL and/or VR. MM >= M.
  189: *> \endverbatim
  190: *>
  191: *> \param[out] M
  192: *> \verbatim
  193: *>          M is INTEGER
  194: *>          The number of columns in the arrays VL and/or VR required to
  195: *>          store the eigenvectors; each selected real eigenvector
  196: *>          occupies one column and each selected complex eigenvector
  197: *>          occupies two columns.
  198: *> \endverbatim
  199: *>
  200: *> \param[out] WORK
  201: *> \verbatim
  202: *>          WORK is DOUBLE PRECISION array, dimension ((N+2)*N)
  203: *> \endverbatim
  204: *>
  205: *> \param[out] IFAILL
  206: *> \verbatim
  207: *>          IFAILL is INTEGER array, dimension (MM)
  208: *>          If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
  209: *>          eigenvector in the i-th column of VL (corresponding to the
  210: *>          eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
  211: *>          eigenvector converged satisfactorily. If the i-th and (i+1)th
  212: *>          columns of VL hold a complex eigenvector, then IFAILL(i) and
  213: *>          IFAILL(i+1) are set to the same value.
  214: *>          If SIDE = 'R', IFAILL is not referenced.
  215: *> \endverbatim
  216: *>
  217: *> \param[out] IFAILR
  218: *> \verbatim
  219: *>          IFAILR is INTEGER array, dimension (MM)
  220: *>          If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
  221: *>          eigenvector in the i-th column of VR (corresponding to the
  222: *>          eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
  223: *>          eigenvector converged satisfactorily. If the i-th and (i+1)th
  224: *>          columns of VR hold a complex eigenvector, then IFAILR(i) and
  225: *>          IFAILR(i+1) are set to the same value.
  226: *>          If SIDE = 'L', IFAILR is not referenced.
  227: *> \endverbatim
  228: *>
  229: *> \param[out] INFO
  230: *> \verbatim
  231: *>          INFO is INTEGER
  232: *>          = 0:  successful exit
  233: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  234: *>          > 0:  if INFO = i, i is the number of eigenvectors which
  235: *>                failed to converge; see IFAILL and IFAILR for further
  236: *>                details.
  237: *> \endverbatim
  238: *
  239: *  Authors:
  240: *  ========
  241: *
  242: *> \author Univ. of Tennessee 
  243: *> \author Univ. of California Berkeley 
  244: *> \author Univ. of Colorado Denver 
  245: *> \author NAG Ltd. 
  246: *
  247: *> \date November 2013
  248: *
  249: *> \ingroup doubleOTHERcomputational
  250: *
  251: *> \par Further Details:
  252: *  =====================
  253: *>
  254: *> \verbatim
  255: *>
  256: *>  Each eigenvector is normalized so that the element of largest
  257: *>  magnitude has magnitude 1; here the magnitude of a complex number
  258: *>  (x,y) is taken to be |x|+|y|.
  259: *> \endverbatim
  260: *>
  261: *  =====================================================================
  262:       SUBROUTINE DHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI,
  263:      $                   VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL,
  264:      $                   IFAILR, INFO )
  265: *
  266: *  -- LAPACK computational routine (version 3.5.0) --
  267: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  268: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  269: *     November 2013
  270: *
  271: *     .. Scalar Arguments ..
  272:       CHARACTER          EIGSRC, INITV, SIDE
  273:       INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
  274: *     ..
  275: *     .. Array Arguments ..
  276:       LOGICAL            SELECT( * )
  277:       INTEGER            IFAILL( * ), IFAILR( * )
  278:       DOUBLE PRECISION   H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
  279:      $                   WI( * ), WORK( * ), WR( * )
  280: *     ..
  281: *
  282: *  =====================================================================
  283: *
  284: *     .. Parameters ..
  285:       DOUBLE PRECISION   ZERO, ONE
  286:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  287: *     ..
  288: *     .. Local Scalars ..
  289:       LOGICAL            BOTHV, FROMQR, LEFTV, NOINIT, PAIR, RIGHTV
  290:       INTEGER            I, IINFO, K, KL, KLN, KR, KSI, KSR, LDWORK
  291:       DOUBLE PRECISION   BIGNUM, EPS3, HNORM, SMLNUM, ULP, UNFL, WKI,
  292:      $                   WKR
  293: *     ..
  294: *     .. External Functions ..
  295:       LOGICAL            LSAME, DISNAN
  296:       DOUBLE PRECISION   DLAMCH, DLANHS
  297:       EXTERNAL           LSAME, DLAMCH, DLANHS, DISNAN
  298: *     ..
  299: *     .. External Subroutines ..
  300:       EXTERNAL           DLAEIN, XERBLA
  301: *     ..
  302: *     .. Intrinsic Functions ..
  303:       INTRINSIC          ABS, MAX
  304: *     ..
  305: *     .. Executable Statements ..
  306: *
  307: *     Decode and test the input parameters.
  308: *
  309:       BOTHV = LSAME( SIDE, 'B' )
  310:       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
  311:       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
  312: *
  313:       FROMQR = LSAME( EIGSRC, 'Q' )
  314: *
  315:       NOINIT = LSAME( INITV, 'N' )
  316: *
  317: *     Set M to the number of columns required to store the selected
  318: *     eigenvectors, and standardize the array SELECT.
  319: *
  320:       M = 0
  321:       PAIR = .FALSE.
  322:       DO 10 K = 1, N
  323:          IF( PAIR ) THEN
  324:             PAIR = .FALSE.
  325:             SELECT( K ) = .FALSE.
  326:          ELSE
  327:             IF( WI( K ).EQ.ZERO ) THEN
  328:                IF( SELECT( K ) )
  329:      $            M = M + 1
  330:             ELSE
  331:                PAIR = .TRUE.
  332:                IF( SELECT( K ) .OR. SELECT( K+1 ) ) THEN
  333:                   SELECT( K ) = .TRUE.
  334:                   M = M + 2
  335:                END IF
  336:             END IF
  337:          END IF
  338:    10 CONTINUE
  339: *
  340:       INFO = 0
  341:       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
  342:          INFO = -1
  343:       ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN
  344:          INFO = -2
  345:       ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN
  346:          INFO = -3
  347:       ELSE IF( N.LT.0 ) THEN
  348:          INFO = -5
  349:       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
  350:          INFO = -7
  351:       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
  352:          INFO = -11
  353:       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
  354:          INFO = -13
  355:       ELSE IF( MM.LT.M ) THEN
  356:          INFO = -14
  357:       END IF
  358:       IF( INFO.NE.0 ) THEN
  359:          CALL XERBLA( 'DHSEIN', -INFO )
  360:          RETURN
  361:       END IF
  362: *
  363: *     Quick return if possible.
  364: *
  365:       IF( N.EQ.0 )
  366:      $   RETURN
  367: *
  368: *     Set machine-dependent constants.
  369: *
  370:       UNFL = DLAMCH( 'Safe minimum' )
  371:       ULP = DLAMCH( 'Precision' )
  372:       SMLNUM = UNFL*( N / ULP )
  373:       BIGNUM = ( ONE-ULP ) / SMLNUM
  374: *
  375:       LDWORK = N + 1
  376: *
  377:       KL = 1
  378:       KLN = 0
  379:       IF( FROMQR ) THEN
  380:          KR = 0
  381:       ELSE
  382:          KR = N
  383:       END IF
  384:       KSR = 1
  385: *
  386:       DO 120 K = 1, N
  387:          IF( SELECT( K ) ) THEN
  388: *
  389: *           Compute eigenvector(s) corresponding to W(K).
  390: *
  391:             IF( FROMQR ) THEN
  392: *
  393: *              If affiliation of eigenvalues is known, check whether
  394: *              the matrix splits.
  395: *
  396: *              Determine KL and KR such that 1 <= KL <= K <= KR <= N
  397: *              and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
  398: *              KR = N).
  399: *
  400: *              Then inverse iteration can be performed with the
  401: *              submatrix H(KL:N,KL:N) for a left eigenvector, and with
  402: *              the submatrix H(1:KR,1:KR) for a right eigenvector.
  403: *
  404:                DO 20 I = K, KL + 1, -1
  405:                   IF( H( I, I-1 ).EQ.ZERO )
  406:      $               GO TO 30
  407:    20          CONTINUE
  408:    30          CONTINUE
  409:                KL = I
  410:                IF( K.GT.KR ) THEN
  411:                   DO 40 I = K, N - 1
  412:                      IF( H( I+1, I ).EQ.ZERO )
  413:      $                  GO TO 50
  414:    40             CONTINUE
  415:    50             CONTINUE
  416:                   KR = I
  417:                END IF
  418:             END IF
  419: *
  420:             IF( KL.NE.KLN ) THEN
  421:                KLN = KL
  422: *
  423: *              Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
  424: *              has not ben computed before.
  425: *
  426:                HNORM = DLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, WORK )
  427:                IF( DISNAN( HNORM ) ) THEN
  428:                   INFO = -6
  429:                   RETURN
  430:                ELSE IF( HNORM.GT.ZERO ) THEN
  431:                   EPS3 = HNORM*ULP
  432:                ELSE
  433:                   EPS3 = SMLNUM
  434:                END IF
  435:             END IF
  436: *
  437: *           Perturb eigenvalue if it is close to any previous
  438: *           selected eigenvalues affiliated to the submatrix
  439: *           H(KL:KR,KL:KR). Close roots are modified by EPS3.
  440: *
  441:             WKR = WR( K )
  442:             WKI = WI( K )
  443:    60       CONTINUE
  444:             DO 70 I = K - 1, KL, -1
  445:                IF( SELECT( I ) .AND. ABS( WR( I )-WKR )+
  446:      $             ABS( WI( I )-WKI ).LT.EPS3 ) THEN
  447:                   WKR = WKR + EPS3
  448:                   GO TO 60
  449:                END IF
  450:    70       CONTINUE
  451:             WR( K ) = WKR
  452: *
  453:             PAIR = WKI.NE.ZERO
  454:             IF( PAIR ) THEN
  455:                KSI = KSR + 1
  456:             ELSE
  457:                KSI = KSR
  458:             END IF
  459:             IF( LEFTV ) THEN
  460: *
  461: *              Compute left eigenvector.
  462: *
  463:                CALL DLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
  464:      $                      WKR, WKI, VL( KL, KSR ), VL( KL, KSI ),
  465:      $                      WORK, LDWORK, WORK( N*N+N+1 ), EPS3, SMLNUM,
  466:      $                      BIGNUM, IINFO )
  467:                IF( IINFO.GT.0 ) THEN
  468:                   IF( PAIR ) THEN
  469:                      INFO = INFO + 2
  470:                   ELSE
  471:                      INFO = INFO + 1
  472:                   END IF
  473:                   IFAILL( KSR ) = K
  474:                   IFAILL( KSI ) = K
  475:                ELSE
  476:                   IFAILL( KSR ) = 0
  477:                   IFAILL( KSI ) = 0
  478:                END IF
  479:                DO 80 I = 1, KL - 1
  480:                   VL( I, KSR ) = ZERO
  481:    80          CONTINUE
  482:                IF( PAIR ) THEN
  483:                   DO 90 I = 1, KL - 1
  484:                      VL( I, KSI ) = ZERO
  485:    90             CONTINUE
  486:                END IF
  487:             END IF
  488:             IF( RIGHTV ) THEN
  489: *
  490: *              Compute right eigenvector.
  491: *
  492:                CALL DLAEIN( .TRUE., NOINIT, KR, H, LDH, WKR, WKI,
  493:      $                      VR( 1, KSR ), VR( 1, KSI ), WORK, LDWORK,
  494:      $                      WORK( N*N+N+1 ), EPS3, SMLNUM, BIGNUM,
  495:      $                      IINFO )
  496:                IF( IINFO.GT.0 ) THEN
  497:                   IF( PAIR ) THEN
  498:                      INFO = INFO + 2
  499:                   ELSE
  500:                      INFO = INFO + 1
  501:                   END IF
  502:                   IFAILR( KSR ) = K
  503:                   IFAILR( KSI ) = K
  504:                ELSE
  505:                   IFAILR( KSR ) = 0
  506:                   IFAILR( KSI ) = 0
  507:                END IF
  508:                DO 100 I = KR + 1, N
  509:                   VR( I, KSR ) = ZERO
  510:   100          CONTINUE
  511:                IF( PAIR ) THEN
  512:                   DO 110 I = KR + 1, N
  513:                      VR( I, KSI ) = ZERO
  514:   110             CONTINUE
  515:                END IF
  516:             END IF
  517: *
  518:             IF( PAIR ) THEN
  519:                KSR = KSR + 2
  520:             ELSE
  521:                KSR = KSR + 1
  522:             END IF
  523:          END IF
  524:   120 CONTINUE
  525: *
  526:       RETURN
  527: *
  528: *     End of DHSEIN
  529: *
  530:       END