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