Annotation of rpl/lapack/lapack/dlaneg.f, revision 1.7
1.5 bertrand 1: INTEGER FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R )
1.1 bertrand 2: IMPLICIT NONE
3: *
1.5 bertrand 4: * -- LAPACK auxiliary routine (version 3.2.2) --
1.1 bertrand 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.5 bertrand 7: * June 2010
1.1 bertrand 8: *
9: * .. Scalar Arguments ..
10: INTEGER N, R
11: DOUBLE PRECISION PIVMIN, SIGMA
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION D( * ), LLD( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DLANEG computes the Sturm count, the number of negative pivots
21: * encountered while factoring tridiagonal T - sigma I = L D L^T.
22: * This implementation works directly on the factors without forming
23: * the tridiagonal matrix T. The Sturm count is also the number of
24: * eigenvalues of T less than sigma.
25: *
26: * This routine is called from DLARRB.
27: *
28: * The current routine does not use the PIVMIN parameter but rather
29: * requires IEEE-754 propagation of Infinities and NaNs. This
30: * routine also has no input range restrictions but does require
31: * default exception handling such that x/0 produces Inf when x is
32: * non-zero, and Inf/Inf produces NaN. For more information, see:
33: *
34: * Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
35: * Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
36: * Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624
37: * (Tech report version in LAWN 172 with the same title.)
38: *
39: * Arguments
40: * =========
41: *
42: * N (input) INTEGER
43: * The order of the matrix.
44: *
45: * D (input) DOUBLE PRECISION array, dimension (N)
46: * The N diagonal elements of the diagonal matrix D.
47: *
48: * LLD (input) DOUBLE PRECISION array, dimension (N-1)
49: * The (N-1) elements L(i)*L(i)*D(i).
50: *
51: * SIGMA (input) DOUBLE PRECISION
52: * Shift amount in T - sigma I = L D L^T.
53: *
54: * PIVMIN (input) DOUBLE PRECISION
55: * The minimum pivot in the Sturm sequence. May be used
56: * when zero pivots are encountered on non-IEEE-754
57: * architectures.
58: *
59: * R (input) INTEGER
60: * The twist index for the twisted factorization that is used
61: * for the negcount.
62: *
63: * Further Details
64: * ===============
65: *
66: * Based on contributions by
67: * Osni Marques, LBNL/NERSC, USA
68: * Christof Voemel, University of California, Berkeley, USA
69: * Jason Riedy, University of California, Berkeley, USA
70: *
71: * =====================================================================
72: *
73: * .. Parameters ..
74: DOUBLE PRECISION ZERO, ONE
75: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
76: * Some architectures propagate Infinities and NaNs very slowly, so
77: * the code computes counts in BLKLEN chunks. Then a NaN can
78: * propagate at most BLKLEN columns before being detected. This is
79: * not a general tuning parameter; it needs only to be just large
80: * enough that the overhead is tiny in common cases.
81: INTEGER BLKLEN
82: PARAMETER ( BLKLEN = 128 )
83: * ..
84: * .. Local Scalars ..
85: INTEGER BJ, J, NEG1, NEG2, NEGCNT
86: DOUBLE PRECISION BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
87: LOGICAL SAWNAN
88: * ..
89: * .. Intrinsic Functions ..
90: INTRINSIC MIN, MAX
91: * ..
92: * .. External Functions ..
93: LOGICAL DISNAN
94: EXTERNAL DISNAN
95: * ..
96: * .. Executable Statements ..
97:
98: NEGCNT = 0
99:
100: * I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
101: T = -SIGMA
102: DO 210 BJ = 1, R-1, BLKLEN
103: NEG1 = 0
104: BSAV = T
105: DO 21 J = BJ, MIN(BJ+BLKLEN-1, R-1)
106: DPLUS = D( J ) + T
107: IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
108: TMP = T / DPLUS
109: T = TMP * LLD( J ) - SIGMA
110: 21 CONTINUE
111: SAWNAN = DISNAN( T )
112: * Run a slower version of the above loop if a NaN is detected.
113: * A NaN should occur only with a zero pivot after an infinite
114: * pivot. In that case, substituting 1 for T/DPLUS is the
115: * correct limit.
116: IF( SAWNAN ) THEN
117: NEG1 = 0
118: T = BSAV
119: DO 22 J = BJ, MIN(BJ+BLKLEN-1, R-1)
120: DPLUS = D( J ) + T
121: IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
122: TMP = T / DPLUS
123: IF (DISNAN(TMP)) TMP = ONE
124: T = TMP * LLD(J) - SIGMA
125: 22 CONTINUE
126: END IF
127: NEGCNT = NEGCNT + NEG1
128: 210 CONTINUE
129: *
130: * II) lower part: L D L^T - SIGMA I = U- D- U-^T
131: P = D( N ) - SIGMA
132: DO 230 BJ = N-1, R, -BLKLEN
133: NEG2 = 0
134: BSAV = P
135: DO 23 J = BJ, MAX(BJ-BLKLEN+1, R), -1
136: DMINUS = LLD( J ) + P
137: IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
138: TMP = P / DMINUS
139: P = TMP * D( J ) - SIGMA
140: 23 CONTINUE
141: SAWNAN = DISNAN( P )
142: * As above, run a slower version that substitutes 1 for Inf/Inf.
143: *
144: IF( SAWNAN ) THEN
145: NEG2 = 0
146: P = BSAV
147: DO 24 J = BJ, MAX(BJ-BLKLEN+1, R), -1
148: DMINUS = LLD( J ) + P
149: IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
150: TMP = P / DMINUS
151: IF (DISNAN(TMP)) TMP = ONE
152: P = TMP * D(J) - SIGMA
153: 24 CONTINUE
154: END IF
155: NEGCNT = NEGCNT + NEG2
156: 230 CONTINUE
157: *
158: * III) Twist index
159: * T was shifted by SIGMA initially.
160: GAMMA = (T + SIGMA) + P
161: IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1
162:
163: DLANEG = NEGCNT
164: 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