dlaed6.f

Go to the documentation of this file.

*> \brief \b DLAED6 used by sstedc. Computes one Newton step in solution of the secular equation.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download DLAED6 + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed6.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed6.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed6.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
* 
*       .. Scalar Arguments ..
*       LOGICAL            ORGATI
*       INTEGER            INFO, KNITER
*       DOUBLE PRECISION   FINIT, RHO, TAU
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   D( 3 ), Z( 3 )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] KNITER
*> \verbatim
*>          KNITER is INTEGER
*>               Refer to DLAED4 for its significance.
*> \endverbatim
*>
*> \param[in] ORGATI
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] RHO
*> \verbatim
*>          RHO is DOUBLE PRECISION
*>               Refer to the equation f(x) above.
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (3)
*>               D satisfies d(1) < d(2) < d(3).
*> \endverbatim
*>
*> \param[in] Z
*> \verbatim
*>          Z is DOUBLE PRECISION array, dimension (3)
*>               Each of the elements in z must be positive.
*> \endverbatim
*>
*> \param[in] FINIT
*> \verbatim
*>          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).
*> \endverbatim
*>
*> \param[out] TAU
*> \verbatim
*>          TAU is DOUBLE PRECISION
*>               The root of the equation f(x).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>               = 0: successful exit
*>               > 0: if INFO = 1, failure to converge
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date September 2012
*
*> \ingroup auxOTHERcomputational
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  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.
*> \endverbatim
*
*> \par Contributors:
*  ==================
*>
*>     Ren-Cang Li, Computer Science Division, University of California
*>     at Berkeley, USA
*>
*  =====================================================================
      SUBROUTINE dlaed6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
*
*  -- LAPACK computational routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. 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.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 DLAED6
*
      END