LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlaneg.f
Go to the documentation of this file.
1 *> \brief \b DLANEG 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
9 *> Download DLANEG + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaneg.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaneg.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaneg.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * INTEGER FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R )
22 *
23 * .. Scalar Arguments ..
24 * INTEGER N, R
25 * DOUBLE PRECISION PIVMIN, SIGMA
26 * ..
27 * .. Array Arguments ..
28 * DOUBLE PRECISION D( * ), LLD( * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> DLANEG 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 DLARRB.
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
49 *> non-zero, and Inf/Inf produces NaN. For more information, see:
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
81 *> Shift amount in T - sigma I = L D L^T.
82 *> \endverbatim
83 *>
84 *> \param[in] PIVMIN
85 *> \verbatim
86 *> PIVMIN is DOUBLE PRECISION
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 *> \date September 2012
108 *
109 *> \ingroup auxOTHERauxiliary
110 *
111 *> \par Contributors:
112 * ==================
113 *>
114 *> Osni Marques, LBNL/NERSC, USA \n
115 *> Christof Voemel, University of California, Berkeley, USA \n
116 *> Jason Riedy, University of California, Berkeley, USA \n
117 *>
118 * =====================================================================
119  INTEGER FUNCTION dlaneg( N, D, LLD, SIGMA, PIVMIN, R )
120 *
121 * -- LAPACK auxiliary routine (version 3.4.2) --
122 * -- LAPACK is a software package provided by Univ. of Tennessee, --
123 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
124 * September 2012
125 *
126 * .. Scalar Arguments ..
127  INTEGER n, r
128  DOUBLE PRECISION pivmin, sigma
129 * ..
130 * .. Array Arguments ..
131  DOUBLE PRECISION d( * ), lld( * )
132 * ..
133 *
134 * =====================================================================
135 *
136 * .. Parameters ..
137  DOUBLE PRECISION zero, one
138  parameter( zero = 0.0d0, one = 1.0d0 )
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  DOUBLE PRECISION bsav, dminus, dplus, gamma, p, t, tmp
150  LOGICAL sawnan
151 * ..
152 * .. Intrinsic Functions ..
153  INTRINSIC min, max
154 * ..
155 * .. External Functions ..
156  LOGICAL disnan
157  EXTERNAL disnan
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 = disnan( 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 (disnan(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 = disnan( 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 (disnan(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  dlaneg = negcnt
227  END