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