```      SUBROUTINE SLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
*
*  -- LAPACK routine (version 3.1.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     February 2007
*
*     .. Scalar Arguments ..
LOGICAL            ORGATI
INTEGER            INFO, KNITER
REAL               FINIT, RHO, TAU
*     ..
*     .. Array Arguments ..
REAL               D( 3 ), Z( 3 )
*     ..
*
*  Purpose
*  =======
*
*  SLAED6 computes the positive or negative root (closest to the origin)
*  of
*                   z(1)        z(2)        z(3)
*  f(x) =   rho + --------- + ---------- + ---------
*                  d(1)-x      d(2)-x      d(3)-x
*
*  It is assumed that
*
*        if ORGATI = .true. the root is between d(2) and d(3);
*        otherwise it is between d(1) and d(2)
*
*  This routine will be called by SLAED4 when necessary. In most cases,
*  the root sought is the smallest in magnitude, though it might not be
*  in some extremely rare situations.
*
*  Arguments
*  =========
*
*  KNITER       (input) INTEGER
*               Refer to SLAED4 for its significance.
*
*  ORGATI       (input) LOGICAL
*               If ORGATI is true, the needed root is between d(2) and
*               d(3); otherwise it is between d(1) and d(2).  See
*               SLAED4 for further details.
*
*  RHO          (input) REAL
*               Refer to the equation f(x) above.
*
*  D            (input) REAL array, dimension (3)
*               D satisfies d(1) < d(2) < d(3).
*
*  Z            (input) REAL array, dimension (3)
*               Each of the elements in z must be positive.
*
*  FINIT        (input) REAL
*               The value of f at 0. It is more accurate than the one
*               evaluated inside this routine (if someone wants to do
*               so).
*
*  TAU          (output) REAL
*               The root of the equation f(x).
*
*  INFO         (output) INTEGER
*               = 0: successful exit
*               > 0: if INFO = 1, failure to converge
*
*  Further Details
*  ===============
*
*  30/06/99: Based on contributions by
*     Ren-Cang Li, Computer Science Division, University of California
*     at Berkeley, USA
*
*  10/02/03: This version has a few statements commented out for thread safety
*     (machine parameters are computed on each entry). SJH.
*
*  05/10/06: Modified from a new version of Ren-Cang Li, use
*     Gragg-Thornton-Warner cubic convergent scheme for better stability.
*
*  =====================================================================
*
*     .. Parameters ..
INTEGER            MAXIT
PARAMETER          ( MAXIT = 40 )
REAL               ZERO, ONE, TWO, THREE, FOUR, EIGHT
PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
\$                   THREE = 3.0E0, FOUR = 4.0E0, EIGHT = 8.0E0 )
*     ..
*     .. External Functions ..
REAL               SLAMCH
EXTERNAL           SLAMCH
*     ..
*     .. Local Arrays ..
REAL               DSCALE( 3 ), ZSCALE( 3 )
*     ..
*     .. Local Scalars ..
LOGICAL            SCALE
INTEGER            I, ITER, NITER
REAL               A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F,
\$                   FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1,
\$                   SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4,
\$                   LBD, UBD
*     ..
*     .. Intrinsic Functions ..
INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT
*     ..
*     .. Executable Statements ..
*
INFO = 0
*
IF( ORGATI ) THEN
LBD = D(2)
UBD = D(3)
ELSE
LBD = D(1)
UBD = D(2)
END IF
IF( FINIT .LT. ZERO )THEN
LBD = ZERO
ELSE
UBD = ZERO
END IF
*
NITER = 1
TAU = ZERO
IF( KNITER.EQ.2 ) THEN
IF( ORGATI ) THEN
TEMP = ( D( 3 )-D( 2 ) ) / TWO
C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP )
A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 )
B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 )
ELSE
TEMP = ( D( 1 )-D( 2 ) ) / TWO
C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP )
A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 )
B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 )
END IF
TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
A = A / TEMP
B = B / TEMP
C = C / TEMP
IF( C.EQ.ZERO ) THEN
TAU = B / A
ELSE IF( A.LE.ZERO ) THEN
TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
ELSE
TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
END IF
IF( TAU .LT. LBD .OR. TAU .GT. UBD )
\$      TAU = ( LBD+UBD )/TWO
IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN
TAU = ZERO
ELSE
TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) +
\$                     TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) +
\$                     TAU*Z(3)/( D(3)*( D( 3 )-TAU ) )
IF( TEMP .LE. ZERO )THEN
LBD = TAU
ELSE
UBD = TAU
END IF
IF( ABS( FINIT ).LE.ABS( TEMP ) )
\$         TAU = ZERO
END IF
END IF
*
*     get machine parameters for possible scaling to avoid overflow
*
*     modified by Sven: parameters SMALL1, SMINV1, SMALL2,
*     SMINV2, EPS are not SAVEd anymore between one call to the
*     others but recomputed at each call
*
EPS = SLAMCH( 'Epsilon' )
BASE = SLAMCH( 'Base' )
SMALL1 = BASE**( INT( LOG( SLAMCH( 'SafMin' ) ) / LOG( BASE ) /
\$         THREE ) )
SMINV1 = ONE / SMALL1
SMALL2 = SMALL1*SMALL1
SMINV2 = SMINV1*SMINV1
*
*     Determine if scaling of inputs necessary to avoid overflow
*     when computing 1/TEMP**3
*
IF( ORGATI ) THEN
TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) )
ELSE
TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) )
END IF
SCALE = .FALSE.
IF( TEMP.LE.SMALL1 ) THEN
SCALE = .TRUE.
IF( TEMP.LE.SMALL2 ) THEN
*
*        Scale up by power of radix nearest 1/SAFMIN**(2/3)
*
SCLFAC = SMINV2
SCLINV = SMALL2
ELSE
*
*        Scale up by power of radix nearest 1/SAFMIN**(1/3)
*
SCLFAC = SMINV1
SCLINV = SMALL1
END IF
*
*        Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
*
DO 10 I = 1, 3
DSCALE( I ) = D( I )*SCLFAC
ZSCALE( I ) = Z( I )*SCLFAC
10    CONTINUE
TAU = TAU*SCLFAC
LBD = LBD*SCLFAC
UBD = UBD*SCLFAC
ELSE
*
*        Copy D and Z to DSCALE and ZSCALE
*
DO 20 I = 1, 3
DSCALE( I ) = D( I )
ZSCALE( I ) = Z( I )
20    CONTINUE
END IF
*
FC = ZERO
DF = ZERO
DDF = ZERO
DO 30 I = 1, 3
TEMP = ONE / ( DSCALE( I )-TAU )
TEMP1 = ZSCALE( I )*TEMP
TEMP2 = TEMP1*TEMP
TEMP3 = TEMP2*TEMP
FC = FC + TEMP1 / DSCALE( I )
DF = DF + TEMP2
DDF = DDF + TEMP3
30 CONTINUE
F = FINIT + TAU*FC
*
IF( ABS( F ).LE.ZERO )
\$   GO TO 60
IF( F .LE. ZERO )THEN
LBD = TAU
ELSE
UBD = TAU
END IF
*
*        Iteration begins -- Use Gragg-Thornton-Warner cubic convergent
*                            scheme
*
*     It is not hard to see that
*
*           1) Iterations will go up monotonically
*              if FINIT < 0;
*
*           2) Iterations will go down monotonically
*              if FINIT > 0.
*
ITER = NITER + 1
*
DO 50 NITER = ITER, MAXIT
*
IF( ORGATI ) THEN
TEMP1 = DSCALE( 2 ) - TAU
TEMP2 = DSCALE( 3 ) - TAU
ELSE
TEMP1 = DSCALE( 1 ) - TAU
TEMP2 = DSCALE( 2 ) - TAU
END IF
A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF
B = TEMP1*TEMP2*F
C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF
TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
A = A / TEMP
B = B / TEMP
C = C / TEMP
IF( C.EQ.ZERO ) THEN
ETA = B / A
ELSE IF( A.LE.ZERO ) THEN
ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
ELSE
ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
END IF
IF( F*ETA.GE.ZERO ) THEN
ETA = -F / DF
END IF
*
TAU = TAU + ETA
IF( TAU .LT. LBD .OR. TAU .GT. UBD )
\$      TAU = ( LBD + UBD )/TWO
*
FC = ZERO
ERRETM = ZERO
DF = ZERO
DDF = ZERO
DO 40 I = 1, 3
TEMP = ONE / ( DSCALE( I )-TAU )
TEMP1 = ZSCALE( I )*TEMP
TEMP2 = TEMP1*TEMP
TEMP3 = TEMP2*TEMP
TEMP4 = TEMP1 / DSCALE( I )
FC = FC + TEMP4
ERRETM = ERRETM + ABS( TEMP4 )
DF = DF + TEMP2
DDF = DDF + TEMP3
40    CONTINUE
F = FINIT + TAU*FC
ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +
\$            ABS( TAU )*DF
IF( ABS( F ).LE.EPS*ERRETM )
\$      GO TO 60
IF( F .LE. ZERO )THEN
LBD = TAU
ELSE
UBD = TAU
END IF
50 CONTINUE
INFO = 1
60 CONTINUE
*
*     Undo scaling
*
IF( SCALE )
\$   TAU = TAU*SCLINV
RETURN
*
*     End of SLAED6
*
END

```