Return to dhsein.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack vers la version 3.3.0.

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