1 SUBROUTINE DLARRD( RANGE, ORDER, N, VL, VU, IL, IU, GERS,
2 $ RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
3 $ M, W, WERR, WL, WU, IBLOCK, INDEXW,
4 $ WORK, IWORK, INFO )
5 *
6 * -- LAPACK auxiliary routine (version 3.3.0) --
7 * -- LAPACK is a software package provided by Univ. of Tennessee, --
8 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9 * November 2010
10 *
11 * .. Scalar Arguments ..
12 CHARACTER ORDER, RANGE
13 INTEGER IL, INFO, IU, M, N, NSPLIT
14 DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU
15 * ..
16 * .. Array Arguments ..
17 INTEGER IBLOCK( * ), INDEXW( * ),
18 $ ISPLIT( * ), IWORK( * )
19 DOUBLE PRECISION D( * ), E( * ), E2( * ),
20 $ GERS( * ), W( * ), WERR( * ), WORK( * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * DLARRD computes the eigenvalues of a symmetric tridiagonal
27 * matrix T to suitable accuracy. This is an auxiliary code to be
28 * called from DSTEMR.
29 * The user may ask for all eigenvalues, all eigenvalues
30 * in the half-open interval (VL, VU], or the IL-th through IU-th
31 * eigenvalues.
32 *
33 * To avoid overflow, the matrix must be scaled so that its
34 * largest element is no greater than overflow**(1/2) *
35 * underflow**(1/4) in absolute value, and for greatest
36 * accuracy, it should not be much smaller than that.
37 *
38 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
39 * Matrix", Report CS41, Computer Science Dept., Stanford
40 * University, July 21, 1966.
41 *
42 * Arguments
43 * =========
44 *
45 * RANGE (input) CHARACTER*1
46 * = 'A': ("All") all eigenvalues will be found.
47 * = 'V': ("Value") all eigenvalues in the half-open interval
48 * (VL, VU] will be found.
49 * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
50 * entire matrix) will be found.
51 *
52 * ORDER (input) CHARACTER*1
53 * = 'B': ("By Block") the eigenvalues will be grouped by
54 * split-off block (see IBLOCK, ISPLIT) and
55 * ordered from smallest to largest within
56 * the block.
57 * = 'E': ("Entire matrix")
58 * the eigenvalues for the entire matrix
59 * will be ordered from smallest to
60 * largest.
61 *
62 * N (input) INTEGER
63 * The order of the tridiagonal matrix T. N >= 0.
64 *
65 * VL (input) DOUBLE PRECISION
66 * VU (input) DOUBLE PRECISION
67 * If RANGE='V', the lower and upper bounds of the interval to
68 * be searched for eigenvalues. Eigenvalues less than or equal
69 * to VL, or greater than VU, will not be returned. VL < VU.
70 * Not referenced if RANGE = 'A' or 'I'.
71 *
72 * IL (input) INTEGER
73 * IU (input) INTEGER
74 * If RANGE='I', the indices (in ascending order) of the
75 * smallest and largest eigenvalues to be returned.
76 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
77 * Not referenced if RANGE = 'A' or 'V'.
78 *
79 * GERS (input) DOUBLE PRECISION array, dimension (2*N)
80 * The N Gerschgorin intervals (the i-th Gerschgorin interval
81 * is (GERS(2*i-1), GERS(2*i)).
82 *
83 * RELTOL (input) DOUBLE PRECISION
84 * The minimum relative width of an interval. When an interval
85 * is narrower than RELTOL times the larger (in
86 * magnitude) endpoint, then it is considered to be
87 * sufficiently small, i.e., converged. Note: this should
88 * always be at least radix*machine epsilon.
89 *
90 * D (input) DOUBLE PRECISION array, dimension (N)
91 * The n diagonal elements of the tridiagonal matrix T.
92 *
93 * E (input) DOUBLE PRECISION array, dimension (N-1)
94 * The (n-1) off-diagonal elements of the tridiagonal matrix T.
95 *
96 * E2 (input) DOUBLE PRECISION array, dimension (N-1)
97 * The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
98 *
99 * PIVMIN (input) DOUBLE PRECISION
100 * The minimum pivot allowed in the Sturm sequence for T.
101 *
102 * NSPLIT (input) INTEGER
103 * The number of diagonal blocks in the matrix T.
104 * 1 <= NSPLIT <= N.
105 *
106 * ISPLIT (input) INTEGER array, dimension (N)
107 * The splitting points, at which T breaks up into submatrices.
108 * The first submatrix consists of rows/columns 1 to ISPLIT(1),
109 * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
110 * etc., and the NSPLIT-th consists of rows/columns
111 * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
112 * (Only the first NSPLIT elements will actually be used, but
113 * since the user cannot know a priori what value NSPLIT will
114 * have, N words must be reserved for ISPLIT.)
115 *
116 * M (output) INTEGER
117 * The actual number of eigenvalues found. 0 <= M <= N.
118 * (See also the description of INFO=2,3.)
119 *
120 * W (output) DOUBLE PRECISION array, dimension (N)
121 * On exit, the first M elements of W will contain the
122 * eigenvalue approximations. DLARRD computes an interval
123 * I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue
124 * approximation is given as the interval midpoint
125 * W(j)= ( a_j + b_j)/2. The corresponding error is bounded by
126 * WERR(j) = abs( a_j - b_j)/2
127 *
128 * WERR (output) DOUBLE PRECISION array, dimension (N)
129 * The error bound on the corresponding eigenvalue approximation
130 * in W.
131 *
132 * WL (output) DOUBLE PRECISION
133 * WU (output) DOUBLE PRECISION
134 * The interval (WL, WU] contains all the wanted eigenvalues.
135 * If RANGE='V', then WL=VL and WU=VU.
136 * If RANGE='A', then WL and WU are the global Gerschgorin bounds
137 * on the spectrum.
138 * If RANGE='I', then WL and WU are computed by DLAEBZ from the
139 * index range specified.
140 *
141 * IBLOCK (output) INTEGER array, dimension (N)
142 * At each row/column j where E(j) is zero or small, the
143 * matrix T is considered to split into a block diagonal
144 * matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which
145 * block (from 1 to the number of blocks) the eigenvalue W(i)
146 * belongs. (DLARRD may use the remaining N-M elements as
147 * workspace.)
148 *
149 * INDEXW (output) INTEGER array, dimension (N)
150 * The indices of the eigenvalues within each block (submatrix);
151 * for example, INDEXW(i)= j and IBLOCK(i)=k imply that the
152 * i-th eigenvalue W(i) is the j-th eigenvalue in block k.
153 *
154 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
155 *
156 * IWORK (workspace) INTEGER array, dimension (3*N)
157 *
158 * INFO (output) INTEGER
159 * = 0: successful exit
160 * < 0: if INFO = -i, the i-th argument had an illegal value
161 * > 0: some or all of the eigenvalues failed to converge or
162 * were not computed:
163 * =1 or 3: Bisection failed to converge for some
164 * eigenvalues; these eigenvalues are flagged by a
165 * negative block number. The effect is that the
166 * eigenvalues may not be as accurate as the
167 * absolute and relative tolerances. This is
168 * generally caused by unexpectedly inaccurate
169 * arithmetic.
170 * =2 or 3: RANGE='I' only: Not all of the eigenvalues
171 * IL:IU were found.
172 * Effect: M < IU+1-IL
173 * Cause: non-monotonic arithmetic, causing the
174 * Sturm sequence to be non-monotonic.
175 * Cure: recalculate, using RANGE='A', and pick
176 * out eigenvalues IL:IU. In some cases,
177 * increasing the PARAMETER "FUDGE" may
178 * make things work.
179 * = 4: RANGE='I', and the Gershgorin interval
180 * initially used was too small. No eigenvalues
181 * were computed.
182 * Probable cause: your machine has sloppy
183 * floating-point arithmetic.
184 * Cure: Increase the PARAMETER "FUDGE",
185 * recompile, and try again.
186 *
187 * Internal Parameters
188 * ===================
189 *
190 * FUDGE DOUBLE PRECISION, default = 2
191 * A "fudge factor" to widen the Gershgorin intervals. Ideally,
192 * a value of 1 should work, but on machines with sloppy
193 * arithmetic, this needs to be larger. The default for
194 * publicly released versions should be large enough to handle
195 * the worst machine around. Note that this has no effect
196 * on accuracy of the solution.
197 *
198 * Based on contributions by
199 * W. Kahan, University of California, Berkeley, USA
200 * Beresford Parlett, University of California, Berkeley, USA
201 * Jim Demmel, University of California, Berkeley, USA
202 * Inderjit Dhillon, University of Texas, Austin, USA
203 * Osni Marques, LBNL/NERSC, USA
204 * Christof Voemel, University of California, Berkeley, USA
205 *
206 * =====================================================================
207 *
208 * .. Parameters ..
209 DOUBLE PRECISION ZERO, ONE, TWO, HALF, FUDGE
210 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0,
211 $ TWO = 2.0D0, HALF = ONE/TWO,
212 $ FUDGE = TWO )
213 INTEGER ALLRNG, VALRNG, INDRNG
214 PARAMETER ( ALLRNG = 1, VALRNG = 2, INDRNG = 3 )
215 * ..
216 * .. Local Scalars ..
217 LOGICAL NCNVRG, TOOFEW
218 INTEGER I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
219 $ IM, IN, IOFF, IOUT, IRANGE, ITMAX, ITMP1,
220 $ ITMP2, IW, IWOFF, J, JBLK, JDISC, JE, JEE, NB,
221 $ NWL, NWU
222 DOUBLE PRECISION ATOLI, EPS, GL, GU, RTOLI, TMP1, TMP2,
223 $ TNORM, UFLOW, WKILL, WLU, WUL
224
225 * ..
226 * .. Local Arrays ..
227 INTEGER IDUMMA( 1 )
228 * ..
229 * .. External Functions ..
230 LOGICAL LSAME
231 INTEGER ILAENV
232 DOUBLE PRECISION DLAMCH
233 EXTERNAL LSAME, ILAENV, DLAMCH
234 * ..
235 * .. External Subroutines ..
236 EXTERNAL DLAEBZ
237 * ..
238 * .. Intrinsic Functions ..
239 INTRINSIC ABS, INT, LOG, MAX, MIN
240 * ..
241 * .. Executable Statements ..
242 *
243 INFO = 0
244 *
245 * Decode RANGE
246 *
247 IF( LSAME( RANGE, 'A' ) ) THEN
248 IRANGE = ALLRNG
249 ELSE IF( LSAME( RANGE, 'V' ) ) THEN
250 IRANGE = VALRNG
251 ELSE IF( LSAME( RANGE, 'I' ) ) THEN
252 IRANGE = INDRNG
253 ELSE
254 IRANGE = 0
255 END IF
256 *
257 * Check for Errors
258 *
259 IF( IRANGE.LE.0 ) THEN
260 INFO = -1
261 ELSE IF( .NOT.(LSAME(ORDER,'B').OR.LSAME(ORDER,'E')) ) THEN
262 INFO = -2
263 ELSE IF( N.LT.0 ) THEN
264 INFO = -3
265 ELSE IF( IRANGE.EQ.VALRNG ) THEN
266 IF( VL.GE.VU )
267 $ INFO = -5
268 ELSE IF( IRANGE.EQ.INDRNG .AND.
269 $ ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) THEN
270 INFO = -6
271 ELSE IF( IRANGE.EQ.INDRNG .AND.
272 $ ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN
273 INFO = -7
274 END IF
275 *
276 IF( INFO.NE.0 ) THEN
277 RETURN
278 END IF
279
280 * Initialize error flags
281 INFO = 0
282 NCNVRG = .FALSE.
283 TOOFEW = .FALSE.
284
285 * Quick return if possible
286 M = 0
287 IF( N.EQ.0 ) RETURN
288
289 * Simplification:
290 IF( IRANGE.EQ.INDRNG .AND. IL.EQ.1 .AND. IU.EQ.N ) IRANGE = 1
291
292 * Get machine constants
293 EPS = DLAMCH( 'P' )
294 UFLOW = DLAMCH( 'U' )
295
296
297 * Special Case when N=1
298 * Treat case of 1x1 matrix for quick return
299 IF( N.EQ.1 ) THEN
300 IF( (IRANGE.EQ.ALLRNG).OR.
301 $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
302 $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
303 M = 1
304 W(1) = D(1)
305 * The computation error of the eigenvalue is zero
306 WERR(1) = ZERO
307 IBLOCK( 1 ) = 1
308 INDEXW( 1 ) = 1
309 ENDIF
310 RETURN
311 END IF
312
313 * NB is the minimum vector length for vector bisection, or 0
314 * if only scalar is to be done.
315 NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 )
316 IF( NB.LE.1 ) NB = 0
317
318 * Find global spectral radius
319 GL = D(1)
320 GU = D(1)
321 DO 5 I = 1,N
322 GL = MIN( GL, GERS( 2*I - 1))
323 GU = MAX( GU, GERS(2*I) )
324 5 CONTINUE
325 * Compute global Gerschgorin bounds and spectral diameter
326 TNORM = MAX( ABS( GL ), ABS( GU ) )
327 GL = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
328 GU = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
329 * [JAN/28/2009] remove the line below since SPDIAM variable not use
330 * SPDIAM = GU - GL
331 * Input arguments for DLAEBZ:
332 * The relative tolerance. An interval (a,b] lies within
333 * "relative tolerance" if b-a < RELTOL*max(|a|,|b|),
334 RTOLI = RELTOL
335 * Set the absolute tolerance for interval convergence to zero to force
336 * interval convergence based on relative size of the interval.
337 * This is dangerous because intervals might not converge when RELTOL is
338 * small. But at least a very small number should be selected so that for
339 * strongly graded matrices, the code can get relatively accurate
340 * eigenvalues.
341 ATOLI = FUDGE*TWO*UFLOW + FUDGE*TWO*PIVMIN
342
343 IF( IRANGE.EQ.INDRNG ) THEN
344
345 * RANGE='I': Compute an interval containing eigenvalues
346 * IL through IU. The initial interval [GL,GU] from the global
347 * Gerschgorin bounds GL and GU is refined by DLAEBZ.
348 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
349 $ LOG( TWO ) ) + 2
350 WORK( N+1 ) = GL
351 WORK( N+2 ) = GL
352 WORK( N+3 ) = GU
353 WORK( N+4 ) = GU
354 WORK( N+5 ) = GL
355 WORK( N+6 ) = GU
356 IWORK( 1 ) = -1
357 IWORK( 2 ) = -1
358 IWORK( 3 ) = N + 1
359 IWORK( 4 ) = N + 1
360 IWORK( 5 ) = IL - 1
361 IWORK( 6 ) = IU
362 *
363 CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN,
364 $ D, E, E2, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
365 $ IWORK, W, IBLOCK, IINFO )
366 IF( IINFO .NE. 0 ) THEN
367 INFO = IINFO
368 RETURN
369 END IF
370 * On exit, output intervals may not be ordered by ascending negcount
371 IF( IWORK( 6 ).EQ.IU ) THEN
372 WL = WORK( N+1 )
373 WLU = WORK( N+3 )
374 NWL = IWORK( 1 )
375 WU = WORK( N+4 )
376 WUL = WORK( N+2 )
377 NWU = IWORK( 4 )
378 ELSE
379 WL = WORK( N+2 )
380 WLU = WORK( N+4 )
381 NWL = IWORK( 2 )
382 WU = WORK( N+3 )
383 WUL = WORK( N+1 )
384 NWU = IWORK( 3 )
385 END IF
386 * On exit, the interval [WL, WLU] contains a value with negcount NWL,
387 * and [WUL, WU] contains a value with negcount NWU.
388 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
389 INFO = 4
390 RETURN
391 END IF
392
393 ELSEIF( IRANGE.EQ.VALRNG ) THEN
394 WL = VL
395 WU = VU
396
397 ELSEIF( IRANGE.EQ.ALLRNG ) THEN
398 WL = GL
399 WU = GU
400 ENDIF
401
402
403
404 * Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
405 * NWL accumulates the number of eigenvalues .le. WL,
406 * NWU accumulates the number of eigenvalues .le. WU
407 M = 0
408 IEND = 0
409 INFO = 0
410 NWL = 0
411 NWU = 0
412 *
413 DO 70 JBLK = 1, NSPLIT
414 IOFF = IEND
415 IBEGIN = IOFF + 1
416 IEND = ISPLIT( JBLK )
417 IN = IEND - IOFF
418 *
419 IF( IN.EQ.1 ) THEN
420 * 1x1 block
421 IF( WL.GE.D( IBEGIN )-PIVMIN )
422 $ NWL = NWL + 1
423 IF( WU.GE.D( IBEGIN )-PIVMIN )
424 $ NWU = NWU + 1
425 IF( IRANGE.EQ.ALLRNG .OR.
426 $ ( WL.LT.D( IBEGIN )-PIVMIN
427 $ .AND. WU.GE. D( IBEGIN )-PIVMIN ) ) THEN
428 M = M + 1
429 W( M ) = D( IBEGIN )
430 WERR(M) = ZERO
431 * The gap for a single block doesn't matter for the later
432 * algorithm and is assigned an arbitrary large value
433 IBLOCK( M ) = JBLK
434 INDEXW( M ) = 1
435 END IF
436
437 * Disabled 2x2 case because of a failure on the following matrix
438 * RANGE = 'I', IL = IU = 4
439 * Original Tridiagonal, d = [
440 * -0.150102010615740E+00
441 * -0.849897989384260E+00
442 * -0.128208148052635E-15
443 * 0.128257718286320E-15
444 * ];
445 * e = [
446 * -0.357171383266986E+00
447 * -0.180411241501588E-15
448 * -0.175152352710251E-15
449 * ];
450 *
451 * ELSE IF( IN.EQ.2 ) THEN
452 ** 2x2 block
453 * DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
454 * TMP1 = HALF*(D(IBEGIN)+D(IEND))
455 * L1 = TMP1 - DISC
456 * IF( WL.GE. L1-PIVMIN )
457 * $ NWL = NWL + 1
458 * IF( WU.GE. L1-PIVMIN )
459 * $ NWU = NWU + 1
460 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
461 * $ L1-PIVMIN ) ) THEN
462 * M = M + 1
463 * W( M ) = L1
464 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
465 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
466 * IBLOCK( M ) = JBLK
467 * INDEXW( M ) = 1
468 * ENDIF
469 * L2 = TMP1 + DISC
470 * IF( WL.GE. L2-PIVMIN )
471 * $ NWL = NWL + 1
472 * IF( WU.GE. L2-PIVMIN )
473 * $ NWU = NWU + 1
474 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
475 * $ L2-PIVMIN ) ) THEN
476 * M = M + 1
477 * W( M ) = L2
478 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
479 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
480 * IBLOCK( M ) = JBLK
481 * INDEXW( M ) = 2
482 * ENDIF
483 ELSE
484 * General Case - block of size IN >= 2
485 * Compute local Gerschgorin interval and use it as the initial
486 * interval for DLAEBZ
487 GU = D( IBEGIN )
488 GL = D( IBEGIN )
489 TMP1 = ZERO
490
491 DO 40 J = IBEGIN, IEND
492 GL = MIN( GL, GERS( 2*J - 1))
493 GU = MAX( GU, GERS(2*J) )
494 40 CONTINUE
495 * [JAN/28/2009]
496 * change SPDIAM by TNORM in lines 2 and 3 thereafter
497 * line 1: remove computation of SPDIAM (not useful anymore)
498 * SPDIAM = GU - GL
499 * GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
500 * GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
501 GL = GL - FUDGE*TNORM*EPS*IN - FUDGE*PIVMIN
502 GU = GU + FUDGE*TNORM*EPS*IN + FUDGE*PIVMIN
503 *
504 IF( IRANGE.GT.1 ) THEN
505 IF( GU.LT.WL ) THEN
506 * the local block contains none of the wanted eigenvalues
507 NWL = NWL + IN
508 NWU = NWU + IN
509 GO TO 70
510 END IF
511 * refine search interval if possible, only range (WL,WU] matters
512 GL = MAX( GL, WL )
513 GU = MIN( GU, WU )
514 IF( GL.GE.GU )
515 $ GO TO 70
516 END IF
517
518 * Find negcount of initial interval boundaries GL and GU
519 WORK( N+1 ) = GL
520 WORK( N+IN+1 ) = GU
521 CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
522 $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
523 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
524 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
525 IF( IINFO .NE. 0 ) THEN
526 INFO = IINFO
527 RETURN
528 END IF
529 *
530 NWL = NWL + IWORK( 1 )
531 NWU = NWU + IWORK( IN+1 )
532 IWOFF = M - IWORK( 1 )
533
534 * Compute Eigenvalues
535 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
536 $ LOG( TWO ) ) + 2
537 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
538 $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
539 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
540 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
541 IF( IINFO .NE. 0 ) THEN
542 INFO = IINFO
543 RETURN
544 END IF
545 *
546 * Copy eigenvalues into W and IBLOCK
547 * Use -JBLK for block number for unconverged eigenvalues.
548 * Loop over the number of output intervals from DLAEBZ
549 DO 60 J = 1, IOUT
550 * eigenvalue approximation is middle point of interval
551 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
552 * semi length of error interval
553 TMP2 = HALF*ABS( WORK( J+N )-WORK( J+IN+N ) )
554 IF( J.GT.IOUT-IINFO ) THEN
555 * Flag non-convergence.
556 NCNVRG = .TRUE.
557 IB = -JBLK
558 ELSE
559 IB = JBLK
560 END IF
561 DO 50 JE = IWORK( J ) + 1 + IWOFF,
562 $ IWORK( J+IN ) + IWOFF
563 W( JE ) = TMP1
564 WERR( JE ) = TMP2
565 INDEXW( JE ) = JE - IWOFF
566 IBLOCK( JE ) = IB
567 50 CONTINUE
568 60 CONTINUE
569 *
570 M = M + IM
571 END IF
572 70 CONTINUE
573
574 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
575 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
576 IF( IRANGE.EQ.INDRNG ) THEN
577 IDISCL = IL - 1 - NWL
578 IDISCU = NWU - IU
579 *
580 IF( IDISCL.GT.0 ) THEN
581 IM = 0
582 DO 80 JE = 1, M
583 * Remove some of the smallest eigenvalues from the left so that
584 * at the end IDISCL =0. Move all eigenvalues up to the left.
585 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
586 IDISCL = IDISCL - 1
587 ELSE
588 IM = IM + 1
589 W( IM ) = W( JE )
590 WERR( IM ) = WERR( JE )
591 INDEXW( IM ) = INDEXW( JE )
592 IBLOCK( IM ) = IBLOCK( JE )
593 END IF
594 80 CONTINUE
595 M = IM
596 END IF
597 IF( IDISCU.GT.0 ) THEN
598 * Remove some of the largest eigenvalues from the right so that
599 * at the end IDISCU =0. Move all eigenvalues up to the left.
600 IM=M+1
601 DO 81 JE = M, 1, -1
602 IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
603 IDISCU = IDISCU - 1
604 ELSE
605 IM = IM - 1
606 W( IM ) = W( JE )
607 WERR( IM ) = WERR( JE )
608 INDEXW( IM ) = INDEXW( JE )
609 IBLOCK( IM ) = IBLOCK( JE )
610 END IF
611 81 CONTINUE
612 JEE = 0
613 DO 82 JE = IM, M
614 JEE = JEE + 1
615 W( JEE ) = W( JE )
616 WERR( JEE ) = WERR( JE )
617 INDEXW( JEE ) = INDEXW( JE )
618 IBLOCK( JEE ) = IBLOCK( JE )
619 82 CONTINUE
620 M = M-IM+1
621 END IF
622
623 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
624 * Code to deal with effects of bad arithmetic. (If N(w) is
625 * monotone non-decreasing, this should never happen.)
626 * Some low eigenvalues to be discarded are not in (WL,WLU],
627 * or high eigenvalues to be discarded are not in (WUL,WU]
628 * so just kill off the smallest IDISCL/largest IDISCU
629 * eigenvalues, by marking the corresponding IBLOCK = 0
630 IF( IDISCL.GT.0 ) THEN
631 WKILL = WU
632 DO 100 JDISC = 1, IDISCL
633 IW = 0
634 DO 90 JE = 1, M
635 IF( IBLOCK( JE ).NE.0 .AND.
636 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
637 IW = JE
638 WKILL = W( JE )
639 END IF
640 90 CONTINUE
641 IBLOCK( IW ) = 0
642 100 CONTINUE
643 END IF
644 IF( IDISCU.GT.0 ) THEN
645 WKILL = WL
646 DO 120 JDISC = 1, IDISCU
647 IW = 0
648 DO 110 JE = 1, M
649 IF( IBLOCK( JE ).NE.0 .AND.
650 $ ( W( JE ).GE.WKILL .OR. IW.EQ.0 ) ) THEN
651 IW = JE
652 WKILL = W( JE )
653 END IF
654 110 CONTINUE
655 IBLOCK( IW ) = 0
656 120 CONTINUE
657 END IF
658 * Now erase all eigenvalues with IBLOCK set to zero
659 IM = 0
660 DO 130 JE = 1, M
661 IF( IBLOCK( JE ).NE.0 ) THEN
662 IM = IM + 1
663 W( IM ) = W( JE )
664 WERR( IM ) = WERR( JE )
665 INDEXW( IM ) = INDEXW( JE )
666 IBLOCK( IM ) = IBLOCK( JE )
667 END IF
668 130 CONTINUE
669 M = IM
670 END IF
671 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
672 TOOFEW = .TRUE.
673 END IF
674 END IF
675 *
676 IF(( IRANGE.EQ.ALLRNG .AND. M.NE.N ).OR.
677 $ ( IRANGE.EQ.INDRNG .AND. M.NE.IU-IL+1 ) ) THEN
678 TOOFEW = .TRUE.
679 END IF
680
681 * If ORDER='B', do nothing the eigenvalues are already sorted by
682 * block.
683 * If ORDER='E', sort the eigenvalues from smallest to largest
684
685 IF( LSAME(ORDER,'E') .AND. NSPLIT.GT.1 ) THEN
686 DO 150 JE = 1, M - 1
687 IE = 0
688 TMP1 = W( JE )
689 DO 140 J = JE + 1, M
690 IF( W( J ).LT.TMP1 ) THEN
691 IE = J
692 TMP1 = W( J )
693 END IF
694 140 CONTINUE
695 IF( IE.NE.0 ) THEN
696 TMP2 = WERR( IE )
697 ITMP1 = IBLOCK( IE )
698 ITMP2 = INDEXW( IE )
699 W( IE ) = W( JE )
700 WERR( IE ) = WERR( JE )
701 IBLOCK( IE ) = IBLOCK( JE )
702 INDEXW( IE ) = INDEXW( JE )
703 W( JE ) = TMP1
704 WERR( JE ) = TMP2
705 IBLOCK( JE ) = ITMP1
706 INDEXW( JE ) = ITMP2
707 END IF
708 150 CONTINUE
709 END IF
710 *
711 INFO = 0
712 IF( NCNVRG )
713 $ INFO = INFO + 1
714 IF( TOOFEW )
715 $ INFO = INFO + 2
716 RETURN
717 *
718 * End of DLARRD
719 *
720 END
2 $ RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
3 $ M, W, WERR, WL, WU, IBLOCK, INDEXW,
4 $ WORK, IWORK, INFO )
5 *
6 * -- LAPACK auxiliary routine (version 3.3.0) --
7 * -- LAPACK is a software package provided by Univ. of Tennessee, --
8 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9 * November 2010
10 *
11 * .. Scalar Arguments ..
12 CHARACTER ORDER, RANGE
13 INTEGER IL, INFO, IU, M, N, NSPLIT
14 DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU
15 * ..
16 * .. Array Arguments ..
17 INTEGER IBLOCK( * ), INDEXW( * ),
18 $ ISPLIT( * ), IWORK( * )
19 DOUBLE PRECISION D( * ), E( * ), E2( * ),
20 $ GERS( * ), W( * ), WERR( * ), WORK( * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * DLARRD computes the eigenvalues of a symmetric tridiagonal
27 * matrix T to suitable accuracy. This is an auxiliary code to be
28 * called from DSTEMR.
29 * The user may ask for all eigenvalues, all eigenvalues
30 * in the half-open interval (VL, VU], or the IL-th through IU-th
31 * eigenvalues.
32 *
33 * To avoid overflow, the matrix must be scaled so that its
34 * largest element is no greater than overflow**(1/2) *
35 * underflow**(1/4) in absolute value, and for greatest
36 * accuracy, it should not be much smaller than that.
37 *
38 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
39 * Matrix", Report CS41, Computer Science Dept., Stanford
40 * University, July 21, 1966.
41 *
42 * Arguments
43 * =========
44 *
45 * RANGE (input) CHARACTER*1
46 * = 'A': ("All") all eigenvalues will be found.
47 * = 'V': ("Value") all eigenvalues in the half-open interval
48 * (VL, VU] will be found.
49 * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
50 * entire matrix) will be found.
51 *
52 * ORDER (input) CHARACTER*1
53 * = 'B': ("By Block") the eigenvalues will be grouped by
54 * split-off block (see IBLOCK, ISPLIT) and
55 * ordered from smallest to largest within
56 * the block.
57 * = 'E': ("Entire matrix")
58 * the eigenvalues for the entire matrix
59 * will be ordered from smallest to
60 * largest.
61 *
62 * N (input) INTEGER
63 * The order of the tridiagonal matrix T. N >= 0.
64 *
65 * VL (input) DOUBLE PRECISION
66 * VU (input) DOUBLE PRECISION
67 * If RANGE='V', the lower and upper bounds of the interval to
68 * be searched for eigenvalues. Eigenvalues less than or equal
69 * to VL, or greater than VU, will not be returned. VL < VU.
70 * Not referenced if RANGE = 'A' or 'I'.
71 *
72 * IL (input) INTEGER
73 * IU (input) INTEGER
74 * If RANGE='I', the indices (in ascending order) of the
75 * smallest and largest eigenvalues to be returned.
76 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
77 * Not referenced if RANGE = 'A' or 'V'.
78 *
79 * GERS (input) DOUBLE PRECISION array, dimension (2*N)
80 * The N Gerschgorin intervals (the i-th Gerschgorin interval
81 * is (GERS(2*i-1), GERS(2*i)).
82 *
83 * RELTOL (input) DOUBLE PRECISION
84 * The minimum relative width of an interval. When an interval
85 * is narrower than RELTOL times the larger (in
86 * magnitude) endpoint, then it is considered to be
87 * sufficiently small, i.e., converged. Note: this should
88 * always be at least radix*machine epsilon.
89 *
90 * D (input) DOUBLE PRECISION array, dimension (N)
91 * The n diagonal elements of the tridiagonal matrix T.
92 *
93 * E (input) DOUBLE PRECISION array, dimension (N-1)
94 * The (n-1) off-diagonal elements of the tridiagonal matrix T.
95 *
96 * E2 (input) DOUBLE PRECISION array, dimension (N-1)
97 * The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
98 *
99 * PIVMIN (input) DOUBLE PRECISION
100 * The minimum pivot allowed in the Sturm sequence for T.
101 *
102 * NSPLIT (input) INTEGER
103 * The number of diagonal blocks in the matrix T.
104 * 1 <= NSPLIT <= N.
105 *
106 * ISPLIT (input) INTEGER array, dimension (N)
107 * The splitting points, at which T breaks up into submatrices.
108 * The first submatrix consists of rows/columns 1 to ISPLIT(1),
109 * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
110 * etc., and the NSPLIT-th consists of rows/columns
111 * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
112 * (Only the first NSPLIT elements will actually be used, but
113 * since the user cannot know a priori what value NSPLIT will
114 * have, N words must be reserved for ISPLIT.)
115 *
116 * M (output) INTEGER
117 * The actual number of eigenvalues found. 0 <= M <= N.
118 * (See also the description of INFO=2,3.)
119 *
120 * W (output) DOUBLE PRECISION array, dimension (N)
121 * On exit, the first M elements of W will contain the
122 * eigenvalue approximations. DLARRD computes an interval
123 * I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue
124 * approximation is given as the interval midpoint
125 * W(j)= ( a_j + b_j)/2. The corresponding error is bounded by
126 * WERR(j) = abs( a_j - b_j)/2
127 *
128 * WERR (output) DOUBLE PRECISION array, dimension (N)
129 * The error bound on the corresponding eigenvalue approximation
130 * in W.
131 *
132 * WL (output) DOUBLE PRECISION
133 * WU (output) DOUBLE PRECISION
134 * The interval (WL, WU] contains all the wanted eigenvalues.
135 * If RANGE='V', then WL=VL and WU=VU.
136 * If RANGE='A', then WL and WU are the global Gerschgorin bounds
137 * on the spectrum.
138 * If RANGE='I', then WL and WU are computed by DLAEBZ from the
139 * index range specified.
140 *
141 * IBLOCK (output) INTEGER array, dimension (N)
142 * At each row/column j where E(j) is zero or small, the
143 * matrix T is considered to split into a block diagonal
144 * matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which
145 * block (from 1 to the number of blocks) the eigenvalue W(i)
146 * belongs. (DLARRD may use the remaining N-M elements as
147 * workspace.)
148 *
149 * INDEXW (output) INTEGER array, dimension (N)
150 * The indices of the eigenvalues within each block (submatrix);
151 * for example, INDEXW(i)= j and IBLOCK(i)=k imply that the
152 * i-th eigenvalue W(i) is the j-th eigenvalue in block k.
153 *
154 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
155 *
156 * IWORK (workspace) INTEGER array, dimension (3*N)
157 *
158 * INFO (output) INTEGER
159 * = 0: successful exit
160 * < 0: if INFO = -i, the i-th argument had an illegal value
161 * > 0: some or all of the eigenvalues failed to converge or
162 * were not computed:
163 * =1 or 3: Bisection failed to converge for some
164 * eigenvalues; these eigenvalues are flagged by a
165 * negative block number. The effect is that the
166 * eigenvalues may not be as accurate as the
167 * absolute and relative tolerances. This is
168 * generally caused by unexpectedly inaccurate
169 * arithmetic.
170 * =2 or 3: RANGE='I' only: Not all of the eigenvalues
171 * IL:IU were found.
172 * Effect: M < IU+1-IL
173 * Cause: non-monotonic arithmetic, causing the
174 * Sturm sequence to be non-monotonic.
175 * Cure: recalculate, using RANGE='A', and pick
176 * out eigenvalues IL:IU. In some cases,
177 * increasing the PARAMETER "FUDGE" may
178 * make things work.
179 * = 4: RANGE='I', and the Gershgorin interval
180 * initially used was too small. No eigenvalues
181 * were computed.
182 * Probable cause: your machine has sloppy
183 * floating-point arithmetic.
184 * Cure: Increase the PARAMETER "FUDGE",
185 * recompile, and try again.
186 *
187 * Internal Parameters
188 * ===================
189 *
190 * FUDGE DOUBLE PRECISION, default = 2
191 * A "fudge factor" to widen the Gershgorin intervals. Ideally,
192 * a value of 1 should work, but on machines with sloppy
193 * arithmetic, this needs to be larger. The default for
194 * publicly released versions should be large enough to handle
195 * the worst machine around. Note that this has no effect
196 * on accuracy of the solution.
197 *
198 * Based on contributions by
199 * W. Kahan, University of California, Berkeley, USA
200 * Beresford Parlett, University of California, Berkeley, USA
201 * Jim Demmel, University of California, Berkeley, USA
202 * Inderjit Dhillon, University of Texas, Austin, USA
203 * Osni Marques, LBNL/NERSC, USA
204 * Christof Voemel, University of California, Berkeley, USA
205 *
206 * =====================================================================
207 *
208 * .. Parameters ..
209 DOUBLE PRECISION ZERO, ONE, TWO, HALF, FUDGE
210 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0,
211 $ TWO = 2.0D0, HALF = ONE/TWO,
212 $ FUDGE = TWO )
213 INTEGER ALLRNG, VALRNG, INDRNG
214 PARAMETER ( ALLRNG = 1, VALRNG = 2, INDRNG = 3 )
215 * ..
216 * .. Local Scalars ..
217 LOGICAL NCNVRG, TOOFEW
218 INTEGER I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
219 $ IM, IN, IOFF, IOUT, IRANGE, ITMAX, ITMP1,
220 $ ITMP2, IW, IWOFF, J, JBLK, JDISC, JE, JEE, NB,
221 $ NWL, NWU
222 DOUBLE PRECISION ATOLI, EPS, GL, GU, RTOLI, TMP1, TMP2,
223 $ TNORM, UFLOW, WKILL, WLU, WUL
224
225 * ..
226 * .. Local Arrays ..
227 INTEGER IDUMMA( 1 )
228 * ..
229 * .. External Functions ..
230 LOGICAL LSAME
231 INTEGER ILAENV
232 DOUBLE PRECISION DLAMCH
233 EXTERNAL LSAME, ILAENV, DLAMCH
234 * ..
235 * .. External Subroutines ..
236 EXTERNAL DLAEBZ
237 * ..
238 * .. Intrinsic Functions ..
239 INTRINSIC ABS, INT, LOG, MAX, MIN
240 * ..
241 * .. Executable Statements ..
242 *
243 INFO = 0
244 *
245 * Decode RANGE
246 *
247 IF( LSAME( RANGE, 'A' ) ) THEN
248 IRANGE = ALLRNG
249 ELSE IF( LSAME( RANGE, 'V' ) ) THEN
250 IRANGE = VALRNG
251 ELSE IF( LSAME( RANGE, 'I' ) ) THEN
252 IRANGE = INDRNG
253 ELSE
254 IRANGE = 0
255 END IF
256 *
257 * Check for Errors
258 *
259 IF( IRANGE.LE.0 ) THEN
260 INFO = -1
261 ELSE IF( .NOT.(LSAME(ORDER,'B').OR.LSAME(ORDER,'E')) ) THEN
262 INFO = -2
263 ELSE IF( N.LT.0 ) THEN
264 INFO = -3
265 ELSE IF( IRANGE.EQ.VALRNG ) THEN
266 IF( VL.GE.VU )
267 $ INFO = -5
268 ELSE IF( IRANGE.EQ.INDRNG .AND.
269 $ ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) THEN
270 INFO = -6
271 ELSE IF( IRANGE.EQ.INDRNG .AND.
272 $ ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN
273 INFO = -7
274 END IF
275 *
276 IF( INFO.NE.0 ) THEN
277 RETURN
278 END IF
279
280 * Initialize error flags
281 INFO = 0
282 NCNVRG = .FALSE.
283 TOOFEW = .FALSE.
284
285 * Quick return if possible
286 M = 0
287 IF( N.EQ.0 ) RETURN
288
289 * Simplification:
290 IF( IRANGE.EQ.INDRNG .AND. IL.EQ.1 .AND. IU.EQ.N ) IRANGE = 1
291
292 * Get machine constants
293 EPS = DLAMCH( 'P' )
294 UFLOW = DLAMCH( 'U' )
295
296
297 * Special Case when N=1
298 * Treat case of 1x1 matrix for quick return
299 IF( N.EQ.1 ) THEN
300 IF( (IRANGE.EQ.ALLRNG).OR.
301 $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
302 $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
303 M = 1
304 W(1) = D(1)
305 * The computation error of the eigenvalue is zero
306 WERR(1) = ZERO
307 IBLOCK( 1 ) = 1
308 INDEXW( 1 ) = 1
309 ENDIF
310 RETURN
311 END IF
312
313 * NB is the minimum vector length for vector bisection, or 0
314 * if only scalar is to be done.
315 NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 )
316 IF( NB.LE.1 ) NB = 0
317
318 * Find global spectral radius
319 GL = D(1)
320 GU = D(1)
321 DO 5 I = 1,N
322 GL = MIN( GL, GERS( 2*I - 1))
323 GU = MAX( GU, GERS(2*I) )
324 5 CONTINUE
325 * Compute global Gerschgorin bounds and spectral diameter
326 TNORM = MAX( ABS( GL ), ABS( GU ) )
327 GL = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
328 GU = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
329 * [JAN/28/2009] remove the line below since SPDIAM variable not use
330 * SPDIAM = GU - GL
331 * Input arguments for DLAEBZ:
332 * The relative tolerance. An interval (a,b] lies within
333 * "relative tolerance" if b-a < RELTOL*max(|a|,|b|),
334 RTOLI = RELTOL
335 * Set the absolute tolerance for interval convergence to zero to force
336 * interval convergence based on relative size of the interval.
337 * This is dangerous because intervals might not converge when RELTOL is
338 * small. But at least a very small number should be selected so that for
339 * strongly graded matrices, the code can get relatively accurate
340 * eigenvalues.
341 ATOLI = FUDGE*TWO*UFLOW + FUDGE*TWO*PIVMIN
342
343 IF( IRANGE.EQ.INDRNG ) THEN
344
345 * RANGE='I': Compute an interval containing eigenvalues
346 * IL through IU. The initial interval [GL,GU] from the global
347 * Gerschgorin bounds GL and GU is refined by DLAEBZ.
348 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
349 $ LOG( TWO ) ) + 2
350 WORK( N+1 ) = GL
351 WORK( N+2 ) = GL
352 WORK( N+3 ) = GU
353 WORK( N+4 ) = GU
354 WORK( N+5 ) = GL
355 WORK( N+6 ) = GU
356 IWORK( 1 ) = -1
357 IWORK( 2 ) = -1
358 IWORK( 3 ) = N + 1
359 IWORK( 4 ) = N + 1
360 IWORK( 5 ) = IL - 1
361 IWORK( 6 ) = IU
362 *
363 CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN,
364 $ D, E, E2, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
365 $ IWORK, W, IBLOCK, IINFO )
366 IF( IINFO .NE. 0 ) THEN
367 INFO = IINFO
368 RETURN
369 END IF
370 * On exit, output intervals may not be ordered by ascending negcount
371 IF( IWORK( 6 ).EQ.IU ) THEN
372 WL = WORK( N+1 )
373 WLU = WORK( N+3 )
374 NWL = IWORK( 1 )
375 WU = WORK( N+4 )
376 WUL = WORK( N+2 )
377 NWU = IWORK( 4 )
378 ELSE
379 WL = WORK( N+2 )
380 WLU = WORK( N+4 )
381 NWL = IWORK( 2 )
382 WU = WORK( N+3 )
383 WUL = WORK( N+1 )
384 NWU = IWORK( 3 )
385 END IF
386 * On exit, the interval [WL, WLU] contains a value with negcount NWL,
387 * and [WUL, WU] contains a value with negcount NWU.
388 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
389 INFO = 4
390 RETURN
391 END IF
392
393 ELSEIF( IRANGE.EQ.VALRNG ) THEN
394 WL = VL
395 WU = VU
396
397 ELSEIF( IRANGE.EQ.ALLRNG ) THEN
398 WL = GL
399 WU = GU
400 ENDIF
401
402
403
404 * Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
405 * NWL accumulates the number of eigenvalues .le. WL,
406 * NWU accumulates the number of eigenvalues .le. WU
407 M = 0
408 IEND = 0
409 INFO = 0
410 NWL = 0
411 NWU = 0
412 *
413 DO 70 JBLK = 1, NSPLIT
414 IOFF = IEND
415 IBEGIN = IOFF + 1
416 IEND = ISPLIT( JBLK )
417 IN = IEND - IOFF
418 *
419 IF( IN.EQ.1 ) THEN
420 * 1x1 block
421 IF( WL.GE.D( IBEGIN )-PIVMIN )
422 $ NWL = NWL + 1
423 IF( WU.GE.D( IBEGIN )-PIVMIN )
424 $ NWU = NWU + 1
425 IF( IRANGE.EQ.ALLRNG .OR.
426 $ ( WL.LT.D( IBEGIN )-PIVMIN
427 $ .AND. WU.GE. D( IBEGIN )-PIVMIN ) ) THEN
428 M = M + 1
429 W( M ) = D( IBEGIN )
430 WERR(M) = ZERO
431 * The gap for a single block doesn't matter for the later
432 * algorithm and is assigned an arbitrary large value
433 IBLOCK( M ) = JBLK
434 INDEXW( M ) = 1
435 END IF
436
437 * Disabled 2x2 case because of a failure on the following matrix
438 * RANGE = 'I', IL = IU = 4
439 * Original Tridiagonal, d = [
440 * -0.150102010615740E+00
441 * -0.849897989384260E+00
442 * -0.128208148052635E-15
443 * 0.128257718286320E-15
444 * ];
445 * e = [
446 * -0.357171383266986E+00
447 * -0.180411241501588E-15
448 * -0.175152352710251E-15
449 * ];
450 *
451 * ELSE IF( IN.EQ.2 ) THEN
452 ** 2x2 block
453 * DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
454 * TMP1 = HALF*(D(IBEGIN)+D(IEND))
455 * L1 = TMP1 - DISC
456 * IF( WL.GE. L1-PIVMIN )
457 * $ NWL = NWL + 1
458 * IF( WU.GE. L1-PIVMIN )
459 * $ NWU = NWU + 1
460 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
461 * $ L1-PIVMIN ) ) THEN
462 * M = M + 1
463 * W( M ) = L1
464 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
465 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
466 * IBLOCK( M ) = JBLK
467 * INDEXW( M ) = 1
468 * ENDIF
469 * L2 = TMP1 + DISC
470 * IF( WL.GE. L2-PIVMIN )
471 * $ NWL = NWL + 1
472 * IF( WU.GE. L2-PIVMIN )
473 * $ NWU = NWU + 1
474 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
475 * $ L2-PIVMIN ) ) THEN
476 * M = M + 1
477 * W( M ) = L2
478 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
479 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
480 * IBLOCK( M ) = JBLK
481 * INDEXW( M ) = 2
482 * ENDIF
483 ELSE
484 * General Case - block of size IN >= 2
485 * Compute local Gerschgorin interval and use it as the initial
486 * interval for DLAEBZ
487 GU = D( IBEGIN )
488 GL = D( IBEGIN )
489 TMP1 = ZERO
490
491 DO 40 J = IBEGIN, IEND
492 GL = MIN( GL, GERS( 2*J - 1))
493 GU = MAX( GU, GERS(2*J) )
494 40 CONTINUE
495 * [JAN/28/2009]
496 * change SPDIAM by TNORM in lines 2 and 3 thereafter
497 * line 1: remove computation of SPDIAM (not useful anymore)
498 * SPDIAM = GU - GL
499 * GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
500 * GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
501 GL = GL - FUDGE*TNORM*EPS*IN - FUDGE*PIVMIN
502 GU = GU + FUDGE*TNORM*EPS*IN + FUDGE*PIVMIN
503 *
504 IF( IRANGE.GT.1 ) THEN
505 IF( GU.LT.WL ) THEN
506 * the local block contains none of the wanted eigenvalues
507 NWL = NWL + IN
508 NWU = NWU + IN
509 GO TO 70
510 END IF
511 * refine search interval if possible, only range (WL,WU] matters
512 GL = MAX( GL, WL )
513 GU = MIN( GU, WU )
514 IF( GL.GE.GU )
515 $ GO TO 70
516 END IF
517
518 * Find negcount of initial interval boundaries GL and GU
519 WORK( N+1 ) = GL
520 WORK( N+IN+1 ) = GU
521 CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
522 $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
523 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
524 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
525 IF( IINFO .NE. 0 ) THEN
526 INFO = IINFO
527 RETURN
528 END IF
529 *
530 NWL = NWL + IWORK( 1 )
531 NWU = NWU + IWORK( IN+1 )
532 IWOFF = M - IWORK( 1 )
533
534 * Compute Eigenvalues
535 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
536 $ LOG( TWO ) ) + 2
537 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
538 $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
539 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
540 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
541 IF( IINFO .NE. 0 ) THEN
542 INFO = IINFO
543 RETURN
544 END IF
545 *
546 * Copy eigenvalues into W and IBLOCK
547 * Use -JBLK for block number for unconverged eigenvalues.
548 * Loop over the number of output intervals from DLAEBZ
549 DO 60 J = 1, IOUT
550 * eigenvalue approximation is middle point of interval
551 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
552 * semi length of error interval
553 TMP2 = HALF*ABS( WORK( J+N )-WORK( J+IN+N ) )
554 IF( J.GT.IOUT-IINFO ) THEN
555 * Flag non-convergence.
556 NCNVRG = .TRUE.
557 IB = -JBLK
558 ELSE
559 IB = JBLK
560 END IF
561 DO 50 JE = IWORK( J ) + 1 + IWOFF,
562 $ IWORK( J+IN ) + IWOFF
563 W( JE ) = TMP1
564 WERR( JE ) = TMP2
565 INDEXW( JE ) = JE - IWOFF
566 IBLOCK( JE ) = IB
567 50 CONTINUE
568 60 CONTINUE
569 *
570 M = M + IM
571 END IF
572 70 CONTINUE
573
574 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
575 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
576 IF( IRANGE.EQ.INDRNG ) THEN
577 IDISCL = IL - 1 - NWL
578 IDISCU = NWU - IU
579 *
580 IF( IDISCL.GT.0 ) THEN
581 IM = 0
582 DO 80 JE = 1, M
583 * Remove some of the smallest eigenvalues from the left so that
584 * at the end IDISCL =0. Move all eigenvalues up to the left.
585 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
586 IDISCL = IDISCL - 1
587 ELSE
588 IM = IM + 1
589 W( IM ) = W( JE )
590 WERR( IM ) = WERR( JE )
591 INDEXW( IM ) = INDEXW( JE )
592 IBLOCK( IM ) = IBLOCK( JE )
593 END IF
594 80 CONTINUE
595 M = IM
596 END IF
597 IF( IDISCU.GT.0 ) THEN
598 * Remove some of the largest eigenvalues from the right so that
599 * at the end IDISCU =0. Move all eigenvalues up to the left.
600 IM=M+1
601 DO 81 JE = M, 1, -1
602 IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
603 IDISCU = IDISCU - 1
604 ELSE
605 IM = IM - 1
606 W( IM ) = W( JE )
607 WERR( IM ) = WERR( JE )
608 INDEXW( IM ) = INDEXW( JE )
609 IBLOCK( IM ) = IBLOCK( JE )
610 END IF
611 81 CONTINUE
612 JEE = 0
613 DO 82 JE = IM, M
614 JEE = JEE + 1
615 W( JEE ) = W( JE )
616 WERR( JEE ) = WERR( JE )
617 INDEXW( JEE ) = INDEXW( JE )
618 IBLOCK( JEE ) = IBLOCK( JE )
619 82 CONTINUE
620 M = M-IM+1
621 END IF
622
623 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
624 * Code to deal with effects of bad arithmetic. (If N(w) is
625 * monotone non-decreasing, this should never happen.)
626 * Some low eigenvalues to be discarded are not in (WL,WLU],
627 * or high eigenvalues to be discarded are not in (WUL,WU]
628 * so just kill off the smallest IDISCL/largest IDISCU
629 * eigenvalues, by marking the corresponding IBLOCK = 0
630 IF( IDISCL.GT.0 ) THEN
631 WKILL = WU
632 DO 100 JDISC = 1, IDISCL
633 IW = 0
634 DO 90 JE = 1, M
635 IF( IBLOCK( JE ).NE.0 .AND.
636 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
637 IW = JE
638 WKILL = W( JE )
639 END IF
640 90 CONTINUE
641 IBLOCK( IW ) = 0
642 100 CONTINUE
643 END IF
644 IF( IDISCU.GT.0 ) THEN
645 WKILL = WL
646 DO 120 JDISC = 1, IDISCU
647 IW = 0
648 DO 110 JE = 1, M
649 IF( IBLOCK( JE ).NE.0 .AND.
650 $ ( W( JE ).GE.WKILL .OR. IW.EQ.0 ) ) THEN
651 IW = JE
652 WKILL = W( JE )
653 END IF
654 110 CONTINUE
655 IBLOCK( IW ) = 0
656 120 CONTINUE
657 END IF
658 * Now erase all eigenvalues with IBLOCK set to zero
659 IM = 0
660 DO 130 JE = 1, M
661 IF( IBLOCK( JE ).NE.0 ) THEN
662 IM = IM + 1
663 W( IM ) = W( JE )
664 WERR( IM ) = WERR( JE )
665 INDEXW( IM ) = INDEXW( JE )
666 IBLOCK( IM ) = IBLOCK( JE )
667 END IF
668 130 CONTINUE
669 M = IM
670 END IF
671 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
672 TOOFEW = .TRUE.
673 END IF
674 END IF
675 *
676 IF(( IRANGE.EQ.ALLRNG .AND. M.NE.N ).OR.
677 $ ( IRANGE.EQ.INDRNG .AND. M.NE.IU-IL+1 ) ) THEN
678 TOOFEW = .TRUE.
679 END IF
680
681 * If ORDER='B', do nothing the eigenvalues are already sorted by
682 * block.
683 * If ORDER='E', sort the eigenvalues from smallest to largest
684
685 IF( LSAME(ORDER,'E') .AND. NSPLIT.GT.1 ) THEN
686 DO 150 JE = 1, M - 1
687 IE = 0
688 TMP1 = W( JE )
689 DO 140 J = JE + 1, M
690 IF( W( J ).LT.TMP1 ) THEN
691 IE = J
692 TMP1 = W( J )
693 END IF
694 140 CONTINUE
695 IF( IE.NE.0 ) THEN
696 TMP2 = WERR( IE )
697 ITMP1 = IBLOCK( IE )
698 ITMP2 = INDEXW( IE )
699 W( IE ) = W( JE )
700 WERR( IE ) = WERR( JE )
701 IBLOCK( IE ) = IBLOCK( JE )
702 INDEXW( IE ) = INDEXW( JE )
703 W( JE ) = TMP1
704 WERR( JE ) = TMP2
705 IBLOCK( JE ) = ITMP1
706 INDEXW( JE ) = ITMP2
707 END IF
708 150 CONTINUE
709 END IF
710 *
711 INFO = 0
712 IF( NCNVRG )
713 $ INFO = INFO + 1
714 IF( TOOFEW )
715 $ INFO = INFO + 2
716 RETURN
717 *
718 * End of DLARRD
719 *
720 END