Annotation of rpl/lapack/lapack/dlarrb.f, revision 1.20
1.11 bertrand 1: *> \brief \b DLARRB provides limited bisection to locate eigenvalues for more accuracy.
1.8 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.15 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.8 bertrand 7: *
8: *> \htmlonly
1.15 bertrand 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">
1.8 bertrand 15: *> [TXT]</a>
1.15 bertrand 16: *> \endhtmlonly
1.8 bertrand 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 )
1.15 bertrand 24: *
1.8 bertrand 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: * ..
1.15 bertrand 34: *
1.8 bertrand 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
1.19 bertrand 94: *> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
1.8 bertrand 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
1.17 bertrand 110: *> estimates of the eigenvalues of L D L^T indexed IFIRST through
1.8 bertrand 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
1.19 bertrand 120: *> eigenvalues I and I+1. Note that if IFIRST = ILAST
1.8 bertrand 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: *
1.15 bertrand 176: *> \author Univ. of Tennessee
177: *> \author Univ. of California Berkeley
178: *> \author Univ. of Colorado Denver
179: *> \author NAG Ltd.
1.8 bertrand 180: *
1.15 bertrand 181: *> \ingroup OTHERauxiliary
1.8 bertrand 182: *
183: *> \par Contributors:
184: * ==================
185: *>
186: *> Beresford Parlett, University of California, Berkeley, USA \n
187: *> Jim Demmel, University of California, Berkeley, USA \n
188: *> Inderjit Dhillon, University of Texas, Austin, USA \n
189: *> Osni Marques, LBNL/NERSC, USA \n
190: *> Christof Voemel, University of California, Berkeley, USA
191: *
192: * =====================================================================
1.1 bertrand 193: SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
194: $ RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
195: $ PIVMIN, SPDIAM, TWIST, INFO )
196: *
1.20 ! bertrand 197: * -- LAPACK auxiliary routine --
1.1 bertrand 198: * -- LAPACK is a software package provided by Univ. of Tennessee, --
199: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
200: *
201: * .. Scalar Arguments ..
202: INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
203: DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPDIAM
204: * ..
205: * .. Array Arguments ..
206: INTEGER IWORK( * )
207: DOUBLE PRECISION D( * ), LLD( * ), W( * ),
208: $ WERR( * ), WGAP( * ), WORK( * )
209: * ..
210: *
211: * =====================================================================
212: *
213: * .. Parameters ..
214: DOUBLE PRECISION ZERO, TWO, HALF
215: PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0,
216: $ HALF = 0.5D0 )
217: INTEGER MAXITR
218: * ..
219: * .. Local Scalars ..
220: INTEGER I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
221: $ OLNINT, PREV, R
222: DOUBLE PRECISION BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
223: $ RGAP, RIGHT, TMP, WIDTH
224: * ..
225: * .. External Functions ..
226: INTEGER DLANEG
227: EXTERNAL DLANEG
228: *
229: * ..
230: * .. Intrinsic Functions ..
231: INTRINSIC ABS, MAX, MIN
232: * ..
233: * .. Executable Statements ..
234: *
235: INFO = 0
236: *
1.17 bertrand 237: * Quick return if possible
238: *
239: IF( N.LE.0 ) THEN
240: RETURN
241: END IF
242: *
1.1 bertrand 243: MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
244: $ LOG( TWO ) ) + 2
245: MNWDTH = TWO * PIVMIN
246: *
247: R = TWIST
248: IF((R.LT.1).OR.(R.GT.N)) R = N
249: *
250: * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
251: * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
252: * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
253: * for an unconverged interval is set to the index of the next unconverged
254: * interval, and is -1 or 0 for a converged interval. Thus a linked
255: * list of unconverged intervals is set up.
256: *
257: I1 = IFIRST
258: * The number of unconverged intervals
259: NINT = 0
260: * The last unconverged interval found
261: PREV = 0
262:
263: RGAP = WGAP( I1-OFFSET )
264: DO 75 I = I1, ILAST
265: K = 2*I
266: II = I - OFFSET
267: LEFT = W( II ) - WERR( II )
268: RIGHT = W( II ) + WERR( II )
269: LGAP = RGAP
270: RGAP = WGAP( II )
271: GAP = MIN( LGAP, RGAP )
272:
273: * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
274: * Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
275: *
276: * Do while( NEGCNT(LEFT).GT.I-1 )
277: *
278: BACK = WERR( II )
279: 20 CONTINUE
280: NEGCNT = DLANEG( N, D, LLD, LEFT, PIVMIN, R )
281: IF( NEGCNT.GT.I-1 ) THEN
282: LEFT = LEFT - BACK
283: BACK = TWO*BACK
284: GO TO 20
285: END IF
286: *
287: * Do while( NEGCNT(RIGHT).LT.I )
288: * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
289: *
290: BACK = WERR( II )
291: 50 CONTINUE
292:
293: NEGCNT = DLANEG( N, D, LLD, RIGHT, PIVMIN, R )
294: IF( NEGCNT.LT.I ) THEN
295: RIGHT = RIGHT + BACK
296: BACK = TWO*BACK
297: GO TO 50
298: END IF
299: WIDTH = HALF*ABS( LEFT - RIGHT )
300: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
301: CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
302: IF( WIDTH.LE.CVRGD .OR. WIDTH.LE.MNWDTH ) THEN
303: * This interval has already converged and does not need refinement.
304: * (Note that the gaps might change through refining the
305: * eigenvalues, however, they can only get bigger.)
306: * Remove it from the list.
307: IWORK( K-1 ) = -1
308: * Make sure that I1 always points to the first unconverged interval
309: IF((I.EQ.I1).AND.(I.LT.ILAST)) I1 = I + 1
310: IF((PREV.GE.I1).AND.(I.LE.ILAST)) IWORK( 2*PREV-1 ) = I + 1
311: ELSE
312: * unconverged interval found
313: PREV = I
314: NINT = NINT + 1
315: IWORK( K-1 ) = I + 1
316: IWORK( K ) = NEGCNT
317: END IF
318: WORK( K-1 ) = LEFT
319: WORK( K ) = RIGHT
320: 75 CONTINUE
321:
322: *
323: * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
324: * and while (ITER.LT.MAXITR)
325: *
326: ITER = 0
327: 80 CONTINUE
328: PREV = I1 - 1
329: I = I1
330: OLNINT = NINT
331:
332: DO 100 IP = 1, OLNINT
333: K = 2*I
334: II = I - OFFSET
335: RGAP = WGAP( II )
336: LGAP = RGAP
337: IF(II.GT.1) LGAP = WGAP( II-1 )
338: GAP = MIN( LGAP, RGAP )
339: NEXT = IWORK( K-1 )
340: LEFT = WORK( K-1 )
341: RIGHT = WORK( K )
342: MID = HALF*( LEFT + RIGHT )
343:
344: * semiwidth of interval
345: WIDTH = RIGHT - MID
346: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
347: CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
348: IF( ( WIDTH.LE.CVRGD ) .OR. ( WIDTH.LE.MNWDTH ).OR.
349: $ ( ITER.EQ.MAXITR ) )THEN
350: * reduce number of unconverged intervals
351: NINT = NINT - 1
352: * Mark interval as converged.
353: IWORK( K-1 ) = 0
354: IF( I1.EQ.I ) THEN
355: I1 = NEXT
356: ELSE
357: * Prev holds the last unconverged interval previously examined
358: IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
359: END IF
360: I = NEXT
361: GO TO 100
362: END IF
363: PREV = I
364: *
365: * Perform one bisection step
366: *
367: NEGCNT = DLANEG( N, D, LLD, MID, PIVMIN, R )
368: IF( NEGCNT.LE.I-1 ) THEN
369: WORK( K-1 ) = MID
370: ELSE
371: WORK( K ) = MID
372: END IF
373: I = NEXT
374: 100 CONTINUE
375: ITER = ITER + 1
376: * do another loop if there are still unconverged intervals
377: * However, in the last iteration, all intervals are accepted
378: * since this is the best we can do.
379: IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
380: *
381: *
382: * At this point, all the intervals have converged
383: DO 110 I = IFIRST, ILAST
384: K = 2*I
385: II = I - OFFSET
386: * All intervals marked by '0' have been refined.
387: IF( IWORK( K-1 ).EQ.0 ) THEN
388: W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
389: WERR( II ) = WORK( K ) - W( II )
390: END IF
391: 110 CONTINUE
392: *
393: DO 111 I = IFIRST+1, ILAST
394: K = 2*I
395: II = I - OFFSET
396: WGAP( II-1 ) = MAX( ZERO,
397: $ W(II) - WERR (II) - W( II-1 ) - WERR( II-1 ))
398: 111 CONTINUE
399:
400: RETURN
401: *
402: * End of DLARRB
403: *
404: 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