1: *> \brief \b DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLARRF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrf.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND,
22: * W, WGAP, WERR,
23: * SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
24: * DPLUS, LPLUS, WORK, INFO )
25: *
26: * .. Scalar Arguments ..
27: * INTEGER CLSTRT, CLEND, INFO, N
28: * DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
29: * ..
30: * .. Array Arguments ..
31: * DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ),
32: * $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
33: * ..
34: *
35: *
36: *> \par Purpose:
37: * =============
38: *>
39: *> \verbatim
40: *>
41: *> Given the initial representation L D L^T and its cluster of close
42: *> eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
43: *> W( CLEND ), DLARRF finds a new relatively robust representation
44: *> L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
45: *> eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
46: *> \endverbatim
47: *
48: * Arguments:
49: * ==========
50: *
51: *> \param[in] N
52: *> \verbatim
53: *> N is INTEGER
54: *> The order of the matrix (subblock, if the matrix split).
55: *> \endverbatim
56: *>
57: *> \param[in] D
58: *> \verbatim
59: *> D is DOUBLE PRECISION array, dimension (N)
60: *> The N diagonal elements of the diagonal matrix D.
61: *> \endverbatim
62: *>
63: *> \param[in] L
64: *> \verbatim
65: *> L is DOUBLE PRECISION array, dimension (N-1)
66: *> The (N-1) subdiagonal elements of the unit bidiagonal
67: *> matrix L.
68: *> \endverbatim
69: *>
70: *> \param[in] LD
71: *> \verbatim
72: *> LD is DOUBLE PRECISION array, dimension (N-1)
73: *> The (N-1) elements L(i)*D(i).
74: *> \endverbatim
75: *>
76: *> \param[in] CLSTRT
77: *> \verbatim
78: *> CLSTRT is INTEGER
79: *> The index of the first eigenvalue in the cluster.
80: *> \endverbatim
81: *>
82: *> \param[in] CLEND
83: *> \verbatim
84: *> CLEND is INTEGER
85: *> The index of the last eigenvalue in the cluster.
86: *> \endverbatim
87: *>
88: *> \param[in] W
89: *> \verbatim
90: *> W is DOUBLE PRECISION array, dimension
91: *> dimension is >= (CLEND-CLSTRT+1)
92: *> The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
93: *> W( CLSTRT ) through W( CLEND ) form the cluster of relatively
94: *> close eigenalues.
95: *> \endverbatim
96: *>
97: *> \param[in,out] WGAP
98: *> \verbatim
99: *> WGAP is DOUBLE PRECISION array, dimension
100: *> dimension is >= (CLEND-CLSTRT+1)
101: *> The separation from the right neighbor eigenvalue in W.
102: *> \endverbatim
103: *>
104: *> \param[in] WERR
105: *> \verbatim
106: *> WERR is DOUBLE PRECISION array, dimension
107: *> dimension is >= (CLEND-CLSTRT+1)
108: *> WERR contain the semiwidth of the uncertainty
109: *> interval of the corresponding eigenvalue APPROXIMATION in W
110: *> \endverbatim
111: *>
112: *> \param[in] SPDIAM
113: *> \verbatim
114: *> SPDIAM is DOUBLE PRECISION
115: *> estimate of the spectral diameter obtained from the
116: *> Gerschgorin intervals
117: *> \endverbatim
118: *>
119: *> \param[in] CLGAPL
120: *> \verbatim
121: *> CLGAPL is DOUBLE PRECISION
122: *> \endverbatim
123: *>
124: *> \param[in] CLGAPR
125: *> \verbatim
126: *> CLGAPR is DOUBLE PRECISION
127: *> absolute gap on each end of the cluster.
128: *> Set by the calling routine to protect against shifts too close
129: *> to eigenvalues outside the cluster.
130: *> \endverbatim
131: *>
132: *> \param[in] PIVMIN
133: *> \verbatim
134: *> PIVMIN is DOUBLE PRECISION
135: *> The minimum pivot allowed in the Sturm sequence.
136: *> \endverbatim
137: *>
138: *> \param[out] SIGMA
139: *> \verbatim
140: *> SIGMA is DOUBLE PRECISION
141: *> The shift used to form L(+) D(+) L(+)^T.
142: *> \endverbatim
143: *>
144: *> \param[out] DPLUS
145: *> \verbatim
146: *> DPLUS is DOUBLE PRECISION array, dimension (N)
147: *> The N diagonal elements of the diagonal matrix D(+).
148: *> \endverbatim
149: *>
150: *> \param[out] LPLUS
151: *> \verbatim
152: *> LPLUS is DOUBLE PRECISION array, dimension (N-1)
153: *> The first (N-1) elements of LPLUS contain the subdiagonal
154: *> elements of the unit bidiagonal matrix L(+).
155: *> \endverbatim
156: *>
157: *> \param[out] WORK
158: *> \verbatim
159: *> WORK is DOUBLE PRECISION array, dimension (2*N)
160: *> Workspace.
161: *> \endverbatim
162: *>
163: *> \param[out] INFO
164: *> \verbatim
165: *> INFO is INTEGER
166: *> Signals processing OK (=0) or failure (=1)
167: *> \endverbatim
168: *
169: * Authors:
170: * ========
171: *
172: *> \author Univ. of Tennessee
173: *> \author Univ. of California Berkeley
174: *> \author Univ. of Colorado Denver
175: *> \author NAG Ltd.
176: *
177: *> \date June 2016
178: *
179: *> \ingroup OTHERauxiliary
180: *
181: *> \par Contributors:
182: * ==================
183: *>
184: *> Beresford Parlett, University of California, Berkeley, USA \n
185: *> Jim Demmel, University of California, Berkeley, USA \n
186: *> Inderjit Dhillon, University of Texas, Austin, USA \n
187: *> Osni Marques, LBNL/NERSC, USA \n
188: *> Christof Voemel, University of California, Berkeley, USA
189: *
190: * =====================================================================
191: SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND,
192: $ W, WGAP, WERR,
193: $ SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
194: $ DPLUS, LPLUS, WORK, INFO )
195: *
196: * -- LAPACK auxiliary routine (version 3.7.1) --
197: * -- LAPACK is a software package provided by Univ. of Tennessee, --
198: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
199: * June 2016
200: *
201: * .. Scalar Arguments ..
202: INTEGER CLSTRT, CLEND, INFO, N
203: DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
204: * ..
205: * .. Array Arguments ..
206: DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ),
207: $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
208: * ..
209: *
210: * =====================================================================
211: *
212: * .. Parameters ..
213: DOUBLE PRECISION FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO
214: PARAMETER ( ONE = 1.0D0, TWO = 2.0D0, FOUR = 4.0D0,
215: $ QUART = 0.25D0,
216: $ MAXGROWTH1 = 8.D0,
217: $ MAXGROWTH2 = 8.D0 )
218: * ..
219: * .. Local Scalars ..
220: LOGICAL DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1
221: INTEGER I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT
222: PARAMETER ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 )
223: DOUBLE PRECISION AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL,
224: $ FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA,
225: $ MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX,
226: $ RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2
227: * ..
228: * .. External Functions ..
229: LOGICAL DISNAN
230: DOUBLE PRECISION DLAMCH
231: EXTERNAL DISNAN, DLAMCH
232: * ..
233: * .. External Subroutines ..
234: EXTERNAL DCOPY
235: * ..
236: * .. Intrinsic Functions ..
237: INTRINSIC ABS
238: * ..
239: * .. Executable Statements ..
240: *
241: INFO = 0
242: *
243: * Quick return if possible
244: *
245: IF( N.LE.0 ) THEN
246: RETURN
247: END IF
248: *
249: FACT = DBLE(2**KTRYMAX)
250: EPS = DLAMCH( 'Precision' )
251: SHIFT = 0
252: FORCER = .FALSE.
253:
254:
255: * Note that we cannot guarantee that for any of the shifts tried,
256: * the factorization has a small or even moderate element growth.
257: * There could be Ritz values at both ends of the cluster and despite
258: * backing off, there are examples where all factorizations tried
259: * (in IEEE mode, allowing zero pivots & infinities) have INFINITE
260: * element growth.
261: * For this reason, we should use PIVMIN in this subroutine so that at
262: * least the L D L^T factorization exists. It can be checked afterwards
263: * whether the element growth caused bad residuals/orthogonality.
264:
265: * Decide whether the code should accept the best among all
266: * representations despite large element growth or signal INFO=1
267: * Setting NOFAIL to .FALSE. for quick fix for bug 113
268: NOFAIL = .FALSE.
269: *
270:
271: * Compute the average gap length of the cluster
272: CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT)
273: AVGAP = CLWDTH / DBLE(CLEND-CLSTRT)
274: MINGAP = MIN(CLGAPL, CLGAPR)
275: * Initial values for shifts to both ends of cluster
276: LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT )
277: RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND )
278:
279: * Use a small fudge to make sure that we really shift to the outside
280: LSIGMA = LSIGMA - ABS(LSIGMA)* FOUR * EPS
281: RSIGMA = RSIGMA + ABS(RSIGMA)* FOUR * EPS
282:
283: * Compute upper bounds for how much to back off the initial shifts
284: LDMAX = QUART * MINGAP + TWO * PIVMIN
285: RDMAX = QUART * MINGAP + TWO * PIVMIN
286:
287: LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT
288: RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT
289: *
290: * Initialize the record of the best representation found
291: *
292: S = DLAMCH( 'S' )
293: SMLGROWTH = ONE / S
294: FAIL = DBLE(N-1)*MINGAP/(SPDIAM*EPS)
295: FAIL2 = DBLE(N-1)*MINGAP/(SPDIAM*SQRT(EPS))
296: BESTSHIFT = LSIGMA
297: *
298: * while (KTRY <= KTRYMAX)
299: KTRY = 0
300: GROWTHBOUND = MAXGROWTH1*SPDIAM
301:
302: 5 CONTINUE
303: SAWNAN1 = .FALSE.
304: SAWNAN2 = .FALSE.
305: * Ensure that we do not back off too much of the initial shifts
306: LDELTA = MIN(LDMAX,LDELTA)
307: RDELTA = MIN(RDMAX,RDELTA)
308:
309: * Compute the element growth when shifting to both ends of the cluster
310: * accept the shift if there is no element growth at one of the two ends
311:
312: * Left end
313: S = -LSIGMA
314: DPLUS( 1 ) = D( 1 ) + S
315: IF(ABS(DPLUS(1)).LT.PIVMIN) THEN
316: DPLUS(1) = -PIVMIN
317: * Need to set SAWNAN1 because refined RRR test should not be used
318: * in this case
319: SAWNAN1 = .TRUE.
320: ENDIF
321: MAX1 = ABS( DPLUS( 1 ) )
322: DO 6 I = 1, N - 1
323: LPLUS( I ) = LD( I ) / DPLUS( I )
324: S = S*LPLUS( I )*L( I ) - LSIGMA
325: DPLUS( I+1 ) = D( I+1 ) + S
326: IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN
327: DPLUS(I+1) = -PIVMIN
328: * Need to set SAWNAN1 because refined RRR test should not be used
329: * in this case
330: SAWNAN1 = .TRUE.
331: ENDIF
332: MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) )
333: 6 CONTINUE
334: SAWNAN1 = SAWNAN1 .OR. DISNAN( MAX1 )
335:
336: IF( FORCER .OR.
337: $ (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN
338: SIGMA = LSIGMA
339: SHIFT = SLEFT
340: GOTO 100
341: ENDIF
342:
343: * Right end
344: S = -RSIGMA
345: WORK( 1 ) = D( 1 ) + S
346: IF(ABS(WORK(1)).LT.PIVMIN) THEN
347: WORK(1) = -PIVMIN
348: * Need to set SAWNAN2 because refined RRR test should not be used
349: * in this case
350: SAWNAN2 = .TRUE.
351: ENDIF
352: MAX2 = ABS( WORK( 1 ) )
353: DO 7 I = 1, N - 1
354: WORK( N+I ) = LD( I ) / WORK( I )
355: S = S*WORK( N+I )*L( I ) - RSIGMA
356: WORK( I+1 ) = D( I+1 ) + S
357: IF(ABS(WORK(I+1)).LT.PIVMIN) THEN
358: WORK(I+1) = -PIVMIN
359: * Need to set SAWNAN2 because refined RRR test should not be used
360: * in this case
361: SAWNAN2 = .TRUE.
362: ENDIF
363: MAX2 = MAX( MAX2,ABS(WORK(I+1)) )
364: 7 CONTINUE
365: SAWNAN2 = SAWNAN2 .OR. DISNAN( MAX2 )
366:
367: IF( FORCER .OR.
368: $ (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN
369: SIGMA = RSIGMA
370: SHIFT = SRIGHT
371: GOTO 100
372: ENDIF
373: * If we are at this point, both shifts led to too much element growth
374:
375: * Record the better of the two shifts (provided it didn't lead to NaN)
376: IF(SAWNAN1.AND.SAWNAN2) THEN
377: * both MAX1 and MAX2 are NaN
378: GOTO 50
379: ELSE
380: IF( .NOT.SAWNAN1 ) THEN
381: INDX = 1
382: IF(MAX1.LE.SMLGROWTH) THEN
383: SMLGROWTH = MAX1
384: BESTSHIFT = LSIGMA
385: ENDIF
386: ENDIF
387: IF( .NOT.SAWNAN2 ) THEN
388: IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2
389: IF(MAX2.LE.SMLGROWTH) THEN
390: SMLGROWTH = MAX2
391: BESTSHIFT = RSIGMA
392: ENDIF
393: ENDIF
394: ENDIF
395:
396: * If we are here, both the left and the right shift led to
397: * element growth. If the element growth is moderate, then
398: * we may still accept the representation, if it passes a
399: * refined test for RRR. This test supposes that no NaN occurred.
400: * Moreover, we use the refined RRR test only for isolated clusters.
401: IF((CLWDTH.LT.MINGAP/DBLE(128)) .AND.
402: $ (MIN(MAX1,MAX2).LT.FAIL2)
403: $ .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN
404: DORRR1 = .TRUE.
405: ELSE
406: DORRR1 = .FALSE.
407: ENDIF
408: TRYRRR1 = .TRUE.
409: IF( TRYRRR1 .AND. DORRR1 ) THEN
410: IF(INDX.EQ.1) THEN
411: TMP = ABS( DPLUS( N ) )
412: ZNM2 = ONE
413: PROD = ONE
414: OLDP = ONE
415: DO 15 I = N-1, 1, -1
416: IF( PROD .LE. EPS ) THEN
417: PROD =
418: $ ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP
419: ELSE
420: PROD = PROD*ABS(WORK(N+I))
421: END IF
422: OLDP = PROD
423: ZNM2 = ZNM2 + PROD**2
424: TMP = MAX( TMP, ABS( DPLUS( I ) * PROD ))
425: 15 CONTINUE
426: RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) )
427: IF (RRR1.LE.MAXGROWTH2) THEN
428: SIGMA = LSIGMA
429: SHIFT = SLEFT
430: GOTO 100
431: ENDIF
432: ELSE IF(INDX.EQ.2) THEN
433: TMP = ABS( WORK( N ) )
434: ZNM2 = ONE
435: PROD = ONE
436: OLDP = ONE
437: DO 16 I = N-1, 1, -1
438: IF( PROD .LE. EPS ) THEN
439: PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP
440: ELSE
441: PROD = PROD*ABS(LPLUS(I))
442: END IF
443: OLDP = PROD
444: ZNM2 = ZNM2 + PROD**2
445: TMP = MAX( TMP, ABS( WORK( I ) * PROD ))
446: 16 CONTINUE
447: RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) )
448: IF (RRR2.LE.MAXGROWTH2) THEN
449: SIGMA = RSIGMA
450: SHIFT = SRIGHT
451: GOTO 100
452: ENDIF
453: END IF
454: ENDIF
455:
456: 50 CONTINUE
457:
458: IF (KTRY.LT.KTRYMAX) THEN
459: * If we are here, both shifts failed also the RRR test.
460: * Back off to the outside
461: LSIGMA = MAX( LSIGMA - LDELTA,
462: $ LSIGMA - LDMAX)
463: RSIGMA = MIN( RSIGMA + RDELTA,
464: $ RSIGMA + RDMAX )
465: LDELTA = TWO * LDELTA
466: RDELTA = TWO * RDELTA
467: KTRY = KTRY + 1
468: GOTO 5
469: ELSE
470: * None of the representations investigated satisfied our
471: * criteria. Take the best one we found.
472: IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN
473: LSIGMA = BESTSHIFT
474: RSIGMA = BESTSHIFT
475: FORCER = .TRUE.
476: GOTO 5
477: ELSE
478: INFO = 1
479: RETURN
480: ENDIF
481: END IF
482:
483: 100 CONTINUE
484: IF (SHIFT.EQ.SLEFT) THEN
485: ELSEIF (SHIFT.EQ.SRIGHT) THEN
486: * store new L and D back into DPLUS, LPLUS
487: CALL DCOPY( N, WORK, 1, DPLUS, 1 )
488: CALL DCOPY( N-1, WORK(N+1), 1, LPLUS, 1 )
489: ENDIF
490:
491: RETURN
492: *
493: * End of DLARRF
494: *
495: 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