Annotation of rpl/lapack/lapack/dlarrd.f, revision 1.18
1.12 bertrand 1: *> \brief \b DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.17 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.17 bertrand 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">
1.9 bertrand 15: *> [TXT]</a>
1.17 bertrand 16: *> \endhtmlonly
1.9 bertrand 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 )
1.17 bertrand 25: *
1.9 bertrand 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: * ..
1.17 bertrand 37: *
1.9 bertrand 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
1.15 bertrand 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'.
1.9 bertrand 99: *> \endverbatim
100: *>
101: *> \param[in] VU
102: *> \verbatim
103: *> VU is DOUBLE PRECISION
1.15 bertrand 104: *> If RANGE='V', the upper bound of the interval to
1.9 bertrand 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
1.15 bertrand 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'.
1.9 bertrand 117: *> \endverbatim
118: *>
119: *> \param[in] IU
120: *> \verbatim
121: *> IU is INTEGER
1.15 bertrand 122: *> If RANGE='I', the index of the
123: *> largest eigenvalue to be returned.
1.9 bertrand 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: *
1.17 bertrand 317: *> \author Univ. of Tennessee
318: *> \author Univ. of California Berkeley
319: *> \author Univ. of Colorado Denver
320: *> \author NAG Ltd.
1.9 bertrand 321: *
1.15 bertrand 322: *> \date June 2016
1.9 bertrand 323: *
1.17 bertrand 324: *> \ingroup OTHERauxiliary
1.9 bertrand 325: *
326: * =====================================================================
1.1 bertrand 327: SUBROUTINE DLARRD( RANGE, ORDER, N, VL, VU, IL, IU, GERS,
328: $ RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
329: $ M, W, WERR, WL, WU, IBLOCK, INDEXW,
330: $ WORK, IWORK, INFO )
331: *
1.17 bertrand 332: * -- LAPACK auxiliary routine (version 3.7.0) --
1.1 bertrand 333: * -- LAPACK is a software package provided by Univ. of Tennessee, --
334: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.15 bertrand 335: * June 2016
1.1 bertrand 336: *
337: * .. Scalar Arguments ..
338: CHARACTER ORDER, RANGE
339: INTEGER IL, INFO, IU, M, N, NSPLIT
340: DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU
341: * ..
342: * .. Array Arguments ..
343: INTEGER IBLOCK( * ), INDEXW( * ),
344: $ ISPLIT( * ), IWORK( * )
345: DOUBLE PRECISION D( * ), E( * ), E2( * ),
346: $ GERS( * ), W( * ), WERR( * ), WORK( * )
347: * ..
348: *
349: * =====================================================================
350: *
351: * .. Parameters ..
352: DOUBLE PRECISION ZERO, ONE, TWO, HALF, FUDGE
353: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0,
354: $ TWO = 2.0D0, HALF = ONE/TWO,
355: $ FUDGE = TWO )
356: INTEGER ALLRNG, VALRNG, INDRNG
357: PARAMETER ( ALLRNG = 1, VALRNG = 2, INDRNG = 3 )
358: * ..
359: * .. Local Scalars ..
360: LOGICAL NCNVRG, TOOFEW
361: INTEGER I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
362: $ IM, IN, IOFF, IOUT, IRANGE, ITMAX, ITMP1,
363: $ ITMP2, IW, IWOFF, J, JBLK, JDISC, JE, JEE, NB,
364: $ NWL, NWU
365: DOUBLE PRECISION ATOLI, EPS, GL, GU, RTOLI, TMP1, TMP2,
366: $ TNORM, UFLOW, WKILL, WLU, WUL
367:
368: * ..
369: * .. Local Arrays ..
370: INTEGER IDUMMA( 1 )
371: * ..
372: * .. External Functions ..
373: LOGICAL LSAME
374: INTEGER ILAENV
375: DOUBLE PRECISION DLAMCH
376: EXTERNAL LSAME, ILAENV, DLAMCH
377: * ..
378: * .. External Subroutines ..
379: EXTERNAL DLAEBZ
380: * ..
381: * .. Intrinsic Functions ..
382: INTRINSIC ABS, INT, LOG, MAX, MIN
383: * ..
384: * .. Executable Statements ..
385: *
386: INFO = 0
387: *
388: * Decode RANGE
389: *
390: IF( LSAME( RANGE, 'A' ) ) THEN
391: IRANGE = ALLRNG
392: ELSE IF( LSAME( RANGE, 'V' ) ) THEN
393: IRANGE = VALRNG
394: ELSE IF( LSAME( RANGE, 'I' ) ) THEN
395: IRANGE = INDRNG
396: ELSE
397: IRANGE = 0
398: END IF
399: *
400: * Check for Errors
401: *
402: IF( IRANGE.LE.0 ) THEN
403: INFO = -1
404: ELSE IF( .NOT.(LSAME(ORDER,'B').OR.LSAME(ORDER,'E')) ) THEN
405: INFO = -2
406: ELSE IF( N.LT.0 ) THEN
407: INFO = -3
408: ELSE IF( IRANGE.EQ.VALRNG ) THEN
409: IF( VL.GE.VU )
410: $ INFO = -5
411: ELSE IF( IRANGE.EQ.INDRNG .AND.
412: $ ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) THEN
413: INFO = -6
414: ELSE IF( IRANGE.EQ.INDRNG .AND.
415: $ ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN
416: INFO = -7
417: END IF
418: *
419: IF( INFO.NE.0 ) THEN
420: RETURN
421: END IF
422:
423: * Initialize error flags
424: INFO = 0
425: NCNVRG = .FALSE.
426: TOOFEW = .FALSE.
427:
428: * Quick return if possible
429: M = 0
430: IF( N.EQ.0 ) RETURN
431:
432: * Simplification:
433: IF( IRANGE.EQ.INDRNG .AND. IL.EQ.1 .AND. IU.EQ.N ) IRANGE = 1
434:
435: * Get machine constants
436: EPS = DLAMCH( 'P' )
437: UFLOW = DLAMCH( 'U' )
438:
439:
440: * Special Case when N=1
441: * Treat case of 1x1 matrix for quick return
442: IF( N.EQ.1 ) THEN
443: IF( (IRANGE.EQ.ALLRNG).OR.
444: $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
445: $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
446: M = 1
447: W(1) = D(1)
448: * The computation error of the eigenvalue is zero
449: WERR(1) = ZERO
450: IBLOCK( 1 ) = 1
451: INDEXW( 1 ) = 1
452: ENDIF
453: RETURN
454: END IF
455:
456: * NB is the minimum vector length for vector bisection, or 0
457: * if only scalar is to be done.
458: NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 )
459: IF( NB.LE.1 ) NB = 0
460:
461: * Find global spectral radius
462: GL = D(1)
463: GU = D(1)
464: DO 5 I = 1,N
465: GL = MIN( GL, GERS( 2*I - 1))
466: GU = MAX( GU, GERS(2*I) )
467: 5 CONTINUE
468: * Compute global Gerschgorin bounds and spectral diameter
469: TNORM = MAX( ABS( GL ), ABS( GU ) )
470: GL = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
471: GU = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
472: * [JAN/28/2009] remove the line below since SPDIAM variable not use
473: * SPDIAM = GU - GL
474: * Input arguments for DLAEBZ:
475: * The relative tolerance. An interval (a,b] lies within
476: * "relative tolerance" if b-a < RELTOL*max(|a|,|b|),
477: RTOLI = RELTOL
478: * Set the absolute tolerance for interval convergence to zero to force
479: * interval convergence based on relative size of the interval.
480: * This is dangerous because intervals might not converge when RELTOL is
481: * small. But at least a very small number should be selected so that for
482: * strongly graded matrices, the code can get relatively accurate
483: * eigenvalues.
484: ATOLI = FUDGE*TWO*UFLOW + FUDGE*TWO*PIVMIN
485:
486: IF( IRANGE.EQ.INDRNG ) THEN
487:
488: * RANGE='I': Compute an interval containing eigenvalues
489: * IL through IU. The initial interval [GL,GU] from the global
490: * Gerschgorin bounds GL and GU is refined by DLAEBZ.
491: ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
492: $ LOG( TWO ) ) + 2
493: WORK( N+1 ) = GL
494: WORK( N+2 ) = GL
495: WORK( N+3 ) = GU
496: WORK( N+4 ) = GU
497: WORK( N+5 ) = GL
498: WORK( N+6 ) = GU
499: IWORK( 1 ) = -1
500: IWORK( 2 ) = -1
501: IWORK( 3 ) = N + 1
502: IWORK( 4 ) = N + 1
503: IWORK( 5 ) = IL - 1
504: IWORK( 6 ) = IU
505: *
506: CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN,
507: $ D, E, E2, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
508: $ IWORK, W, IBLOCK, IINFO )
509: IF( IINFO .NE. 0 ) THEN
510: INFO = IINFO
511: RETURN
512: END IF
513: * On exit, output intervals may not be ordered by ascending negcount
514: IF( IWORK( 6 ).EQ.IU ) THEN
515: WL = WORK( N+1 )
516: WLU = WORK( N+3 )
517: NWL = IWORK( 1 )
518: WU = WORK( N+4 )
519: WUL = WORK( N+2 )
520: NWU = IWORK( 4 )
521: ELSE
522: WL = WORK( N+2 )
523: WLU = WORK( N+4 )
524: NWL = IWORK( 2 )
525: WU = WORK( N+3 )
526: WUL = WORK( N+1 )
527: NWU = IWORK( 3 )
528: END IF
529: * On exit, the interval [WL, WLU] contains a value with negcount NWL,
530: * and [WUL, WU] contains a value with negcount NWU.
531: IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
532: INFO = 4
533: RETURN
534: END IF
535:
536: ELSEIF( IRANGE.EQ.VALRNG ) THEN
537: WL = VL
538: WU = VU
539:
540: ELSEIF( IRANGE.EQ.ALLRNG ) THEN
541: WL = GL
542: WU = GU
543: ENDIF
544:
545:
546:
547: * Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
548: * NWL accumulates the number of eigenvalues .le. WL,
549: * NWU accumulates the number of eigenvalues .le. WU
550: M = 0
551: IEND = 0
552: INFO = 0
553: NWL = 0
554: NWU = 0
555: *
556: DO 70 JBLK = 1, NSPLIT
557: IOFF = IEND
558: IBEGIN = IOFF + 1
559: IEND = ISPLIT( JBLK )
560: IN = IEND - IOFF
561: *
562: IF( IN.EQ.1 ) THEN
563: * 1x1 block
564: IF( WL.GE.D( IBEGIN )-PIVMIN )
565: $ NWL = NWL + 1
566: IF( WU.GE.D( IBEGIN )-PIVMIN )
567: $ NWU = NWU + 1
568: IF( IRANGE.EQ.ALLRNG .OR.
569: $ ( WL.LT.D( IBEGIN )-PIVMIN
570: $ .AND. WU.GE. D( IBEGIN )-PIVMIN ) ) THEN
571: M = M + 1
572: W( M ) = D( IBEGIN )
573: WERR(M) = ZERO
574: * The gap for a single block doesn't matter for the later
575: * algorithm and is assigned an arbitrary large value
576: IBLOCK( M ) = JBLK
577: INDEXW( M ) = 1
578: END IF
579:
580: * Disabled 2x2 case because of a failure on the following matrix
581: * RANGE = 'I', IL = IU = 4
582: * Original Tridiagonal, d = [
583: * -0.150102010615740E+00
584: * -0.849897989384260E+00
585: * -0.128208148052635E-15
586: * 0.128257718286320E-15
587: * ];
588: * e = [
589: * -0.357171383266986E+00
590: * -0.180411241501588E-15
591: * -0.175152352710251E-15
592: * ];
593: *
594: * ELSE IF( IN.EQ.2 ) THEN
595: ** 2x2 block
596: * DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
597: * TMP1 = HALF*(D(IBEGIN)+D(IEND))
598: * L1 = TMP1 - DISC
599: * IF( WL.GE. L1-PIVMIN )
600: * $ NWL = NWL + 1
601: * IF( WU.GE. L1-PIVMIN )
602: * $ NWU = NWU + 1
603: * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
604: * $ L1-PIVMIN ) ) THEN
605: * M = M + 1
606: * W( M ) = L1
607: ** The uncertainty of eigenvalues of a 2x2 matrix is very small
608: * WERR( M ) = EPS * ABS( W( M ) ) * TWO
609: * IBLOCK( M ) = JBLK
610: * INDEXW( M ) = 1
611: * ENDIF
612: * L2 = TMP1 + DISC
613: * IF( WL.GE. L2-PIVMIN )
614: * $ NWL = NWL + 1
615: * IF( WU.GE. L2-PIVMIN )
616: * $ NWU = NWU + 1
617: * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
618: * $ L2-PIVMIN ) ) THEN
619: * M = M + 1
620: * W( M ) = L2
621: ** The uncertainty of eigenvalues of a 2x2 matrix is very small
622: * WERR( M ) = EPS * ABS( W( M ) ) * TWO
623: * IBLOCK( M ) = JBLK
624: * INDEXW( M ) = 2
625: * ENDIF
626: ELSE
627: * General Case - block of size IN >= 2
628: * Compute local Gerschgorin interval and use it as the initial
629: * interval for DLAEBZ
630: GU = D( IBEGIN )
631: GL = D( IBEGIN )
632: TMP1 = ZERO
633:
634: DO 40 J = IBEGIN, IEND
635: GL = MIN( GL, GERS( 2*J - 1))
636: GU = MAX( GU, GERS(2*J) )
637: 40 CONTINUE
638: * [JAN/28/2009]
639: * change SPDIAM by TNORM in lines 2 and 3 thereafter
640: * line 1: remove computation of SPDIAM (not useful anymore)
641: * SPDIAM = GU - GL
642: * GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
643: * GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
644: GL = GL - FUDGE*TNORM*EPS*IN - FUDGE*PIVMIN
645: GU = GU + FUDGE*TNORM*EPS*IN + FUDGE*PIVMIN
646: *
647: IF( IRANGE.GT.1 ) THEN
648: IF( GU.LT.WL ) THEN
649: * the local block contains none of the wanted eigenvalues
650: NWL = NWL + IN
651: NWU = NWU + IN
652: GO TO 70
653: END IF
654: * refine search interval if possible, only range (WL,WU] matters
655: GL = MAX( GL, WL )
656: GU = MIN( GU, WU )
657: IF( GL.GE.GU )
658: $ GO TO 70
659: END IF
660:
661: * Find negcount of initial interval boundaries GL and GU
662: WORK( N+1 ) = GL
663: WORK( N+IN+1 ) = GU
664: CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
665: $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
666: $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
667: $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
668: IF( IINFO .NE. 0 ) THEN
669: INFO = IINFO
670: RETURN
671: END IF
672: *
673: NWL = NWL + IWORK( 1 )
674: NWU = NWU + IWORK( IN+1 )
675: IWOFF = M - IWORK( 1 )
676:
677: * Compute Eigenvalues
678: ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
679: $ LOG( TWO ) ) + 2
680: CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
681: $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
682: $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
683: $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
684: IF( IINFO .NE. 0 ) THEN
685: INFO = IINFO
686: RETURN
687: END IF
688: *
689: * Copy eigenvalues into W and IBLOCK
690: * Use -JBLK for block number for unconverged eigenvalues.
691: * Loop over the number of output intervals from DLAEBZ
692: DO 60 J = 1, IOUT
693: * eigenvalue approximation is middle point of interval
694: TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
695: * semi length of error interval
696: TMP2 = HALF*ABS( WORK( J+N )-WORK( J+IN+N ) )
697: IF( J.GT.IOUT-IINFO ) THEN
698: * Flag non-convergence.
699: NCNVRG = .TRUE.
700: IB = -JBLK
701: ELSE
702: IB = JBLK
703: END IF
704: DO 50 JE = IWORK( J ) + 1 + IWOFF,
705: $ IWORK( J+IN ) + IWOFF
706: W( JE ) = TMP1
707: WERR( JE ) = TMP2
708: INDEXW( JE ) = JE - IWOFF
709: IBLOCK( JE ) = IB
710: 50 CONTINUE
711: 60 CONTINUE
712: *
713: M = M + IM
714: END IF
715: 70 CONTINUE
716:
717: * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
718: * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
719: IF( IRANGE.EQ.INDRNG ) THEN
720: IDISCL = IL - 1 - NWL
721: IDISCU = NWU - IU
722: *
723: IF( IDISCL.GT.0 ) THEN
724: IM = 0
725: DO 80 JE = 1, M
726: * Remove some of the smallest eigenvalues from the left so that
727: * at the end IDISCL =0. Move all eigenvalues up to the left.
728: IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
729: IDISCL = IDISCL - 1
730: ELSE
731: IM = IM + 1
732: W( IM ) = W( JE )
733: WERR( IM ) = WERR( JE )
734: INDEXW( IM ) = INDEXW( JE )
735: IBLOCK( IM ) = IBLOCK( JE )
736: END IF
737: 80 CONTINUE
738: M = IM
739: END IF
740: IF( IDISCU.GT.0 ) THEN
741: * Remove some of the largest eigenvalues from the right so that
742: * at the end IDISCU =0. Move all eigenvalues up to the left.
743: IM=M+1
744: DO 81 JE = M, 1, -1
745: IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
746: IDISCU = IDISCU - 1
747: ELSE
748: IM = IM - 1
749: W( IM ) = W( JE )
750: WERR( IM ) = WERR( JE )
751: INDEXW( IM ) = INDEXW( JE )
752: IBLOCK( IM ) = IBLOCK( JE )
753: END IF
754: 81 CONTINUE
755: JEE = 0
756: DO 82 JE = IM, M
757: JEE = JEE + 1
758: W( JEE ) = W( JE )
759: WERR( JEE ) = WERR( JE )
760: INDEXW( JEE ) = INDEXW( JE )
761: IBLOCK( JEE ) = IBLOCK( JE )
762: 82 CONTINUE
763: M = M-IM+1
764: END IF
765:
766: IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
767: * Code to deal with effects of bad arithmetic. (If N(w) is
768: * monotone non-decreasing, this should never happen.)
769: * Some low eigenvalues to be discarded are not in (WL,WLU],
770: * or high eigenvalues to be discarded are not in (WUL,WU]
771: * so just kill off the smallest IDISCL/largest IDISCU
772: * eigenvalues, by marking the corresponding IBLOCK = 0
773: IF( IDISCL.GT.0 ) THEN
774: WKILL = WU
775: DO 100 JDISC = 1, IDISCL
776: IW = 0
777: DO 90 JE = 1, M
778: IF( IBLOCK( JE ).NE.0 .AND.
779: $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
780: IW = JE
781: WKILL = W( JE )
782: END IF
783: 90 CONTINUE
784: IBLOCK( IW ) = 0
785: 100 CONTINUE
786: END IF
787: IF( IDISCU.GT.0 ) THEN
788: WKILL = WL
789: DO 120 JDISC = 1, IDISCU
790: IW = 0
791: DO 110 JE = 1, M
792: IF( IBLOCK( JE ).NE.0 .AND.
793: $ ( W( JE ).GE.WKILL .OR. IW.EQ.0 ) ) THEN
794: IW = JE
795: WKILL = W( JE )
796: END IF
797: 110 CONTINUE
798: IBLOCK( IW ) = 0
799: 120 CONTINUE
800: END IF
801: * Now erase all eigenvalues with IBLOCK set to zero
802: IM = 0
803: DO 130 JE = 1, M
804: IF( IBLOCK( JE ).NE.0 ) THEN
805: IM = IM + 1
806: W( IM ) = W( JE )
807: WERR( IM ) = WERR( JE )
808: INDEXW( IM ) = INDEXW( JE )
809: IBLOCK( IM ) = IBLOCK( JE )
810: END IF
811: 130 CONTINUE
812: M = IM
813: END IF
814: IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
815: TOOFEW = .TRUE.
816: END IF
817: END IF
818: *
819: IF(( IRANGE.EQ.ALLRNG .AND. M.NE.N ).OR.
820: $ ( IRANGE.EQ.INDRNG .AND. M.NE.IU-IL+1 ) ) THEN
821: TOOFEW = .TRUE.
822: END IF
823:
824: * If ORDER='B', do nothing the eigenvalues are already sorted by
825: * block.
826: * If ORDER='E', sort the eigenvalues from smallest to largest
827:
828: IF( LSAME(ORDER,'E') .AND. NSPLIT.GT.1 ) THEN
829: DO 150 JE = 1, M - 1
830: IE = 0
831: TMP1 = W( JE )
832: DO 140 J = JE + 1, M
833: IF( W( J ).LT.TMP1 ) THEN
834: IE = J
835: TMP1 = W( J )
836: END IF
837: 140 CONTINUE
838: IF( IE.NE.0 ) THEN
839: TMP2 = WERR( IE )
840: ITMP1 = IBLOCK( IE )
841: ITMP2 = INDEXW( IE )
842: W( IE ) = W( JE )
843: WERR( IE ) = WERR( JE )
844: IBLOCK( IE ) = IBLOCK( JE )
845: INDEXW( IE ) = INDEXW( JE )
846: W( JE ) = TMP1
847: WERR( JE ) = TMP2
848: IBLOCK( JE ) = ITMP1
849: INDEXW( JE ) = ITMP2
850: END IF
851: 150 CONTINUE
852: END IF
853: *
854: INFO = 0
855: IF( NCNVRG )
856: $ INFO = INFO + 1
857: IF( TOOFEW )
858: $ INFO = INFO + 2
859: RETURN
860: *
861: * End of DLARRD
862: *
863: 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