Annotation of rpl/lapack/lapack/dlarre.f, revision 1.24
1.14 bertrand 1: *> \brief \b DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues.
1.11 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.19 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.11 bertrand 7: *
8: *> \htmlonly
1.19 bertrand 9: *> Download DLARRE + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarre.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarre.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarre.f">
1.11 bertrand 15: *> [TXT]</a>
1.19 bertrand 16: *> \endhtmlonly
1.11 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
22: * RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
23: * W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
24: * WORK, IWORK, INFO )
1.19 bertrand 25: *
1.11 bertrand 26: * .. Scalar Arguments ..
27: * CHARACTER RANGE
28: * INTEGER IL, INFO, IU, M, N, NSPLIT
29: * DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
30: * ..
31: * .. Array Arguments ..
32: * INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
33: * $ INDEXW( * )
34: * DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
35: * $ W( * ),WERR( * ), WGAP( * ), WORK( * )
36: * ..
1.19 bertrand 37: *
1.11 bertrand 38: *
39: *> \par Purpose:
40: * =============
41: *>
42: *> \verbatim
43: *>
44: *> To find the desired eigenvalues of a given real symmetric
45: *> tridiagonal matrix T, DLARRE sets any "small" off-diagonal
46: *> elements to zero, and for each unreduced block T_i, it finds
47: *> (a) a suitable shift at one end of the block's spectrum,
48: *> (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
49: *> (c) eigenvalues of each L_i D_i L_i^T.
50: *> The representations and eigenvalues found are then used by
51: *> DSTEMR to compute the eigenvectors of T.
52: *> The accuracy varies depending on whether bisection is used to
53: *> find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to
54: *> conpute all and then discard any unwanted one.
55: *> As an added benefit, DLARRE also outputs the n
56: *> Gerschgorin intervals for the matrices L_i D_i L_i^T.
57: *> \endverbatim
58: *
59: * Arguments:
60: * ==========
61: *
62: *> \param[in] RANGE
63: *> \verbatim
64: *> RANGE is CHARACTER*1
65: *> = 'A': ("All") all eigenvalues will be found.
66: *> = 'V': ("Value") all eigenvalues in the half-open interval
67: *> (VL, VU] will be found.
68: *> = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
69: *> entire matrix) will be found.
70: *> \endverbatim
71: *>
72: *> \param[in] N
73: *> \verbatim
74: *> N is INTEGER
75: *> The order of the matrix. N > 0.
76: *> \endverbatim
77: *>
78: *> \param[in,out] VL
79: *> \verbatim
80: *> VL is DOUBLE PRECISION
1.17 bertrand 81: *> If RANGE='V', the lower bound for the eigenvalues.
82: *> Eigenvalues less than or equal to VL, or greater than VU,
83: *> will not be returned. VL < VU.
84: *> If RANGE='I' or ='A', DLARRE computes bounds on the desired
85: *> part of the spectrum.
1.11 bertrand 86: *> \endverbatim
87: *>
88: *> \param[in,out] VU
89: *> \verbatim
90: *> VU is DOUBLE PRECISION
1.17 bertrand 91: *> If RANGE='V', the upper bound for the eigenvalues.
1.11 bertrand 92: *> Eigenvalues less than or equal to VL, or greater than VU,
93: *> will not be returned. VL < VU.
94: *> If RANGE='I' or ='A', DLARRE computes bounds on the desired
95: *> part of the spectrum.
96: *> \endverbatim
97: *>
98: *> \param[in] IL
99: *> \verbatim
100: *> IL is INTEGER
1.17 bertrand 101: *> If RANGE='I', the index of the
102: *> smallest eigenvalue to be returned.
103: *> 1 <= IL <= IU <= N.
1.11 bertrand 104: *> \endverbatim
105: *>
106: *> \param[in] IU
107: *> \verbatim
108: *> IU is INTEGER
1.17 bertrand 109: *> If RANGE='I', the index of the
110: *> largest eigenvalue to be returned.
1.11 bertrand 111: *> 1 <= IL <= IU <= N.
112: *> \endverbatim
113: *>
114: *> \param[in,out] D
115: *> \verbatim
116: *> D is DOUBLE PRECISION array, dimension (N)
117: *> On entry, the N diagonal elements of the tridiagonal
118: *> matrix T.
119: *> On exit, the N diagonal elements of the diagonal
120: *> matrices D_i.
121: *> \endverbatim
122: *>
123: *> \param[in,out] E
124: *> \verbatim
125: *> E is DOUBLE PRECISION array, dimension (N)
126: *> On entry, the first (N-1) entries contain the subdiagonal
127: *> elements of the tridiagonal matrix T; E(N) need not be set.
128: *> On exit, E contains the subdiagonal elements of the unit
129: *> bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
130: *> 1 <= I <= NSPLIT, contain the base points sigma_i on output.
131: *> \endverbatim
132: *>
133: *> \param[in,out] E2
134: *> \verbatim
135: *> E2 is DOUBLE PRECISION array, dimension (N)
136: *> On entry, the first (N-1) entries contain the SQUARES of the
137: *> subdiagonal elements of the tridiagonal matrix T;
138: *> E2(N) need not be set.
139: *> On exit, the entries E2( ISPLIT( I ) ),
140: *> 1 <= I <= NSPLIT, have been set to zero
141: *> \endverbatim
142: *>
143: *> \param[in] RTOL1
144: *> \verbatim
145: *> RTOL1 is DOUBLE PRECISION
146: *> \endverbatim
147: *>
148: *> \param[in] RTOL2
149: *> \verbatim
150: *> RTOL2 is DOUBLE PRECISION
151: *> Parameters for bisection.
152: *> An interval [LEFT,RIGHT] has converged if
1.23 bertrand 153: *> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
1.11 bertrand 154: *> \endverbatim
155: *>
156: *> \param[in] SPLTOL
157: *> \verbatim
158: *> SPLTOL is DOUBLE PRECISION
159: *> The threshold for splitting.
160: *> \endverbatim
161: *>
162: *> \param[out] NSPLIT
163: *> \verbatim
164: *> NSPLIT is INTEGER
165: *> The number of blocks T splits into. 1 <= NSPLIT <= N.
166: *> \endverbatim
167: *>
168: *> \param[out] ISPLIT
169: *> \verbatim
170: *> ISPLIT is INTEGER array, dimension (N)
171: *> The splitting points, at which T breaks up into blocks.
172: *> The first block 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: *> \endverbatim
177: *>
178: *> \param[out] M
179: *> \verbatim
180: *> M is INTEGER
181: *> The total number of eigenvalues (of all L_i D_i L_i^T)
182: *> found.
183: *> \endverbatim
184: *>
185: *> \param[out] W
186: *> \verbatim
187: *> W is DOUBLE PRECISION array, dimension (N)
188: *> The first M elements contain the eigenvalues. The
189: *> eigenvalues of each of the blocks, L_i D_i L_i^T, are
190: *> sorted in ascending order ( DLARRE may use the
191: *> remaining N-M elements as workspace).
192: *> \endverbatim
193: *>
194: *> \param[out] WERR
195: *> \verbatim
196: *> WERR is DOUBLE PRECISION array, dimension (N)
197: *> The error bound on the corresponding eigenvalue in W.
198: *> \endverbatim
199: *>
200: *> \param[out] WGAP
201: *> \verbatim
202: *> WGAP is DOUBLE PRECISION array, dimension (N)
203: *> The separation from the right neighbor eigenvalue in W.
204: *> The gap is only with respect to the eigenvalues of the same block
205: *> as each block has its own representation tree.
206: *> Exception: at the right end of a block we store the left gap
207: *> \endverbatim
208: *>
209: *> \param[out] IBLOCK
210: *> \verbatim
211: *> IBLOCK is INTEGER array, dimension (N)
212: *> The indices of the blocks (submatrices) associated with the
213: *> corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
214: *> W(i) belongs to the first block from the top, =2 if W(i)
215: *> belongs to the second block, etc.
216: *> \endverbatim
217: *>
218: *> \param[out] INDEXW
219: *> \verbatim
220: *> INDEXW is INTEGER array, dimension (N)
221: *> The indices of the eigenvalues within each block (submatrix);
222: *> for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
223: *> i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
224: *> \endverbatim
225: *>
226: *> \param[out] GERS
227: *> \verbatim
228: *> GERS is DOUBLE PRECISION array, dimension (2*N)
229: *> The N Gerschgorin intervals (the i-th Gerschgorin interval
230: *> is (GERS(2*i-1), GERS(2*i)).
231: *> \endverbatim
232: *>
233: *> \param[out] PIVMIN
234: *> \verbatim
235: *> PIVMIN is DOUBLE PRECISION
236: *> The minimum pivot in the Sturm sequence for T.
237: *> \endverbatim
238: *>
239: *> \param[out] WORK
240: *> \verbatim
241: *> WORK is DOUBLE PRECISION array, dimension (6*N)
242: *> Workspace.
243: *> \endverbatim
244: *>
245: *> \param[out] IWORK
246: *> \verbatim
247: *> IWORK is INTEGER array, dimension (5*N)
248: *> Workspace.
249: *> \endverbatim
250: *>
251: *> \param[out] INFO
252: *> \verbatim
253: *> INFO is INTEGER
254: *> = 0: successful exit
1.17 bertrand 255: *> > 0: A problem occurred in DLARRE.
1.11 bertrand 256: *> < 0: One of the called subroutines signaled an internal problem.
257: *> Needs inspection of the corresponding parameter IINFO
258: *> for further information.
259: *>
260: *> =-1: Problem in DLARRD.
261: *> = 2: No base representation could be found in MAXTRY iterations.
262: *> Increasing MAXTRY and recompilation might be a remedy.
263: *> =-3: Problem in DLARRB when computing the refined root
264: *> representation for DLASQ2.
265: *> =-4: Problem in DLARRB when preforming bisection on the
266: *> desired part of the spectrum.
267: *> =-5: Problem in DLASQ2.
268: *> =-6: Problem in DLASQ2.
269: *> \endverbatim
270: *
271: * Authors:
272: * ========
273: *
1.19 bertrand 274: *> \author Univ. of Tennessee
275: *> \author Univ. of California Berkeley
276: *> \author Univ. of Colorado Denver
277: *> \author NAG Ltd.
1.11 bertrand 278: *
1.19 bertrand 279: *> \ingroup OTHERauxiliary
1.11 bertrand 280: *
281: *> \par Further Details:
282: * =====================
283: *>
284: *> \verbatim
285: *>
286: *> The base representations are required to suffer very little
287: *> element growth and consequently define all their eigenvalues to
288: *> high relative accuracy.
289: *> \endverbatim
290: *
291: *> \par Contributors:
292: * ==================
293: *>
294: *> Beresford Parlett, University of California, Berkeley, USA \n
295: *> Jim Demmel, University of California, Berkeley, USA \n
296: *> Inderjit Dhillon, University of Texas, Austin, USA \n
297: *> Osni Marques, LBNL/NERSC, USA \n
298: *> Christof Voemel, University of California, Berkeley, USA \n
299: *>
300: * =====================================================================
1.1 bertrand 301: SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
302: $ RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
303: $ W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
304: $ WORK, IWORK, INFO )
305: *
1.24 ! bertrand 306: * -- LAPACK auxiliary routine --
1.1 bertrand 307: * -- LAPACK is a software package provided by Univ. of Tennessee, --
308: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
309: *
310: * .. Scalar Arguments ..
311: CHARACTER RANGE
312: INTEGER IL, INFO, IU, M, N, NSPLIT
313: DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
314: * ..
315: * .. Array Arguments ..
316: INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
317: $ INDEXW( * )
318: DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
319: $ W( * ),WERR( * ), WGAP( * ), WORK( * )
320: * ..
321: *
322: * =====================================================================
323: *
324: * .. Parameters ..
325: DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
326: $ MAXGROWTH, ONE, PERT, TWO, ZERO
327: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0,
328: $ TWO = 2.0D0, FOUR=4.0D0,
329: $ HNDRD = 100.0D0,
330: $ PERT = 8.0D0,
331: $ HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF,
332: $ MAXGROWTH = 64.0D0, FUDGE = 2.0D0 )
333: INTEGER MAXTRY, ALLRNG, INDRNG, VALRNG
334: PARAMETER ( MAXTRY = 6, ALLRNG = 1, INDRNG = 2,
335: $ VALRNG = 3 )
336: * ..
337: * .. Local Scalars ..
338: LOGICAL FORCEB, NOREP, USEDQD
339: INTEGER CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
340: $ IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
341: $ WBEGIN, WEND
342: DOUBLE PRECISION AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
343: $ EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL,
344: $ RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM,
345: $ TAU, TMP, TMP1
346:
347:
348: * ..
349: * .. Local Arrays ..
350: INTEGER ISEED( 4 )
351: * ..
352: * .. External Functions ..
353: LOGICAL LSAME
354: DOUBLE PRECISION DLAMCH
355: EXTERNAL DLAMCH, LSAME
356:
357: * ..
358: * .. External Subroutines ..
359: EXTERNAL DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD,
1.21 bertrand 360: $ DLASQ2, DLARRK
1.1 bertrand 361: * ..
362: * .. Intrinsic Functions ..
363: INTRINSIC ABS, MAX, MIN
364:
365: * ..
366: * .. Executable Statements ..
367: *
368:
369: INFO = 0
1.21 bertrand 370: *
371: * Quick return if possible
372: *
373: IF( N.LE.0 ) THEN
374: RETURN
375: END IF
1.1 bertrand 376: *
377: * Decode RANGE
378: *
379: IF( LSAME( RANGE, 'A' ) ) THEN
380: IRANGE = ALLRNG
381: ELSE IF( LSAME( RANGE, 'V' ) ) THEN
382: IRANGE = VALRNG
383: ELSE IF( LSAME( RANGE, 'I' ) ) THEN
384: IRANGE = INDRNG
385: END IF
386:
387: M = 0
388:
389: * Get machine constants
390: SAFMIN = DLAMCH( 'S' )
391: EPS = DLAMCH( 'P' )
392:
393: * Set parameters
394: RTL = SQRT(EPS)
395: BSRTOL = SQRT(EPS)
396:
397: * Treat case of 1x1 matrix for quick return
398: IF( N.EQ.1 ) THEN
399: IF( (IRANGE.EQ.ALLRNG).OR.
400: $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
401: $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
402: M = 1
403: W(1) = D(1)
404: * The computation error of the eigenvalue is zero
405: WERR(1) = ZERO
406: WGAP(1) = ZERO
407: IBLOCK( 1 ) = 1
408: INDEXW( 1 ) = 1
409: GERS(1) = D( 1 )
410: GERS(2) = D( 1 )
411: ENDIF
412: * store the shift for the initial RRR, which is zero in this case
413: E(1) = ZERO
414: RETURN
415: END IF
416:
417: * General case: tridiagonal matrix of order > 1
418: *
419: * Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
420: * Compute maximum off-diagonal entry and pivmin.
421: GL = D(1)
422: GU = D(1)
423: EOLD = ZERO
424: EMAX = ZERO
425: E(N) = ZERO
426: DO 5 I = 1,N
427: WERR(I) = ZERO
428: WGAP(I) = ZERO
429: EABS = ABS( E(I) )
430: IF( EABS .GE. EMAX ) THEN
431: EMAX = EABS
432: END IF
433: TMP1 = EABS + EOLD
434: GERS( 2*I-1) = D(I) - TMP1
435: GL = MIN( GL, GERS( 2*I - 1))
436: GERS( 2*I ) = D(I) + TMP1
437: GU = MAX( GU, GERS(2*I) )
438: EOLD = EABS
439: 5 CONTINUE
440: * The minimum pivot allowed in the Sturm sequence for T
441: PIVMIN = SAFMIN * MAX( ONE, EMAX**2 )
442: * Compute spectral diameter. The Gerschgorin bounds give an
443: * estimate that is wrong by at most a factor of SQRT(2)
444: SPDIAM = GU - GL
445:
446: * Compute splitting points
447: CALL DLARRA( N, D, E, E2, SPLTOL, SPDIAM,
448: $ NSPLIT, ISPLIT, IINFO )
449:
450: * Can force use of bisection instead of faster DQDS.
451: * Option left in the code for future multisection work.
452: FORCEB = .FALSE.
453:
454: * Initialize USEDQD, DQDS should be used for ALLRNG unless someone
455: * explicitly wants bisection.
456: USEDQD = (( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB))
457:
458: IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN
459: * Set interval [VL,VU] that contains all eigenvalues
460: VL = GL
461: VU = GU
462: ELSE
463: * We call DLARRD to find crude approximations to the eigenvalues
464: * in the desired range. In case IRANGE = INDRNG, we also obtain the
465: * interval (VL,VU] that contains all the wanted eigenvalues.
466: * An interval [LEFT,RIGHT] has converged if
467: * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
468: * DLARRD needs a WORK of size 4*N, IWORK of size 3*N
469: CALL DLARRD( RANGE, 'B', N, VL, VU, IL, IU, GERS,
470: $ BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
471: $ MM, W, WERR, VL, VU, IBLOCK, INDEXW,
472: $ WORK, IWORK, IINFO )
473: IF( IINFO.NE.0 ) THEN
474: INFO = -1
475: RETURN
476: ENDIF
477: * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
478: DO 14 I = MM+1,N
479: W( I ) = ZERO
480: WERR( I ) = ZERO
481: IBLOCK( I ) = 0
482: INDEXW( I ) = 0
483: 14 CONTINUE
484: END IF
485:
486:
487: ***
488: * Loop over unreduced blocks
489: IBEGIN = 1
490: WBEGIN = 1
491: DO 170 JBLK = 1, NSPLIT
492: IEND = ISPLIT( JBLK )
493: IN = IEND - IBEGIN + 1
494:
495: * 1 X 1 block
496: IF( IN.EQ.1 ) THEN
497: IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND.
498: $ ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) )
499: $ .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK))
500: $ ) THEN
501: M = M + 1
502: W( M ) = D( IBEGIN )
503: WERR(M) = ZERO
504: * The gap for a single block doesn't matter for the later
505: * algorithm and is assigned an arbitrary large value
506: WGAP(M) = ZERO
507: IBLOCK( M ) = JBLK
508: INDEXW( M ) = 1
509: WBEGIN = WBEGIN + 1
510: ENDIF
511: * E( IEND ) holds the shift for the initial RRR
512: E( IEND ) = ZERO
513: IBEGIN = IEND + 1
514: GO TO 170
515: END IF
516: *
517: * Blocks of size larger than 1x1
518: *
519: * E( IEND ) will hold the shift for the initial RRR, for now set it =0
520: E( IEND ) = ZERO
521: *
522: * Find local outer bounds GL,GU for the block
523: GL = D(IBEGIN)
524: GU = D(IBEGIN)
525: DO 15 I = IBEGIN , IEND
526: GL = MIN( GERS( 2*I-1 ), GL )
527: GU = MAX( GERS( 2*I ), GU )
528: 15 CONTINUE
529: SPDIAM = GU - GL
530:
531: IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN
532: * Count the number of eigenvalues in the current block.
533: MB = 0
534: DO 20 I = WBEGIN,MM
535: IF( IBLOCK(I).EQ.JBLK ) THEN
536: MB = MB+1
537: ELSE
538: GOTO 21
539: ENDIF
540: 20 CONTINUE
541: 21 CONTINUE
542:
543: IF( MB.EQ.0) THEN
544: * No eigenvalue in the current block lies in the desired range
545: * E( IEND ) holds the shift for the initial RRR
546: E( IEND ) = ZERO
547: IBEGIN = IEND + 1
548: GO TO 170
549: ELSE
550:
551: * Decide whether dqds or bisection is more efficient
552: USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) )
553: WEND = WBEGIN + MB - 1
554: * Calculate gaps for the current block
555: * In later stages, when representations for individual
556: * eigenvalues are different, we use SIGMA = E( IEND ).
557: SIGMA = ZERO
558: DO 30 I = WBEGIN, WEND - 1
559: WGAP( I ) = MAX( ZERO,
560: $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
561: 30 CONTINUE
562: WGAP( WEND ) = MAX( ZERO,
563: $ VU - SIGMA - (W( WEND )+WERR( WEND )))
564: * Find local index of the first and last desired evalue.
565: INDL = INDEXW(WBEGIN)
566: INDU = INDEXW( WEND )
567: ENDIF
568: ENDIF
569: IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN
570: * Case of DQDS
571: * Find approximations to the extremal eigenvalues of the block
572: CALL DLARRK( IN, 1, GL, GU, D(IBEGIN),
573: $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
574: IF( IINFO.NE.0 ) THEN
575: INFO = -1
576: RETURN
577: ENDIF
578: ISLEFT = MAX(GL, TMP - TMP1
579: $ - HNDRD * EPS* ABS(TMP - TMP1))
580:
581: CALL DLARRK( IN, IN, GL, GU, D(IBEGIN),
582: $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
583: IF( IINFO.NE.0 ) THEN
584: INFO = -1
585: RETURN
586: ENDIF
587: ISRGHT = MIN(GU, TMP + TMP1
588: $ + HNDRD * EPS * ABS(TMP + TMP1))
589: * Improve the estimate of the spectral diameter
590: SPDIAM = ISRGHT - ISLEFT
591: ELSE
592: * Case of bisection
593: * Find approximations to the wanted extremal eigenvalues
594: ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN)
595: $ - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) ))
596: ISRGHT = MIN(GU,W(WEND) + WERR(WEND)
597: $ + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND)))
598: ENDIF
599:
600:
601: * Decide whether the base representation for the current block
602: * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
603: * should be on the left or the right end of the current block.
604: * The strategy is to shift to the end which is "more populated"
605: * Furthermore, decide whether to use DQDS for the computation of
606: * the eigenvalue approximations at the end of DLARRE or bisection.
607: * dqds is chosen if all eigenvalues are desired or the number of
608: * eigenvalues to be computed is large compared to the blocksize.
609: IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
610: * If all the eigenvalues have to be computed, we use dqd
611: USEDQD = .TRUE.
612: * INDL is the local index of the first eigenvalue to compute
613: INDL = 1
614: INDU = IN
615: * MB = number of eigenvalues to compute
616: MB = IN
617: WEND = WBEGIN + MB - 1
618: * Define 1/4 and 3/4 points of the spectrum
619: S1 = ISLEFT + FOURTH * SPDIAM
620: S2 = ISRGHT - FOURTH * SPDIAM
621: ELSE
622: * DLARRD has computed IBLOCK and INDEXW for each eigenvalue
623: * approximation.
624: * choose sigma
625: IF( USEDQD ) THEN
626: S1 = ISLEFT + FOURTH * SPDIAM
627: S2 = ISRGHT - FOURTH * SPDIAM
628: ELSE
629: TMP = MIN(ISRGHT,VU) - MAX(ISLEFT,VL)
630: S1 = MAX(ISLEFT,VL) + FOURTH * TMP
631: S2 = MIN(ISRGHT,VU) - FOURTH * TMP
632: ENDIF
633: ENDIF
634:
635: * Compute the negcount at the 1/4 and 3/4 points
636: IF(MB.GT.1) THEN
637: CALL DLARRC( 'T', IN, S1, S2, D(IBEGIN),
638: $ E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO)
639: ENDIF
640:
641: IF(MB.EQ.1) THEN
642: SIGMA = GL
643: SGNDEF = ONE
644: ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN
645: IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
646: SIGMA = MAX(ISLEFT,GL)
647: ELSEIF( USEDQD ) THEN
648: * use Gerschgorin bound as shift to get pos def matrix
649: * for dqds
650: SIGMA = ISLEFT
651: ELSE
652: * use approximation of the first desired eigenvalue of the
653: * block as shift
654: SIGMA = MAX(ISLEFT,VL)
655: ENDIF
656: SGNDEF = ONE
657: ELSE
658: IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
659: SIGMA = MIN(ISRGHT,GU)
660: ELSEIF( USEDQD ) THEN
661: * use Gerschgorin bound as shift to get neg def matrix
662: * for dqds
663: SIGMA = ISRGHT
664: ELSE
665: * use approximation of the first desired eigenvalue of the
666: * block as shift
667: SIGMA = MIN(ISRGHT,VU)
668: ENDIF
669: SGNDEF = -ONE
670: ENDIF
671:
672:
673: * An initial SIGMA has been chosen that will be used for computing
674: * T - SIGMA I = L D L^T
675: * Define the increment TAU of the shift in case the initial shift
676: * needs to be refined to obtain a factorization with not too much
677: * element growth.
678: IF( USEDQD ) THEN
679: * The initial SIGMA was to the outer end of the spectrum
680: * the matrix is definite and we need not retreat.
681: TAU = SPDIAM*EPS*N + TWO*PIVMIN
1.10 bertrand 682: TAU = MAX( TAU,TWO*EPS*ABS(SIGMA) )
1.1 bertrand 683: ELSE
684: IF(MB.GT.1) THEN
685: CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
686: AVGAP = ABS(CLWDTH / DBLE(WEND-WBEGIN))
687: IF( SGNDEF.EQ.ONE ) THEN
688: TAU = HALF*MAX(WGAP(WBEGIN),AVGAP)
689: TAU = MAX(TAU,WERR(WBEGIN))
690: ELSE
691: TAU = HALF*MAX(WGAP(WEND-1),AVGAP)
692: TAU = MAX(TAU,WERR(WEND))
693: ENDIF
694: ELSE
695: TAU = WERR(WBEGIN)
696: ENDIF
697: ENDIF
698: *
699: DO 80 IDUM = 1, MAXTRY
700: * Compute L D L^T factorization of tridiagonal matrix T - sigma I.
701: * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
702: * pivots in WORK(2*IN+1:3*IN)
703: DPIVOT = D( IBEGIN ) - SIGMA
704: WORK( 1 ) = DPIVOT
705: DMAX = ABS( WORK(1) )
706: J = IBEGIN
707: DO 70 I = 1, IN - 1
708: WORK( 2*IN+I ) = ONE / WORK( I )
709: TMP = E( J )*WORK( 2*IN+I )
710: WORK( IN+I ) = TMP
711: DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J )
712: WORK( I+1 ) = DPIVOT
713: DMAX = MAX( DMAX, ABS(DPIVOT) )
714: J = J + 1
715: 70 CONTINUE
716: * check for element growth
717: IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN
718: NOREP = .TRUE.
719: ELSE
720: NOREP = .FALSE.
721: ENDIF
722: IF( USEDQD .AND. .NOT.NOREP ) THEN
723: * Ensure the definiteness of the representation
724: * All entries of D (of L D L^T) must have the same sign
725: DO 71 I = 1, IN
726: TMP = SGNDEF*WORK( I )
727: IF( TMP.LT.ZERO ) NOREP = .TRUE.
728: 71 CONTINUE
729: ENDIF
730: IF(NOREP) THEN
731: * Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
732: * shift which makes the matrix definite. So we should end up
733: * here really only in the case of IRANGE = VALRNG or INDRNG.
734: IF( IDUM.EQ.MAXTRY-1 ) THEN
735: IF( SGNDEF.EQ.ONE ) THEN
736: * The fudged Gerschgorin shift should succeed
737: SIGMA =
738: $ GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN
739: ELSE
740: SIGMA =
741: $ GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN
742: END IF
743: ELSE
744: SIGMA = SIGMA - SGNDEF * TAU
745: TAU = TWO * TAU
746: END IF
747: ELSE
748: * an initial RRR is found
749: GO TO 83
750: END IF
751: 80 CONTINUE
752: * if the program reaches this point, no base representation could be
753: * found in MAXTRY iterations.
754: INFO = 2
755: RETURN
756:
757: 83 CONTINUE
758: * At this point, we have found an initial base representation
759: * T - SIGMA I = L D L^T with not too much element growth.
760: * Store the shift.
761: E( IEND ) = SIGMA
762: * Store D and L.
763: CALL DCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
764: CALL DCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )
765:
766:
767: IF(MB.GT.1 ) THEN
768: *
769: * Perturb each entry of the base representation by a small
770: * (but random) relative amount to overcome difficulties with
771: * glued matrices.
772: *
773: DO 122 I = 1, 4
774: ISEED( I ) = 1
775: 122 CONTINUE
776:
777: CALL DLARNV(2, ISEED, 2*IN-1, WORK(1))
778: DO 125 I = 1,IN-1
779: D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
780: E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
781: 125 CONTINUE
782: D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN))
783: *
784: ENDIF
785: *
786: * Don't update the Gerschgorin intervals because keeping track
787: * of the updates would be too much work in DLARRV.
788: * We update W instead and use it to locate the proper Gerschgorin
789: * intervals.
790:
791: * Compute the required eigenvalues of L D L' by bisection or dqds
792: IF ( .NOT.USEDQD ) THEN
793: * If DLARRD has been used, shift the eigenvalue approximations
794: * according to their representation. This is necessary for
795: * a uniform DLARRV since dqds computes eigenvalues of the
796: * shifted representation. In DLARRV, W will always hold the
797: * UNshifted eigenvalue approximation.
798: DO 134 J=WBEGIN,WEND
799: W(J) = W(J) - SIGMA
800: WERR(J) = WERR(J) + ABS(W(J)) * EPS
801: 134 CONTINUE
802: * call DLARRB to reduce eigenvalue error of the approximations
803: * from DLARRD
804: DO 135 I = IBEGIN, IEND-1
805: WORK( I ) = D( I ) * E( I )**2
806: 135 CONTINUE
807: * use bisection to find EV from INDL to INDU
808: CALL DLARRB(IN, D(IBEGIN), WORK(IBEGIN),
809: $ INDL, INDU, RTOL1, RTOL2, INDL-1,
810: $ W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN),
811: $ WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM,
812: $ IN, IINFO )
813: IF( IINFO .NE. 0 ) THEN
814: INFO = -4
815: RETURN
816: END IF
817: * DLARRB computes all gaps correctly except for the last one
818: * Record distance to VU/GU
819: WGAP( WEND ) = MAX( ZERO,
820: $ ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) )
821: DO 138 I = INDL, INDU
822: M = M + 1
823: IBLOCK(M) = JBLK
824: INDEXW(M) = I
825: 138 CONTINUE
826: ELSE
827: * Call dqds to get all eigs (and then possibly delete unwanted
828: * eigenvalues).
829: * Note that dqds finds the eigenvalues of the L D L^T representation
830: * of T to high relative accuracy. High relative accuracy
831: * might be lost when the shift of the RRR is subtracted to obtain
832: * the eigenvalues of T. However, T is not guaranteed to define its
833: * eigenvalues to high relative accuracy anyway.
834: * Set RTOL to the order of the tolerance used in DLASQ2
835: * This is an ESTIMATED error, the worst case bound is 4*N*EPS
836: * which is usually too large and requires unnecessary work to be
837: * done by bisection when computing the eigenvectors
838: RTOL = LOG(DBLE(IN)) * FOUR * EPS
839: J = IBEGIN
840: DO 140 I = 1, IN - 1
841: WORK( 2*I-1 ) = ABS( D( J ) )
842: WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
843: J = J + 1
844: 140 CONTINUE
845: WORK( 2*IN-1 ) = ABS( D( IEND ) )
846: WORK( 2*IN ) = ZERO
847: CALL DLASQ2( IN, WORK, IINFO )
848: IF( IINFO .NE. 0 ) THEN
849: * If IINFO = -5 then an index is part of a tight cluster
850: * and should be changed. The index is in IWORK(1) and the
851: * gap is in WORK(N+1)
852: INFO = -5
853: RETURN
854: ELSE
855: * Test that all eigenvalues are positive as expected
856: DO 149 I = 1, IN
857: IF( WORK( I ).LT.ZERO ) THEN
858: INFO = -6
859: RETURN
860: ENDIF
861: 149 CONTINUE
862: END IF
863: IF( SGNDEF.GT.ZERO ) THEN
864: DO 150 I = INDL, INDU
865: M = M + 1
866: W( M ) = WORK( IN-I+1 )
867: IBLOCK( M ) = JBLK
868: INDEXW( M ) = I
869: 150 CONTINUE
870: ELSE
871: DO 160 I = INDL, INDU
872: M = M + 1
873: W( M ) = -WORK( I )
874: IBLOCK( M ) = JBLK
875: INDEXW( M ) = I
876: 160 CONTINUE
877: END IF
878:
879: DO 165 I = M - MB + 1, M
880: * the value of RTOL below should be the tolerance in DLASQ2
881: WERR( I ) = RTOL * ABS( W(I) )
882: 165 CONTINUE
883: DO 166 I = M - MB + 1, M - 1
884: * compute the right gap between the intervals
885: WGAP( I ) = MAX( ZERO,
886: $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
887: 166 CONTINUE
888: WGAP( M ) = MAX( ZERO,
889: $ ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) )
890: END IF
891: * proceed with next block
892: IBEGIN = IEND + 1
893: WBEGIN = WEND + 1
894: 170 CONTINUE
895: *
896:
897: RETURN
898: *
1.24 ! bertrand 899: * End of DLARRE
1.1 bertrand 900: *
901: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarre.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack