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