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 137 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: