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