Annotation of rpl/lapack/lapack/dlarrb.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
! 2: $ RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
! 3: $ PIVMIN, SPDIAM, TWIST, INFO )
! 4: *
! 5: * -- LAPACK auxiliary routine (version 3.2) --
! 6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 8: * November 2006
! 9: *
! 10: * .. Scalar Arguments ..
! 11: INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
! 12: DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPDIAM
! 13: * ..
! 14: * .. Array Arguments ..
! 15: INTEGER IWORK( * )
! 16: DOUBLE PRECISION D( * ), LLD( * ), W( * ),
! 17: $ WERR( * ), WGAP( * ), WORK( * )
! 18: * ..
! 19: *
! 20: * Purpose
! 21: * =======
! 22: *
! 23: * Given the relatively robust representation(RRR) L D L^T, DLARRB
! 24: * does "limited" bisection to refine the eigenvalues of L D L^T,
! 25: * W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
! 26: * guesses for these eigenvalues are input in W, the corresponding estimate
! 27: * of the error in these guesses and their gaps are input in WERR
! 28: * and WGAP, respectively. During bisection, intervals
! 29: * [left, right] are maintained by storing their mid-points and
! 30: * semi-widths in the arrays W and WERR respectively.
! 31: *
! 32: * Arguments
! 33: * =========
! 34: *
! 35: * N (input) INTEGER
! 36: * The order of the matrix.
! 37: *
! 38: * D (input) DOUBLE PRECISION array, dimension (N)
! 39: * The N diagonal elements of the diagonal matrix D.
! 40: *
! 41: * LLD (input) DOUBLE PRECISION array, dimension (N-1)
! 42: * The (N-1) elements L(i)*L(i)*D(i).
! 43: *
! 44: * IFIRST (input) INTEGER
! 45: * The index of the first eigenvalue to be computed.
! 46: *
! 47: * ILAST (input) INTEGER
! 48: * The index of the last eigenvalue to be computed.
! 49: *
! 50: * RTOL1 (input) DOUBLE PRECISION
! 51: * RTOL2 (input) DOUBLE PRECISION
! 52: * Tolerance for the convergence of the bisection intervals.
! 53: * An interval [LEFT,RIGHT] has converged if
! 54: * RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
! 55: * where GAP is the (estimated) distance to the nearest
! 56: * eigenvalue.
! 57: *
! 58: * OFFSET (input) INTEGER
! 59: * Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
! 60: * through ILAST-OFFSET elements of these arrays are to be used.
! 61: *
! 62: * W (input/output) DOUBLE PRECISION array, dimension (N)
! 63: * On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
! 64: * estimates of the eigenvalues of L D L^T indexed IFIRST throug
! 65: * ILAST.
! 66: * On output, these estimates are refined.
! 67: *
! 68: * WGAP (input/output) DOUBLE PRECISION array, dimension (N-1)
! 69: * On input, the (estimated) gaps between consecutive
! 70: * eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
! 71: * eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST
! 72: * then WGAP(IFIRST-OFFSET) must be set to ZERO.
! 73: * On output, these gaps are refined.
! 74: *
! 75: * WERR (input/output) DOUBLE PRECISION array, dimension (N)
! 76: * On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
! 77: * the errors in the estimates of the corresponding elements in W.
! 78: * On output, these errors are refined.
! 79: *
! 80: * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
! 81: * Workspace.
! 82: *
! 83: * IWORK (workspace) INTEGER array, dimension (2*N)
! 84: * Workspace.
! 85: *
! 86: * PIVMIN (input) DOUBLE PRECISION
! 87: * The minimum pivot in the Sturm sequence.
! 88: *
! 89: * SPDIAM (input) DOUBLE PRECISION
! 90: * The spectral diameter of the matrix.
! 91: *
! 92: * TWIST (input) INTEGER
! 93: * The twist index for the twisted factorization that is used
! 94: * for the negcount.
! 95: * TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
! 96: * TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
! 97: * TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
! 98: *
! 99: * INFO (output) INTEGER
! 100: * Error flag.
! 101: *
! 102: * Further Details
! 103: * ===============
! 104: *
! 105: * Based on contributions by
! 106: * Beresford Parlett, University of California, Berkeley, USA
! 107: * Jim Demmel, University of California, Berkeley, USA
! 108: * Inderjit Dhillon, University of Texas, Austin, USA
! 109: * Osni Marques, LBNL/NERSC, USA
! 110: * Christof Voemel, University of California, Berkeley, USA
! 111: *
! 112: * =====================================================================
! 113: *
! 114: * .. Parameters ..
! 115: DOUBLE PRECISION ZERO, TWO, HALF
! 116: PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0,
! 117: $ HALF = 0.5D0 )
! 118: INTEGER MAXITR
! 119: * ..
! 120: * .. Local Scalars ..
! 121: INTEGER I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
! 122: $ OLNINT, PREV, R
! 123: DOUBLE PRECISION BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
! 124: $ RGAP, RIGHT, TMP, WIDTH
! 125: * ..
! 126: * .. External Functions ..
! 127: INTEGER DLANEG
! 128: EXTERNAL DLANEG
! 129: *
! 130: * ..
! 131: * .. Intrinsic Functions ..
! 132: INTRINSIC ABS, MAX, MIN
! 133: * ..
! 134: * .. Executable Statements ..
! 135: *
! 136: INFO = 0
! 137: *
! 138: MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
! 139: $ LOG( TWO ) ) + 2
! 140: MNWDTH = TWO * PIVMIN
! 141: *
! 142: R = TWIST
! 143: IF((R.LT.1).OR.(R.GT.N)) R = N
! 144: *
! 145: * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
! 146: * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
! 147: * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
! 148: * for an unconverged interval is set to the index of the next unconverged
! 149: * interval, and is -1 or 0 for a converged interval. Thus a linked
! 150: * list of unconverged intervals is set up.
! 151: *
! 152: I1 = IFIRST
! 153: * The number of unconverged intervals
! 154: NINT = 0
! 155: * The last unconverged interval found
! 156: PREV = 0
! 157:
! 158: RGAP = WGAP( I1-OFFSET )
! 159: DO 75 I = I1, ILAST
! 160: K = 2*I
! 161: II = I - OFFSET
! 162: LEFT = W( II ) - WERR( II )
! 163: RIGHT = W( II ) + WERR( II )
! 164: LGAP = RGAP
! 165: RGAP = WGAP( II )
! 166: GAP = MIN( LGAP, RGAP )
! 167:
! 168: * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
! 169: * Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
! 170: *
! 171: * Do while( NEGCNT(LEFT).GT.I-1 )
! 172: *
! 173: BACK = WERR( II )
! 174: 20 CONTINUE
! 175: NEGCNT = DLANEG( N, D, LLD, LEFT, PIVMIN, R )
! 176: IF( NEGCNT.GT.I-1 ) THEN
! 177: LEFT = LEFT - BACK
! 178: BACK = TWO*BACK
! 179: GO TO 20
! 180: END IF
! 181: *
! 182: * Do while( NEGCNT(RIGHT).LT.I )
! 183: * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
! 184: *
! 185: BACK = WERR( II )
! 186: 50 CONTINUE
! 187:
! 188: NEGCNT = DLANEG( N, D, LLD, RIGHT, PIVMIN, R )
! 189: IF( NEGCNT.LT.I ) THEN
! 190: RIGHT = RIGHT + BACK
! 191: BACK = TWO*BACK
! 192: GO TO 50
! 193: END IF
! 194: WIDTH = HALF*ABS( LEFT - RIGHT )
! 195: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
! 196: CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
! 197: IF( WIDTH.LE.CVRGD .OR. WIDTH.LE.MNWDTH ) THEN
! 198: * This interval has already converged and does not need refinement.
! 199: * (Note that the gaps might change through refining the
! 200: * eigenvalues, however, they can only get bigger.)
! 201: * Remove it from the list.
! 202: IWORK( K-1 ) = -1
! 203: * Make sure that I1 always points to the first unconverged interval
! 204: IF((I.EQ.I1).AND.(I.LT.ILAST)) I1 = I + 1
! 205: IF((PREV.GE.I1).AND.(I.LE.ILAST)) IWORK( 2*PREV-1 ) = I + 1
! 206: ELSE
! 207: * unconverged interval found
! 208: PREV = I
! 209: NINT = NINT + 1
! 210: IWORK( K-1 ) = I + 1
! 211: IWORK( K ) = NEGCNT
! 212: END IF
! 213: WORK( K-1 ) = LEFT
! 214: WORK( K ) = RIGHT
! 215: 75 CONTINUE
! 216:
! 217: *
! 218: * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
! 219: * and while (ITER.LT.MAXITR)
! 220: *
! 221: ITER = 0
! 222: 80 CONTINUE
! 223: PREV = I1 - 1
! 224: I = I1
! 225: OLNINT = NINT
! 226:
! 227: DO 100 IP = 1, OLNINT
! 228: K = 2*I
! 229: II = I - OFFSET
! 230: RGAP = WGAP( II )
! 231: LGAP = RGAP
! 232: IF(II.GT.1) LGAP = WGAP( II-1 )
! 233: GAP = MIN( LGAP, RGAP )
! 234: NEXT = IWORK( K-1 )
! 235: LEFT = WORK( K-1 )
! 236: RIGHT = WORK( K )
! 237: MID = HALF*( LEFT + RIGHT )
! 238:
! 239: * semiwidth of interval
! 240: WIDTH = RIGHT - MID
! 241: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
! 242: CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
! 243: IF( ( WIDTH.LE.CVRGD ) .OR. ( WIDTH.LE.MNWDTH ).OR.
! 244: $ ( ITER.EQ.MAXITR ) )THEN
! 245: * reduce number of unconverged intervals
! 246: NINT = NINT - 1
! 247: * Mark interval as converged.
! 248: IWORK( K-1 ) = 0
! 249: IF( I1.EQ.I ) THEN
! 250: I1 = NEXT
! 251: ELSE
! 252: * Prev holds the last unconverged interval previously examined
! 253: IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
! 254: END IF
! 255: I = NEXT
! 256: GO TO 100
! 257: END IF
! 258: PREV = I
! 259: *
! 260: * Perform one bisection step
! 261: *
! 262: NEGCNT = DLANEG( N, D, LLD, MID, PIVMIN, R )
! 263: IF( NEGCNT.LE.I-1 ) THEN
! 264: WORK( K-1 ) = MID
! 265: ELSE
! 266: WORK( K ) = MID
! 267: END IF
! 268: I = NEXT
! 269: 100 CONTINUE
! 270: ITER = ITER + 1
! 271: * do another loop if there are still unconverged intervals
! 272: * However, in the last iteration, all intervals are accepted
! 273: * since this is the best we can do.
! 274: IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
! 275: *
! 276: *
! 277: * At this point, all the intervals have converged
! 278: DO 110 I = IFIRST, ILAST
! 279: K = 2*I
! 280: II = I - OFFSET
! 281: * All intervals marked by '0' have been refined.
! 282: IF( IWORK( K-1 ).EQ.0 ) THEN
! 283: W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
! 284: WERR( II ) = WORK( K ) - W( II )
! 285: END IF
! 286: 110 CONTINUE
! 287: *
! 288: DO 111 I = IFIRST+1, ILAST
! 289: K = 2*I
! 290: II = I - OFFSET
! 291: WGAP( II-1 ) = MAX( ZERO,
! 292: $ W(II) - WERR (II) - W( II-1 ) - WERR( II-1 ))
! 293: 111 CONTINUE
! 294:
! 295: RETURN
! 296: *
! 297: * End of DLARRB
! 298: *
! 299: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarrb.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack