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.2) --
8: * -- LAPACK is a software package provided by Univ. of Tennessee, --
9: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
10: * November 2006
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
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: ELSE
536: IF(MB.GT.1) THEN
537: CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
538: AVGAP = ABS(CLWDTH / DBLE(WEND-WBEGIN))
539: IF( SGNDEF.EQ.ONE ) THEN
540: TAU = HALF*MAX(WGAP(WBEGIN),AVGAP)
541: TAU = MAX(TAU,WERR(WBEGIN))
542: ELSE
543: TAU = HALF*MAX(WGAP(WEND-1),AVGAP)
544: TAU = MAX(TAU,WERR(WEND))
545: ENDIF
546: ELSE
547: TAU = WERR(WBEGIN)
548: ENDIF
549: ENDIF
550: *
551: DO 80 IDUM = 1, MAXTRY
552: * Compute L D L^T factorization of tridiagonal matrix T - sigma I.
553: * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
554: * pivots in WORK(2*IN+1:3*IN)
555: DPIVOT = D( IBEGIN ) - SIGMA
556: WORK( 1 ) = DPIVOT
557: DMAX = ABS( WORK(1) )
558: J = IBEGIN
559: DO 70 I = 1, IN - 1
560: WORK( 2*IN+I ) = ONE / WORK( I )
561: TMP = E( J )*WORK( 2*IN+I )
562: WORK( IN+I ) = TMP
563: DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J )
564: WORK( I+1 ) = DPIVOT
565: DMAX = MAX( DMAX, ABS(DPIVOT) )
566: J = J + 1
567: 70 CONTINUE
568: * check for element growth
569: IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN
570: NOREP = .TRUE.
571: ELSE
572: NOREP = .FALSE.
573: ENDIF
574: IF( USEDQD .AND. .NOT.NOREP ) THEN
575: * Ensure the definiteness of the representation
576: * All entries of D (of L D L^T) must have the same sign
577: DO 71 I = 1, IN
578: TMP = SGNDEF*WORK( I )
579: IF( TMP.LT.ZERO ) NOREP = .TRUE.
580: 71 CONTINUE
581: ENDIF
582: IF(NOREP) THEN
583: * Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
584: * shift which makes the matrix definite. So we should end up
585: * here really only in the case of IRANGE = VALRNG or INDRNG.
586: IF( IDUM.EQ.MAXTRY-1 ) THEN
587: IF( SGNDEF.EQ.ONE ) THEN
588: * The fudged Gerschgorin shift should succeed
589: SIGMA =
590: $ GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN
591: ELSE
592: SIGMA =
593: $ GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN
594: END IF
595: ELSE
596: SIGMA = SIGMA - SGNDEF * TAU
597: TAU = TWO * TAU
598: END IF
599: ELSE
600: * an initial RRR is found
601: GO TO 83
602: END IF
603: 80 CONTINUE
604: * if the program reaches this point, no base representation could be
605: * found in MAXTRY iterations.
606: INFO = 2
607: RETURN
608:
609: 83 CONTINUE
610: * At this point, we have found an initial base representation
611: * T - SIGMA I = L D L^T with not too much element growth.
612: * Store the shift.
613: E( IEND ) = SIGMA
614: * Store D and L.
615: CALL DCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
616: CALL DCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )
617:
618:
619: IF(MB.GT.1 ) THEN
620: *
621: * Perturb each entry of the base representation by a small
622: * (but random) relative amount to overcome difficulties with
623: * glued matrices.
624: *
625: DO 122 I = 1, 4
626: ISEED( I ) = 1
627: 122 CONTINUE
628:
629: CALL DLARNV(2, ISEED, 2*IN-1, WORK(1))
630: DO 125 I = 1,IN-1
631: D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
632: E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
633: 125 CONTINUE
634: D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN))
635: *
636: ENDIF
637: *
638: * Don't update the Gerschgorin intervals because keeping track
639: * of the updates would be too much work in DLARRV.
640: * We update W instead and use it to locate the proper Gerschgorin
641: * intervals.
642:
643: * Compute the required eigenvalues of L D L' by bisection or dqds
644: IF ( .NOT.USEDQD ) THEN
645: * If DLARRD has been used, shift the eigenvalue approximations
646: * according to their representation. This is necessary for
647: * a uniform DLARRV since dqds computes eigenvalues of the
648: * shifted representation. In DLARRV, W will always hold the
649: * UNshifted eigenvalue approximation.
650: DO 134 J=WBEGIN,WEND
651: W(J) = W(J) - SIGMA
652: WERR(J) = WERR(J) + ABS(W(J)) * EPS
653: 134 CONTINUE
654: * call DLARRB to reduce eigenvalue error of the approximations
655: * from DLARRD
656: DO 135 I = IBEGIN, IEND-1
657: WORK( I ) = D( I ) * E( I )**2
658: 135 CONTINUE
659: * use bisection to find EV from INDL to INDU
660: CALL DLARRB(IN, D(IBEGIN), WORK(IBEGIN),
661: $ INDL, INDU, RTOL1, RTOL2, INDL-1,
662: $ W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN),
663: $ WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM,
664: $ IN, IINFO )
665: IF( IINFO .NE. 0 ) THEN
666: INFO = -4
667: RETURN
668: END IF
669: * DLARRB computes all gaps correctly except for the last one
670: * Record distance to VU/GU
671: WGAP( WEND ) = MAX( ZERO,
672: $ ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) )
673: DO 138 I = INDL, INDU
674: M = M + 1
675: IBLOCK(M) = JBLK
676: INDEXW(M) = I
677: 138 CONTINUE
678: ELSE
679: * Call dqds to get all eigs (and then possibly delete unwanted
680: * eigenvalues).
681: * Note that dqds finds the eigenvalues of the L D L^T representation
682: * of T to high relative accuracy. High relative accuracy
683: * might be lost when the shift of the RRR is subtracted to obtain
684: * the eigenvalues of T. However, T is not guaranteed to define its
685: * eigenvalues to high relative accuracy anyway.
686: * Set RTOL to the order of the tolerance used in DLASQ2
687: * This is an ESTIMATED error, the worst case bound is 4*N*EPS
688: * which is usually too large and requires unnecessary work to be
689: * done by bisection when computing the eigenvectors
690: RTOL = LOG(DBLE(IN)) * FOUR * EPS
691: J = IBEGIN
692: DO 140 I = 1, IN - 1
693: WORK( 2*I-1 ) = ABS( D( J ) )
694: WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
695: J = J + 1
696: 140 CONTINUE
697: WORK( 2*IN-1 ) = ABS( D( IEND ) )
698: WORK( 2*IN ) = ZERO
699: CALL DLASQ2( IN, WORK, IINFO )
700: IF( IINFO .NE. 0 ) THEN
701: * If IINFO = -5 then an index is part of a tight cluster
702: * and should be changed. The index is in IWORK(1) and the
703: * gap is in WORK(N+1)
704: INFO = -5
705: RETURN
706: ELSE
707: * Test that all eigenvalues are positive as expected
708: DO 149 I = 1, IN
709: IF( WORK( I ).LT.ZERO ) THEN
710: INFO = -6
711: RETURN
712: ENDIF
713: 149 CONTINUE
714: END IF
715: IF( SGNDEF.GT.ZERO ) THEN
716: DO 150 I = INDL, INDU
717: M = M + 1
718: W( M ) = WORK( IN-I+1 )
719: IBLOCK( M ) = JBLK
720: INDEXW( M ) = I
721: 150 CONTINUE
722: ELSE
723: DO 160 I = INDL, INDU
724: M = M + 1
725: W( M ) = -WORK( I )
726: IBLOCK( M ) = JBLK
727: INDEXW( M ) = I
728: 160 CONTINUE
729: END IF
730:
731: DO 165 I = M - MB + 1, M
732: * the value of RTOL below should be the tolerance in DLASQ2
733: WERR( I ) = RTOL * ABS( W(I) )
734: 165 CONTINUE
735: DO 166 I = M - MB + 1, M - 1
736: * compute the right gap between the intervals
737: WGAP( I ) = MAX( ZERO,
738: $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
739: 166 CONTINUE
740: WGAP( M ) = MAX( ZERO,
741: $ ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) )
742: END IF
743: * proceed with next block
744: IBEGIN = IEND + 1
745: WBEGIN = WEND + 1
746: 170 CONTINUE
747: *
748:
749: RETURN
750: *
751: * end of DLARRE
752: *
753: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarre.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack