 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ slaneg()

 integer function slaneg ( integer N, real, dimension( * ) D, real, dimension( * ) LLD, real SIGMA, real PIVMIN, integer R )

SLANEG computes the Sturm count.

Purpose:
SLANEG computes the Sturm count, the number of negative pivots
encountered while factoring tridiagonal T - sigma I = L D L^T.
This implementation works directly on the factors without forming
the tridiagonal matrix T.  The Sturm count is also the number of
eigenvalues of T less than sigma.

This routine is called from SLARRB.

The current routine does not use the PIVMIN parameter but rather
requires IEEE-754 propagation of Infinities and NaNs.  This
routine also has no input range restrictions but does require
default exception handling such that x/0 produces Inf when x is

Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
(Tech report version in LAWN 172 with the same title.)
Parameters
 [in] N N is INTEGER The order of the matrix. [in] D D is REAL array, dimension (N) The N diagonal elements of the diagonal matrix D. [in] LLD LLD is REAL array, dimension (N-1) The (N-1) elements L(i)*L(i)*D(i). [in] SIGMA SIGMA is REAL Shift amount in T - sigma I = L D L^T. [in] PIVMIN PIVMIN is REAL The minimum pivot in the Sturm sequence. May be used when zero pivots are encountered on non-IEEE-754 architectures. [in] R R is INTEGER The twist index for the twisted factorization that is used for the negcount.
Contributors:
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA
Jason Riedy, University of California, Berkeley, USA

Definition at line 117 of file slaneg.f.

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
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
Here is the call graph for this function: