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