Annotation of rpl/lapack/lapack/dlarrf.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLARRF
! 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 splitted).
! 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 November 2011
! 178: *
! 179: *> \ingroup auxOTHERauxiliary
! 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: * =====================================================================
1.1 bertrand 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: *
1.9 ! bertrand 196: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 197: * -- LAPACK is a software package provided by Univ. of Tennessee, --
198: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 199: * November 2011
! 200: *
1.1 bertrand 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: $ ZERO
215: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
216: $ FOUR = 4.0D0, QUART = 0.25D0,
217: $ MAXGROWTH1 = 8.D0,
218: $ MAXGROWTH2 = 8.D0 )
219: * ..
220: * .. Local Scalars ..
221: LOGICAL DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1
222: INTEGER I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT
223: PARAMETER ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 )
224: DOUBLE PRECISION AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL,
225: $ FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA,
226: $ MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX,
227: $ RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2
228: * ..
229: * .. External Functions ..
230: LOGICAL DISNAN
231: DOUBLE PRECISION DLAMCH
232: EXTERNAL DISNAN, DLAMCH
233: * ..
234: * .. External Subroutines ..
235: EXTERNAL DCOPY
236: * ..
237: * .. Intrinsic Functions ..
238: INTRINSIC ABS
239: * ..
240: * .. Executable Statements ..
241: *
242: INFO = 0
243: FACT = DBLE(2**KTRYMAX)
244: EPS = DLAMCH( 'Precision' )
245: SHIFT = 0
246: FORCER = .FALSE.
247:
248:
249: * Note that we cannot guarantee that for any of the shifts tried,
250: * the factorization has a small or even moderate element growth.
251: * There could be Ritz values at both ends of the cluster and despite
252: * backing off, there are examples where all factorizations tried
253: * (in IEEE mode, allowing zero pivots & infinities) have INFINITE
254: * element growth.
255: * For this reason, we should use PIVMIN in this subroutine so that at
256: * least the L D L^T factorization exists. It can be checked afterwards
257: * whether the element growth caused bad residuals/orthogonality.
258:
259: * Decide whether the code should accept the best among all
260: * representations despite large element growth or signal INFO=1
261: NOFAIL = .TRUE.
262: *
263:
264: * Compute the average gap length of the cluster
265: CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT)
266: AVGAP = CLWDTH / DBLE(CLEND-CLSTRT)
267: MINGAP = MIN(CLGAPL, CLGAPR)
268: * Initial values for shifts to both ends of cluster
269: LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT )
270: RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND )
271:
272: * Use a small fudge to make sure that we really shift to the outside
273: LSIGMA = LSIGMA - ABS(LSIGMA)* FOUR * EPS
274: RSIGMA = RSIGMA + ABS(RSIGMA)* FOUR * EPS
275:
276: * Compute upper bounds for how much to back off the initial shifts
277: LDMAX = QUART * MINGAP + TWO * PIVMIN
278: RDMAX = QUART * MINGAP + TWO * PIVMIN
279:
280: LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT
281: RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT
282: *
283: * Initialize the record of the best representation found
284: *
285: S = DLAMCH( 'S' )
286: SMLGROWTH = ONE / S
287: FAIL = DBLE(N-1)*MINGAP/(SPDIAM*EPS)
288: FAIL2 = DBLE(N-1)*MINGAP/(SPDIAM*SQRT(EPS))
289: BESTSHIFT = LSIGMA
290: *
291: * while (KTRY <= KTRYMAX)
292: KTRY = 0
293: GROWTHBOUND = MAXGROWTH1*SPDIAM
294:
295: 5 CONTINUE
296: SAWNAN1 = .FALSE.
297: SAWNAN2 = .FALSE.
298: * Ensure that we do not back off too much of the initial shifts
299: LDELTA = MIN(LDMAX,LDELTA)
300: RDELTA = MIN(RDMAX,RDELTA)
301:
302: * Compute the element growth when shifting to both ends of the cluster
303: * accept the shift if there is no element growth at one of the two ends
304:
305: * Left end
306: S = -LSIGMA
307: DPLUS( 1 ) = D( 1 ) + S
308: IF(ABS(DPLUS(1)).LT.PIVMIN) THEN
309: DPLUS(1) = -PIVMIN
310: * Need to set SAWNAN1 because refined RRR test should not be used
311: * in this case
312: SAWNAN1 = .TRUE.
313: ENDIF
314: MAX1 = ABS( DPLUS( 1 ) )
315: DO 6 I = 1, N - 1
316: LPLUS( I ) = LD( I ) / DPLUS( I )
317: S = S*LPLUS( I )*L( I ) - LSIGMA
318: DPLUS( I+1 ) = D( I+1 ) + S
319: IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN
320: DPLUS(I+1) = -PIVMIN
321: * Need to set SAWNAN1 because refined RRR test should not be used
322: * in this case
323: SAWNAN1 = .TRUE.
324: ENDIF
325: MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) )
326: 6 CONTINUE
327: SAWNAN1 = SAWNAN1 .OR. DISNAN( MAX1 )
328:
329: IF( FORCER .OR.
330: $ (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN
331: SIGMA = LSIGMA
332: SHIFT = SLEFT
333: GOTO 100
334: ENDIF
335:
336: * Right end
337: S = -RSIGMA
338: WORK( 1 ) = D( 1 ) + S
339: IF(ABS(WORK(1)).LT.PIVMIN) THEN
340: WORK(1) = -PIVMIN
341: * Need to set SAWNAN2 because refined RRR test should not be used
342: * in this case
343: SAWNAN2 = .TRUE.
344: ENDIF
345: MAX2 = ABS( WORK( 1 ) )
346: DO 7 I = 1, N - 1
347: WORK( N+I ) = LD( I ) / WORK( I )
348: S = S*WORK( N+I )*L( I ) - RSIGMA
349: WORK( I+1 ) = D( I+1 ) + S
350: IF(ABS(WORK(I+1)).LT.PIVMIN) THEN
351: WORK(I+1) = -PIVMIN
352: * Need to set SAWNAN2 because refined RRR test should not be used
353: * in this case
354: SAWNAN2 = .TRUE.
355: ENDIF
356: MAX2 = MAX( MAX2,ABS(WORK(I+1)) )
357: 7 CONTINUE
358: SAWNAN2 = SAWNAN2 .OR. DISNAN( MAX2 )
359:
360: IF( FORCER .OR.
361: $ (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN
362: SIGMA = RSIGMA
363: SHIFT = SRIGHT
364: GOTO 100
365: ENDIF
366: * If we are at this point, both shifts led to too much element growth
367:
368: * Record the better of the two shifts (provided it didn't lead to NaN)
369: IF(SAWNAN1.AND.SAWNAN2) THEN
370: * both MAX1 and MAX2 are NaN
371: GOTO 50
372: ELSE
373: IF( .NOT.SAWNAN1 ) THEN
374: INDX = 1
375: IF(MAX1.LE.SMLGROWTH) THEN
376: SMLGROWTH = MAX1
377: BESTSHIFT = LSIGMA
378: ENDIF
379: ENDIF
380: IF( .NOT.SAWNAN2 ) THEN
381: IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2
382: IF(MAX2.LE.SMLGROWTH) THEN
383: SMLGROWTH = MAX2
384: BESTSHIFT = RSIGMA
385: ENDIF
386: ENDIF
387: ENDIF
388:
389: * If we are here, both the left and the right shift led to
390: * element growth. If the element growth is moderate, then
391: * we may still accept the representation, if it passes a
392: * refined test for RRR. This test supposes that no NaN occurred.
393: * Moreover, we use the refined RRR test only for isolated clusters.
394: IF((CLWDTH.LT.MINGAP/DBLE(128)) .AND.
395: $ (MIN(MAX1,MAX2).LT.FAIL2)
396: $ .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN
397: DORRR1 = .TRUE.
398: ELSE
399: DORRR1 = .FALSE.
400: ENDIF
401: TRYRRR1 = .TRUE.
402: IF( TRYRRR1 .AND. DORRR1 ) THEN
403: IF(INDX.EQ.1) THEN
404: TMP = ABS( DPLUS( N ) )
405: ZNM2 = ONE
406: PROD = ONE
407: OLDP = ONE
408: DO 15 I = N-1, 1, -1
409: IF( PROD .LE. EPS ) THEN
410: PROD =
411: $ ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP
412: ELSE
413: PROD = PROD*ABS(WORK(N+I))
414: END IF
415: OLDP = PROD
416: ZNM2 = ZNM2 + PROD**2
417: TMP = MAX( TMP, ABS( DPLUS( I ) * PROD ))
418: 15 CONTINUE
419: RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) )
420: IF (RRR1.LE.MAXGROWTH2) THEN
421: SIGMA = LSIGMA
422: SHIFT = SLEFT
423: GOTO 100
424: ENDIF
425: ELSE IF(INDX.EQ.2) THEN
426: TMP = ABS( WORK( N ) )
427: ZNM2 = ONE
428: PROD = ONE
429: OLDP = ONE
430: DO 16 I = N-1, 1, -1
431: IF( PROD .LE. EPS ) THEN
432: PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP
433: ELSE
434: PROD = PROD*ABS(LPLUS(I))
435: END IF
436: OLDP = PROD
437: ZNM2 = ZNM2 + PROD**2
438: TMP = MAX( TMP, ABS( WORK( I ) * PROD ))
439: 16 CONTINUE
440: RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) )
441: IF (RRR2.LE.MAXGROWTH2) THEN
442: SIGMA = RSIGMA
443: SHIFT = SRIGHT
444: GOTO 100
445: ENDIF
446: END IF
447: ENDIF
448:
449: 50 CONTINUE
450:
451: IF (KTRY.LT.KTRYMAX) THEN
452: * If we are here, both shifts failed also the RRR test.
453: * Back off to the outside
454: LSIGMA = MAX( LSIGMA - LDELTA,
455: $ LSIGMA - LDMAX)
456: RSIGMA = MIN( RSIGMA + RDELTA,
457: $ RSIGMA + RDMAX )
458: LDELTA = TWO * LDELTA
459: RDELTA = TWO * RDELTA
460: KTRY = KTRY + 1
461: GOTO 5
462: ELSE
463: * None of the representations investigated satisfied our
464: * criteria. Take the best one we found.
465: IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN
466: LSIGMA = BESTSHIFT
467: RSIGMA = BESTSHIFT
468: FORCER = .TRUE.
469: GOTO 5
470: ELSE
471: INFO = 1
472: RETURN
473: ENDIF
474: END IF
475:
476: 100 CONTINUE
477: IF (SHIFT.EQ.SLEFT) THEN
478: ELSEIF (SHIFT.EQ.SRIGHT) THEN
479: * store new L and D back into DPLUS, LPLUS
480: CALL DCOPY( N, WORK, 1, DPLUS, 1 )
481: CALL DCOPY( N-1, WORK(N+1), 1, LPLUS, 1 )
482: ENDIF
483:
484: RETURN
485: *
486: * End of DLARRF
487: *
488: 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