subroutine aptnewt (a, np, xinit, xmin, xmax, nitmax, tol,
     &                    xroot, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTNEWT
c
c     call aptnewt (a, np, xinit, xmin, xmax, nitmax, tol, xroot, nerr)
c
c     Version:  aptnewt  Updated    1997 September 3 13:20.
c               aptnewt  Originated 1997 July 30 15:40.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To find a root xroot of the polynomial equation in x,
c               P(x) = SUM {a(n)*x**n}, (n = 0, np), with an initial guess of
c               xinit, a search range from xmin to xmax, allowing no more then
c               nitmax Newtonian iterations, and a relative precision of tol.
c               Flag nerr, if non-zero, indicates an input error, or a failure
c               to find a root.
c
c     Method:   A modified Newtonian iteration method is used, with an initial
c               x of xinit.  The n'th iterate of x is found as follows: 
c               x(n) = x(n-1) - P{x(n-1)} / P'{x(n-1)}, where P'(x) is the
c               derivative of P(x) with respect to x.  A root has been found
c               when the relative error in calculating P(x(n)) is no more than
c               tol.  A maximum of nitmax iterations are allowed.  If P'(x) is
c               zero, or if x gets stuck at xmin or xmax, a new trial value of
c               x is randomly selected between xmin and xmax.
c
c     Note:     Roots at a minimum or maximum of P(x) can not be found by this
c               method.  Try finding the roots of P'(x), and seeing if P(x) is
c               zero at any of those roots.  Repeat for higher derivatives,
c               if necessary.
c               If an interval (xmin, xmax) has been found for which P(xmin)
c               and P(xmax) have opposite signs, a good initial guess xinit is
c               (xmin * P(xmax) - xmax * P(xmin)) / (P(xmax) - P(xmin)).
c               Up to np roots may exist between xmin and xmax.
c
c     Input:    a, np, xinit, xmin, xmax, nitmax, tol.
c
c     Output:   xroot, nerr
c
c     Glossary:
c
c     a         Input    The array of coefficients of the polynomial equation
c                          in x of degree np:
c                          P(x) = SUM {a(n) * x**n} (n = 0, np).
c                          Must be dimensioned a(0:np) or bigger.
c
c     nitmax    Input    The maximum number of iterations to be tried.
c                          Should be at least -3.32 log (tol), or about 50
c                          for 64-bit floating point numbers.
c
c     nerr      Output   Indicates no input errors were found, and a root was
c                          found, if zero.  Otherwise:
c                           1  np was not positive.
c                           2  xmin was not less than xmax.
c                           3  xinit was not between xmin and xmax.
c                           4  nitmax was not positive.
c                           5  no convergence after nitmax iterations.
c
c     np        Input    The largest exponent of x in the polynomial equation
c                          P(x), for which the coefficient is non-zero.
c                          Must be positive.  If 4, use aptquar.  If 3, use
c                          aptcubs.  If 2, use aptqrts or aptqrtv.
c                          If 1, xroot = -a(0) / a(1).
c
c     tol       Input    Numerical iteration convergence control.
c                          If zero, a value of 1.e-11 is used instead.
c                          If using 64-bit floating point numbers, a value
c                          between 1.e-11 and 1.e-14 is recommended.
c
c     xinit     Input    The initial value of x to try.  Must be between
c                          xmin and xmax.
c
c     xmax      Input    The maximum allowed value of x to try.
c                          Must be more than xmin.
c
c     xmin      Input    The minimum allowed value of x to try.
c                          Must be less than xmax.
c
c     xroot     Output   A root of the polynomial P(x), for which the following
c                          conditions are satisfied:
c                          The initial guess is xinit;
c                          the root is between xmin and xmax;
c                          no more than nitmax iterations are required;
c                          the relative error in calculating P(xroot) is no
c                          more than tol.
c                          If no root was found, -1.e99 is returned.
c                          Up to np roots may exist between xmin and xmax.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

      dimension a (0:1)               ! The coefficients of x**n.

cbugc***DEBUG begins.
cbug 9901 format (/ 'aptnewt using Newton''s iteration to find a root',
cbug     &  ' of a polynomial.')
cbug 9902 format ('tol=',1pe22.14,'  np=',i5,'  nitmax=',i6)
cbug 9903 format ('xinit,min,max=',1p3e22.14)
cbug 9904 format (i5,2x,1pe22.14)
cbug      write ( 3, 9901)
cbug      write ( 3, 9902) tol, np, nitmax
cbug      write ( 3, 9903) xinit, xmin, xmax
cbug      write ( 3, 9904) (n, a(n), n = 0, np)
cbugc***DEBUG ends.

c.... Initialize.

      nerr  = 0
      xroot = -1.e99
      niter = 0
      tolx  = tol
      if (tolx .eq. 0.0) tolx = 1.e-11

c.... Test for input errors.

      if (np .le. 0) then
        nerr = 1
        go to 210
      endif

      xerr = 2.0 * tolx * (abs (xmax) + abs (xmin))
      if ((xmin - xmax) .gt. xerr) then
        nerr = 2
        go to 210
      endif

      xerr1 = 2.0 * tolx * (abs (xinit) + abs (xmin))
      xerr2 = 2.0 * tolx * (abs (xinit) + abs (xmax))
      if (((xmin - xinit) .gt. xerr1)  .or.
     &    ((xinit - xmax) .gt. xerr2)) then
        nerr = 3
        go to 210
      endif

      if (nitmax .le. 0) then
        nerr = 4
        go to 210
      endif

c.... Initialize iterate.

      xold = xinit
      xnew = xinit

c.... Iterate using Newtonian iteration.

  110 niter = niter + 1

c.... Find value of polynomial.

      poly  = a(0)
      perr  = 0.0                     ! Error in polynomial values.
      xperr = 0.0
      xpow  = 1.0

      do 120 n = 1, np
        xperr = abs (xold) * (xperr + tolx * abs (xpow))
        xpow  = xpow * xold
        perr  = perr + abs (a(n)) * (xperr + tolx * abs (xpow))
        poly  = poly + a(n) * xpow
  120 continue

cbugcbugc***DEBUG begins.
cbugcbug 9911 format ('n=',i5,'  x,p,e=',1p3e20.12)
cbugcbug      write ( 3, 9911) niter, xold, poly, perr
cbugcbugc***DEBUG ends.

c.... Use iterate if polynomial value is less than allowed error.

      if (abs (poly) .le. perr) then
        xroot = xold
cbugcbugc***DEBUG begins.
cbugcbug 9912 format ('CONVERGED.')
cbugcbug      write ( 3, 9912)
cbugcbugc***DEBUG ends.
        go to 210
      endif

c.... Find value of slope of polynomial.

      slope = a(1)
      serr  = 0.0                     ! Error in slope value.
      xperr = 0.0
      xpow  = 1.0

      do 130 n = 2, np
        xperr = abs (xold) * (xperr + tolx * abs (xpow))
        xpow  = xpow * xold
        serr  = serr + n * abs (a(n)) * (xperr + tolx * abs (xpow))
        slope = slope + n * a(n) * xpow
  130 continue

c.... Pick an x at random in interval if slope is zero.

      if (abs (slope) .le. serr) then
cbugcbugc***DEBUG begins.
cbugcbug 9913 format ('PICKING NEW XINT.           SLOPE =',1pe22.14)
cbugcbug      write ( 3, 9913) slope
cbugcbugc***DEBUG ends.
        xnew = xmin + ranf () * (xmax - xmin)
        go to 200
      endif

c.... Find next iterate.

      xnew = xold - poly / slope
      if (xnew .lt. xmin) xnew = xmin
      if (xnew .gt. xmax) xnew = xmax

c.... Pick an x at random in interval if stuck at xmin or xmax.

      if ((xold .eq. xmin) .and. (xnew .eq. xmin)) then
cbugcbugc***DEBUG begins.
cbugcbug 9914 format ('STUCK AT XMIN.              SLOPE =',1pe22.14)
cbugcbug      write ( 3, 9914) slope
cbugcbugc***DEBUG ends.
        xnew = xmin + ranf () * (xmax - xmin)
      endif

      if ((xold .eq. xmax) .and. (xnew .eq. xmax)) then
cbugcbugc***DEBUG begins.
cbugcbug 9915 format ('STUCK AT XMAX.              SLOPE =',1pe22.14)
cbugcbug      write ( 3, 9915) slope
cbugcbugc***DEBUG ends.
        xnew = xmin + ranf () * (xmax - xmin)
      endif

 200  xold = xnew

      if (niter .le. nitmax) go to 110

cbugcbugc***DEBUG begins.
cbugcbug 9916 format ('TOO MANY ITERATIONS.')
cbugcbug      write ( 3, 9916)
cbugcbugc***DEBUG ends.
      nerr = 5

  210 continue
cbugc***DEBUG begins.
cbug 9940 format (/ 'aptnewt results.  nerr=',i3,'  nit=',i5,
cbug     &  '  xroot=',1pe22.14)
cbug      write ( 3, 9940) nerr, niter, xroot
cbugc***DEBUG ends.

      return

c.... End of subroutine aptnewt.      (+1 line.)
      end

UCRL-WEB-209832