1 INTEGER FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R )
2 IMPLICIT NONE
3 *
4 * -- LAPACK auxiliary routine (version 3.2.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * June 2010
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
2 IMPLICIT NONE
3 *
4 * -- LAPACK auxiliary routine (version 3.2.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * June 2010
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