LAPACK  3.8.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.

Download SLANEG + dependencies [TGZ] [ZIP] [TXT]

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
 non-zero, and Inf/Inf produces NaN.  For more information, see:

   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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016
Contributors:
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA
Jason Riedy, University of California, Berkeley, USA

Definition at line 120 of file slaneg.f.

120 *
121 * -- LAPACK auxiliary routine (version 3.7.0) --
122 * -- LAPACK is a software package provided by Univ. of Tennessee, --
123 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
124 * December 2016
125 *
126 * .. Scalar Arguments ..
127  INTEGER n, r
128  REAL pivmin, sigma
129 * ..
130 * .. Array Arguments ..
131  REAL d( * ), lld( * )
132 * ..
133 *
134 * =====================================================================
135 *
136 * .. Parameters ..
137  REAL zero, one
138  parameter( zero = 0.0e0, one = 1.0e0 )
139 * Some architectures propagate Infinities and NaNs very slowly, so
140 * the code computes counts in BLKLEN chunks. Then a NaN can
141 * propagate at most BLKLEN columns before being detected. This is
142 * not a general tuning parameter; it needs only to be just large
143 * enough that the overhead is tiny in common cases.
144  INTEGER blklen
145  parameter( blklen = 128 )
146 * ..
147 * .. Local Scalars ..
148  INTEGER bj, j, neg1, neg2, negcnt
149  REAL bsav, dminus, dplus, gamma, p, t, tmp
150  LOGICAL sawnan
151 * ..
152 * .. Intrinsic Functions ..
153  INTRINSIC min, max
154 * ..
155 * .. External Functions ..
156  LOGICAL sisnan
157  EXTERNAL sisnan
158 * ..
159 * .. Executable Statements ..
160 
161  negcnt = 0
162 
163 * I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
164  t = -sigma
165  DO 210 bj = 1, r-1, blklen
166  neg1 = 0
167  bsav = t
168  DO 21 j = bj, min(bj+blklen-1, r-1)
169  dplus = d( j ) + t
170  IF( dplus.LT.zero ) neg1 = neg1 + 1
171  tmp = t / dplus
172  t = tmp * lld( j ) - sigma
173  21 CONTINUE
174  sawnan = sisnan( t )
175 * Run a slower version of the above loop if a NaN is detected.
176 * A NaN should occur only with a zero pivot after an infinite
177 * pivot. In that case, substituting 1 for T/DPLUS is the
178 * correct limit.
179  IF( sawnan ) THEN
180  neg1 = 0
181  t = bsav
182  DO 22 j = bj, min(bj+blklen-1, r-1)
183  dplus = d( j ) + t
184  IF( dplus.LT.zero ) neg1 = neg1 + 1
185  tmp = t / dplus
186  IF (sisnan(tmp)) tmp = one
187  t = tmp * lld(j) - sigma
188  22 CONTINUE
189  END IF
190  negcnt = negcnt + neg1
191  210 CONTINUE
192 *
193 * II) lower part: L D L^T - SIGMA I = U- D- U-^T
194  p = d( n ) - sigma
195  DO 230 bj = n-1, r, -blklen
196  neg2 = 0
197  bsav = p
198  DO 23 j = bj, max(bj-blklen+1, r), -1
199  dminus = lld( j ) + p
200  IF( dminus.LT.zero ) neg2 = neg2 + 1
201  tmp = p / dminus
202  p = tmp * d( j ) - sigma
203  23 CONTINUE
204  sawnan = sisnan( p )
205 * As above, run a slower version that substitutes 1 for Inf/Inf.
206 *
207  IF( sawnan ) THEN
208  neg2 = 0
209  p = bsav
210  DO 24 j = bj, max(bj-blklen+1, r), -1
211  dminus = lld( j ) + p
212  IF( dminus.LT.zero ) neg2 = neg2 + 1
213  tmp = p / dminus
214  IF (sisnan(tmp)) tmp = one
215  p = tmp * d(j) - sigma
216  24 CONTINUE
217  END IF
218  negcnt = negcnt + neg2
219  230 CONTINUE
220 *
221 * III) Twist index
222 * T was shifted by SIGMA initially.
223  gamma = (t + sigma) + p
224  IF( gamma.LT.zero ) negcnt = negcnt+1
225 
226  slaneg = negcnt
logical function sisnan(SIN)
SISNAN tests input for NaN.
Definition: sisnan.f:61
integer function slaneg(N, D, LLD, SIGMA, PIVMIN, R)
SLANEG computes the Sturm count.
Definition: slaneg.f:120