Annotation of rpl/lapack/lapack/dlarrb.f, revision 1.12
1.11 bertrand 1: *> \brief \b DLARRB provides limited bisection to locate eigenvalues for more accuracy.
1.8 bertrand 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: *
1.11 bertrand 181: *> \date September 2012
1.8 bertrand 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.11 bertrand 199: * -- LAPACK auxiliary routine (version 3.4.2) --
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.11 bertrand 202: * September 2012
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