1 SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST,
2 $ RTOL, OFFSET, W, WERR, WORK, IWORK,
3 $ PIVMIN, SPDIAM, INFO )
4 *
5 * -- LAPACK auxiliary routine (version 3.2.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 * June 2010
9 *
10 * .. Scalar Arguments ..
11 INTEGER IFIRST, ILAST, INFO, N, OFFSET
12 DOUBLE PRECISION PIVMIN, RTOL, SPDIAM
13 * ..
14 * .. Array Arguments ..
15 INTEGER IWORK( * )
16 DOUBLE PRECISION D( * ), E2( * ), W( * ),
17 $ WERR( * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * Given the initial eigenvalue approximations of T, DLARRJ
24 * does bisection to refine the eigenvalues of T,
25 * W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
26 * guesses for these eigenvalues are input in W, the corresponding estimate
27 * of the error in these guesses in WERR. During bisection, intervals
28 * [left, right] are maintained by storing their mid-points and
29 * semi-widths in the arrays W and WERR respectively.
30 *
31 * Arguments
32 * =========
33 *
34 * N (input) INTEGER
35 * The order of the matrix.
36 *
37 * D (input) DOUBLE PRECISION array, dimension (N)
38 * The N diagonal elements of T.
39 *
40 * E2 (input) DOUBLE PRECISION array, dimension (N-1)
41 * The Squares of the (N-1) subdiagonal elements of T.
42 *
43 * IFIRST (input) INTEGER
44 * The index of the first eigenvalue to be computed.
45 *
46 * ILAST (input) INTEGER
47 * The index of the last eigenvalue to be computed.
48 *
49 * RTOL (input) DOUBLE PRECISION
50 * Tolerance for the convergence of the bisection intervals.
51 * An interval [LEFT,RIGHT] has converged if
52 * RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|).
53 *
54 * OFFSET (input) INTEGER
55 * Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET
56 * through ILAST-OFFSET elements of these arrays are to be used.
57 *
58 * W (input/output) DOUBLE PRECISION array, dimension (N)
59 * On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
60 * estimates of the eigenvalues of L D L^T indexed IFIRST through
61 * ILAST.
62 * On output, these estimates are refined.
63 *
64 * WERR (input/output) DOUBLE PRECISION array, dimension (N)
65 * On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
66 * the errors in the estimates of the corresponding elements in W.
67 * On output, these errors are refined.
68 *
69 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
70 * Workspace.
71 *
72 * IWORK (workspace) INTEGER array, dimension (2*N)
73 * Workspace.
74 *
75 * PIVMIN (input) DOUBLE PRECISION
76 * The minimum pivot in the Sturm sequence for T.
77 *
78 * SPDIAM (input) DOUBLE PRECISION
79 * The spectral diameter of T.
80 *
81 * INFO (output) INTEGER
82 * Error flag.
83 *
84 * Further Details
85 * ===============
86 *
87 * Based on contributions by
88 * Beresford Parlett, University of California, Berkeley, USA
89 * Jim Demmel, University of California, Berkeley, USA
90 * Inderjit Dhillon, University of Texas, Austin, USA
91 * Osni Marques, LBNL/NERSC, USA
92 * Christof Voemel, University of California, Berkeley, USA
93 *
94 * =====================================================================
95 *
96 * .. Parameters ..
97 DOUBLE PRECISION ZERO, ONE, TWO, HALF
98 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
99 $ HALF = 0.5D0 )
100 INTEGER MAXITR
101 * ..
102 * .. Local Scalars ..
103 INTEGER CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
104 $ OLNINT, P, PREV, SAVI1
105 DOUBLE PRECISION DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH
106 *
107 * ..
108 * .. Intrinsic Functions ..
109 INTRINSIC ABS, MAX
110 * ..
111 * .. Executable Statements ..
112 *
113 INFO = 0
114 *
115 MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
116 $ LOG( TWO ) ) + 2
117 *
118 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
119 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
120 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
121 * for an unconverged interval is set to the index of the next unconverged
122 * interval, and is -1 or 0 for a converged interval. Thus a linked
123 * list of unconverged intervals is set up.
124 *
125
126 I1 = IFIRST
127 I2 = ILAST
128 * The number of unconverged intervals
129 NINT = 0
130 * The last unconverged interval found
131 PREV = 0
132 DO 75 I = I1, I2
133 K = 2*I
134 II = I - OFFSET
135 LEFT = W( II ) - WERR( II )
136 MID = W(II)
137 RIGHT = W( II ) + WERR( II )
138 WIDTH = RIGHT - MID
139 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
140
141 * The following test prevents the test of converged intervals
142 IF( WIDTH.LT.RTOL*TMP ) THEN
143 * This interval has already converged and does not need refinement.
144 * (Note that the gaps might change through refining the
145 * eigenvalues, however, they can only get bigger.)
146 * Remove it from the list.
147 IWORK( K-1 ) = -1
148 * Make sure that I1 always points to the first unconverged interval
149 IF((I.EQ.I1).AND.(I.LT.I2)) I1 = I + 1
150 IF((PREV.GE.I1).AND.(I.LE.I2)) IWORK( 2*PREV-1 ) = I + 1
151 ELSE
152 * unconverged interval found
153 PREV = I
154 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
155 *
156 * Do while( CNT(LEFT).GT.I-1 )
157 *
158 FAC = ONE
159 20 CONTINUE
160 CNT = 0
161 S = LEFT
162 DPLUS = D( 1 ) - S
163 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
164 DO 30 J = 2, N
165 DPLUS = D( J ) - S - E2( J-1 )/DPLUS
166 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
167 30 CONTINUE
168 IF( CNT.GT.I-1 ) THEN
169 LEFT = LEFT - WERR( II )*FAC
170 FAC = TWO*FAC
171 GO TO 20
172 END IF
173 *
174 * Do while( CNT(RIGHT).LT.I )
175 *
176 FAC = ONE
177 50 CONTINUE
178 CNT = 0
179 S = RIGHT
180 DPLUS = D( 1 ) - S
181 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
182 DO 60 J = 2, N
183 DPLUS = D( J ) - S - E2( J-1 )/DPLUS
184 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
185 60 CONTINUE
186 IF( CNT.LT.I ) THEN
187 RIGHT = RIGHT + WERR( II )*FAC
188 FAC = TWO*FAC
189 GO TO 50
190 END IF
191 NINT = NINT + 1
192 IWORK( K-1 ) = I + 1
193 IWORK( K ) = CNT
194 END IF
195 WORK( K-1 ) = LEFT
196 WORK( K ) = RIGHT
197 75 CONTINUE
198
199
200 SAVI1 = I1
201 *
202 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
203 * and while (ITER.LT.MAXITR)
204 *
205 ITER = 0
206 80 CONTINUE
207 PREV = I1 - 1
208 I = I1
209 OLNINT = NINT
210
211 DO 100 P = 1, OLNINT
212 K = 2*I
213 II = I - OFFSET
214 NEXT = IWORK( K-1 )
215 LEFT = WORK( K-1 )
216 RIGHT = WORK( K )
217 MID = HALF*( LEFT + RIGHT )
218
219 * semiwidth of interval
220 WIDTH = RIGHT - MID
221 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
222
223 IF( ( WIDTH.LT.RTOL*TMP ) .OR.
224 $ (ITER.EQ.MAXITR) )THEN
225 * reduce number of unconverged intervals
226 NINT = NINT - 1
227 * Mark interval as converged.
228 IWORK( K-1 ) = 0
229 IF( I1.EQ.I ) THEN
230 I1 = NEXT
231 ELSE
232 * Prev holds the last unconverged interval previously examined
233 IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
234 END IF
235 I = NEXT
236 GO TO 100
237 END IF
238 PREV = I
239 *
240 * Perform one bisection step
241 *
242 CNT = 0
243 S = MID
244 DPLUS = D( 1 ) - S
245 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
246 DO 90 J = 2, N
247 DPLUS = D( J ) - S - E2( J-1 )/DPLUS
248 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
249 90 CONTINUE
250 IF( CNT.LE.I-1 ) THEN
251 WORK( K-1 ) = MID
252 ELSE
253 WORK( K ) = MID
254 END IF
255 I = NEXT
256
257 100 CONTINUE
258 ITER = ITER + 1
259 * do another loop if there are still unconverged intervals
260 * However, in the last iteration, all intervals are accepted
261 * since this is the best we can do.
262 IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
263 *
264 *
265 * At this point, all the intervals have converged
266 DO 110 I = SAVI1, ILAST
267 K = 2*I
268 II = I - OFFSET
269 * All intervals marked by '0' have been refined.
270 IF( IWORK( K-1 ).EQ.0 ) THEN
271 W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
272 WERR( II ) = WORK( K ) - W( II )
273 END IF
274 110 CONTINUE
275 *
276
277 RETURN
278 *
279 * End of DLARRJ
280 *
281 END
2 $ RTOL, OFFSET, W, WERR, WORK, IWORK,
3 $ PIVMIN, SPDIAM, INFO )
4 *
5 * -- LAPACK auxiliary routine (version 3.2.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 * June 2010
9 *
10 * .. Scalar Arguments ..
11 INTEGER IFIRST, ILAST, INFO, N, OFFSET
12 DOUBLE PRECISION PIVMIN, RTOL, SPDIAM
13 * ..
14 * .. Array Arguments ..
15 INTEGER IWORK( * )
16 DOUBLE PRECISION D( * ), E2( * ), W( * ),
17 $ WERR( * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * Given the initial eigenvalue approximations of T, DLARRJ
24 * does bisection to refine the eigenvalues of T,
25 * W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
26 * guesses for these eigenvalues are input in W, the corresponding estimate
27 * of the error in these guesses in WERR. During bisection, intervals
28 * [left, right] are maintained by storing their mid-points and
29 * semi-widths in the arrays W and WERR respectively.
30 *
31 * Arguments
32 * =========
33 *
34 * N (input) INTEGER
35 * The order of the matrix.
36 *
37 * D (input) DOUBLE PRECISION array, dimension (N)
38 * The N diagonal elements of T.
39 *
40 * E2 (input) DOUBLE PRECISION array, dimension (N-1)
41 * The Squares of the (N-1) subdiagonal elements of T.
42 *
43 * IFIRST (input) INTEGER
44 * The index of the first eigenvalue to be computed.
45 *
46 * ILAST (input) INTEGER
47 * The index of the last eigenvalue to be computed.
48 *
49 * RTOL (input) DOUBLE PRECISION
50 * Tolerance for the convergence of the bisection intervals.
51 * An interval [LEFT,RIGHT] has converged if
52 * RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|).
53 *
54 * OFFSET (input) INTEGER
55 * Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET
56 * through ILAST-OFFSET elements of these arrays are to be used.
57 *
58 * W (input/output) DOUBLE PRECISION array, dimension (N)
59 * On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
60 * estimates of the eigenvalues of L D L^T indexed IFIRST through
61 * ILAST.
62 * On output, these estimates are refined.
63 *
64 * WERR (input/output) DOUBLE PRECISION array, dimension (N)
65 * On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
66 * the errors in the estimates of the corresponding elements in W.
67 * On output, these errors are refined.
68 *
69 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
70 * Workspace.
71 *
72 * IWORK (workspace) INTEGER array, dimension (2*N)
73 * Workspace.
74 *
75 * PIVMIN (input) DOUBLE PRECISION
76 * The minimum pivot in the Sturm sequence for T.
77 *
78 * SPDIAM (input) DOUBLE PRECISION
79 * The spectral diameter of T.
80 *
81 * INFO (output) INTEGER
82 * Error flag.
83 *
84 * Further Details
85 * ===============
86 *
87 * Based on contributions by
88 * Beresford Parlett, University of California, Berkeley, USA
89 * Jim Demmel, University of California, Berkeley, USA
90 * Inderjit Dhillon, University of Texas, Austin, USA
91 * Osni Marques, LBNL/NERSC, USA
92 * Christof Voemel, University of California, Berkeley, USA
93 *
94 * =====================================================================
95 *
96 * .. Parameters ..
97 DOUBLE PRECISION ZERO, ONE, TWO, HALF
98 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
99 $ HALF = 0.5D0 )
100 INTEGER MAXITR
101 * ..
102 * .. Local Scalars ..
103 INTEGER CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
104 $ OLNINT, P, PREV, SAVI1
105 DOUBLE PRECISION DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH
106 *
107 * ..
108 * .. Intrinsic Functions ..
109 INTRINSIC ABS, MAX
110 * ..
111 * .. Executable Statements ..
112 *
113 INFO = 0
114 *
115 MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
116 $ LOG( TWO ) ) + 2
117 *
118 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
119 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
120 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
121 * for an unconverged interval is set to the index of the next unconverged
122 * interval, and is -1 or 0 for a converged interval. Thus a linked
123 * list of unconverged intervals is set up.
124 *
125
126 I1 = IFIRST
127 I2 = ILAST
128 * The number of unconverged intervals
129 NINT = 0
130 * The last unconverged interval found
131 PREV = 0
132 DO 75 I = I1, I2
133 K = 2*I
134 II = I - OFFSET
135 LEFT = W( II ) - WERR( II )
136 MID = W(II)
137 RIGHT = W( II ) + WERR( II )
138 WIDTH = RIGHT - MID
139 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
140
141 * The following test prevents the test of converged intervals
142 IF( WIDTH.LT.RTOL*TMP ) THEN
143 * This interval has already converged and does not need refinement.
144 * (Note that the gaps might change through refining the
145 * eigenvalues, however, they can only get bigger.)
146 * Remove it from the list.
147 IWORK( K-1 ) = -1
148 * Make sure that I1 always points to the first unconverged interval
149 IF((I.EQ.I1).AND.(I.LT.I2)) I1 = I + 1
150 IF((PREV.GE.I1).AND.(I.LE.I2)) IWORK( 2*PREV-1 ) = I + 1
151 ELSE
152 * unconverged interval found
153 PREV = I
154 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
155 *
156 * Do while( CNT(LEFT).GT.I-1 )
157 *
158 FAC = ONE
159 20 CONTINUE
160 CNT = 0
161 S = LEFT
162 DPLUS = D( 1 ) - S
163 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
164 DO 30 J = 2, N
165 DPLUS = D( J ) - S - E2( J-1 )/DPLUS
166 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
167 30 CONTINUE
168 IF( CNT.GT.I-1 ) THEN
169 LEFT = LEFT - WERR( II )*FAC
170 FAC = TWO*FAC
171 GO TO 20
172 END IF
173 *
174 * Do while( CNT(RIGHT).LT.I )
175 *
176 FAC = ONE
177 50 CONTINUE
178 CNT = 0
179 S = RIGHT
180 DPLUS = D( 1 ) - S
181 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
182 DO 60 J = 2, N
183 DPLUS = D( J ) - S - E2( J-1 )/DPLUS
184 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
185 60 CONTINUE
186 IF( CNT.LT.I ) THEN
187 RIGHT = RIGHT + WERR( II )*FAC
188 FAC = TWO*FAC
189 GO TO 50
190 END IF
191 NINT = NINT + 1
192 IWORK( K-1 ) = I + 1
193 IWORK( K ) = CNT
194 END IF
195 WORK( K-1 ) = LEFT
196 WORK( K ) = RIGHT
197 75 CONTINUE
198
199
200 SAVI1 = I1
201 *
202 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
203 * and while (ITER.LT.MAXITR)
204 *
205 ITER = 0
206 80 CONTINUE
207 PREV = I1 - 1
208 I = I1
209 OLNINT = NINT
210
211 DO 100 P = 1, OLNINT
212 K = 2*I
213 II = I - OFFSET
214 NEXT = IWORK( K-1 )
215 LEFT = WORK( K-1 )
216 RIGHT = WORK( K )
217 MID = HALF*( LEFT + RIGHT )
218
219 * semiwidth of interval
220 WIDTH = RIGHT - MID
221 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
222
223 IF( ( WIDTH.LT.RTOL*TMP ) .OR.
224 $ (ITER.EQ.MAXITR) )THEN
225 * reduce number of unconverged intervals
226 NINT = NINT - 1
227 * Mark interval as converged.
228 IWORK( K-1 ) = 0
229 IF( I1.EQ.I ) THEN
230 I1 = NEXT
231 ELSE
232 * Prev holds the last unconverged interval previously examined
233 IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
234 END IF
235 I = NEXT
236 GO TO 100
237 END IF
238 PREV = I
239 *
240 * Perform one bisection step
241 *
242 CNT = 0
243 S = MID
244 DPLUS = D( 1 ) - S
245 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
246 DO 90 J = 2, N
247 DPLUS = D( J ) - S - E2( J-1 )/DPLUS
248 IF( DPLUS.LT.ZERO ) CNT = CNT + 1
249 90 CONTINUE
250 IF( CNT.LE.I-1 ) THEN
251 WORK( K-1 ) = MID
252 ELSE
253 WORK( K ) = MID
254 END IF
255 I = NEXT
256
257 100 CONTINUE
258 ITER = ITER + 1
259 * do another loop if there are still unconverged intervals
260 * However, in the last iteration, all intervals are accepted
261 * since this is the best we can do.
262 IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
263 *
264 *
265 * At this point, all the intervals have converged
266 DO 110 I = SAVI1, ILAST
267 K = 2*I
268 II = I - OFFSET
269 * All intervals marked by '0' have been refined.
270 IF( IWORK( K-1 ).EQ.0 ) THEN
271 W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
272 WERR( II ) = WORK( K ) - W( II )
273 END IF
274 110 CONTINUE
275 *
276
277 RETURN
278 *
279 * End of DLARRJ
280 *
281 END