3 integer slaneg_(integer *n, real *d__, real *lld, real *sigma, real *pivmin,
6 /* System generated locals */
7 integer ret_val, i__1, i__2, i__3, i__4;
15 real bsav, gamma, dplus;
18 extern logical sisnan_(real *);
22 /* -- LAPACK auxiliary routine (version 3.1) -- */
23 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
26 /* .. Scalar Arguments .. */
28 /* .. Array Arguments .. */
34 /* SLANEG computes the Sturm count, the number of negative pivots */
35 /* encountered while factoring tridiagonal T - sigma I = L D L^T. */
36 /* This implementation works directly on the factors without forming */
37 /* the tridiagonal matrix T. The Sturm count is also the number of */
38 /* eigenvalues of T less than sigma. */
40 /* This routine is called from SLARRB. */
42 /* The current routine does not use the PIVMIN parameter but rather */
43 /* requires IEEE-754 propagation of Infinities and NaNs. This */
44 /* routine also has no input range restrictions but does require */
45 /* default exception handling such that x/0 produces Inf when x is */
46 /* non-zero, and Inf/Inf produces NaN. For more information, see: */
48 /* Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in */
49 /* Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on */
50 /* Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624 */
51 /* (Tech report version in LAWN 172 with the same title.) */
56 /* N (input) INTEGER */
57 /* The order of the matrix. */
59 /* D (input) REAL array, dimension (N) */
60 /* The N diagonal elements of the diagonal matrix D. */
62 /* LLD (input) REAL array, dimension (N-1) */
63 /* The (N-1) elements L(i)*L(i)*D(i). */
65 /* SIGMA (input) REAL */
66 /* Shift amount in T - sigma I = L D L^T. */
68 /* PIVMIN (input) REAL */
69 /* The minimum pivot in the Sturm sequence. May be used */
70 /* when zero pivots are encountered on non-IEEE-754 */
73 /* R (input) INTEGER */
74 /* The twist index for the twisted factorization that is used */
75 /* for the negcount. */
80 /* Based on contributions by */
81 /* Osni Marques, LBNL/NERSC, USA */
82 /* Christof Voemel, University of California, Berkeley, USA */
83 /* Jason Riedy, University of California, Berkeley, USA */
85 /* ===================================================================== */
87 /* .. Parameters .. */
88 /* Some architectures propagate Infinities and NaNs very slowly, so */
89 /* the code computes counts in BLKLEN chunks. Then a NaN can */
90 /* propagate at most BLKLEN columns before being detected. This is */
91 /* not a general tuning parameter; it needs only to be just large */
92 /* enough that the overhead is tiny in common cases. */
94 /* .. Local Scalars .. */
96 /* .. Intrinsic Functions .. */
98 /* .. External Functions .. */
100 /* .. Executable Statements .. */
101 /* Parameter adjustments */
107 /* I) upper part: L D L^T - SIGMA I = L+ D+ L+^T */
110 for (bj = 1; bj <= i__1; bj += 128) {
114 i__3 = bj + 127, i__4 = *r__ - 1;
115 i__2 = min(i__3,i__4);
116 for (j = bj; j <= i__2; ++j) {
122 t = tmp * lld[j] - *sigma;
125 sawnan = sisnan_(&t);
126 /* Run a slower version of the above loop if a NaN is detected. */
127 /* A NaN should occur only with a zero pivot after an infinite */
128 /* pivot. In that case, substituting 1 for T/DPLUS is the */
134 i__3 = bj + 127, i__4 = *r__ - 1;
135 i__2 = min(i__3,i__4);
136 for (j = bj; j <= i__2; ++j) {
145 t = tmp * lld[j] - *sigma;
153 /* II) lower part: L D L^T - SIGMA I = U- D- U-^T */
154 p = d__[*n] - *sigma;
156 for (bj = *n - 1; bj >= i__1; bj += -128) {
161 i__2 = max(i__3,*r__);
162 for (j = bj; j >= i__2; --j) {
168 p = tmp * d__[j] - *sigma;
171 sawnan = sisnan_(&p);
172 /* As above, run a slower version that substitutes 1 for Inf/Inf. */
179 i__2 = max(i__3,*r__);
180 for (j = bj; j >= i__2; --j) {
189 p = tmp * d__[j] - *sigma;
197 /* III) Twist index */
198 /* T was shifted by SIGMA initially. */
199 gamma = t + *sigma + p;