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: *> \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 DLANEG( 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: DOUBLE PRECISION PIVMIN, SIGMA
126: * ..
127: * .. Array Arguments ..
128: DOUBLE PRECISION D( * ), LLD( * )
129: * ..
130: *
131: * =====================================================================
132: *
133: * .. Parameters ..
134: DOUBLE PRECISION ZERO, ONE
135: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
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: DOUBLE PRECISION BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
147: LOGICAL SAWNAN
148: * ..
149: * .. Intrinsic Functions ..
150: INTRINSIC MIN, MAX
151: * ..
152: * .. External Functions ..
153: LOGICAL DISNAN
154: EXTERNAL DISNAN
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 = DISNAN( 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 (DISNAN(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 = DISNAN( 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 (DISNAN(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: DLANEG = NEGCNT
224: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaneg.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack