Return to dlarrk.f CVS log

Up to [local] / rpl / lapack / lapack

En route pour la 4.0.15 !

    1:       SUBROUTINE DLARRK( N, IW, GL, GU,
    2:      $                    D, E2, PIVMIN, RELTOL, W, WERR, INFO)
    3:       IMPLICIT NONE
    4: *
    5: *  -- LAPACK auxiliary 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:       INTEGER   INFO, IW, N
   12:       DOUBLE PRECISION    PIVMIN, RELTOL, GL, GU, W, WERR
   13: *     ..
   14: *     .. Array Arguments ..
   15:       DOUBLE PRECISION   D( * ), E2( * )
   16: *     ..
   17: *
   18: *  Purpose
   19: *  =======
   20: *
   21: *  DLARRK computes one eigenvalue of a symmetric tridiagonal
   22: *  matrix T to suitable accuracy. This is an auxiliary code to be
   23: *  called from DSTEMR.
   24: *
   25: *  To avoid overflow, the matrix must be scaled so that its
   26: *  largest element is no greater than overflow**(1/2) *
   27: *  underflow**(1/4) in absolute value, and for greatest
   28: *  accuracy, it should not be much smaller than that.
   29: *
   30: *  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
   31: *  Matrix", Report CS41, Computer Science Dept., Stanford
   32: *  University, July 21, 1966.
   33: *
   34: *  Arguments
   35: *  =========
   36: *
   37: *  N       (input) INTEGER
   38: *          The order of the tridiagonal matrix T.  N >= 0.
   39: *
   40: *  IW      (input) INTEGER
   41: *          The index of the eigenvalues to be returned.
   42: *
   43: *  GL      (input) DOUBLE PRECISION
   44: *  GU      (input) DOUBLE PRECISION
   45: *          An upper and a lower bound on the eigenvalue.
   46: *
   47: *  D       (input) DOUBLE PRECISION array, dimension (N)
   48: *          The n diagonal elements of the tridiagonal matrix T.
   49: *
   50: *  E2      (input) DOUBLE PRECISION array, dimension (N-1)
   51: *          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
   52: *
   53: *  PIVMIN  (input) DOUBLE PRECISION
   54: *          The minimum pivot allowed in the Sturm sequence for T.
   55: *
   56: *  RELTOL  (input) DOUBLE PRECISION
   57: *          The minimum relative width of an interval.  When an interval
   58: *          is narrower than RELTOL times the larger (in
   59: *          magnitude) endpoint, then it is considered to be
   60: *          sufficiently small, i.e., converged.  Note: this should
   61: *          always be at least radix*machine epsilon.
   62: *
   63: *  W       (output) DOUBLE PRECISION
   64: *
   65: *  WERR    (output) DOUBLE PRECISION
   66: *          The error bound on the corresponding eigenvalue approximation
   67: *          in W.
   68: *
   69: *  INFO    (output) INTEGER
   70: *          = 0:       Eigenvalue converged
   71: *          = -1:      Eigenvalue did NOT converge
   72: *
   73: *  Internal Parameters
   74: *  ===================
   75: *
   76: *  FUDGE   DOUBLE PRECISION, default = 2
   77: *          A "fudge factor" to widen the Gershgorin intervals.
   78: *
   79: *  =====================================================================
   80: *
   81: *     .. Parameters ..
   82:       DOUBLE PRECISION   FUDGE, HALF, TWO, ZERO
   83:       PARAMETER          ( HALF = 0.5D0, TWO = 2.0D0,
   84:      $                     FUDGE = TWO, ZERO = 0.0D0 )
   85: *     ..
   86: *     .. Local Scalars ..
   87:       INTEGER   I, IT, ITMAX, NEGCNT
   88:       DOUBLE PRECISION   ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
   89:      $                   TMP2, TNORM
   90: *     ..
   91: *     .. External Functions ..
   92:       DOUBLE PRECISION   DLAMCH
   93:       EXTERNAL   DLAMCH
   94: *     ..
   95: *     .. Intrinsic Functions ..
   96:       INTRINSIC          ABS, INT, LOG, MAX
   97: *     ..
   98: *     .. Executable Statements ..
   99: *
  100: *     Get machine constants
  101:       EPS = DLAMCH( 'P' )
  102: 
  103:       TNORM = MAX( ABS( GL ), ABS( GU ) )
  104:       RTOLI = RELTOL
  105:       ATOLI = FUDGE*TWO*PIVMIN
  106: 
  107:       ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
  108:      $           LOG( TWO ) ) + 2
  109: 
  110:       INFO = -1
  111: 
  112:       LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
  113:       RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
  114:       IT = 0
  115: 
  116:  10   CONTINUE
  117: *
  118: *     Check if interval converged or maximum number of iterations reached
  119: *
  120:       TMP1 = ABS( RIGHT - LEFT )
  121:       TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
  122:       IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
  123:          INFO = 0
  124:          GOTO 30
  125:       ENDIF
  126:       IF(IT.GT.ITMAX)
  127:      $   GOTO 30
  128: 
  129: *
  130: *     Count number of negative pivots for mid-point
  131: *
  132:       IT = IT + 1
  133:       MID = HALF * (LEFT + RIGHT)
  134:       NEGCNT = 0
  135:       TMP1 = D( 1 ) - MID
  136:       IF( ABS( TMP1 ).LT.PIVMIN )
  137:      $   TMP1 = -PIVMIN
  138:       IF( TMP1.LE.ZERO )
  139:      $   NEGCNT = NEGCNT + 1
  140: *
  141:       DO 20 I = 2, N
  142:          TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
  143:          IF( ABS( TMP1 ).LT.PIVMIN )
  144:      $      TMP1 = -PIVMIN
  145:          IF( TMP1.LE.ZERO )
  146:      $      NEGCNT = NEGCNT + 1
  147:  20   CONTINUE
  148: 
  149:       IF(NEGCNT.GE.IW) THEN
  150:          RIGHT = MID
  151:       ELSE
  152:          LEFT = MID
  153:       ENDIF
  154:       GOTO 10
  155: 
  156:  30   CONTINUE
  157: *
  158: *     Converged or maximum number of iterations reached
  159: *
  160:       W = HALF * (LEFT + RIGHT)
  161:       WERR = HALF * ABS( RIGHT - LEFT )
  162: 
  163:       RETURN
  164: *
  165: *     End of DLARRK
  166: *
  167:       END