Annotation of rpl/lapack/lapack/dlarrb.f, revision 1.18
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
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
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
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: *
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.17 bertrand 181: *> \date June 2017
1.8 bertrand 182: *
1.15 bertrand 183: *> \ingroup OTHERauxiliary
1.8 bertrand 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.17 bertrand 199: * -- LAPACK auxiliary routine (version 3.7.1) --
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.17 bertrand 202: * June 2017
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: *
1.17 bertrand 240: * Quick return if possible
241: *
242: IF( N.LE.0 ) THEN
243: RETURN
244: END IF
245: *
1.1 bertrand 246: MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
247: $ LOG( TWO ) ) + 2
248: MNWDTH = TWO * PIVMIN
249: *
250: R = TWIST
251: IF((R.LT.1).OR.(R.GT.N)) R = N
252: *
253: * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
254: * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
255: * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
256: * for an unconverged interval is set to the index of the next unconverged
257: * interval, and is -1 or 0 for a converged interval. Thus a linked
258: * list of unconverged intervals is set up.
259: *
260: I1 = IFIRST
261: * The number of unconverged intervals
262: NINT = 0
263: * The last unconverged interval found
264: PREV = 0
265:
266: RGAP = WGAP( I1-OFFSET )
267: DO 75 I = I1, ILAST
268: K = 2*I
269: II = I - OFFSET
270: LEFT = W( II ) - WERR( II )
271: RIGHT = W( II ) + WERR( II )
272: LGAP = RGAP
273: RGAP = WGAP( II )
274: GAP = MIN( LGAP, RGAP )
275:
276: * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
277: * Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
278: *
279: * Do while( NEGCNT(LEFT).GT.I-1 )
280: *
281: BACK = WERR( II )
282: 20 CONTINUE
283: NEGCNT = DLANEG( N, D, LLD, LEFT, PIVMIN, R )
284: IF( NEGCNT.GT.I-1 ) THEN
285: LEFT = LEFT - BACK
286: BACK = TWO*BACK
287: GO TO 20
288: END IF
289: *
290: * Do while( NEGCNT(RIGHT).LT.I )
291: * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
292: *
293: BACK = WERR( II )
294: 50 CONTINUE
295:
296: NEGCNT = DLANEG( N, D, LLD, RIGHT, PIVMIN, R )
297: IF( NEGCNT.LT.I ) THEN
298: RIGHT = RIGHT + BACK
299: BACK = TWO*BACK
300: GO TO 50
301: END IF
302: WIDTH = HALF*ABS( LEFT - RIGHT )
303: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
304: CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
305: IF( WIDTH.LE.CVRGD .OR. WIDTH.LE.MNWDTH ) THEN
306: * This interval has already converged and does not need refinement.
307: * (Note that the gaps might change through refining the
308: * eigenvalues, however, they can only get bigger.)
309: * Remove it from the list.
310: IWORK( K-1 ) = -1
311: * Make sure that I1 always points to the first unconverged interval
312: IF((I.EQ.I1).AND.(I.LT.ILAST)) I1 = I + 1
313: IF((PREV.GE.I1).AND.(I.LE.ILAST)) IWORK( 2*PREV-1 ) = I + 1
314: ELSE
315: * unconverged interval found
316: PREV = I
317: NINT = NINT + 1
318: IWORK( K-1 ) = I + 1
319: IWORK( K ) = NEGCNT
320: END IF
321: WORK( K-1 ) = LEFT
322: WORK( K ) = RIGHT
323: 75 CONTINUE
324:
325: *
326: * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
327: * and while (ITER.LT.MAXITR)
328: *
329: ITER = 0
330: 80 CONTINUE
331: PREV = I1 - 1
332: I = I1
333: OLNINT = NINT
334:
335: DO 100 IP = 1, OLNINT
336: K = 2*I
337: II = I - OFFSET
338: RGAP = WGAP( II )
339: LGAP = RGAP
340: IF(II.GT.1) LGAP = WGAP( II-1 )
341: GAP = MIN( LGAP, RGAP )
342: NEXT = IWORK( K-1 )
343: LEFT = WORK( K-1 )
344: RIGHT = WORK( K )
345: MID = HALF*( LEFT + RIGHT )
346:
347: * semiwidth of interval
348: WIDTH = RIGHT - MID
349: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
350: CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
351: IF( ( WIDTH.LE.CVRGD ) .OR. ( WIDTH.LE.MNWDTH ).OR.
352: $ ( ITER.EQ.MAXITR ) )THEN
353: * reduce number of unconverged intervals
354: NINT = NINT - 1
355: * Mark interval as converged.
356: IWORK( K-1 ) = 0
357: IF( I1.EQ.I ) THEN
358: I1 = NEXT
359: ELSE
360: * Prev holds the last unconverged interval previously examined
361: IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
362: END IF
363: I = NEXT
364: GO TO 100
365: END IF
366: PREV = I
367: *
368: * Perform one bisection step
369: *
370: NEGCNT = DLANEG( N, D, LLD, MID, PIVMIN, R )
371: IF( NEGCNT.LE.I-1 ) THEN
372: WORK( K-1 ) = MID
373: ELSE
374: WORK( K ) = MID
375: END IF
376: I = NEXT
377: 100 CONTINUE
378: ITER = ITER + 1
379: * do another loop if there are still unconverged intervals
380: * However, in the last iteration, all intervals are accepted
381: * since this is the best we can do.
382: IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
383: *
384: *
385: * At this point, all the intervals have converged
386: DO 110 I = IFIRST, ILAST
387: K = 2*I
388: II = I - OFFSET
389: * All intervals marked by '0' have been refined.
390: IF( IWORK( K-1 ).EQ.0 ) THEN
391: W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
392: WERR( II ) = WORK( K ) - W( II )
393: END IF
394: 110 CONTINUE
395: *
396: DO 111 I = IFIRST+1, ILAST
397: K = 2*I
398: II = I - OFFSET
399: WGAP( II-1 ) = MAX( ZERO,
400: $ W(II) - WERR (II) - W( II-1 ) - WERR( II-1 ))
401: 111 CONTINUE
402:
403: RETURN
404: *
405: * End of DLARRB
406: *
407: 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