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.2.1) --
7: * -- LAPACK is a software package provided by Univ. of Tennessee, --
8: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9: * -- April 2009 --
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
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
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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack