Annotation of rpl/lapack/lapack/dlarrj.f, revision 1.4
1.1 bertrand 1: SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST,
2: $ RTOL, OFFSET, W, WERR, WORK, IWORK,
3: $ PIVMIN, SPDIAM, 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
12: DOUBLE PRECISION PIVMIN, RTOL, SPDIAM
13: * ..
14: * .. Array Arguments ..
15: INTEGER IWORK( * )
16: DOUBLE PRECISION D( * ), E2( * ), W( * ),
17: $ WERR( * ), WORK( * )
18: * ..
19: *
20: * Purpose
21: * =======
22: *
23: * Given the initial eigenvalue approximations of T, DLARRJ
24: * does bisection to refine the eigenvalues of 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 in WERR. During bisection, intervals
28: * [left, right] are maintained by storing their mid-points and
29: * semi-widths in the arrays W and WERR respectively.
30: *
31: * Arguments
32: * =========
33: *
34: * N (input) INTEGER
35: * The order of the matrix.
36: *
37: * D (input) DOUBLE PRECISION array, dimension (N)
38: * The N diagonal elements of T.
39: *
40: * E2 (input) DOUBLE PRECISION array, dimension (N-1)
41: * The Squares of the (N-1) subdiagonal elements of T.
42: *
43: * IFIRST (input) INTEGER
44: * The index of the first eigenvalue to be computed.
45: *
46: * ILAST (input) INTEGER
47: * The index of the last eigenvalue to be computed.
48: *
49: * RTOL (input) DOUBLE PRECISION
50: * Tolerance for the convergence of the bisection intervals.
51: * An interval [LEFT,RIGHT] has converged if
52: * RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|).
53: *
54: * OFFSET (input) INTEGER
55: * Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET
56: * through ILAST-OFFSET elements of these arrays are to be used.
57: *
58: * W (input/output) DOUBLE PRECISION array, dimension (N)
59: * On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
60: * estimates of the eigenvalues of L D L^T indexed IFIRST through
61: * ILAST.
62: * On output, these estimates are refined.
63: *
64: * WERR (input/output) DOUBLE PRECISION array, dimension (N)
65: * On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
66: * the errors in the estimates of the corresponding elements in W.
67: * On output, these errors are refined.
68: *
69: * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
70: * Workspace.
71: *
72: * IWORK (workspace) INTEGER array, dimension (2*N)
73: * Workspace.
74: *
75: * PIVMIN (input) DOUBLE PRECISION
76: * The minimum pivot in the Sturm sequence for T.
77: *
78: * SPDIAM (input) DOUBLE PRECISION
79: * The spectral diameter of T.
80: *
81: * INFO (output) INTEGER
82: * Error flag.
83: *
84: * Further Details
85: * ===============
86: *
87: * Based on contributions by
88: * Beresford Parlett, University of California, Berkeley, USA
89: * Jim Demmel, University of California, Berkeley, USA
90: * Inderjit Dhillon, University of Texas, Austin, USA
91: * Osni Marques, LBNL/NERSC, USA
92: * Christof Voemel, University of California, Berkeley, USA
93: *
94: * =====================================================================
95: *
96: * .. Parameters ..
97: DOUBLE PRECISION ZERO, ONE, TWO, HALF
98: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
99: $ HALF = 0.5D0 )
100: INTEGER MAXITR
101: * ..
102: * .. Local Scalars ..
103: INTEGER CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
104: $ OLNINT, P, PREV, SAVI1
105: DOUBLE PRECISION DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH
106: *
107: * ..
108: * .. Intrinsic Functions ..
109: INTRINSIC ABS, MAX
110: * ..
111: * .. Executable Statements ..
112: *
113: INFO = 0
114: *
115: MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
116: $ LOG( TWO ) ) + 2
117: *
118: * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
119: * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
120: * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
121: * for an unconverged interval is set to the index of the next unconverged
122: * interval, and is -1 or 0 for a converged interval. Thus a linked
123: * list of unconverged intervals is set up.
124: *
125:
126: I1 = IFIRST
127: I2 = ILAST
128: * The number of unconverged intervals
129: NINT = 0
130: * The last unconverged interval found
131: PREV = 0
132: DO 75 I = I1, I2
133: K = 2*I
134: II = I - OFFSET
135: LEFT = W( II ) - WERR( II )
136: MID = W(II)
137: RIGHT = W( II ) + WERR( II )
138: WIDTH = RIGHT - MID
139: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
140:
141: * The following test prevents the test of converged intervals
142: IF( WIDTH.LT.RTOL*TMP ) THEN
143: * This interval has already converged and does not need refinement.
144: * (Note that the gaps might change through refining the
145: * eigenvalues, however, they can only get bigger.)
146: * Remove it from the list.
147: IWORK( K-1 ) = -1
148: * Make sure that I1 always points to the first unconverged interval
149: IF((I.EQ.I1).AND.(I.LT.I2)) I1 = I + 1
150: IF((PREV.GE.I1).AND.(I.LE.I2)) IWORK( 2*PREV-1 ) = I + 1
151: ELSE
152: * unconverged interval found
153: PREV = I
154: * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
155: *
156: * Do while( CNT(LEFT).GT.I-1 )
157: *
158: FAC = ONE
159: 20 CONTINUE
160: CNT = 0
161: S = LEFT
162: DPLUS = D( 1 ) - S
163: IF( DPLUS.LT.ZERO ) CNT = CNT + 1
164: DO 30 J = 2, N
165: DPLUS = D( J ) - S - E2( J-1 )/DPLUS
166: IF( DPLUS.LT.ZERO ) CNT = CNT + 1
167: 30 CONTINUE
168: IF( CNT.GT.I-1 ) THEN
169: LEFT = LEFT - WERR( II )*FAC
170: FAC = TWO*FAC
171: GO TO 20
172: END IF
173: *
174: * Do while( CNT(RIGHT).LT.I )
175: *
176: FAC = ONE
177: 50 CONTINUE
178: CNT = 0
179: S = RIGHT
180: DPLUS = D( 1 ) - S
181: IF( DPLUS.LT.ZERO ) CNT = CNT + 1
182: DO 60 J = 2, N
183: DPLUS = D( J ) - S - E2( J-1 )/DPLUS
184: IF( DPLUS.LT.ZERO ) CNT = CNT + 1
185: 60 CONTINUE
186: IF( CNT.LT.I ) THEN
187: RIGHT = RIGHT + WERR( II )*FAC
188: FAC = TWO*FAC
189: GO TO 50
190: END IF
191: NINT = NINT + 1
192: IWORK( K-1 ) = I + 1
193: IWORK( K ) = CNT
194: END IF
195: WORK( K-1 ) = LEFT
196: WORK( K ) = RIGHT
197: 75 CONTINUE
198:
199:
200: SAVI1 = I1
201: *
202: * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
203: * and while (ITER.LT.MAXITR)
204: *
205: ITER = 0
206: 80 CONTINUE
207: PREV = I1 - 1
208: I = I1
209: OLNINT = NINT
210:
211: DO 100 P = 1, OLNINT
212: K = 2*I
213: II = I - OFFSET
214: NEXT = IWORK( K-1 )
215: LEFT = WORK( K-1 )
216: RIGHT = WORK( K )
217: MID = HALF*( LEFT + RIGHT )
218:
219: * semiwidth of interval
220: WIDTH = RIGHT - MID
221: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
222:
223: IF( ( WIDTH.LT.RTOL*TMP ) .OR.
224: $ (ITER.EQ.MAXITR) )THEN
225: * reduce number of unconverged intervals
226: NINT = NINT - 1
227: * Mark interval as converged.
228: IWORK( K-1 ) = 0
229: IF( I1.EQ.I ) THEN
230: I1 = NEXT
231: ELSE
232: * Prev holds the last unconverged interval previously examined
233: IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
234: END IF
235: I = NEXT
236: GO TO 100
237: END IF
238: PREV = I
239: *
240: * Perform one bisection step
241: *
242: CNT = 0
243: S = MID
244: DPLUS = D( 1 ) - S
245: IF( DPLUS.LT.ZERO ) CNT = CNT + 1
246: DO 90 J = 2, N
247: DPLUS = D( J ) - S - E2( J-1 )/DPLUS
248: IF( DPLUS.LT.ZERO ) CNT = CNT + 1
249: 90 CONTINUE
250: IF( CNT.LE.I-1 ) THEN
251: WORK( K-1 ) = MID
252: ELSE
253: WORK( K ) = MID
254: END IF
255: I = NEXT
256:
257: 100 CONTINUE
258: ITER = ITER + 1
259: * do another loop if there are still unconverged intervals
260: * However, in the last iteration, all intervals are accepted
261: * since this is the best we can do.
262: IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
263: *
264: *
265: * At this point, all the intervals have converged
266: DO 110 I = SAVI1, ILAST
267: K = 2*I
268: II = I - OFFSET
269: * All intervals marked by '0' have been refined.
270: IF( IWORK( K-1 ).EQ.0 ) THEN
271: W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
272: WERR( II ) = WORK( K ) - W( II )
273: END IF
274: 110 CONTINUE
275: *
276:
277: RETURN
278: *
279: * End of DLARRJ
280: *
281: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarrj.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack