1*> \brief \b SLANEG computes the Sturm count.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SLANEG + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaneg.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaneg.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaneg.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       INTEGER FUNCTION SLANEG( N, D, LLD, SIGMA, PIVMIN, R )
22*
23*       .. Scalar Arguments ..
24*       INTEGER            N, R
25*       REAL               PIVMIN, SIGMA
26*       ..
27*       .. Array Arguments ..
28*       REAL               D( * ), LLD( * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> SLANEG computes the Sturm count, the number of negative pivots
38*> encountered while factoring tridiagonal T - sigma I = L D L^T.
39*> This implementation works directly on the factors without forming
40*> the tridiagonal matrix T.  The Sturm count is also the number of
41*> eigenvalues of T less than sigma.
42*>
43*> This routine is called from SLARRB.
44*>
45*> The current routine does not use the PIVMIN parameter but rather
46*> requires IEEE-754 propagation of Infinities and NaNs.  This
47*> routine also has no input range restrictions but does require
48*> default exception handling such that x/0 produces Inf when x is
49*> non-zero, and Inf/Inf produces NaN.  For more information, see:
50*>
51*>   Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
52*>   Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
53*>   Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
54*>   (Tech report version in LAWN 172 with the same title.)
55*> \endverbatim
56*
57*  Arguments:
58*  ==========
59*
60*> \param[in] N
61*> \verbatim
62*>          N is INTEGER
63*>          The order of the matrix.
64*> \endverbatim
65*>
66*> \param[in] D
67*> \verbatim
68*>          D is REAL array, dimension (N)
69*>          The N diagonal elements of the diagonal matrix D.
70*> \endverbatim
71*>
72*> \param[in] LLD
73*> \verbatim
74*>          LLD is REAL array, dimension (N-1)
75*>          The (N-1) elements L(i)*L(i)*D(i).
76*> \endverbatim
77*>
78*> \param[in] SIGMA
79*> \verbatim
80*>          SIGMA is REAL
81*>          Shift amount in T - sigma I = L D L^T.
82*> \endverbatim
83*>
84*> \param[in] PIVMIN
85*> \verbatim
86*>          PIVMIN is REAL
87*>          The minimum pivot in the Sturm sequence.  May be used
88*>          when zero pivots are encountered on non-IEEE-754
89*>          architectures.
90*> \endverbatim
91*>
92*> \param[in] R
93*> \verbatim
94*>          R is INTEGER
95*>          The twist index for the twisted factorization that is used
96*>          for the negcount.
97*> \endverbatim
98*
99*  Authors:
100*  ========
101*
102*> \author Univ. of Tennessee
103*> \author Univ. of California Berkeley
104*> \author Univ. of Colorado Denver
105*> \author NAG Ltd.
106*
107*> \ingroup OTHERauxiliary
108*
109*> \par Contributors:
110*  ==================
111*>
112*>     Osni Marques, LBNL/NERSC, USA \n
113*>     Christof Voemel, University of California, Berkeley, USA \n
114*>     Jason Riedy, University of California, Berkeley, USA \n
115*>
116*  =====================================================================
117      INTEGER FUNCTION SLANEG( N, D, LLD, SIGMA, PIVMIN, R )
118*
119*  -- LAPACK auxiliary routine --
120*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
121*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
122*
123*     .. Scalar Arguments ..
124      INTEGER            N, R
125      REAL               PIVMIN, SIGMA
126*     ..
127*     .. Array Arguments ..
128      REAL               D( * ), LLD( * )
129*     ..
130*
131*  =====================================================================
132*
133*     .. Parameters ..
134      REAL               ZERO, ONE
135      PARAMETER        ( ZERO = 0.0E0, ONE = 1.0E0 )
136*     Some architectures propagate Infinities and NaNs very slowly, so
137*     the code computes counts in BLKLEN chunks.  Then a NaN can
138*     propagate at most BLKLEN columns before being detected.  This is
139*     not a general tuning parameter; it needs only to be just large
140*     enough that the overhead is tiny in common cases.
141      INTEGER BLKLEN
142      PARAMETER ( BLKLEN = 128 )
143*     ..
144*     .. Local Scalars ..
145      INTEGER            BJ, J, NEG1, NEG2, NEGCNT
146      REAL               BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
147      LOGICAL SAWNAN
148*     ..
149*     .. Intrinsic Functions ..
150      INTRINSIC MIN, MAX
151*     ..
152*     .. External Functions ..
153      LOGICAL SISNAN
154      EXTERNAL SISNAN
155*     ..
156*     .. Executable Statements ..
157
158      NEGCNT = 0
159
160*     I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
161      T = -SIGMA
162      DO 210 BJ = 1, R-1, BLKLEN
163         NEG1 = 0
164         BSAV = T
165         DO 21 J = BJ, MIN(BJ+BLKLEN-1, R-1)
166            DPLUS = D( J ) + T
167            IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
168            TMP = T / DPLUS
169            T = TMP * LLD( J ) - SIGMA
170 21      CONTINUE
171         SAWNAN = SISNAN( T )
172*     Run a slower version of the above loop if a NaN is detected.
173*     A NaN should occur only with a zero pivot after an infinite
174*     pivot.  In that case, substituting 1 for T/DPLUS is the
175*     correct limit.
176         IF( SAWNAN ) THEN
177            NEG1 = 0
178            T = BSAV
179            DO 22 J = BJ, MIN(BJ+BLKLEN-1, R-1)
180               DPLUS = D( J ) + T
181               IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
182               TMP = T / DPLUS
183               IF (SISNAN(TMP)) TMP = ONE
184               T = TMP * LLD(J) - SIGMA
185 22         CONTINUE
186         END IF
187         NEGCNT = NEGCNT + NEG1
188 210  CONTINUE
189*
190*     II) lower part: L D L^T - SIGMA I = U- D- U-^T
191      P = D( N ) - SIGMA
192      DO 230 BJ = N-1, R, -BLKLEN
193         NEG2 = 0
194         BSAV = P
195         DO 23 J = BJ, MAX(BJ-BLKLEN+1, R), -1
196            DMINUS = LLD( J ) + P
197            IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
198            TMP = P / DMINUS
199            P = TMP * D( J ) - SIGMA
200 23      CONTINUE
201         SAWNAN = SISNAN( P )
202*     As above, run a slower version that substitutes 1 for Inf/Inf.
203*
204         IF( SAWNAN ) THEN
205            NEG2 = 0
206            P = BSAV
207            DO 24 J = BJ, MAX(BJ-BLKLEN+1, R), -1
208               DMINUS = LLD( J ) + P
209               IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
210               TMP = P / DMINUS
211               IF (SISNAN(TMP)) TMP = ONE
212               P = TMP * D(J) - SIGMA
213 24         CONTINUE
214         END IF
215         NEGCNT = NEGCNT + NEG2
216 230  CONTINUE
217*
218*     III) Twist index
219*       T was shifted by SIGMA initially.
220      GAMMA = (T + SIGMA) + P
221      IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1
222
223      SLANEG = NEGCNT
224      END
225