dlaed6()

subroutine dlaed6	(	integer	kniter,
		logical	orgati,
		double precision	rho,
		double precision, dimension( 3 )	d,
		double precision, dimension( 3 )	z,
		double precision	finit,
		double precision	tau,
		integer	info
	)

DLAED6 used by DSTEDC. Computes one Newton step in solution of the secular equation.

Download DLAED6 + dependencies [TGZ] [ZIP] [TXT]

Purpose:

 DLAED6 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 DLAED4 when necessary. In most cases,
 the root sought is the smallest in magnitude, though it might not be
 in some extremely rare situations.

Parameters

[in]	KNITER	KNITER is INTEGER Refer to DLAED4 for its significance.
[in]	ORGATI	ORGATI is 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 DLAED4 for further details.
[in]	RHO	RHO is DOUBLE PRECISION Refer to the equation f(x) above.
[in]	D	D is DOUBLE PRECISION array, dimension (3) D satisfies d(1) < d(2) < d(3).
[in]	Z	Z is DOUBLE PRECISION array, dimension (3) Each of the elements in z must be positive.
[in]	FINIT	FINIT is DOUBLE PRECISION The value of f at 0. It is more accurate than the one evaluated inside this routine (if someone wants to do so).
[out]	TAU	TAU is DOUBLE PRECISION The root of the equation f(x).
[out]	INFO	INFO is INTEGER = 0: successful exit > 0: if INFO = 1, failure to converge

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

  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.

Contributors:

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

Definition at line 139 of file dlaed6.f.

*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      LOGICAL            ORGATI
      INTEGER            INFO, KNITER
      DOUBLE PRECISION   FINIT, RHO, TAU
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   D( 3 ), Z( 3 )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            MAXIT
      parameter( maxit = 40 )
      DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, EIGHT
      parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
     $                   three = 3.0d0, four = 4.0d0, eight = 8.0d0 )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           dlamch
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   DSCALE( 3 ), ZSCALE( 3 )
*     ..
*     .. Local Scalars ..
      LOGICAL            SCALE
      INTEGER            I, ITER, NITER
      DOUBLE PRECISION   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 = dlamch( 'Epsilon' )
      base = dlamch( 'Base' )
      small1 = base**( int( log( dlamch( '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
            IF ( ( dscale( i )-tau ).NE.zero ) THEN
               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
            ELSE
               GO TO 60
            END IF
   40    CONTINUE
         f = finit + tau*fc
         erretm = eight*( abs( finit )+abs( tau )*erretm ) +
     $            abs( tau )*df
         IF( ( abs( f ).LE.four*eps*erretm ) .OR.
     $      ( (ubd-lbd).LE.four*eps*abs(tau) )  )
     $      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 DLAED6
*

Here is the caller graph for this function: