001:       FUNCTION SLANEG( N, D, LLD, SIGMA, PIVMIN, R )
002:       IMPLICIT NONE
003:       INTEGER SLANEG
004: *
005: *  -- LAPACK auxiliary routine (version 3.2) --
006: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *
010: *     .. Scalar Arguments ..
011:       INTEGER            N, R
012:       REAL               PIVMIN, SIGMA
013: *     ..
014: *     .. Array Arguments ..
015:       REAL               D( * ), LLD( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  SLANEG computes the Sturm count, the number of negative pivots
022: *  encountered while factoring tridiagonal T - sigma I = L D L^T.
023: *  This implementation works directly on the factors without forming
024: *  the tridiagonal matrix T.  The Sturm count is also the number of
025: *  eigenvalues of T less than sigma.
026: *
027: *  This routine is called from SLARRB.
028: *
029: *  The current routine does not use the PIVMIN parameter but rather
030: *  requires IEEE-754 propagation of Infinities and NaNs.  This
031: *  routine also has no input range restrictions but does require
032: *  default exception handling such that x/0 produces Inf when x is
033: *  non-zero, and Inf/Inf produces NaN.  For more information, see:
034: *
035: *    Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
036: *    Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
037: *    Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
038: *    (Tech report version in LAWN 172 with the same title.)
039: *
040: *  Arguments
041: *  =========
042: *
043: *  N       (input) INTEGER
044: *          The order of the matrix.
045: *
046: *  D       (input) REAL             array, dimension (N)
047: *          The N diagonal elements of the diagonal matrix D.
048: *
049: *  LLD     (input) REAL             array, dimension (N-1)
050: *          The (N-1) elements L(i)*L(i)*D(i).
051: *
052: *  SIGMA   (input) REAL            
053: *          Shift amount in T - sigma I = L D L^T.
054: *
055: *  PIVMIN  (input) REAL            
056: *          The minimum pivot in the Sturm sequence.  May be used
057: *          when zero pivots are encountered on non-IEEE-754
058: *          architectures.
059: *
060: *  R       (input) INTEGER
061: *          The twist index for the twisted factorization that is used
062: *          for the negcount.
063: *
064: *  Further Details
065: *  ===============
066: *
067: *  Based on contributions by
068: *     Osni Marques, LBNL/NERSC, USA
069: *     Christof Voemel, University of California, Berkeley, USA
070: *     Jason Riedy, University of California, Berkeley, USA
071: *
072: *  =====================================================================
073: *
074: *     .. Parameters ..
075:       REAL               ZERO, ONE
076:       PARAMETER        ( ZERO = 0.0E0, ONE = 1.0E0 )
077: *     Some architectures propagate Infinities and NaNs very slowly, so
078: *     the code computes counts in BLKLEN chunks.  Then a NaN can
079: *     propagate at most BLKLEN columns before being detected.  This is
080: *     not a general tuning parameter; it needs only to be just large
081: *     enough that the overhead is tiny in common cases.
082:       INTEGER BLKLEN
083:       PARAMETER ( BLKLEN = 128 )
084: *     ..
085: *     .. Local Scalars ..
086:       INTEGER            BJ, J, NEG1, NEG2, NEGCNT
087:       REAL               BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
088:       LOGICAL SAWNAN
089: *     ..
090: *     .. Intrinsic Functions ..
091:       INTRINSIC MIN, MAX
092: *     ..
093: *     .. External Functions ..
094:       LOGICAL SISNAN
095:       EXTERNAL SISNAN
096: *     ..
097: *     .. Executable Statements ..
098: 
099:       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 = SISNAN( 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 (SISNAN(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 = SISNAN( 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 (SISNAN(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:       SLANEG = NEGCNT
165:       END
166: