Return to dlaneg.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R )
    2:       IMPLICIT NONE
    3:       INTEGER DLANEG
    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            N, R
   12:       DOUBLE PRECISION   PIVMIN, SIGMA
   13: *     ..
   14: *     .. Array Arguments ..
   15:       DOUBLE PRECISION   D( * ), LLD( * )
   16: *     ..
   17: *
   18: *  Purpose
   19: *  =======
   20: *
   21: *  DLANEG computes the Sturm count, the number of negative pivots
   22: *  encountered while factoring tridiagonal T - sigma I = L D L^T.
   23: *  This implementation works directly on the factors without forming
   24: *  the tridiagonal matrix T.  The Sturm count is also the number of
   25: *  eigenvalues of T less than sigma.
   26: *
   27: *  This routine is called from DLARRB.
   28: *
   29: *  The current routine does not use the PIVMIN parameter but rather
   30: *  requires IEEE-754 propagation of Infinities and NaNs.  This
   31: *  routine also has no input range restrictions but does require
   32: *  default exception handling such that x/0 produces Inf when x is
   33: *  non-zero, and Inf/Inf produces NaN.  For more information, see:
   34: *
   35: *    Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
   36: *    Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
   37: *    Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
   38: *    (Tech report version in LAWN 172 with the same title.)
   39: *
   40: *  Arguments
   41: *  =========
   42: *
   43: *  N       (input) INTEGER
   44: *          The order of the matrix.
   45: *
   46: *  D       (input) DOUBLE PRECISION array, dimension (N)
   47: *          The N diagonal elements of the diagonal matrix D.
   48: *
   49: *  LLD     (input) DOUBLE PRECISION array, dimension (N-1)
   50: *          The (N-1) elements L(i)*L(i)*D(i).
   51: *
   52: *  SIGMA   (input) DOUBLE PRECISION
   53: *          Shift amount in T - sigma I = L D L^T.
   54: *
   55: *  PIVMIN  (input) DOUBLE PRECISION
   56: *          The minimum pivot in the Sturm sequence.  May be used
   57: *          when zero pivots are encountered on non-IEEE-754
   58: *          architectures.
   59: *
   60: *  R       (input) INTEGER
   61: *          The twist index for the twisted factorization that is used
   62: *          for the negcount.
   63: *
   64: *  Further Details
   65: *  ===============
   66: *
   67: *  Based on contributions by
   68: *     Osni Marques, LBNL/NERSC, USA
   69: *     Christof Voemel, University of California, Berkeley, USA
   70: *     Jason Riedy, University of California, Berkeley, USA
   71: *
   72: *  =====================================================================
   73: *
   74: *     .. Parameters ..
   75:       DOUBLE PRECISION   ZERO, ONE
   76:       PARAMETER        ( ZERO = 0.0D0, ONE = 1.0D0 )
   77: *     Some architectures propagate Infinities and NaNs very slowly, so
   78: *     the code computes counts in BLKLEN chunks.  Then a NaN can
   79: *     propagate at most BLKLEN columns before being detected.  This is
   80: *     not a general tuning parameter; it needs only to be just large
   81: *     enough that the overhead is tiny in common cases.
   82:       INTEGER BLKLEN
   83:       PARAMETER ( BLKLEN = 128 )
   84: *     ..
   85: *     .. Local Scalars ..
   86:       INTEGER            BJ, J, NEG1, NEG2, NEGCNT
   87:       DOUBLE PRECISION   BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
   88:       LOGICAL SAWNAN
   89: *     ..
   90: *     .. Intrinsic Functions ..
   91:       INTRINSIC MIN, MAX
   92: *     ..
   93: *     .. External Functions ..
   94:       LOGICAL DISNAN
   95:       EXTERNAL DISNAN
   96: *     ..
   97: *     .. Executable Statements ..
   98: 
   99:       NEGCNT = 0
  100: 
  101: *     I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
  102:       T = -SIGMA
  103:       DO 210 BJ = 1, R-1, BLKLEN
  104:          NEG1 = 0
  105:          BSAV = T
  106:          DO 21 J = BJ, MIN(BJ+BLKLEN-1, R-1)
  107:             DPLUS = D( J ) + T
  108:             IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
  109:             TMP = T / DPLUS
  110:             T = TMP * LLD( J ) - SIGMA
  111:  21      CONTINUE
  112:          SAWNAN = DISNAN( T )
  113: *     Run a slower version of the above loop if a NaN is detected.
  114: *     A NaN should occur only with a zero pivot after an infinite
  115: *     pivot.  In that case, substituting 1 for T/DPLUS is the
  116: *     correct limit.
  117:          IF( SAWNAN ) THEN
  118:             NEG1 = 0
  119:             T = BSAV
  120:             DO 22 J = BJ, MIN(BJ+BLKLEN-1, R-1)
  121:                DPLUS = D( J ) + T
  122:                IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
  123:                TMP = T / DPLUS
  124:                IF (DISNAN(TMP)) TMP = ONE
  125:                T = TMP * LLD(J) - SIGMA
  126:  22         CONTINUE
  127:          END IF
  128:          NEGCNT = NEGCNT + NEG1
  129:  210  CONTINUE
  130: *
  131: *     II) lower part: L D L^T - SIGMA I = U- D- U-^T
  132:       P = D( N ) - SIGMA
  133:       DO 230 BJ = N-1, R, -BLKLEN
  134:          NEG2 = 0
  135:          BSAV = P
  136:          DO 23 J = BJ, MAX(BJ-BLKLEN+1, R), -1
  137:             DMINUS = LLD( J ) + P
  138:             IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
  139:             TMP = P / DMINUS
  140:             P = TMP * D( J ) - SIGMA
  141:  23      CONTINUE
  142:          SAWNAN = DISNAN( P )
  143: *     As above, run a slower version that substitutes 1 for Inf/Inf.
  144: *
  145:          IF( SAWNAN ) THEN
  146:             NEG2 = 0
  147:             P = BSAV
  148:             DO 24 J = BJ, MAX(BJ-BLKLEN+1, R), -1
  149:                DMINUS = LLD( J ) + P
  150:                IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
  151:                TMP = P / DMINUS
  152:                IF (DISNAN(TMP)) TMP = ONE
  153:                P = TMP * D(J) - SIGMA
  154:  24         CONTINUE
  155:          END IF
  156:          NEGCNT = NEGCNT + NEG2
  157:  230  CONTINUE
  158: *
  159: *     III) Twist index
  160: *       T was shifted by SIGMA initially.
  161:       GAMMA = (T + SIGMA) + P
  162:       IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1
  163: 
  164:       DLANEG = NEGCNT
  165:       END