Annotation of rpl/lapack/lapack/dlarrb.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b DLARRB
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLARRB + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrb.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrb.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrb.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
! 22: * RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
! 23: * PIVMIN, SPDIAM, TWIST, INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
! 27: * DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPDIAM
! 28: * ..
! 29: * .. Array Arguments ..
! 30: * INTEGER IWORK( * )
! 31: * DOUBLE PRECISION D( * ), LLD( * ), W( * ),
! 32: * $ WERR( * ), WGAP( * ), WORK( * )
! 33: * ..
! 34: *
! 35: *
! 36: *> \par Purpose:
! 37: * =============
! 38: *>
! 39: *> \verbatim
! 40: *>
! 41: *> Given the relatively robust representation(RRR) L D L^T, DLARRB
! 42: *> does "limited" bisection to refine the eigenvalues of L D L^T,
! 43: *> W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
! 44: *> guesses for these eigenvalues are input in W, the corresponding estimate
! 45: *> of the error in these guesses and their gaps are input in WERR
! 46: *> and WGAP, respectively. During bisection, intervals
! 47: *> [left, right] are maintained by storing their mid-points and
! 48: *> semi-widths in the arrays W and WERR respectively.
! 49: *> \endverbatim
! 50: *
! 51: * Arguments:
! 52: * ==========
! 53: *
! 54: *> \param[in] N
! 55: *> \verbatim
! 56: *> N is INTEGER
! 57: *> The order of the matrix.
! 58: *> \endverbatim
! 59: *>
! 60: *> \param[in] D
! 61: *> \verbatim
! 62: *> D is DOUBLE PRECISION array, dimension (N)
! 63: *> The N diagonal elements of the diagonal matrix D.
! 64: *> \endverbatim
! 65: *>
! 66: *> \param[in] LLD
! 67: *> \verbatim
! 68: *> LLD is DOUBLE PRECISION array, dimension (N-1)
! 69: *> The (N-1) elements L(i)*L(i)*D(i).
! 70: *> \endverbatim
! 71: *>
! 72: *> \param[in] IFIRST
! 73: *> \verbatim
! 74: *> IFIRST is INTEGER
! 75: *> The index of the first eigenvalue to be computed.
! 76: *> \endverbatim
! 77: *>
! 78: *> \param[in] ILAST
! 79: *> \verbatim
! 80: *> ILAST is INTEGER
! 81: *> The index of the last eigenvalue to be computed.
! 82: *> \endverbatim
! 83: *>
! 84: *> \param[in] RTOL1
! 85: *> \verbatim
! 86: *> RTOL1 is DOUBLE PRECISION
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in] RTOL2
! 90: *> \verbatim
! 91: *> RTOL2 is DOUBLE PRECISION
! 92: *> Tolerance for the convergence of the bisection intervals.
! 93: *> An interval [LEFT,RIGHT] has converged if
! 94: *> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
! 95: *> where GAP is the (estimated) distance to the nearest
! 96: *> eigenvalue.
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] OFFSET
! 100: *> \verbatim
! 101: *> OFFSET is INTEGER
! 102: *> Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
! 103: *> through ILAST-OFFSET elements of these arrays are to be used.
! 104: *> \endverbatim
! 105: *>
! 106: *> \param[in,out] W
! 107: *> \verbatim
! 108: *> W is DOUBLE PRECISION array, dimension (N)
! 109: *> On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
! 110: *> estimates of the eigenvalues of L D L^T indexed IFIRST throug
! 111: *> ILAST.
! 112: *> On output, these estimates are refined.
! 113: *> \endverbatim
! 114: *>
! 115: *> \param[in,out] WGAP
! 116: *> \verbatim
! 117: *> WGAP is DOUBLE PRECISION array, dimension (N-1)
! 118: *> On input, the (estimated) gaps between consecutive
! 119: *> eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
! 120: *> eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST
! 121: *> then WGAP(IFIRST-OFFSET) must be set to ZERO.
! 122: *> On output, these gaps are refined.
! 123: *> \endverbatim
! 124: *>
! 125: *> \param[in,out] WERR
! 126: *> \verbatim
! 127: *> WERR is DOUBLE PRECISION array, dimension (N)
! 128: *> On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
! 129: *> the errors in the estimates of the corresponding elements in W.
! 130: *> On output, these errors are refined.
! 131: *> \endverbatim
! 132: *>
! 133: *> \param[out] WORK
! 134: *> \verbatim
! 135: *> WORK is DOUBLE PRECISION array, dimension (2*N)
! 136: *> Workspace.
! 137: *> \endverbatim
! 138: *>
! 139: *> \param[out] IWORK
! 140: *> \verbatim
! 141: *> IWORK is INTEGER array, dimension (2*N)
! 142: *> Workspace.
! 143: *> \endverbatim
! 144: *>
! 145: *> \param[in] PIVMIN
! 146: *> \verbatim
! 147: *> PIVMIN is DOUBLE PRECISION
! 148: *> The minimum pivot in the Sturm sequence.
! 149: *> \endverbatim
! 150: *>
! 151: *> \param[in] SPDIAM
! 152: *> \verbatim
! 153: *> SPDIAM is DOUBLE PRECISION
! 154: *> The spectral diameter of the matrix.
! 155: *> \endverbatim
! 156: *>
! 157: *> \param[in] TWIST
! 158: *> \verbatim
! 159: *> TWIST is INTEGER
! 160: *> The twist index for the twisted factorization that is used
! 161: *> for the negcount.
! 162: *> TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
! 163: *> TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
! 164: *> TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
! 165: *> \endverbatim
! 166: *>
! 167: *> \param[out] INFO
! 168: *> \verbatim
! 169: *> INFO is INTEGER
! 170: *> Error flag.
! 171: *> \endverbatim
! 172: *
! 173: * Authors:
! 174: * ========
! 175: *
! 176: *> \author Univ. of Tennessee
! 177: *> \author Univ. of California Berkeley
! 178: *> \author Univ. of Colorado Denver
! 179: *> \author NAG Ltd.
! 180: *
! 181: *> \date November 2011
! 182: *
! 183: *> \ingroup auxOTHERauxiliary
! 184: *
! 185: *> \par Contributors:
! 186: * ==================
! 187: *>
! 188: *> Beresford Parlett, University of California, Berkeley, USA \n
! 189: *> Jim Demmel, University of California, Berkeley, USA \n
! 190: *> Inderjit Dhillon, University of Texas, Austin, USA \n
! 191: *> Osni Marques, LBNL/NERSC, USA \n
! 192: *> Christof Voemel, University of California, Berkeley, USA
! 193: *
! 194: * =====================================================================
1.1 bertrand 195: SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
196: $ RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
197: $ PIVMIN, SPDIAM, TWIST, INFO )
198: *
1.8 ! bertrand 199: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 200: * -- LAPACK is a software package provided by Univ. of Tennessee, --
201: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 202: * November 2011
1.1 bertrand 203: *
204: * .. Scalar Arguments ..
205: INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
206: DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPDIAM
207: * ..
208: * .. Array Arguments ..
209: INTEGER IWORK( * )
210: DOUBLE PRECISION D( * ), LLD( * ), W( * ),
211: $ WERR( * ), WGAP( * ), WORK( * )
212: * ..
213: *
214: * =====================================================================
215: *
216: * .. Parameters ..
217: DOUBLE PRECISION ZERO, TWO, HALF
218: PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0,
219: $ HALF = 0.5D0 )
220: INTEGER MAXITR
221: * ..
222: * .. Local Scalars ..
223: INTEGER I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
224: $ OLNINT, PREV, R
225: DOUBLE PRECISION BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
226: $ RGAP, RIGHT, TMP, WIDTH
227: * ..
228: * .. External Functions ..
229: INTEGER DLANEG
230: EXTERNAL DLANEG
231: *
232: * ..
233: * .. Intrinsic Functions ..
234: INTRINSIC ABS, MAX, MIN
235: * ..
236: * .. Executable Statements ..
237: *
238: INFO = 0
239: *
240: MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
241: $ LOG( TWO ) ) + 2
242: MNWDTH = TWO * PIVMIN
243: *
244: R = TWIST
245: IF((R.LT.1).OR.(R.GT.N)) R = N
246: *
247: * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
248: * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
249: * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
250: * for an unconverged interval is set to the index of the next unconverged
251: * interval, and is -1 or 0 for a converged interval. Thus a linked
252: * list of unconverged intervals is set up.
253: *
254: I1 = IFIRST
255: * The number of unconverged intervals
256: NINT = 0
257: * The last unconverged interval found
258: PREV = 0
259:
260: RGAP = WGAP( I1-OFFSET )
261: DO 75 I = I1, ILAST
262: K = 2*I
263: II = I - OFFSET
264: LEFT = W( II ) - WERR( II )
265: RIGHT = W( II ) + WERR( II )
266: LGAP = RGAP
267: RGAP = WGAP( II )
268: GAP = MIN( LGAP, RGAP )
269:
270: * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
271: * Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
272: *
273: * Do while( NEGCNT(LEFT).GT.I-1 )
274: *
275: BACK = WERR( II )
276: 20 CONTINUE
277: NEGCNT = DLANEG( N, D, LLD, LEFT, PIVMIN, R )
278: IF( NEGCNT.GT.I-1 ) THEN
279: LEFT = LEFT - BACK
280: BACK = TWO*BACK
281: GO TO 20
282: END IF
283: *
284: * Do while( NEGCNT(RIGHT).LT.I )
285: * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
286: *
287: BACK = WERR( II )
288: 50 CONTINUE
289:
290: NEGCNT = DLANEG( N, D, LLD, RIGHT, PIVMIN, R )
291: IF( NEGCNT.LT.I ) THEN
292: RIGHT = RIGHT + BACK
293: BACK = TWO*BACK
294: GO TO 50
295: END IF
296: WIDTH = HALF*ABS( LEFT - RIGHT )
297: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
298: CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
299: IF( WIDTH.LE.CVRGD .OR. WIDTH.LE.MNWDTH ) THEN
300: * This interval has already converged and does not need refinement.
301: * (Note that the gaps might change through refining the
302: * eigenvalues, however, they can only get bigger.)
303: * Remove it from the list.
304: IWORK( K-1 ) = -1
305: * Make sure that I1 always points to the first unconverged interval
306: IF((I.EQ.I1).AND.(I.LT.ILAST)) I1 = I + 1
307: IF((PREV.GE.I1).AND.(I.LE.ILAST)) IWORK( 2*PREV-1 ) = I + 1
308: ELSE
309: * unconverged interval found
310: PREV = I
311: NINT = NINT + 1
312: IWORK( K-1 ) = I + 1
313: IWORK( K ) = NEGCNT
314: END IF
315: WORK( K-1 ) = LEFT
316: WORK( K ) = RIGHT
317: 75 CONTINUE
318:
319: *
320: * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
321: * and while (ITER.LT.MAXITR)
322: *
323: ITER = 0
324: 80 CONTINUE
325: PREV = I1 - 1
326: I = I1
327: OLNINT = NINT
328:
329: DO 100 IP = 1, OLNINT
330: K = 2*I
331: II = I - OFFSET
332: RGAP = WGAP( II )
333: LGAP = RGAP
334: IF(II.GT.1) LGAP = WGAP( II-1 )
335: GAP = MIN( LGAP, RGAP )
336: NEXT = IWORK( K-1 )
337: LEFT = WORK( K-1 )
338: RIGHT = WORK( K )
339: MID = HALF*( LEFT + RIGHT )
340:
341: * semiwidth of interval
342: WIDTH = RIGHT - MID
343: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
344: CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
345: IF( ( WIDTH.LE.CVRGD ) .OR. ( WIDTH.LE.MNWDTH ).OR.
346: $ ( ITER.EQ.MAXITR ) )THEN
347: * reduce number of unconverged intervals
348: NINT = NINT - 1
349: * Mark interval as converged.
350: IWORK( K-1 ) = 0
351: IF( I1.EQ.I ) THEN
352: I1 = NEXT
353: ELSE
354: * Prev holds the last unconverged interval previously examined
355: IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
356: END IF
357: I = NEXT
358: GO TO 100
359: END IF
360: PREV = I
361: *
362: * Perform one bisection step
363: *
364: NEGCNT = DLANEG( N, D, LLD, MID, PIVMIN, R )
365: IF( NEGCNT.LE.I-1 ) THEN
366: WORK( K-1 ) = MID
367: ELSE
368: WORK( K ) = MID
369: END IF
370: I = NEXT
371: 100 CONTINUE
372: ITER = ITER + 1
373: * do another loop if there are still unconverged intervals
374: * However, in the last iteration, all intervals are accepted
375: * since this is the best we can do.
376: IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
377: *
378: *
379: * At this point, all the intervals have converged
380: DO 110 I = IFIRST, ILAST
381: K = 2*I
382: II = I - OFFSET
383: * All intervals marked by '0' have been refined.
384: IF( IWORK( K-1 ).EQ.0 ) THEN
385: W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
386: WERR( II ) = WORK( K ) - W( II )
387: END IF
388: 110 CONTINUE
389: *
390: DO 111 I = IFIRST+1, ILAST
391: K = 2*I
392: II = I - OFFSET
393: WGAP( II-1 ) = MAX( ZERO,
394: $ W(II) - WERR (II) - W( II-1 ) - WERR( II-1 ))
395: 111 CONTINUE
396:
397: RETURN
398: *
399: * End of DLARRB
400: *
401: 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