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