001:       FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R )
002:       IMPLICIT NONE
003:       INTEGER DLANEG
004: *
005: *  -- LAPACK auxiliary routine (version 3.2) --
006: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            N, R
011:       DOUBLE PRECISION   PIVMIN, SIGMA
012: *     ..
013: *     .. Array Arguments ..
014:       DOUBLE PRECISION   D( * ), LLD( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DLANEG computes the Sturm count, the number of negative pivots
021: *  encountered while factoring tridiagonal T - sigma I = L D L^T.
022: *  This implementation works directly on the factors without forming
023: *  the tridiagonal matrix T.  The Sturm count is also the number of
024: *  eigenvalues of T less than sigma.
025: *
026: *  This routine is called from DLARRB.
027: *
028: *  The current routine does not use the PIVMIN parameter but rather
029: *  requires IEEE-754 propagation of Infinities and NaNs.  This
030: *  routine also has no input range restrictions but does require
031: *  default exception handling such that x/0 produces Inf when x is
032: *  non-zero, and Inf/Inf produces NaN.  For more information, see:
033: *
034: *    Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
035: *    Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
036: *    Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
037: *    (Tech report version in LAWN 172 with the same title.)
038: *
039: *  Arguments
040: *  =========
041: *
042: *  N       (input) INTEGER
043: *          The order of the matrix.
044: *
045: *  D       (input) DOUBLE PRECISION array, dimension (N)
046: *          The N diagonal elements of the diagonal matrix D.
047: *
048: *  LLD     (input) DOUBLE PRECISION array, dimension (N-1)
049: *          The (N-1) elements L(i)*L(i)*D(i).
050: *
051: *  SIGMA   (input) DOUBLE PRECISION
052: *          Shift amount in T - sigma I = L D L^T.
053: *
054: *  PIVMIN  (input) DOUBLE PRECISION
055: *          The minimum pivot in the Sturm sequence.  May be used
056: *          when zero pivots are encountered on non-IEEE-754
057: *          architectures.
058: *
059: *  R       (input) INTEGER
060: *          The twist index for the twisted factorization that is used
061: *          for the negcount.
062: *
063: *  Further Details
064: *  ===============
065: *
066: *  Based on contributions by
067: *     Osni Marques, LBNL/NERSC, USA
068: *     Christof Voemel, University of California, Berkeley, USA
069: *     Jason Riedy, University of California, Berkeley, USA
070: *
071: *  =====================================================================
072: *
073: *     .. Parameters ..
074:       DOUBLE PRECISION   ZERO, ONE
075:       PARAMETER        ( ZERO = 0.0D0, ONE = 1.0D0 )
076: *     Some architectures propagate Infinities and NaNs very slowly, so
077: *     the code computes counts in BLKLEN chunks.  Then a NaN can
078: *     propagate at most BLKLEN columns before being detected.  This is
079: *     not a general tuning parameter; it needs only to be just large
080: *     enough that the overhead is tiny in common cases.
081:       INTEGER BLKLEN
082:       PARAMETER ( BLKLEN = 128 )
083: *     ..
084: *     .. Local Scalars ..
085:       INTEGER            BJ, J, NEG1, NEG2, NEGCNT
086:       DOUBLE PRECISION   BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
087:       LOGICAL SAWNAN
088: *     ..
089: *     .. Intrinsic Functions ..
090:       INTRINSIC MIN, MAX
091: *     ..
092: *     .. External Functions ..
093:       LOGICAL DISNAN
094:       EXTERNAL DISNAN
095: *     ..
096: *     .. Executable Statements ..
097: 
098:       NEGCNT = 0
099: 
100: *     I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
101:       T = -SIGMA
102:       DO 210 BJ = 1, R-1, BLKLEN
103:          NEG1 = 0
104:          BSAV = T
105:          DO 21 J = BJ, MIN(BJ+BLKLEN-1, R-1)
106:             DPLUS = D( J ) + T
107:             IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
108:             TMP = T / DPLUS
109:             T = TMP * LLD( J ) - SIGMA
110:  21      CONTINUE
111:          SAWNAN = DISNAN( T )
112: *     Run a slower version of the above loop if a NaN is detected.
113: *     A NaN should occur only with a zero pivot after an infinite
114: *     pivot.  In that case, substituting 1 for T/DPLUS is the
115: *     correct limit.
116:          IF( SAWNAN ) THEN
117:             NEG1 = 0
118:             T = BSAV
119:             DO 22 J = BJ, MIN(BJ+BLKLEN-1, R-1)
120:                DPLUS = D( J ) + T
121:                IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
122:                TMP = T / DPLUS
123:                IF (DISNAN(TMP)) TMP = ONE
124:                T = TMP * LLD(J) - SIGMA
125:  22         CONTINUE
126:          END IF
127:          NEGCNT = NEGCNT + NEG1
128:  210  CONTINUE
129: *
130: *     II) lower part: L D L^T - SIGMA I = U- D- U-^T
131:       P = D( N ) - SIGMA
132:       DO 230 BJ = N-1, R, -BLKLEN
133:          NEG2 = 0
134:          BSAV = P
135:          DO 23 J = BJ, MAX(BJ-BLKLEN+1, R), -1
136:             DMINUS = LLD( J ) + P
137:             IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
138:             TMP = P / DMINUS
139:             P = TMP * D( J ) - SIGMA
140:  23      CONTINUE
141:          SAWNAN = DISNAN( P )
142: *     As above, run a slower version that substitutes 1 for Inf/Inf.
143: *
144:          IF( SAWNAN ) THEN
145:             NEG2 = 0
146:             P = BSAV
147:             DO 24 J = BJ, MAX(BJ-BLKLEN+1, R), -1
148:                DMINUS = LLD( J ) + P
149:                IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
150:                TMP = P / DMINUS
151:                IF (DISNAN(TMP)) TMP = ONE
152:                P = TMP * D(J) - SIGMA
153:  24         CONTINUE
154:          END IF
155:          NEGCNT = NEGCNT + NEG2
156:  230  CONTINUE
157: *
158: *     III) Twist index
159: *       T was shifted by SIGMA initially.
160:       GAMMA = (T + SIGMA) + P
161:       IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1
162: 
163:       DLANEG = NEGCNT
164:       END
165: