LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
slaneg.f
Go to the documentation of this file.
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
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
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
integer function slaneg(N, D, LLD, SIGMA, PIVMIN, R)
SLANEG computes the Sturm count.
Definition: slaneg.f:118
logical function sisnan(SIN)
SISNAN tests input for NaN.
Definition: sisnan.f:59