1 SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
2 $ RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
3 $ W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
4 $ WORK, IWORK, INFO )
5 IMPLICIT NONE
6 *
7 * -- LAPACK auxiliary routine (version 3.3.1) --
8 * -- LAPACK is a software package provided by Univ. of Tennessee, --
9 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
10 * -- April 2011 --
11 *
12 * .. Scalar Arguments ..
13 CHARACTER RANGE
14 INTEGER IL, INFO, IU, M, N, NSPLIT
15 DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
16 * ..
17 * .. Array Arguments ..
18 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
19 $ INDEXW( * )
20 DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
21 $ W( * ),WERR( * ), WGAP( * ), WORK( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * To find the desired eigenvalues of a given real symmetric
28 * tridiagonal matrix T, DLARRE sets any "small" off-diagonal
29 * elements to zero, and for each unreduced block T_i, it finds
30 * (a) a suitable shift at one end of the block's spectrum,
31 * (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
32 * (c) eigenvalues of each L_i D_i L_i^T.
33 * The representations and eigenvalues found are then used by
34 * DSTEMR to compute the eigenvectors of T.
35 * The accuracy varies depending on whether bisection is used to
36 * find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to
37 * conpute all and then discard any unwanted one.
38 * As an added benefit, DLARRE also outputs the n
39 * Gerschgorin intervals for the matrices L_i D_i L_i^T.
40 *
41 * Arguments
42 * =========
43 *
44 * RANGE (input) CHARACTER*1
45 * = 'A': ("All") all eigenvalues will be found.
46 * = 'V': ("Value") all eigenvalues in the half-open interval
47 * (VL, VU] will be found.
48 * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
49 * entire matrix) will be found.
50 *
51 * N (input) INTEGER
52 * The order of the matrix. N > 0.
53 *
54 * VL (input/output) DOUBLE PRECISION
55 * VU (input/output) DOUBLE PRECISION
56 * If RANGE='V', the lower and upper bounds for the eigenvalues.
57 * Eigenvalues less than or equal to VL, or greater than VU,
58 * will not be returned. VL < VU.
59 * If RANGE='I' or ='A', DLARRE computes bounds on the desired
60 * part of the spectrum.
61 *
62 * IL (input) INTEGER
63 * IU (input) INTEGER
64 * If RANGE='I', the indices (in ascending order) of the
65 * smallest and largest eigenvalues to be returned.
66 * 1 <= IL <= IU <= N.
67 *
68 * D (input/output) DOUBLE PRECISION array, dimension (N)
69 * On entry, the N diagonal elements of the tridiagonal
70 * matrix T.
71 * On exit, the N diagonal elements of the diagonal
72 * matrices D_i.
73 *
74 * E (input/output) DOUBLE PRECISION array, dimension (N)
75 * On entry, the first (N-1) entries contain the subdiagonal
76 * elements of the tridiagonal matrix T; E(N) need not be set.
77 * On exit, E contains the subdiagonal elements of the unit
78 * bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
79 * 1 <= I <= NSPLIT, contain the base points sigma_i on output.
80 *
81 * E2 (input/output) DOUBLE PRECISION array, dimension (N)
82 * On entry, the first (N-1) entries contain the SQUARES of the
83 * subdiagonal elements of the tridiagonal matrix T;
84 * E2(N) need not be set.
85 * On exit, the entries E2( ISPLIT( I ) ),
86 * 1 <= I <= NSPLIT, have been set to zero
87 *
88 * RTOL1 (input) DOUBLE PRECISION
89 * RTOL2 (input) DOUBLE PRECISION
90 * Parameters for bisection.
91 * An interval [LEFT,RIGHT] has converged if
92 * RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
93 *
94 * SPLTOL (input) DOUBLE PRECISION
95 * The threshold for splitting.
96 *
97 * NSPLIT (output) INTEGER
98 * The number of blocks T splits into. 1 <= NSPLIT <= N.
99 *
100 * ISPLIT (output) INTEGER array, dimension (N)
101 * The splitting points, at which T breaks up into blocks.
102 * The first block consists of rows/columns 1 to ISPLIT(1),
103 * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
104 * etc., and the NSPLIT-th consists of rows/columns
105 * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
106 *
107 * M (output) INTEGER
108 * The total number of eigenvalues (of all L_i D_i L_i^T)
109 * found.
110 *
111 * W (output) DOUBLE PRECISION array, dimension (N)
112 * The first M elements contain the eigenvalues. The
113 * eigenvalues of each of the blocks, L_i D_i L_i^T, are
114 * sorted in ascending order ( DLARRE may use the
115 * remaining N-M elements as workspace).
116 *
117 * WERR (output) DOUBLE PRECISION array, dimension (N)
118 * The error bound on the corresponding eigenvalue in W.
119 *
120 * WGAP (output) DOUBLE PRECISION array, dimension (N)
121 * The separation from the right neighbor eigenvalue in W.
122 * The gap is only with respect to the eigenvalues of the same block
123 * as each block has its own representation tree.
124 * Exception: at the right end of a block we store the left gap
125 *
126 * IBLOCK (output) INTEGER array, dimension (N)
127 * The indices of the blocks (submatrices) associated with the
128 * corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
129 * W(i) belongs to the first block from the top, =2 if W(i)
130 * belongs to the second block, etc.
131 *
132 * INDEXW (output) INTEGER array, dimension (N)
133 * The indices of the eigenvalues within each block (submatrix);
134 * for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
135 * i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
136 *
137 * GERS (output) DOUBLE PRECISION array, dimension (2*N)
138 * The N Gerschgorin intervals (the i-th Gerschgorin interval
139 * is (GERS(2*i-1), GERS(2*i)).
140 *
141 * PIVMIN (output) DOUBLE PRECISION
142 * The minimum pivot in the Sturm sequence for T.
143 *
144 * WORK (workspace) DOUBLE PRECISION array, dimension (6*N)
145 * Workspace.
146 *
147 * IWORK (workspace) INTEGER array, dimension (5*N)
148 * Workspace.
149 *
150 * INFO (output) INTEGER
151 * = 0: successful exit
152 * > 0: A problem occured in DLARRE.
153 * < 0: One of the called subroutines signaled an internal problem.
154 * Needs inspection of the corresponding parameter IINFO
155 * for further information.
156 *
157 * =-1: Problem in DLARRD.
158 * = 2: No base representation could be found in MAXTRY iterations.
159 * Increasing MAXTRY and recompilation might be a remedy.
160 * =-3: Problem in DLARRB when computing the refined root
161 * representation for DLASQ2.
162 * =-4: Problem in DLARRB when preforming bisection on the
163 * desired part of the spectrum.
164 * =-5: Problem in DLASQ2.
165 * =-6: Problem in DLASQ2.
166 *
167 * Further Details
168 * The base representations are required to suffer very little
169 * element growth and consequently define all their eigenvalues to
170 * high relative accuracy.
171 * ===============
172 *
173 * Based on contributions by
174 * Beresford Parlett, University of California, Berkeley, USA
175 * Jim Demmel, University of California, Berkeley, USA
176 * Inderjit Dhillon, University of Texas, Austin, USA
177 * Osni Marques, LBNL/NERSC, USA
178 * Christof Voemel, University of California, Berkeley, USA
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183 DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
184 $ MAXGROWTH, ONE, PERT, TWO, ZERO
185 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0,
186 $ TWO = 2.0D0, FOUR=4.0D0,
187 $ HNDRD = 100.0D0,
188 $ PERT = 8.0D0,
189 $ HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF,
190 $ MAXGROWTH = 64.0D0, FUDGE = 2.0D0 )
191 INTEGER MAXTRY, ALLRNG, INDRNG, VALRNG
192 PARAMETER ( MAXTRY = 6, ALLRNG = 1, INDRNG = 2,
193 $ VALRNG = 3 )
194 * ..
195 * .. Local Scalars ..
196 LOGICAL FORCEB, NOREP, USEDQD
197 INTEGER CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
198 $ IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
199 $ WBEGIN, WEND
200 DOUBLE PRECISION AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
201 $ EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL,
202 $ RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM,
203 $ TAU, TMP, TMP1
204
205
206 * ..
207 * .. Local Arrays ..
208 INTEGER ISEED( 4 )
209 * ..
210 * .. External Functions ..
211 LOGICAL LSAME
212 DOUBLE PRECISION DLAMCH
213 EXTERNAL DLAMCH, LSAME
214
215 * ..
216 * .. External Subroutines ..
217 EXTERNAL DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD,
218 $ DLASQ2
219 * ..
220 * .. Intrinsic Functions ..
221 INTRINSIC ABS, MAX, MIN
222
223 * ..
224 * .. Executable Statements ..
225 *
226
227 INFO = 0
228
229 *
230 * Decode RANGE
231 *
232 IF( LSAME( RANGE, 'A' ) ) THEN
233 IRANGE = ALLRNG
234 ELSE IF( LSAME( RANGE, 'V' ) ) THEN
235 IRANGE = VALRNG
236 ELSE IF( LSAME( RANGE, 'I' ) ) THEN
237 IRANGE = INDRNG
238 END IF
239
240 M = 0
241
242 * Get machine constants
243 SAFMIN = DLAMCH( 'S' )
244 EPS = DLAMCH( 'P' )
245
246 * Set parameters
247 RTL = SQRT(EPS)
248 BSRTOL = SQRT(EPS)
249
250 * Treat case of 1x1 matrix for quick return
251 IF( N.EQ.1 ) THEN
252 IF( (IRANGE.EQ.ALLRNG).OR.
253 $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
254 $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
255 M = 1
256 W(1) = D(1)
257 * The computation error of the eigenvalue is zero
258 WERR(1) = ZERO
259 WGAP(1) = ZERO
260 IBLOCK( 1 ) = 1
261 INDEXW( 1 ) = 1
262 GERS(1) = D( 1 )
263 GERS(2) = D( 1 )
264 ENDIF
265 * store the shift for the initial RRR, which is zero in this case
266 E(1) = ZERO
267 RETURN
268 END IF
269
270 * General case: tridiagonal matrix of order > 1
271 *
272 * Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
273 * Compute maximum off-diagonal entry and pivmin.
274 GL = D(1)
275 GU = D(1)
276 EOLD = ZERO
277 EMAX = ZERO
278 E(N) = ZERO
279 DO 5 I = 1,N
280 WERR(I) = ZERO
281 WGAP(I) = ZERO
282 EABS = ABS( E(I) )
283 IF( EABS .GE. EMAX ) THEN
284 EMAX = EABS
285 END IF
286 TMP1 = EABS + EOLD
287 GERS( 2*I-1) = D(I) - TMP1
288 GL = MIN( GL, GERS( 2*I - 1))
289 GERS( 2*I ) = D(I) + TMP1
290 GU = MAX( GU, GERS(2*I) )
291 EOLD = EABS
292 5 CONTINUE
293 * The minimum pivot allowed in the Sturm sequence for T
294 PIVMIN = SAFMIN * MAX( ONE, EMAX**2 )
295 * Compute spectral diameter. The Gerschgorin bounds give an
296 * estimate that is wrong by at most a factor of SQRT(2)
297 SPDIAM = GU - GL
298
299 * Compute splitting points
300 CALL DLARRA( N, D, E, E2, SPLTOL, SPDIAM,
301 $ NSPLIT, ISPLIT, IINFO )
302
303 * Can force use of bisection instead of faster DQDS.
304 * Option left in the code for future multisection work.
305 FORCEB = .FALSE.
306
307 * Initialize USEDQD, DQDS should be used for ALLRNG unless someone
308 * explicitly wants bisection.
309 USEDQD = (( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB))
310
311 IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN
312 * Set interval [VL,VU] that contains all eigenvalues
313 VL = GL
314 VU = GU
315 ELSE
316 * We call DLARRD to find crude approximations to the eigenvalues
317 * in the desired range. In case IRANGE = INDRNG, we also obtain the
318 * interval (VL,VU] that contains all the wanted eigenvalues.
319 * An interval [LEFT,RIGHT] has converged if
320 * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
321 * DLARRD needs a WORK of size 4*N, IWORK of size 3*N
322 CALL DLARRD( RANGE, 'B', N, VL, VU, IL, IU, GERS,
323 $ BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
324 $ MM, W, WERR, VL, VU, IBLOCK, INDEXW,
325 $ WORK, IWORK, IINFO )
326 IF( IINFO.NE.0 ) THEN
327 INFO = -1
328 RETURN
329 ENDIF
330 * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
331 DO 14 I = MM+1,N
332 W( I ) = ZERO
333 WERR( I ) = ZERO
334 IBLOCK( I ) = 0
335 INDEXW( I ) = 0
336 14 CONTINUE
337 END IF
338
339
340 ***
341 * Loop over unreduced blocks
342 IBEGIN = 1
343 WBEGIN = 1
344 DO 170 JBLK = 1, NSPLIT
345 IEND = ISPLIT( JBLK )
346 IN = IEND - IBEGIN + 1
347
348 * 1 X 1 block
349 IF( IN.EQ.1 ) THEN
350 IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND.
351 $ ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) )
352 $ .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK))
353 $ ) THEN
354 M = M + 1
355 W( M ) = D( IBEGIN )
356 WERR(M) = ZERO
357 * The gap for a single block doesn't matter for the later
358 * algorithm and is assigned an arbitrary large value
359 WGAP(M) = ZERO
360 IBLOCK( M ) = JBLK
361 INDEXW( M ) = 1
362 WBEGIN = WBEGIN + 1
363 ENDIF
364 * E( IEND ) holds the shift for the initial RRR
365 E( IEND ) = ZERO
366 IBEGIN = IEND + 1
367 GO TO 170
368 END IF
369 *
370 * Blocks of size larger than 1x1
371 *
372 * E( IEND ) will hold the shift for the initial RRR, for now set it =0
373 E( IEND ) = ZERO
374 *
375 * Find local outer bounds GL,GU for the block
376 GL = D(IBEGIN)
377 GU = D(IBEGIN)
378 DO 15 I = IBEGIN , IEND
379 GL = MIN( GERS( 2*I-1 ), GL )
380 GU = MAX( GERS( 2*I ), GU )
381 15 CONTINUE
382 SPDIAM = GU - GL
383
384 IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN
385 * Count the number of eigenvalues in the current block.
386 MB = 0
387 DO 20 I = WBEGIN,MM
388 IF( IBLOCK(I).EQ.JBLK ) THEN
389 MB = MB+1
390 ELSE
391 GOTO 21
392 ENDIF
393 20 CONTINUE
394 21 CONTINUE
395
396 IF( MB.EQ.0) THEN
397 * No eigenvalue in the current block lies in the desired range
398 * E( IEND ) holds the shift for the initial RRR
399 E( IEND ) = ZERO
400 IBEGIN = IEND + 1
401 GO TO 170
402 ELSE
403
404 * Decide whether dqds or bisection is more efficient
405 USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) )
406 WEND = WBEGIN + MB - 1
407 * Calculate gaps for the current block
408 * In later stages, when representations for individual
409 * eigenvalues are different, we use SIGMA = E( IEND ).
410 SIGMA = ZERO
411 DO 30 I = WBEGIN, WEND - 1
412 WGAP( I ) = MAX( ZERO,
413 $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
414 30 CONTINUE
415 WGAP( WEND ) = MAX( ZERO,
416 $ VU - SIGMA - (W( WEND )+WERR( WEND )))
417 * Find local index of the first and last desired evalue.
418 INDL = INDEXW(WBEGIN)
419 INDU = INDEXW( WEND )
420 ENDIF
421 ENDIF
422 IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN
423 * Case of DQDS
424 * Find approximations to the extremal eigenvalues of the block
425 CALL DLARRK( IN, 1, GL, GU, D(IBEGIN),
426 $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
427 IF( IINFO.NE.0 ) THEN
428 INFO = -1
429 RETURN
430 ENDIF
431 ISLEFT = MAX(GL, TMP - TMP1
432 $ - HNDRD * EPS* ABS(TMP - TMP1))
433
434 CALL DLARRK( IN, IN, GL, GU, D(IBEGIN),
435 $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
436 IF( IINFO.NE.0 ) THEN
437 INFO = -1
438 RETURN
439 ENDIF
440 ISRGHT = MIN(GU, TMP + TMP1
441 $ + HNDRD * EPS * ABS(TMP + TMP1))
442 * Improve the estimate of the spectral diameter
443 SPDIAM = ISRGHT - ISLEFT
444 ELSE
445 * Case of bisection
446 * Find approximations to the wanted extremal eigenvalues
447 ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN)
448 $ - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) ))
449 ISRGHT = MIN(GU,W(WEND) + WERR(WEND)
450 $ + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND)))
451 ENDIF
452
453
454 * Decide whether the base representation for the current block
455 * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
456 * should be on the left or the right end of the current block.
457 * The strategy is to shift to the end which is "more populated"
458 * Furthermore, decide whether to use DQDS for the computation of
459 * the eigenvalue approximations at the end of DLARRE or bisection.
460 * dqds is chosen if all eigenvalues are desired or the number of
461 * eigenvalues to be computed is large compared to the blocksize.
462 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
463 * If all the eigenvalues have to be computed, we use dqd
464 USEDQD = .TRUE.
465 * INDL is the local index of the first eigenvalue to compute
466 INDL = 1
467 INDU = IN
468 * MB = number of eigenvalues to compute
469 MB = IN
470 WEND = WBEGIN + MB - 1
471 * Define 1/4 and 3/4 points of the spectrum
472 S1 = ISLEFT + FOURTH * SPDIAM
473 S2 = ISRGHT - FOURTH * SPDIAM
474 ELSE
475 * DLARRD has computed IBLOCK and INDEXW for each eigenvalue
476 * approximation.
477 * choose sigma
478 IF( USEDQD ) THEN
479 S1 = ISLEFT + FOURTH * SPDIAM
480 S2 = ISRGHT - FOURTH * SPDIAM
481 ELSE
482 TMP = MIN(ISRGHT,VU) - MAX(ISLEFT,VL)
483 S1 = MAX(ISLEFT,VL) + FOURTH * TMP
484 S2 = MIN(ISRGHT,VU) - FOURTH * TMP
485 ENDIF
486 ENDIF
487
488 * Compute the negcount at the 1/4 and 3/4 points
489 IF(MB.GT.1) THEN
490 CALL DLARRC( 'T', IN, S1, S2, D(IBEGIN),
491 $ E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO)
492 ENDIF
493
494 IF(MB.EQ.1) THEN
495 SIGMA = GL
496 SGNDEF = ONE
497 ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN
498 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
499 SIGMA = MAX(ISLEFT,GL)
500 ELSEIF( USEDQD ) THEN
501 * use Gerschgorin bound as shift to get pos def matrix
502 * for dqds
503 SIGMA = ISLEFT
504 ELSE
505 * use approximation of the first desired eigenvalue of the
506 * block as shift
507 SIGMA = MAX(ISLEFT,VL)
508 ENDIF
509 SGNDEF = ONE
510 ELSE
511 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
512 SIGMA = MIN(ISRGHT,GU)
513 ELSEIF( USEDQD ) THEN
514 * use Gerschgorin bound as shift to get neg def matrix
515 * for dqds
516 SIGMA = ISRGHT
517 ELSE
518 * use approximation of the first desired eigenvalue of the
519 * block as shift
520 SIGMA = MIN(ISRGHT,VU)
521 ENDIF
522 SGNDEF = -ONE
523 ENDIF
524
525
526 * An initial SIGMA has been chosen that will be used for computing
527 * T - SIGMA I = L D L^T
528 * Define the increment TAU of the shift in case the initial shift
529 * needs to be refined to obtain a factorization with not too much
530 * element growth.
531 IF( USEDQD ) THEN
532 * The initial SIGMA was to the outer end of the spectrum
533 * the matrix is definite and we need not retreat.
534 TAU = SPDIAM*EPS*N + TWO*PIVMIN
535 TAU = MAX( TAU,TWO*EPS*ABS(SIGMA) )
536 ELSE
537 IF(MB.GT.1) THEN
538 CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
539 AVGAP = ABS(CLWDTH / DBLE(WEND-WBEGIN))
540 IF( SGNDEF.EQ.ONE ) THEN
541 TAU = HALF*MAX(WGAP(WBEGIN),AVGAP)
542 TAU = MAX(TAU,WERR(WBEGIN))
543 ELSE
544 TAU = HALF*MAX(WGAP(WEND-1),AVGAP)
545 TAU = MAX(TAU,WERR(WEND))
546 ENDIF
547 ELSE
548 TAU = WERR(WBEGIN)
549 ENDIF
550 ENDIF
551 *
552 DO 80 IDUM = 1, MAXTRY
553 * Compute L D L^T factorization of tridiagonal matrix T - sigma I.
554 * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
555 * pivots in WORK(2*IN+1:3*IN)
556 DPIVOT = D( IBEGIN ) - SIGMA
557 WORK( 1 ) = DPIVOT
558 DMAX = ABS( WORK(1) )
559 J = IBEGIN
560 DO 70 I = 1, IN - 1
561 WORK( 2*IN+I ) = ONE / WORK( I )
562 TMP = E( J )*WORK( 2*IN+I )
563 WORK( IN+I ) = TMP
564 DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J )
565 WORK( I+1 ) = DPIVOT
566 DMAX = MAX( DMAX, ABS(DPIVOT) )
567 J = J + 1
568 70 CONTINUE
569 * check for element growth
570 IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN
571 NOREP = .TRUE.
572 ELSE
573 NOREP = .FALSE.
574 ENDIF
575 IF( USEDQD .AND. .NOT.NOREP ) THEN
576 * Ensure the definiteness of the representation
577 * All entries of D (of L D L^T) must have the same sign
578 DO 71 I = 1, IN
579 TMP = SGNDEF*WORK( I )
580 IF( TMP.LT.ZERO ) NOREP = .TRUE.
581 71 CONTINUE
582 ENDIF
583 IF(NOREP) THEN
584 * Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
585 * shift which makes the matrix definite. So we should end up
586 * here really only in the case of IRANGE = VALRNG or INDRNG.
587 IF( IDUM.EQ.MAXTRY-1 ) THEN
588 IF( SGNDEF.EQ.ONE ) THEN
589 * The fudged Gerschgorin shift should succeed
590 SIGMA =
591 $ GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN
592 ELSE
593 SIGMA =
594 $ GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN
595 END IF
596 ELSE
597 SIGMA = SIGMA - SGNDEF * TAU
598 TAU = TWO * TAU
599 END IF
600 ELSE
601 * an initial RRR is found
602 GO TO 83
603 END IF
604 80 CONTINUE
605 * if the program reaches this point, no base representation could be
606 * found in MAXTRY iterations.
607 INFO = 2
608 RETURN
609
610 83 CONTINUE
611 * At this point, we have found an initial base representation
612 * T - SIGMA I = L D L^T with not too much element growth.
613 * Store the shift.
614 E( IEND ) = SIGMA
615 * Store D and L.
616 CALL DCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
617 CALL DCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )
618
619
620 IF(MB.GT.1 ) THEN
621 *
622 * Perturb each entry of the base representation by a small
623 * (but random) relative amount to overcome difficulties with
624 * glued matrices.
625 *
626 DO 122 I = 1, 4
627 ISEED( I ) = 1
628 122 CONTINUE
629
630 CALL DLARNV(2, ISEED, 2*IN-1, WORK(1))
631 DO 125 I = 1,IN-1
632 D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
633 E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
634 125 CONTINUE
635 D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN))
636 *
637 ENDIF
638 *
639 * Don't update the Gerschgorin intervals because keeping track
640 * of the updates would be too much work in DLARRV.
641 * We update W instead and use it to locate the proper Gerschgorin
642 * intervals.
643
644 * Compute the required eigenvalues of L D L' by bisection or dqds
645 IF ( .NOT.USEDQD ) THEN
646 * If DLARRD has been used, shift the eigenvalue approximations
647 * according to their representation. This is necessary for
648 * a uniform DLARRV since dqds computes eigenvalues of the
649 * shifted representation. In DLARRV, W will always hold the
650 * UNshifted eigenvalue approximation.
651 DO 134 J=WBEGIN,WEND
652 W(J) = W(J) - SIGMA
653 WERR(J) = WERR(J) + ABS(W(J)) * EPS
654 134 CONTINUE
655 * call DLARRB to reduce eigenvalue error of the approximations
656 * from DLARRD
657 DO 135 I = IBEGIN, IEND-1
658 WORK( I ) = D( I ) * E( I )**2
659 135 CONTINUE
660 * use bisection to find EV from INDL to INDU
661 CALL DLARRB(IN, D(IBEGIN), WORK(IBEGIN),
662 $ INDL, INDU, RTOL1, RTOL2, INDL-1,
663 $ W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN),
664 $ WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM,
665 $ IN, IINFO )
666 IF( IINFO .NE. 0 ) THEN
667 INFO = -4
668 RETURN
669 END IF
670 * DLARRB computes all gaps correctly except for the last one
671 * Record distance to VU/GU
672 WGAP( WEND ) = MAX( ZERO,
673 $ ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) )
674 DO 138 I = INDL, INDU
675 M = M + 1
676 IBLOCK(M) = JBLK
677 INDEXW(M) = I
678 138 CONTINUE
679 ELSE
680 * Call dqds to get all eigs (and then possibly delete unwanted
681 * eigenvalues).
682 * Note that dqds finds the eigenvalues of the L D L^T representation
683 * of T to high relative accuracy. High relative accuracy
684 * might be lost when the shift of the RRR is subtracted to obtain
685 * the eigenvalues of T. However, T is not guaranteed to define its
686 * eigenvalues to high relative accuracy anyway.
687 * Set RTOL to the order of the tolerance used in DLASQ2
688 * This is an ESTIMATED error, the worst case bound is 4*N*EPS
689 * which is usually too large and requires unnecessary work to be
690 * done by bisection when computing the eigenvectors
691 RTOL = LOG(DBLE(IN)) * FOUR * EPS
692 J = IBEGIN
693 DO 140 I = 1, IN - 1
694 WORK( 2*I-1 ) = ABS( D( J ) )
695 WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
696 J = J + 1
697 140 CONTINUE
698 WORK( 2*IN-1 ) = ABS( D( IEND ) )
699 WORK( 2*IN ) = ZERO
700 CALL DLASQ2( IN, WORK, IINFO )
701 IF( IINFO .NE. 0 ) THEN
702 * If IINFO = -5 then an index is part of a tight cluster
703 * and should be changed. The index is in IWORK(1) and the
704 * gap is in WORK(N+1)
705 INFO = -5
706 RETURN
707 ELSE
708 * Test that all eigenvalues are positive as expected
709 DO 149 I = 1, IN
710 IF( WORK( I ).LT.ZERO ) THEN
711 INFO = -6
712 RETURN
713 ENDIF
714 149 CONTINUE
715 END IF
716 IF( SGNDEF.GT.ZERO ) THEN
717 DO 150 I = INDL, INDU
718 M = M + 1
719 W( M ) = WORK( IN-I+1 )
720 IBLOCK( M ) = JBLK
721 INDEXW( M ) = I
722 150 CONTINUE
723 ELSE
724 DO 160 I = INDL, INDU
725 M = M + 1
726 W( M ) = -WORK( I )
727 IBLOCK( M ) = JBLK
728 INDEXW( M ) = I
729 160 CONTINUE
730 END IF
731
732 DO 165 I = M - MB + 1, M
733 * the value of RTOL below should be the tolerance in DLASQ2
734 WERR( I ) = RTOL * ABS( W(I) )
735 165 CONTINUE
736 DO 166 I = M - MB + 1, M - 1
737 * compute the right gap between the intervals
738 WGAP( I ) = MAX( ZERO,
739 $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
740 166 CONTINUE
741 WGAP( M ) = MAX( ZERO,
742 $ ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) )
743 END IF
744 * proceed with next block
745 IBEGIN = IEND + 1
746 WBEGIN = WEND + 1
747 170 CONTINUE
748 *
749
750 RETURN
751 *
752 * end of DLARRE
753 *
754 END
2 $ RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
3 $ W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
4 $ WORK, IWORK, INFO )
5 IMPLICIT NONE
6 *
7 * -- LAPACK auxiliary routine (version 3.3.1) --
8 * -- LAPACK is a software package provided by Univ. of Tennessee, --
9 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
10 * -- April 2011 --
11 *
12 * .. Scalar Arguments ..
13 CHARACTER RANGE
14 INTEGER IL, INFO, IU, M, N, NSPLIT
15 DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
16 * ..
17 * .. Array Arguments ..
18 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
19 $ INDEXW( * )
20 DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
21 $ W( * ),WERR( * ), WGAP( * ), WORK( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * To find the desired eigenvalues of a given real symmetric
28 * tridiagonal matrix T, DLARRE sets any "small" off-diagonal
29 * elements to zero, and for each unreduced block T_i, it finds
30 * (a) a suitable shift at one end of the block's spectrum,
31 * (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
32 * (c) eigenvalues of each L_i D_i L_i^T.
33 * The representations and eigenvalues found are then used by
34 * DSTEMR to compute the eigenvectors of T.
35 * The accuracy varies depending on whether bisection is used to
36 * find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to
37 * conpute all and then discard any unwanted one.
38 * As an added benefit, DLARRE also outputs the n
39 * Gerschgorin intervals for the matrices L_i D_i L_i^T.
40 *
41 * Arguments
42 * =========
43 *
44 * RANGE (input) CHARACTER*1
45 * = 'A': ("All") all eigenvalues will be found.
46 * = 'V': ("Value") all eigenvalues in the half-open interval
47 * (VL, VU] will be found.
48 * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
49 * entire matrix) will be found.
50 *
51 * N (input) INTEGER
52 * The order of the matrix. N > 0.
53 *
54 * VL (input/output) DOUBLE PRECISION
55 * VU (input/output) DOUBLE PRECISION
56 * If RANGE='V', the lower and upper bounds for the eigenvalues.
57 * Eigenvalues less than or equal to VL, or greater than VU,
58 * will not be returned. VL < VU.
59 * If RANGE='I' or ='A', DLARRE computes bounds on the desired
60 * part of the spectrum.
61 *
62 * IL (input) INTEGER
63 * IU (input) INTEGER
64 * If RANGE='I', the indices (in ascending order) of the
65 * smallest and largest eigenvalues to be returned.
66 * 1 <= IL <= IU <= N.
67 *
68 * D (input/output) DOUBLE PRECISION array, dimension (N)
69 * On entry, the N diagonal elements of the tridiagonal
70 * matrix T.
71 * On exit, the N diagonal elements of the diagonal
72 * matrices D_i.
73 *
74 * E (input/output) DOUBLE PRECISION array, dimension (N)
75 * On entry, the first (N-1) entries contain the subdiagonal
76 * elements of the tridiagonal matrix T; E(N) need not be set.
77 * On exit, E contains the subdiagonal elements of the unit
78 * bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
79 * 1 <= I <= NSPLIT, contain the base points sigma_i on output.
80 *
81 * E2 (input/output) DOUBLE PRECISION array, dimension (N)
82 * On entry, the first (N-1) entries contain the SQUARES of the
83 * subdiagonal elements of the tridiagonal matrix T;
84 * E2(N) need not be set.
85 * On exit, the entries E2( ISPLIT( I ) ),
86 * 1 <= I <= NSPLIT, have been set to zero
87 *
88 * RTOL1 (input) DOUBLE PRECISION
89 * RTOL2 (input) DOUBLE PRECISION
90 * Parameters for bisection.
91 * An interval [LEFT,RIGHT] has converged if
92 * RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
93 *
94 * SPLTOL (input) DOUBLE PRECISION
95 * The threshold for splitting.
96 *
97 * NSPLIT (output) INTEGER
98 * The number of blocks T splits into. 1 <= NSPLIT <= N.
99 *
100 * ISPLIT (output) INTEGER array, dimension (N)
101 * The splitting points, at which T breaks up into blocks.
102 * The first block consists of rows/columns 1 to ISPLIT(1),
103 * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
104 * etc., and the NSPLIT-th consists of rows/columns
105 * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
106 *
107 * M (output) INTEGER
108 * The total number of eigenvalues (of all L_i D_i L_i^T)
109 * found.
110 *
111 * W (output) DOUBLE PRECISION array, dimension (N)
112 * The first M elements contain the eigenvalues. The
113 * eigenvalues of each of the blocks, L_i D_i L_i^T, are
114 * sorted in ascending order ( DLARRE may use the
115 * remaining N-M elements as workspace).
116 *
117 * WERR (output) DOUBLE PRECISION array, dimension (N)
118 * The error bound on the corresponding eigenvalue in W.
119 *
120 * WGAP (output) DOUBLE PRECISION array, dimension (N)
121 * The separation from the right neighbor eigenvalue in W.
122 * The gap is only with respect to the eigenvalues of the same block
123 * as each block has its own representation tree.
124 * Exception: at the right end of a block we store the left gap
125 *
126 * IBLOCK (output) INTEGER array, dimension (N)
127 * The indices of the blocks (submatrices) associated with the
128 * corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
129 * W(i) belongs to the first block from the top, =2 if W(i)
130 * belongs to the second block, etc.
131 *
132 * INDEXW (output) INTEGER array, dimension (N)
133 * The indices of the eigenvalues within each block (submatrix);
134 * for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
135 * i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
136 *
137 * GERS (output) DOUBLE PRECISION array, dimension (2*N)
138 * The N Gerschgorin intervals (the i-th Gerschgorin interval
139 * is (GERS(2*i-1), GERS(2*i)).
140 *
141 * PIVMIN (output) DOUBLE PRECISION
142 * The minimum pivot in the Sturm sequence for T.
143 *
144 * WORK (workspace) DOUBLE PRECISION array, dimension (6*N)
145 * Workspace.
146 *
147 * IWORK (workspace) INTEGER array, dimension (5*N)
148 * Workspace.
149 *
150 * INFO (output) INTEGER
151 * = 0: successful exit
152 * > 0: A problem occured in DLARRE.
153 * < 0: One of the called subroutines signaled an internal problem.
154 * Needs inspection of the corresponding parameter IINFO
155 * for further information.
156 *
157 * =-1: Problem in DLARRD.
158 * = 2: No base representation could be found in MAXTRY iterations.
159 * Increasing MAXTRY and recompilation might be a remedy.
160 * =-3: Problem in DLARRB when computing the refined root
161 * representation for DLASQ2.
162 * =-4: Problem in DLARRB when preforming bisection on the
163 * desired part of the spectrum.
164 * =-5: Problem in DLASQ2.
165 * =-6: Problem in DLASQ2.
166 *
167 * Further Details
168 * The base representations are required to suffer very little
169 * element growth and consequently define all their eigenvalues to
170 * high relative accuracy.
171 * ===============
172 *
173 * Based on contributions by
174 * Beresford Parlett, University of California, Berkeley, USA
175 * Jim Demmel, University of California, Berkeley, USA
176 * Inderjit Dhillon, University of Texas, Austin, USA
177 * Osni Marques, LBNL/NERSC, USA
178 * Christof Voemel, University of California, Berkeley, USA
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183 DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
184 $ MAXGROWTH, ONE, PERT, TWO, ZERO
185 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0,
186 $ TWO = 2.0D0, FOUR=4.0D0,
187 $ HNDRD = 100.0D0,
188 $ PERT = 8.0D0,
189 $ HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF,
190 $ MAXGROWTH = 64.0D0, FUDGE = 2.0D0 )
191 INTEGER MAXTRY, ALLRNG, INDRNG, VALRNG
192 PARAMETER ( MAXTRY = 6, ALLRNG = 1, INDRNG = 2,
193 $ VALRNG = 3 )
194 * ..
195 * .. Local Scalars ..
196 LOGICAL FORCEB, NOREP, USEDQD
197 INTEGER CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
198 $ IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
199 $ WBEGIN, WEND
200 DOUBLE PRECISION AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
201 $ EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL,
202 $ RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM,
203 $ TAU, TMP, TMP1
204
205
206 * ..
207 * .. Local Arrays ..
208 INTEGER ISEED( 4 )
209 * ..
210 * .. External Functions ..
211 LOGICAL LSAME
212 DOUBLE PRECISION DLAMCH
213 EXTERNAL DLAMCH, LSAME
214
215 * ..
216 * .. External Subroutines ..
217 EXTERNAL DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD,
218 $ DLASQ2
219 * ..
220 * .. Intrinsic Functions ..
221 INTRINSIC ABS, MAX, MIN
222
223 * ..
224 * .. Executable Statements ..
225 *
226
227 INFO = 0
228
229 *
230 * Decode RANGE
231 *
232 IF( LSAME( RANGE, 'A' ) ) THEN
233 IRANGE = ALLRNG
234 ELSE IF( LSAME( RANGE, 'V' ) ) THEN
235 IRANGE = VALRNG
236 ELSE IF( LSAME( RANGE, 'I' ) ) THEN
237 IRANGE = INDRNG
238 END IF
239
240 M = 0
241
242 * Get machine constants
243 SAFMIN = DLAMCH( 'S' )
244 EPS = DLAMCH( 'P' )
245
246 * Set parameters
247 RTL = SQRT(EPS)
248 BSRTOL = SQRT(EPS)
249
250 * Treat case of 1x1 matrix for quick return
251 IF( N.EQ.1 ) THEN
252 IF( (IRANGE.EQ.ALLRNG).OR.
253 $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
254 $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
255 M = 1
256 W(1) = D(1)
257 * The computation error of the eigenvalue is zero
258 WERR(1) = ZERO
259 WGAP(1) = ZERO
260 IBLOCK( 1 ) = 1
261 INDEXW( 1 ) = 1
262 GERS(1) = D( 1 )
263 GERS(2) = D( 1 )
264 ENDIF
265 * store the shift for the initial RRR, which is zero in this case
266 E(1) = ZERO
267 RETURN
268 END IF
269
270 * General case: tridiagonal matrix of order > 1
271 *
272 * Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
273 * Compute maximum off-diagonal entry and pivmin.
274 GL = D(1)
275 GU = D(1)
276 EOLD = ZERO
277 EMAX = ZERO
278 E(N) = ZERO
279 DO 5 I = 1,N
280 WERR(I) = ZERO
281 WGAP(I) = ZERO
282 EABS = ABS( E(I) )
283 IF( EABS .GE. EMAX ) THEN
284 EMAX = EABS
285 END IF
286 TMP1 = EABS + EOLD
287 GERS( 2*I-1) = D(I) - TMP1
288 GL = MIN( GL, GERS( 2*I - 1))
289 GERS( 2*I ) = D(I) + TMP1
290 GU = MAX( GU, GERS(2*I) )
291 EOLD = EABS
292 5 CONTINUE
293 * The minimum pivot allowed in the Sturm sequence for T
294 PIVMIN = SAFMIN * MAX( ONE, EMAX**2 )
295 * Compute spectral diameter. The Gerschgorin bounds give an
296 * estimate that is wrong by at most a factor of SQRT(2)
297 SPDIAM = GU - GL
298
299 * Compute splitting points
300 CALL DLARRA( N, D, E, E2, SPLTOL, SPDIAM,
301 $ NSPLIT, ISPLIT, IINFO )
302
303 * Can force use of bisection instead of faster DQDS.
304 * Option left in the code for future multisection work.
305 FORCEB = .FALSE.
306
307 * Initialize USEDQD, DQDS should be used for ALLRNG unless someone
308 * explicitly wants bisection.
309 USEDQD = (( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB))
310
311 IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN
312 * Set interval [VL,VU] that contains all eigenvalues
313 VL = GL
314 VU = GU
315 ELSE
316 * We call DLARRD to find crude approximations to the eigenvalues
317 * in the desired range. In case IRANGE = INDRNG, we also obtain the
318 * interval (VL,VU] that contains all the wanted eigenvalues.
319 * An interval [LEFT,RIGHT] has converged if
320 * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
321 * DLARRD needs a WORK of size 4*N, IWORK of size 3*N
322 CALL DLARRD( RANGE, 'B', N, VL, VU, IL, IU, GERS,
323 $ BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
324 $ MM, W, WERR, VL, VU, IBLOCK, INDEXW,
325 $ WORK, IWORK, IINFO )
326 IF( IINFO.NE.0 ) THEN
327 INFO = -1
328 RETURN
329 ENDIF
330 * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
331 DO 14 I = MM+1,N
332 W( I ) = ZERO
333 WERR( I ) = ZERO
334 IBLOCK( I ) = 0
335 INDEXW( I ) = 0
336 14 CONTINUE
337 END IF
338
339
340 ***
341 * Loop over unreduced blocks
342 IBEGIN = 1
343 WBEGIN = 1
344 DO 170 JBLK = 1, NSPLIT
345 IEND = ISPLIT( JBLK )
346 IN = IEND - IBEGIN + 1
347
348 * 1 X 1 block
349 IF( IN.EQ.1 ) THEN
350 IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND.
351 $ ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) )
352 $ .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK))
353 $ ) THEN
354 M = M + 1
355 W( M ) = D( IBEGIN )
356 WERR(M) = ZERO
357 * The gap for a single block doesn't matter for the later
358 * algorithm and is assigned an arbitrary large value
359 WGAP(M) = ZERO
360 IBLOCK( M ) = JBLK
361 INDEXW( M ) = 1
362 WBEGIN = WBEGIN + 1
363 ENDIF
364 * E( IEND ) holds the shift for the initial RRR
365 E( IEND ) = ZERO
366 IBEGIN = IEND + 1
367 GO TO 170
368 END IF
369 *
370 * Blocks of size larger than 1x1
371 *
372 * E( IEND ) will hold the shift for the initial RRR, for now set it =0
373 E( IEND ) = ZERO
374 *
375 * Find local outer bounds GL,GU for the block
376 GL = D(IBEGIN)
377 GU = D(IBEGIN)
378 DO 15 I = IBEGIN , IEND
379 GL = MIN( GERS( 2*I-1 ), GL )
380 GU = MAX( GERS( 2*I ), GU )
381 15 CONTINUE
382 SPDIAM = GU - GL
383
384 IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN
385 * Count the number of eigenvalues in the current block.
386 MB = 0
387 DO 20 I = WBEGIN,MM
388 IF( IBLOCK(I).EQ.JBLK ) THEN
389 MB = MB+1
390 ELSE
391 GOTO 21
392 ENDIF
393 20 CONTINUE
394 21 CONTINUE
395
396 IF( MB.EQ.0) THEN
397 * No eigenvalue in the current block lies in the desired range
398 * E( IEND ) holds the shift for the initial RRR
399 E( IEND ) = ZERO
400 IBEGIN = IEND + 1
401 GO TO 170
402 ELSE
403
404 * Decide whether dqds or bisection is more efficient
405 USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) )
406 WEND = WBEGIN + MB - 1
407 * Calculate gaps for the current block
408 * In later stages, when representations for individual
409 * eigenvalues are different, we use SIGMA = E( IEND ).
410 SIGMA = ZERO
411 DO 30 I = WBEGIN, WEND - 1
412 WGAP( I ) = MAX( ZERO,
413 $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
414 30 CONTINUE
415 WGAP( WEND ) = MAX( ZERO,
416 $ VU - SIGMA - (W( WEND )+WERR( WEND )))
417 * Find local index of the first and last desired evalue.
418 INDL = INDEXW(WBEGIN)
419 INDU = INDEXW( WEND )
420 ENDIF
421 ENDIF
422 IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN
423 * Case of DQDS
424 * Find approximations to the extremal eigenvalues of the block
425 CALL DLARRK( IN, 1, GL, GU, D(IBEGIN),
426 $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
427 IF( IINFO.NE.0 ) THEN
428 INFO = -1
429 RETURN
430 ENDIF
431 ISLEFT = MAX(GL, TMP - TMP1
432 $ - HNDRD * EPS* ABS(TMP - TMP1))
433
434 CALL DLARRK( IN, IN, GL, GU, D(IBEGIN),
435 $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
436 IF( IINFO.NE.0 ) THEN
437 INFO = -1
438 RETURN
439 ENDIF
440 ISRGHT = MIN(GU, TMP + TMP1
441 $ + HNDRD * EPS * ABS(TMP + TMP1))
442 * Improve the estimate of the spectral diameter
443 SPDIAM = ISRGHT - ISLEFT
444 ELSE
445 * Case of bisection
446 * Find approximations to the wanted extremal eigenvalues
447 ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN)
448 $ - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) ))
449 ISRGHT = MIN(GU,W(WEND) + WERR(WEND)
450 $ + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND)))
451 ENDIF
452
453
454 * Decide whether the base representation for the current block
455 * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
456 * should be on the left or the right end of the current block.
457 * The strategy is to shift to the end which is "more populated"
458 * Furthermore, decide whether to use DQDS for the computation of
459 * the eigenvalue approximations at the end of DLARRE or bisection.
460 * dqds is chosen if all eigenvalues are desired or the number of
461 * eigenvalues to be computed is large compared to the blocksize.
462 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
463 * If all the eigenvalues have to be computed, we use dqd
464 USEDQD = .TRUE.
465 * INDL is the local index of the first eigenvalue to compute
466 INDL = 1
467 INDU = IN
468 * MB = number of eigenvalues to compute
469 MB = IN
470 WEND = WBEGIN + MB - 1
471 * Define 1/4 and 3/4 points of the spectrum
472 S1 = ISLEFT + FOURTH * SPDIAM
473 S2 = ISRGHT - FOURTH * SPDIAM
474 ELSE
475 * DLARRD has computed IBLOCK and INDEXW for each eigenvalue
476 * approximation.
477 * choose sigma
478 IF( USEDQD ) THEN
479 S1 = ISLEFT + FOURTH * SPDIAM
480 S2 = ISRGHT - FOURTH * SPDIAM
481 ELSE
482 TMP = MIN(ISRGHT,VU) - MAX(ISLEFT,VL)
483 S1 = MAX(ISLEFT,VL) + FOURTH * TMP
484 S2 = MIN(ISRGHT,VU) - FOURTH * TMP
485 ENDIF
486 ENDIF
487
488 * Compute the negcount at the 1/4 and 3/4 points
489 IF(MB.GT.1) THEN
490 CALL DLARRC( 'T', IN, S1, S2, D(IBEGIN),
491 $ E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO)
492 ENDIF
493
494 IF(MB.EQ.1) THEN
495 SIGMA = GL
496 SGNDEF = ONE
497 ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN
498 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
499 SIGMA = MAX(ISLEFT,GL)
500 ELSEIF( USEDQD ) THEN
501 * use Gerschgorin bound as shift to get pos def matrix
502 * for dqds
503 SIGMA = ISLEFT
504 ELSE
505 * use approximation of the first desired eigenvalue of the
506 * block as shift
507 SIGMA = MAX(ISLEFT,VL)
508 ENDIF
509 SGNDEF = ONE
510 ELSE
511 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
512 SIGMA = MIN(ISRGHT,GU)
513 ELSEIF( USEDQD ) THEN
514 * use Gerschgorin bound as shift to get neg def matrix
515 * for dqds
516 SIGMA = ISRGHT
517 ELSE
518 * use approximation of the first desired eigenvalue of the
519 * block as shift
520 SIGMA = MIN(ISRGHT,VU)
521 ENDIF
522 SGNDEF = -ONE
523 ENDIF
524
525
526 * An initial SIGMA has been chosen that will be used for computing
527 * T - SIGMA I = L D L^T
528 * Define the increment TAU of the shift in case the initial shift
529 * needs to be refined to obtain a factorization with not too much
530 * element growth.
531 IF( USEDQD ) THEN
532 * The initial SIGMA was to the outer end of the spectrum
533 * the matrix is definite and we need not retreat.
534 TAU = SPDIAM*EPS*N + TWO*PIVMIN
535 TAU = MAX( TAU,TWO*EPS*ABS(SIGMA) )
536 ELSE
537 IF(MB.GT.1) THEN
538 CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
539 AVGAP = ABS(CLWDTH / DBLE(WEND-WBEGIN))
540 IF( SGNDEF.EQ.ONE ) THEN
541 TAU = HALF*MAX(WGAP(WBEGIN),AVGAP)
542 TAU = MAX(TAU,WERR(WBEGIN))
543 ELSE
544 TAU = HALF*MAX(WGAP(WEND-1),AVGAP)
545 TAU = MAX(TAU,WERR(WEND))
546 ENDIF
547 ELSE
548 TAU = WERR(WBEGIN)
549 ENDIF
550 ENDIF
551 *
552 DO 80 IDUM = 1, MAXTRY
553 * Compute L D L^T factorization of tridiagonal matrix T - sigma I.
554 * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
555 * pivots in WORK(2*IN+1:3*IN)
556 DPIVOT = D( IBEGIN ) - SIGMA
557 WORK( 1 ) = DPIVOT
558 DMAX = ABS( WORK(1) )
559 J = IBEGIN
560 DO 70 I = 1, IN - 1
561 WORK( 2*IN+I ) = ONE / WORK( I )
562 TMP = E( J )*WORK( 2*IN+I )
563 WORK( IN+I ) = TMP
564 DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J )
565 WORK( I+1 ) = DPIVOT
566 DMAX = MAX( DMAX, ABS(DPIVOT) )
567 J = J + 1
568 70 CONTINUE
569 * check for element growth
570 IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN
571 NOREP = .TRUE.
572 ELSE
573 NOREP = .FALSE.
574 ENDIF
575 IF( USEDQD .AND. .NOT.NOREP ) THEN
576 * Ensure the definiteness of the representation
577 * All entries of D (of L D L^T) must have the same sign
578 DO 71 I = 1, IN
579 TMP = SGNDEF*WORK( I )
580 IF( TMP.LT.ZERO ) NOREP = .TRUE.
581 71 CONTINUE
582 ENDIF
583 IF(NOREP) THEN
584 * Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
585 * shift which makes the matrix definite. So we should end up
586 * here really only in the case of IRANGE = VALRNG or INDRNG.
587 IF( IDUM.EQ.MAXTRY-1 ) THEN
588 IF( SGNDEF.EQ.ONE ) THEN
589 * The fudged Gerschgorin shift should succeed
590 SIGMA =
591 $ GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN
592 ELSE
593 SIGMA =
594 $ GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN
595 END IF
596 ELSE
597 SIGMA = SIGMA - SGNDEF * TAU
598 TAU = TWO * TAU
599 END IF
600 ELSE
601 * an initial RRR is found
602 GO TO 83
603 END IF
604 80 CONTINUE
605 * if the program reaches this point, no base representation could be
606 * found in MAXTRY iterations.
607 INFO = 2
608 RETURN
609
610 83 CONTINUE
611 * At this point, we have found an initial base representation
612 * T - SIGMA I = L D L^T with not too much element growth.
613 * Store the shift.
614 E( IEND ) = SIGMA
615 * Store D and L.
616 CALL DCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
617 CALL DCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )
618
619
620 IF(MB.GT.1 ) THEN
621 *
622 * Perturb each entry of the base representation by a small
623 * (but random) relative amount to overcome difficulties with
624 * glued matrices.
625 *
626 DO 122 I = 1, 4
627 ISEED( I ) = 1
628 122 CONTINUE
629
630 CALL DLARNV(2, ISEED, 2*IN-1, WORK(1))
631 DO 125 I = 1,IN-1
632 D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
633 E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
634 125 CONTINUE
635 D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN))
636 *
637 ENDIF
638 *
639 * Don't update the Gerschgorin intervals because keeping track
640 * of the updates would be too much work in DLARRV.
641 * We update W instead and use it to locate the proper Gerschgorin
642 * intervals.
643
644 * Compute the required eigenvalues of L D L' by bisection or dqds
645 IF ( .NOT.USEDQD ) THEN
646 * If DLARRD has been used, shift the eigenvalue approximations
647 * according to their representation. This is necessary for
648 * a uniform DLARRV since dqds computes eigenvalues of the
649 * shifted representation. In DLARRV, W will always hold the
650 * UNshifted eigenvalue approximation.
651 DO 134 J=WBEGIN,WEND
652 W(J) = W(J) - SIGMA
653 WERR(J) = WERR(J) + ABS(W(J)) * EPS
654 134 CONTINUE
655 * call DLARRB to reduce eigenvalue error of the approximations
656 * from DLARRD
657 DO 135 I = IBEGIN, IEND-1
658 WORK( I ) = D( I ) * E( I )**2
659 135 CONTINUE
660 * use bisection to find EV from INDL to INDU
661 CALL DLARRB(IN, D(IBEGIN), WORK(IBEGIN),
662 $ INDL, INDU, RTOL1, RTOL2, INDL-1,
663 $ W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN),
664 $ WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM,
665 $ IN, IINFO )
666 IF( IINFO .NE. 0 ) THEN
667 INFO = -4
668 RETURN
669 END IF
670 * DLARRB computes all gaps correctly except for the last one
671 * Record distance to VU/GU
672 WGAP( WEND ) = MAX( ZERO,
673 $ ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) )
674 DO 138 I = INDL, INDU
675 M = M + 1
676 IBLOCK(M) = JBLK
677 INDEXW(M) = I
678 138 CONTINUE
679 ELSE
680 * Call dqds to get all eigs (and then possibly delete unwanted
681 * eigenvalues).
682 * Note that dqds finds the eigenvalues of the L D L^T representation
683 * of T to high relative accuracy. High relative accuracy
684 * might be lost when the shift of the RRR is subtracted to obtain
685 * the eigenvalues of T. However, T is not guaranteed to define its
686 * eigenvalues to high relative accuracy anyway.
687 * Set RTOL to the order of the tolerance used in DLASQ2
688 * This is an ESTIMATED error, the worst case bound is 4*N*EPS
689 * which is usually too large and requires unnecessary work to be
690 * done by bisection when computing the eigenvectors
691 RTOL = LOG(DBLE(IN)) * FOUR * EPS
692 J = IBEGIN
693 DO 140 I = 1, IN - 1
694 WORK( 2*I-1 ) = ABS( D( J ) )
695 WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
696 J = J + 1
697 140 CONTINUE
698 WORK( 2*IN-1 ) = ABS( D( IEND ) )
699 WORK( 2*IN ) = ZERO
700 CALL DLASQ2( IN, WORK, IINFO )
701 IF( IINFO .NE. 0 ) THEN
702 * If IINFO = -5 then an index is part of a tight cluster
703 * and should be changed. The index is in IWORK(1) and the
704 * gap is in WORK(N+1)
705 INFO = -5
706 RETURN
707 ELSE
708 * Test that all eigenvalues are positive as expected
709 DO 149 I = 1, IN
710 IF( WORK( I ).LT.ZERO ) THEN
711 INFO = -6
712 RETURN
713 ENDIF
714 149 CONTINUE
715 END IF
716 IF( SGNDEF.GT.ZERO ) THEN
717 DO 150 I = INDL, INDU
718 M = M + 1
719 W( M ) = WORK( IN-I+1 )
720 IBLOCK( M ) = JBLK
721 INDEXW( M ) = I
722 150 CONTINUE
723 ELSE
724 DO 160 I = INDL, INDU
725 M = M + 1
726 W( M ) = -WORK( I )
727 IBLOCK( M ) = JBLK
728 INDEXW( M ) = I
729 160 CONTINUE
730 END IF
731
732 DO 165 I = M - MB + 1, M
733 * the value of RTOL below should be the tolerance in DLASQ2
734 WERR( I ) = RTOL * ABS( W(I) )
735 165 CONTINUE
736 DO 166 I = M - MB + 1, M - 1
737 * compute the right gap between the intervals
738 WGAP( I ) = MAX( ZERO,
739 $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
740 166 CONTINUE
741 WGAP( M ) = MAX( ZERO,
742 $ ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) )
743 END IF
744 * proceed with next block
745 IBEGIN = IEND + 1
746 WBEGIN = WEND + 1
747 170 CONTINUE
748 *
749
750 RETURN
751 *
752 * end of DLARRE
753 *
754 END