1: SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND,
2: $ W, WGAP, WERR,
3: $ SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
4: $ DPLUS, LPLUS, WORK, INFO )
5: *
6: * -- LAPACK auxiliary routine (version 3.2) --
7: * -- LAPACK is a software package provided by Univ. of Tennessee, --
8: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9: * November 2006
10: **
11: * .. Scalar Arguments ..
12: INTEGER CLSTRT, CLEND, INFO, N
13: DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
14: * ..
15: * .. Array Arguments ..
16: DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ),
17: $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
18: * ..
19: *
20: * Purpose
21: * =======
22: *
23: * Given the initial representation L D L^T and its cluster of close
24: * eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
25: * W( CLEND ), DLARRF finds a new relatively robust representation
26: * L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
27: * eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
28: *
29: * Arguments
30: * =========
31: *
32: * N (input) INTEGER
33: * The order of the matrix (subblock, if the matrix splitted).
34: *
35: * D (input) DOUBLE PRECISION array, dimension (N)
36: * The N diagonal elements of the diagonal matrix D.
37: *
38: * L (input) DOUBLE PRECISION array, dimension (N-1)
39: * The (N-1) subdiagonal elements of the unit bidiagonal
40: * matrix L.
41: *
42: * LD (input) DOUBLE PRECISION array, dimension (N-1)
43: * The (N-1) elements L(i)*D(i).
44: *
45: * CLSTRT (input) INTEGER
46: * The index of the first eigenvalue in the cluster.
47: *
48: * CLEND (input) INTEGER
49: * The index of the last eigenvalue in the cluster.
50: *
51: * W (input) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1)
52: * The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
53: * W( CLSTRT ) through W( CLEND ) form the cluster of relatively
54: * close eigenalues.
55: *
56: * WGAP (input/output) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1)
57: * The separation from the right neighbor eigenvalue in W.
58: *
59: * WERR (input) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1)
60: * WERR contain the semiwidth of the uncertainty
61: * interval of the corresponding eigenvalue APPROXIMATION in W
62: *
63: * SPDIAM (input) estimate of the spectral diameter obtained from the
64: * Gerschgorin intervals
65: *
66: * CLGAPL, CLGAPR (input) absolute gap on each end of the cluster.
67: * Set by the calling routine to protect against shifts too close
68: * to eigenvalues outside the cluster.
69: *
70: * PIVMIN (input) DOUBLE PRECISION
71: * The minimum pivot allowed in the Sturm sequence.
72: *
73: * SIGMA (output) DOUBLE PRECISION
74: * The shift used to form L(+) D(+) L(+)^T.
75: *
76: * DPLUS (output) DOUBLE PRECISION array, dimension (N)
77: * The N diagonal elements of the diagonal matrix D(+).
78: *
79: * LPLUS (output) DOUBLE PRECISION array, dimension (N-1)
80: * The first (N-1) elements of LPLUS contain the subdiagonal
81: * elements of the unit bidiagonal matrix L(+).
82: *
83: * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
84: * Workspace.
85: *
86: * Further Details
87: * ===============
88: *
89: * Based on contributions by
90: * Beresford Parlett, University of California, Berkeley, USA
91: * Jim Demmel, University of California, Berkeley, USA
92: * Inderjit Dhillon, University of Texas, Austin, USA
93: * Osni Marques, LBNL/NERSC, USA
94: * Christof Voemel, University of California, Berkeley, USA
95: *
96: * =====================================================================
97: *
98: * .. Parameters ..
99: DOUBLE PRECISION FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO,
100: $ ZERO
101: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
102: $ FOUR = 4.0D0, QUART = 0.25D0,
103: $ MAXGROWTH1 = 8.D0,
104: $ MAXGROWTH2 = 8.D0 )
105: * ..
106: * .. Local Scalars ..
107: LOGICAL DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1
108: INTEGER I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT
109: PARAMETER ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 )
110: DOUBLE PRECISION AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL,
111: $ FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA,
112: $ MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX,
113: $ RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2
114: * ..
115: * .. External Functions ..
116: LOGICAL DISNAN
117: DOUBLE PRECISION DLAMCH
118: EXTERNAL DISNAN, DLAMCH
119: * ..
120: * .. External Subroutines ..
121: EXTERNAL DCOPY
122: * ..
123: * .. Intrinsic Functions ..
124: INTRINSIC ABS
125: * ..
126: * .. Executable Statements ..
127: *
128: INFO = 0
129: FACT = DBLE(2**KTRYMAX)
130: EPS = DLAMCH( 'Precision' )
131: SHIFT = 0
132: FORCER = .FALSE.
133:
134:
135: * Note that we cannot guarantee that for any of the shifts tried,
136: * the factorization has a small or even moderate element growth.
137: * There could be Ritz values at both ends of the cluster and despite
138: * backing off, there are examples where all factorizations tried
139: * (in IEEE mode, allowing zero pivots & infinities) have INFINITE
140: * element growth.
141: * For this reason, we should use PIVMIN in this subroutine so that at
142: * least the L D L^T factorization exists. It can be checked afterwards
143: * whether the element growth caused bad residuals/orthogonality.
144:
145: * Decide whether the code should accept the best among all
146: * representations despite large element growth or signal INFO=1
147: NOFAIL = .TRUE.
148: *
149:
150: * Compute the average gap length of the cluster
151: CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT)
152: AVGAP = CLWDTH / DBLE(CLEND-CLSTRT)
153: MINGAP = MIN(CLGAPL, CLGAPR)
154: * Initial values for shifts to both ends of cluster
155: LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT )
156: RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND )
157:
158: * Use a small fudge to make sure that we really shift to the outside
159: LSIGMA = LSIGMA - ABS(LSIGMA)* FOUR * EPS
160: RSIGMA = RSIGMA + ABS(RSIGMA)* FOUR * EPS
161:
162: * Compute upper bounds for how much to back off the initial shifts
163: LDMAX = QUART * MINGAP + TWO * PIVMIN
164: RDMAX = QUART * MINGAP + TWO * PIVMIN
165:
166: LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT
167: RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT
168: *
169: * Initialize the record of the best representation found
170: *
171: S = DLAMCH( 'S' )
172: SMLGROWTH = ONE / S
173: FAIL = DBLE(N-1)*MINGAP/(SPDIAM*EPS)
174: FAIL2 = DBLE(N-1)*MINGAP/(SPDIAM*SQRT(EPS))
175: BESTSHIFT = LSIGMA
176: *
177: * while (KTRY <= KTRYMAX)
178: KTRY = 0
179: GROWTHBOUND = MAXGROWTH1*SPDIAM
180:
181: 5 CONTINUE
182: SAWNAN1 = .FALSE.
183: SAWNAN2 = .FALSE.
184: * Ensure that we do not back off too much of the initial shifts
185: LDELTA = MIN(LDMAX,LDELTA)
186: RDELTA = MIN(RDMAX,RDELTA)
187:
188: * Compute the element growth when shifting to both ends of the cluster
189: * accept the shift if there is no element growth at one of the two ends
190:
191: * Left end
192: S = -LSIGMA
193: DPLUS( 1 ) = D( 1 ) + S
194: IF(ABS(DPLUS(1)).LT.PIVMIN) THEN
195: DPLUS(1) = -PIVMIN
196: * Need to set SAWNAN1 because refined RRR test should not be used
197: * in this case
198: SAWNAN1 = .TRUE.
199: ENDIF
200: MAX1 = ABS( DPLUS( 1 ) )
201: DO 6 I = 1, N - 1
202: LPLUS( I ) = LD( I ) / DPLUS( I )
203: S = S*LPLUS( I )*L( I ) - LSIGMA
204: DPLUS( I+1 ) = D( I+1 ) + S
205: IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN
206: DPLUS(I+1) = -PIVMIN
207: * Need to set SAWNAN1 because refined RRR test should not be used
208: * in this case
209: SAWNAN1 = .TRUE.
210: ENDIF
211: MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) )
212: 6 CONTINUE
213: SAWNAN1 = SAWNAN1 .OR. DISNAN( MAX1 )
214:
215: IF( FORCER .OR.
216: $ (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN
217: SIGMA = LSIGMA
218: SHIFT = SLEFT
219: GOTO 100
220: ENDIF
221:
222: * Right end
223: S = -RSIGMA
224: WORK( 1 ) = D( 1 ) + S
225: IF(ABS(WORK(1)).LT.PIVMIN) THEN
226: WORK(1) = -PIVMIN
227: * Need to set SAWNAN2 because refined RRR test should not be used
228: * in this case
229: SAWNAN2 = .TRUE.
230: ENDIF
231: MAX2 = ABS( WORK( 1 ) )
232: DO 7 I = 1, N - 1
233: WORK( N+I ) = LD( I ) / WORK( I )
234: S = S*WORK( N+I )*L( I ) - RSIGMA
235: WORK( I+1 ) = D( I+1 ) + S
236: IF(ABS(WORK(I+1)).LT.PIVMIN) THEN
237: WORK(I+1) = -PIVMIN
238: * Need to set SAWNAN2 because refined RRR test should not be used
239: * in this case
240: SAWNAN2 = .TRUE.
241: ENDIF
242: MAX2 = MAX( MAX2,ABS(WORK(I+1)) )
243: 7 CONTINUE
244: SAWNAN2 = SAWNAN2 .OR. DISNAN( MAX2 )
245:
246: IF( FORCER .OR.
247: $ (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN
248: SIGMA = RSIGMA
249: SHIFT = SRIGHT
250: GOTO 100
251: ENDIF
252: * If we are at this point, both shifts led to too much element growth
253:
254: * Record the better of the two shifts (provided it didn't lead to NaN)
255: IF(SAWNAN1.AND.SAWNAN2) THEN
256: * both MAX1 and MAX2 are NaN
257: GOTO 50
258: ELSE
259: IF( .NOT.SAWNAN1 ) THEN
260: INDX = 1
261: IF(MAX1.LE.SMLGROWTH) THEN
262: SMLGROWTH = MAX1
263: BESTSHIFT = LSIGMA
264: ENDIF
265: ENDIF
266: IF( .NOT.SAWNAN2 ) THEN
267: IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2
268: IF(MAX2.LE.SMLGROWTH) THEN
269: SMLGROWTH = MAX2
270: BESTSHIFT = RSIGMA
271: ENDIF
272: ENDIF
273: ENDIF
274:
275: * If we are here, both the left and the right shift led to
276: * element growth. If the element growth is moderate, then
277: * we may still accept the representation, if it passes a
278: * refined test for RRR. This test supposes that no NaN occurred.
279: * Moreover, we use the refined RRR test only for isolated clusters.
280: IF((CLWDTH.LT.MINGAP/DBLE(128)) .AND.
281: $ (MIN(MAX1,MAX2).LT.FAIL2)
282: $ .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN
283: DORRR1 = .TRUE.
284: ELSE
285: DORRR1 = .FALSE.
286: ENDIF
287: TRYRRR1 = .TRUE.
288: IF( TRYRRR1 .AND. DORRR1 ) THEN
289: IF(INDX.EQ.1) THEN
290: TMP = ABS( DPLUS( N ) )
291: ZNM2 = ONE
292: PROD = ONE
293: OLDP = ONE
294: DO 15 I = N-1, 1, -1
295: IF( PROD .LE. EPS ) THEN
296: PROD =
297: $ ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP
298: ELSE
299: PROD = PROD*ABS(WORK(N+I))
300: END IF
301: OLDP = PROD
302: ZNM2 = ZNM2 + PROD**2
303: TMP = MAX( TMP, ABS( DPLUS( I ) * PROD ))
304: 15 CONTINUE
305: RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) )
306: IF (RRR1.LE.MAXGROWTH2) THEN
307: SIGMA = LSIGMA
308: SHIFT = SLEFT
309: GOTO 100
310: ENDIF
311: ELSE IF(INDX.EQ.2) THEN
312: TMP = ABS( WORK( N ) )
313: ZNM2 = ONE
314: PROD = ONE
315: OLDP = ONE
316: DO 16 I = N-1, 1, -1
317: IF( PROD .LE. EPS ) THEN
318: PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP
319: ELSE
320: PROD = PROD*ABS(LPLUS(I))
321: END IF
322: OLDP = PROD
323: ZNM2 = ZNM2 + PROD**2
324: TMP = MAX( TMP, ABS( WORK( I ) * PROD ))
325: 16 CONTINUE
326: RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) )
327: IF (RRR2.LE.MAXGROWTH2) THEN
328: SIGMA = RSIGMA
329: SHIFT = SRIGHT
330: GOTO 100
331: ENDIF
332: END IF
333: ENDIF
334:
335: 50 CONTINUE
336:
337: IF (KTRY.LT.KTRYMAX) THEN
338: * If we are here, both shifts failed also the RRR test.
339: * Back off to the outside
340: LSIGMA = MAX( LSIGMA - LDELTA,
341: $ LSIGMA - LDMAX)
342: RSIGMA = MIN( RSIGMA + RDELTA,
343: $ RSIGMA + RDMAX )
344: LDELTA = TWO * LDELTA
345: RDELTA = TWO * RDELTA
346: KTRY = KTRY + 1
347: GOTO 5
348: ELSE
349: * None of the representations investigated satisfied our
350: * criteria. Take the best one we found.
351: IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN
352: LSIGMA = BESTSHIFT
353: RSIGMA = BESTSHIFT
354: FORCER = .TRUE.
355: GOTO 5
356: ELSE
357: INFO = 1
358: RETURN
359: ENDIF
360: END IF
361:
362: 100 CONTINUE
363: IF (SHIFT.EQ.SLEFT) THEN
364: ELSEIF (SHIFT.EQ.SRIGHT) THEN
365: * store new L and D back into DPLUS, LPLUS
366: CALL DCOPY( N, WORK, 1, DPLUS, 1 )
367: CALL DCOPY( N-1, WORK(N+1), 1, LPLUS, 1 )
368: ENDIF
369:
370: RETURN
371: *
372: * End of DLARRF
373: *
374: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarrf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack