```      SUBROUTINE SLARRF( N, D, L, LD, CLSTRT, CLEND,
\$                   W, WGAP, WERR,
\$                   SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
\$                   DPLUS, LPLUS, WORK, INFO )
*
*  -- LAPACK auxiliary routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
**
*     .. Scalar Arguments ..
INTEGER            CLSTRT, CLEND, INFO, N
REAL               CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
*     ..
*     .. Array Arguments ..
REAL               D( * ), DPLUS( * ), L( * ), LD( * ),
\$          LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  Given the initial representation L D L^T and its cluster of close
*  eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
*  W( CLEND ), SLARRF finds a new relatively robust representation
*  L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
*  eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix (subblock, if the matrix splitted).
*
*  D       (input) REAL             array, dimension (N)
*          The N diagonal elements of the diagonal matrix D.
*
*  L       (input) REAL             array, dimension (N-1)
*          The (N-1) subdiagonal elements of the unit bidiagonal
*          matrix L.
*
*  LD      (input) REAL             array, dimension (N-1)
*          The (N-1) elements L(i)*D(i).
*
*  CLSTRT  (input) INTEGER
*          The index of the first eigenvalue in the cluster.
*
*  CLEND   (input) INTEGER
*          The index of the last eigenvalue in the cluster.
*
*  W       (input) REAL             array, dimension >=  (CLEND-CLSTRT+1)
*          The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
*          W( CLSTRT ) through W( CLEND ) form the cluster of relatively
*          close eigenalues.
*
*  WGAP    (input/output) REAL             array, dimension >=  (CLEND-CLSTRT+1)
*          The separation from the right neighbor eigenvalue in W.
*
*  WERR    (input) REAL             array, dimension >=  (CLEND-CLSTRT+1)
*          WERR contain the semiwidth of the uncertainty
*          interval of the corresponding eigenvalue APPROXIMATION in W
*
*  SPDIAM (input) estimate of the spectral diameter obtained from the
*          Gerschgorin intervals
*
*  CLGAPL, CLGAPR (input) absolute gap on each end of the cluster.
*          Set by the calling routine to protect against shifts too close
*          to eigenvalues outside the cluster.
*
*  PIVMIN  (input) DOUBLE PRECISION
*          The minimum pivot allowed in the Sturm sequence.
*
*  SIGMA   (output) REAL
*          The shift used to form L(+) D(+) L(+)^T.
*
*  DPLUS   (output) REAL             array, dimension (N)
*          The N diagonal elements of the diagonal matrix D(+).
*
*  LPLUS   (output) REAL             array, dimension (N-1)
*          The first (N-1) elements of LPLUS contain the subdiagonal
*          elements of the unit bidiagonal matrix L(+).
*
*  WORK    (workspace) REAL             array, dimension (2*N)
*          Workspace.
*
*  Further Details
*  ===============
*
*  Based on contributions by
*     Beresford Parlett, University of California, Berkeley, USA
*     Jim Demmel, University of California, Berkeley, USA
*     Inderjit Dhillon, University of Texas, Austin, USA
*     Osni Marques, LBNL/NERSC, USA
*     Christof Voemel, University of California, Berkeley, USA
*
*  =====================================================================
*
*     .. Parameters ..
REAL               FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO,
\$                   ZERO
PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
\$                     FOUR = 4.0E0, QUART = 0.25E0,
\$                     MAXGROWTH1 = 8.E0,
\$                     MAXGROWTH2 = 8.E0 )
*     ..
*     .. Local Scalars ..
LOGICAL   DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1
INTEGER            I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT
PARAMETER          ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 )
REAL               AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL,
\$                   FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA,
\$                   MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX,
\$                   RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2
*     ..
*     .. External Functions ..
LOGICAL SISNAN
REAL               SLAMCH
EXTERNAL           SISNAN, SLAMCH
*     ..
*     .. External Subroutines ..
EXTERNAL           SCOPY
*     ..
*     .. Intrinsic Functions ..
INTRINSIC          ABS
*     ..
*     .. Executable Statements ..
*
INFO = 0
FACT = REAL(2**KTRYMAX)
EPS = SLAMCH( 'Precision' )
SHIFT = 0
FORCER = .FALSE.

*     Note that we cannot guarantee that for any of the shifts tried,
*     the factorization has a small or even moderate element growth.
*     There could be Ritz values at both ends of the cluster and despite
*     backing off, there are examples where all factorizations tried
*     (in IEEE mode, allowing zero pivots & infinities) have INFINITE
*     element growth.
*     For this reason, we should use PIVMIN in this subroutine so that at
*     least the L D L^T factorization exists. It can be checked afterwards
*     whether the element growth caused bad residuals/orthogonality.

*     Decide whether the code should accept the best among all
*     representations despite large element growth or signal INFO=1
NOFAIL = .TRUE.
*

*     Compute the average gap length of the cluster
CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT)
AVGAP = CLWDTH / REAL(CLEND-CLSTRT)
MINGAP = MIN(CLGAPL, CLGAPR)
*     Initial values for shifts to both ends of cluster
LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT )
RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND )

*     Use a small fudge to make sure that we really shift to the outside
LSIGMA = LSIGMA - ABS(LSIGMA)* TWO * EPS
RSIGMA = RSIGMA + ABS(RSIGMA)* TWO * EPS

*     Compute upper bounds for how much to back off the initial shifts
LDMAX = QUART * MINGAP + TWO * PIVMIN
RDMAX = QUART * MINGAP + TWO * PIVMIN

LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT
RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT
*
*     Initialize the record of the best representation found
*
S = SLAMCH( 'S' )
SMLGROWTH = ONE / S
FAIL = REAL(N-1)*MINGAP/(SPDIAM*EPS)
FAIL2 = REAL(N-1)*MINGAP/(SPDIAM*SQRT(EPS))
BESTSHIFT = LSIGMA
*
*     while (KTRY <= KTRYMAX)
KTRY = 0
GROWTHBOUND = MAXGROWTH1*SPDIAM

5    CONTINUE
SAWNAN1 = .FALSE.
SAWNAN2 = .FALSE.
*     Ensure that we do not back off too much of the initial shifts
LDELTA = MIN(LDMAX,LDELTA)
RDELTA = MIN(RDMAX,RDELTA)

*     Compute the element growth when shifting to both ends of the cluster
*     accept the shift if there is no element growth at one of the two ends

*     Left end
S = -LSIGMA
DPLUS( 1 ) = D( 1 ) + S
IF(ABS(DPLUS(1)).LT.PIVMIN) THEN
DPLUS(1) = -PIVMIN
*        Need to set SAWNAN1 because refined RRR test should not be used
*        in this case
SAWNAN1 = .TRUE.
ENDIF
MAX1 = ABS( DPLUS( 1 ) )
DO 6 I = 1, N - 1
LPLUS( I ) = LD( I ) / DPLUS( I )
S = S*LPLUS( I )*L( I ) - LSIGMA
DPLUS( I+1 ) = D( I+1 ) + S
IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN
DPLUS(I+1) = -PIVMIN
*           Need to set SAWNAN1 because refined RRR test should not be used
*           in this case
SAWNAN1 = .TRUE.
ENDIF
MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) )
6    CONTINUE
SAWNAN1 = SAWNAN1 .OR.  SISNAN( MAX1 )

IF( FORCER .OR.
\$   (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN
SIGMA = LSIGMA
SHIFT = SLEFT
GOTO 100
ENDIF

*     Right end
S = -RSIGMA
WORK( 1 ) = D( 1 ) + S
IF(ABS(WORK(1)).LT.PIVMIN) THEN
WORK(1) = -PIVMIN
*        Need to set SAWNAN2 because refined RRR test should not be used
*        in this case
SAWNAN2 = .TRUE.
ENDIF
MAX2 = ABS( WORK( 1 ) )
DO 7 I = 1, N - 1
WORK( N+I ) = LD( I ) / WORK( I )
S = S*WORK( N+I )*L( I ) - RSIGMA
WORK( I+1 ) = D( I+1 ) + S
IF(ABS(WORK(I+1)).LT.PIVMIN) THEN
WORK(I+1) = -PIVMIN
*           Need to set SAWNAN2 because refined RRR test should not be used
*           in this case
SAWNAN2 = .TRUE.
ENDIF
MAX2 = MAX( MAX2,ABS(WORK(I+1)) )
7    CONTINUE
SAWNAN2 = SAWNAN2 .OR.  SISNAN( MAX2 )

IF( FORCER .OR.
\$   (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN
SIGMA = RSIGMA
SHIFT = SRIGHT
GOTO 100
ENDIF
*     If we are at this point, both shifts led to too much element growth

*     Record the better of the two shifts (provided it didn't lead to NaN)
IF(SAWNAN1.AND.SAWNAN2) THEN
*        both MAX1 and MAX2 are NaN
GOTO 50
ELSE
IF( .NOT.SAWNAN1 ) THEN
INDX = 1
IF(MAX1.LE.SMLGROWTH) THEN
SMLGROWTH = MAX1
BESTSHIFT = LSIGMA
ENDIF
ENDIF
IF( .NOT.SAWNAN2 ) THEN
IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2
IF(MAX2.LE.SMLGROWTH) THEN
SMLGROWTH = MAX2
BESTSHIFT = RSIGMA
ENDIF
ENDIF
ENDIF

*     If we are here, both the left and the right shift led to
*     element growth. If the element growth is moderate, then
*     we may still accept the representation, if it passes a
*     refined test for RRR. This test supposes that no NaN occurred.
*     Moreover, we use the refined RRR test only for isolated clusters.
IF((CLWDTH.LT.MINGAP/REAL(128)) .AND.
\$   (MIN(MAX1,MAX2).LT.FAIL2)
\$  .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN
DORRR1 = .TRUE.
ELSE
DORRR1 = .FALSE.
ENDIF
TRYRRR1 = .TRUE.
IF( TRYRRR1 .AND. DORRR1 ) THEN
IF(INDX.EQ.1) THEN
TMP = ABS( DPLUS( N ) )
ZNM2 = ONE
PROD = ONE
OLDP = ONE
DO 15 I = N-1, 1, -1
IF( PROD .LE. EPS ) THEN
PROD =
\$         ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP
ELSE
PROD = PROD*ABS(WORK(N+I))
END IF
OLDP = PROD
ZNM2 = ZNM2 + PROD**2
TMP = MAX( TMP, ABS( DPLUS( I ) * PROD ))
15      CONTINUE
RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) )
IF (RRR1.LE.MAXGROWTH2) THEN
SIGMA = LSIGMA
SHIFT = SLEFT
GOTO 100
ENDIF
ELSE IF(INDX.EQ.2) THEN
TMP = ABS( WORK( N ) )
ZNM2 = ONE
PROD = ONE
OLDP = ONE
DO 16 I = N-1, 1, -1
IF( PROD .LE. EPS ) THEN
PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP
ELSE
PROD = PROD*ABS(LPLUS(I))
END IF
OLDP = PROD
ZNM2 = ZNM2 + PROD**2
TMP = MAX( TMP, ABS( WORK( I ) * PROD ))
16      CONTINUE
RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) )
IF (RRR2.LE.MAXGROWTH2) THEN
SIGMA = RSIGMA
SHIFT = SRIGHT
GOTO 100
ENDIF
END IF
ENDIF

50   CONTINUE

IF (KTRY.LT.KTRYMAX) THEN
*        If we are here, both shifts failed also the RRR test.
*        Back off to the outside
LSIGMA = MAX( LSIGMA - LDELTA,
\$     LSIGMA - LDMAX)
RSIGMA = MIN( RSIGMA + RDELTA,
\$     RSIGMA + RDMAX )
LDELTA = TWO * LDELTA
RDELTA = TWO * RDELTA
KTRY = KTRY + 1
GOTO 5
ELSE
*        None of the representations investigated satisfied our
*        criteria. Take the best one we found.
IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN
LSIGMA = BESTSHIFT
RSIGMA = BESTSHIFT
FORCER = .TRUE.
GOTO 5
ELSE
INFO = 1
RETURN
ENDIF
END IF

100  CONTINUE
IF (SHIFT.EQ.SLEFT) THEN
ELSEIF (SHIFT.EQ.SRIGHT) THEN
*        store new L and D back into DPLUS, LPLUS
CALL SCOPY( N, WORK, 1, DPLUS, 1 )
CALL SCOPY( N-1, WORK(N+1), 1, LPLUS, 1 )
ENDIF

RETURN
*
*     End of SLARRF
*
END

```