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