subroutine aptpolr (xmin, dx, xmax, nitmax, a, na, nxmax, tol,
& nx, xpoly, atype, pa, pb, pc, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPOLR
c
c call aptpolr (xmin, dx, xmax, nitmax, a, na, nxmax, tol,
c nx, xpoly, atype, pa, pb, pc, nerr)
c
c Version: aptpolr Updated 1997 August 25 15:45.
c aptpolr Originated 1997 August 13 12:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To search for the real roots, extrema and inflection points of
c the polynomial equation with real coefficients a(n), n = 0, na:
c P(x) = sum (n=0,na) {a(n) * x**n}, and the values pa, pb and pc
c of P(x), P'(x) and P''(x), respectively, at those points.
c The search starts with the trial values xmin (dx) xmax.
c The number nx of roots, extrema and inflection points found is
c limited to nxmax. The type of event is returned in atype.
c Flag nerr, if non-zero, indicates an input error.
c
c Method: For each trial value of x, the values of P(x), P'(x) and
c (if na is 3 or more) P''(x), are found. If any are zero, the
c value of x is saved in xpoly. If any changes sign from the
c preceding trial value of x, a Newtonian iteration is used at an
c interpolated value of x, with a maximum of nitmax iterations,
c and convergence determined by tol, to try to find a root, an
c extremum or (if na is 3 or more) an inflection point,
c respectively, in the interval between the preceding and
c following trial values. In addition, if the value of P'(x) is
c zero or changes sign, Newtonian iteration is used to try to find
c a root.
c
c Input: xmin, dx, xmax, nitmax, a, na, nxmax, tol.
c
c Output: nx, xpoly, atype, pa, pb, pc, nerr.
c
c Glossary:
c
c a Input The array of coefficients of the polynomial equation
c in x of degree na:
c P(x) = SUM {a(n) * x**n} (n = 0, na).
c Must be dimensioned a(0:na) or larger.
c Final coefficients of zero will be ignored.
c Each leading coefficient of zero represents a root
c at x = 0. Unless the slopes and inflection points
c of the original polynomial equation are of interest,
c leading coefficients of zero should not be specified.
c
c atype Output A type character*8 array, indicating the type of event
c at point x = xpoly.
c Roots, at which P(x) = 0:
c "single": P(x) = 0, P'(x) /= 0.
c May also be an inflection point.
c "dbl max": P(x) = 0, P'(x) = 0, P''(x) < 0.
c "dbl min": P(x) = 0, P'(x) = 0, P''(x) > 0.
c "triple": P(x) = 0, P'(x) = 0, P''(x) = 0.
c May be of higher order than triple.
c Extrema, at which P'(x) = 0:
c "maximum": P(x) /= 0, P'(x) = 0, P''(x) < 0.
c "minimum": P(x) /= 0, P'(x) = 0, P''(x) > 0.
c "min/max": P(x) /= 0, P'(x) = 0, P''(x) = 0.
c Is also an inflection point.
c Inflection points, at which P''(x) = 0:
c "inflectn": P(x) /= 0, P'(x) /= 0, P''(x) = 0,
c and na is 3 or more.
c Note: "/= 0" means "does not equal zero".
c
c dx Input The increment between successive initial trial values
c of x. Must not be negative, and must not be zero
c unless xmin = xmax.
c
c nerr Output Indicates no input errors were found, if zero.
c Otherwise:
c 1 xmax < xmin.
c 2 xmax = xmin and dx < 0.
c 3 xmax > xmin and dx <= 0.
c 4 nitmax <= 0.
c 5 na < 1.
c 6 nxmax < 3* (na - 1).
c 7 nx tried to exceed nxmax.
c
c na Input The largest exponent of x in the polynomial equation
c P(x). Final coefficients of zero will be ignored.
c Must be positive. Inflection points will be searched
c for if na is 3 or more.
c
c nx Output The total number of values of xpoly found at roots,
c extrema and inflection points. If nx = nxmax and
c nerr = 7, some points preceding the last point were
c not saved.
c
c nxmax Input The maximum number of values of xpoly that may be
c returned (the array sizes of xpoly, atype, pa, pb
c and pc).
c The maximum possible numbers of real roots, extrema
c and inflection points are na, na - 1 and na - 2,
c respectively, for a total of 3 * (na - 1), so nxmax
c must be at least 3 * (na - 1). However, because of
c numerical uncertainty, additional false values may be
c returned near real values, especially near higher
c order roots, where not only P(x) is zero, but one or
c more higher derivatives of P(x) are zero.
c If more than nxmax values of xpoly are found,
c nerr = 7, and only nxmax values of xpoly will be
c returned.
c
c nitmax Input The maximum number of iterations to be used for the
c Newtonian iteration method.
c Should be at least -3.32 log (tol), or about 50
c for 64-bit floating point numbers.
c
c pa Output The value of P(xpoly). Must be size nxmax or larger.
c
c pb Output The value of P'(xpoly). Must be size nxmax or larger.
c
c pc Output The value of P''(xpoly). Must be size nxmax or larger.
c
c tol Input Numerical precision control.
c If using 64-bit floating point numbers, a value
c between 1.e-11 and 1.e-14 is recommended.
c Values of xpoly, pa, pb and pc less than the relative
c error in their calculation, based on tol, are
c truncated to zero.
c
c xpoly Output The values of x at which roots, extrema and inflection
c points were found, arranged in ascending order.
c All will either be zero, or between xmin and xmax.
c The type of event is indicated by type character*8
c array atype. Must be size nxmax or larger.
c Note: real roots must be separated by extrema, and
c extrema must be separated by inflection points or
c min/max points.
c
c xmax Input The final and maximum trial value of x.
c Must not be less than xmin.
c
c xmin Input The initial and minimum trial value of x.
c Must not exceed xmax.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
dimension a(0:1) ! Coefficients of polynomial equation.
dimension xpoly(1) ! Root, extrema, inflection point x.
dimension atype(1) ! Root, extrema, inflection point type.
character*8 atype ! Root, extrema, inflection point type.
dimension pa(1) ! Value of P(x) at x.
dimension pb(1) ! Value of P'(x) at x.
dimension pc(1) ! Value of P''(x) at x.
c.... Local variables.
dimension b(0:64) ! Coefficients of 1st derivative equation.
dimension c(0:64) ! Coefficients of 2nd derivative equation.
character*8 aft ! Root, extrema, inflection point type.
c=======================================================================********
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptpolr finding roots of a polynomial. nitmax =',i3 /
cbug & 'xmin,dx,xmax= ',1p3e22.14 /
cbug & 'tol= ',1pe22.14,' nxmax=',i5 )
cbug 9902 format ('n=',i3,' a(n)=',1pe22.14)
cbug write ( 3, 9901) nitmax, xmin, dx, xmax, tol, nxmax
cbug if (na .gt. 0) then
cbug write ( 3, 9902) (n, a(n), n = 0, na)
cbug endif
cbugc***DEBUG ends.
c.... Test for input errors.
nerr = 0
if (xmax .lt. xmin) then
nerr = 1
go to 999
elseif ((xmax .eq. xmin) .and. (dx .lt. 0.0)) then
nerr = 2
go to 999
elseif ((xmax .gt. xmin) .and. (dx .le. 0.0)) then
nerr = 3
go to 999
endif
if (nitmax .le. 0) then
nerr = 4
go to 999
endif
if (na .lt. 1) then
nerr = 5
go to 999
endif
if (nxmax .lt. 3 * (na - 1)) then
nerr = 6
go to 999
endif
c.... Eliminate trailing zero coefficients.
nar = na ! Index of last non-zero coefficient.
110 if (a(nar) .eq. 0.0) then
nar = nar - 1
go to 110
endif
c.... Initialize.
ntota = 0 ! Total # of roots.
ntotb = 0 ! Total # of extrema.
ntotc = 0 ! Total # of inflection points.
nx = 0 ! Total # of roots, extrema, infl. pts.
xpolys = -1.e99 ! Last saved value of xpoly.
do 120 n = 0, nar
a(n) = a(n) ! Coefficients of polynomial.
b(n) = 0.0 ! Coefficients of P'(x).
c(n) = 0.0 ! Coefficients of P''(x).
120 continue
do 130 n = 1, nxmax
xpoly(n) = -1.e99
atype(n) = 'unknown'
130 continue
c.... Check for leading zero coefficient of the polynomial equation.
cbugcbugc***DEBUG begins.
cbugcbug 9971 format (' x (root ) =',1pe20.12)
cbugcbugc***DEBUG ends.
if (a(0) .eq. 0.0) then ! One root is zero.
if (xpolys .ne. 0.0) then ! Found a new root at zero.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'root'
xpoly(nx) = 0.0
xpolys = 0.0
cbugcbugc***DEBUG begins.
cbugcbug xpolya = 0.0
cbugcbug write ( 3, 9971) xpolya
cbugcbugc***DEBUG ends.
endif ! Found a new root at zero.
endif ! Tested for leading zero coefficient.
c.... Find the coefficients of the 1st derivative equation.
nb = nar - 1
if (nb .ge. 0) then
do 160 n = 0, nb
b(n) = (n + 1) * a(n+1)
160 continue
endif
c.... Check for leading zero coefficient of the 1st derivative equation.
cbugcbugc***DEBUG begins.
cbugcbug 9972 format (' x (extremum) =',1pe20.12)
cbugcbugc***DEBUG ends.
if ((nb .ge. 1) .and. (b(0) .eq. 0.0)) then ! Extremum at x = 0.
if (xpolys .ne. 0.0) then ! Found a new extremum at zero.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'extremum'
xpoly(nx) = 0.0
xpolys = 0.0
cbugcbugc***DEBUG begins.
cbugcbug xpolyb = 0.0
cbugcbug write ( 3, 9972) xpolyb
cbugcbugc***DEBUG ends.
endif ! Found a new extremum at zero.
endif ! Tested for leading zero coefficient.
c.... Find the coefficients of the 2nd derivative equation.
nc = nb - 1
if (nc .ge. 0) then ! Test order.
do 170 n = 0, nc
c(n) = (n + 1) * b(n+1)
170 continue
endif ! Order is 0 or more.
cbugcbugc***DEBUG begins.
cbugcbug 9973 format (' n, a, b, c = ',i3,1p3e20.12)
cbugcbug write ( 3, 9973) (n, a(n), b(n), c(n), n = 0, nar)
cbugcbugc***DEBUG ends.
c.... Check for leading zero coefficient of the 2nd derivative equation.
cbugcbugc***DEBUG begins.
cbugcbug 9974 format (' x (inflectn) =',1pe20.12)
cbugcbugc***DEBUG ends.
if ((nc .ge. 1) .and. (c(0) .eq. 0.0)) then ! Inflection at x = 0.
if (xpolys .ne. 0.0) then ! Found a new inflection point at zero.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'inflectn'
xpoly(nx) = 0.0
xpolys = 0.0
cbugcbugc***DEBUG begins.
cbugcbug xpolyc = 0.0
cbugcbug write ( 3, 9974) xpolyc
cbugcbugc***DEBUG ends.
endif ! Found a new inflection point at zero.
endif ! Tested for leading zero coefficient.
c.... Find the number of trial values of x.
if (dx .eq. 0.0) then
ntry = 1
else
ntry = (xmax - xmin) / dx + 2.0
endif
c.... Initialize for the loop over trial values of x.
xtry = xmin ! Trial value of x.
xerr = 0.0 ! Error in calculating xtry.
polda = 0.0 ! Value of polynomial.
poldb = 0.0 ! Value of 1st derivative of polynomial.
poldc = 0.0 ! Value of 2nd derivative of polynomial.
nmod = 1 + (ntry - 1) / 32 ! Index interval between edits.
c.... Loop over the trial values of x.
cbugcbugc***DEBUG begins.
cbugcbug 9975 format (/ ' Trial',15x,'x',19x,'P(x)',9x,'P''(x)',8x,'P''''(x)')
cbugcbug write ( 3, 9975)
cbugcbugc***DEBUG ends.
do 220 nt = 1, ntry ! Loop over trial values of x.
c.... Find the value of the polynomial equation, its 1st derivative
c.... and 2nd derivative at xtry.
ptrya = a(0) ! Value of polynomial.
ptryb = b(0) ! Value of 1st derivative of polynomial.
ptryc = c(0) ! Value of 2nd derivative of polynomial.
perra = 0.0 ! Error in polynomial value.
perrb = 0.0 ! Error in 1st derivative value.
perrc = 0.0 ! Error in 2nd derivative value.
xperr = 0.0 ! Error in x**n, n = 0, nar.
xpow = 1.0 ! Value of x**n, n = 0, nar.
do 210 n = 1, nar ! Loop over remaining coefficients.
xperr = abs (xtry) * (xperr + tol * abs (xpow))
xpow = xpow * xtry
ptrya = ptrya + a(n) * xpow
ptryb = ptryb + b(n) * xpow
ptryc = ptryc + c(n) * xpow
perra = perra + abs (a(n)) * (xperr + tol * abs (xpow))
perrb = perrb + abs (b(n)) * (xperr + tol * abs (xpow))
perrc = perrc + abs (c(n)) * (xperr + tol * abs (xpow))
210 continue ! End of loop over coefficients.
if (abs (ptrya) .le. perra) ptrya = 0.0
if (abs (ptryb) .le. perrb) ptryb = 0.0
if (abs (ptryc) .le. perrc) ptryc = 0.0
cbugcbugc***DEBUG begins.
cbugcbugc.... Display argument and function values up to 32 times.
cbugcbug 9976 format (' x (trial) =',1pe20.12,1p3e13.5)
cbugcbug if ((mod (nt, nmod) .eq. 0) .or.
cbugcbug & (xtry .eq. xmin) .or.
cbugcbug & (xtry .eq. xmax)) then
cbugcbug write ( 3, 9976) xtry, ptrya, ptryb, ptryc
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Find the search limits for this trial x.
xmint = xtry - dx
if (xmint .lt. xmin) xmint = xmin
xmaxt = xtry + dx
if (xmaxt .gt. xmax) xmaxt = xmax
c.... Try to find a root of the polynomial equation if the polynomial
c.... or the slope changes sign in the preceding interval,
cbugcbugc***DEBUG begins.
cbugcbug 9977 format (' x (root) =',1pe20.12)
cbugcbugc***DEBUG ends.
if (ptrya .eq. 0.0) then ! Root at xtry.
xdiffe = tol * (abs (xtry) + abs (xpolys))
if (abs (xtry - xpolys) .gt. xdiffe) then ! New root at xtry.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'root'
xpoly(nx) = xtry
xpolys = xtry
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9977) xtry
cbugcbugc***DEBUG ends.
endif ! Found a new root at xtry.
elseif (ptrya * polda .lt. 0.0) then ! Polynomial changed sign.
call aptnewt (a, nar, xtry, xmint, xmaxt, nitmax, tol,
& xpolya, nerra)
if (nerra .eq. 0) then
xdiffe = tol * (abs (xpolya) + abs (xpolys))
if (abs (xpolya - xpolys) .gt. xdiffe) then ! New root at xpolya.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'root'
xpoly(nx) = xpolya
xpolys = xpolya
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9977) xpolya
cbugcbugc***DEBUG ends.
endif ! Found a new root at xpolya.
endif
if (xpolya .gt. xtry) then ! Try an interpolated x value.
xint = (xmint * ptrya - xtry * polda) / (ptrya -polda)
if (xint .lt. xmint) xint = xmint
if (xint .gt. xtry) xint = xtry
call aptnewt (a, nar, xint, xmint, xtry, nitmax, tol,
& xpolya, nerra)
if (nerra .eq. 0) then
xdiffe = tol * (abs (xpolya) + abs (xpolys))
if (abs (xpolya - xpolys) .gt. xdiffe) then ! New root at xpolya.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'root'
xpoly(nx) = xpolya
xpolys = xpolya
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9977) xpolya
cbugcbugc***DEBUG ends.
endif ! Found a new root at xpolya.
endif
endif ! Tried an interpolated x value.
elseif (ptryb * poldb .le. 0.0) then ! Slope zero or changed sign.
call aptnewt (a, nar, xtry, xmint, xmaxt, nitmax, tol,
& xpolya, nerra)
if (nerra .eq. 0) then
xdiffe = tol * (abs (xpolya) + abs (xpolys))
if (abs (xpolya - xpolys) .gt. xdiffe) then ! New root at xpolya.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'root'
xpoly(nx) = xpolya
xpolys = xpolya
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9977) xpolya
cbugcbugc***DEBUG ends.
endif ! Found new root at xpolya.
endif
if (xpolya .gt. xtry) then ! Try an interpolated x value.
pdiffb = ptryb - poldb
perrb = tol * (abs (ptryb) + abs (poldb))
if (abs (pdiffb) .le. perrb) pdiffb = 0.0
if (pdiffb .ne. 0.0) then
xint = (xmint * ptryb - xtry * poldb) / pdiffb
else
xint = 0.5 * (xmint + xtry)
endif
if (xint .lt. xmint) xint = xmint
if (xint .gt. xtry) xint = xtry
call aptnewt (a, nar, xint, xmint, xtry, nitmax, tol,
& xpolya, nerra)
if (nerra .eq. 0) then
xdiffe = tol * (abs (xpolya) + abs (xpolys))
if (abs (xpolya - xpolys) .gt. xdiffe) then ! New root at xpolya.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'root'
xpoly(nx) = xpolya
xpolys = xpolya
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9977) xpolya
cbugcbugc***DEBUG ends.
endif ! Found new root at xpolya.
endif
endif ! Tried an interpolated x value.
endif ! Tried to find a root.
c.... Try to find an extremum, if the slope changes sign in the
c.... preceding interval.
if (ptryb .eq. 0.0) then
xdiffe = tol * (abs (xtry) + abs (xpolys))
if (abs (xtry - xpolys) .gt. xdiffe) then ! New extremum at xtry.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'extremum'
xpoly(nx) = xtry
xpolys = xtry
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9972) xtry
cbugcbugc***DEBUG ends.
endif ! Found new root at xtrya.
elseif (ptryb * poldb .lt. 0.0) then
call aptnewt (b, nb, xtry, xmint, xmaxt, nitmax, tol,
& xpolyb, nerrb)
if (nerrb .eq. 0) then
xdiffe = tol * (abs (xpolyb) + abs (xpolys))
if (abs (xpolyb - xpolys) .gt. xdiffe) then ! New extr at xpolyb.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'extremum'
xpoly(nx) = xpolyb
xpolys = xpolyb
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9972) xpolyb
cbugcbugc***DEBUG ends.
endif ! Found new extremum at xpolyb.
endif
if (xpolyb .gt. xtry) then ! Try an interpolated x value.
xint = (xmint * ptryb - xtry * poldb) / (ptryb -poldb)
if (xint .lt. xmint) xint = xmint
if (xint .gt. xtry) xint = xtry
call aptnewt (b, nb, xint, xmint, xtry, nitmax, tol,
& xpolyb, nerrb)
if (nerrb .eq. 0) then
xdiffe = tol * (abs (xpolyb) + abs (xpolys))
if (abs (xpolyb - xpolys) .gt. xdiffe) then ! New extr at xpolyb.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'extremum'
xpoly(nx) = xpolyb
xpolys = xpolyb
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9972) xpolyb
cbugcbugc***DEBUG ends.
endif ! Found new extremum at xpolyb.
endif
endif ! Tried an interpolated x value.
endif ! Tried to find an extremum.
c.... Try to find an inflection point, if the 2nd derivative changes
c.... sign in the preceding interval.
if ((nar .ge. 3) .and. (ptryc .eq. 0.0)) then
xdiffe = tol * (abs (xtry) + abs (xpolys))
if (abs (xtry - xpolys) .gt. xdiffe) then ! New infl point at xtry.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'inflectn'
xpoly(nx) = xtry
xpolys = xtry
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9974) xtry
cbugcbugc***DEBUG ends.
endif ! Found new inflection point at xtry.
elseif (ptryc * poldc .lt. 0.0) then
call aptnewt (c, nc, xtry, xmint, xmaxt, nitmax, tol,
& xpolyc, nerrc)
if (nerrc .eq. 0) then
xdiffe = tol * (abs (xpolyc) + abs (xpolys))
if (abs (xpolyc - xpolys) .gt. xdiffe) then ! New infl at xpolyc.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'inflectn'
xpoly(nx) = xpolyc
xpolys = xpolyc
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9974) xpolyc
cbugcbugc***DEBUG ends.
endif ! Found new inflection point at xpolyc.
endif
if (xpolyc .gt. xtry) then ! Try an interpolated x value.
xint = (xmint * ptryc - xtry * poldc) / (ptryc -poldc)
if (xint .lt. xmint) xint = xmint
if (xint .gt. xtry) xint = xtry
call aptnewt (c, nc, xint, xmint, xtry, nitmax, tol,
& xpolyc, nerrc)
if (nerrc .eq. 0) then
xdiffe = tol * (abs (xpolyc) + abs (xpolys))
if (abs (xpolyc - xpolys) .gt. xdiffe) then ! New infl at xpolyc.
nx = nx + 1
if (nx .gt. nxmax) then
nerr = 7
nx = nxmax
endif
atype(nx) = 'inflectn'
xpoly(nx) = xpolyc
xpolys = xpolyc
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9974) xpolyc
cbugcbugc***DEBUG ends.
endif ! Found new inflection point at xpolyc.
endif
endif ! Tried an interpolated x value.
endif ! Tried to find an inflection point.
c.... Save old values of polynomial, 1st and 2nd derivatives.
polda = ptrya
poldb = ptryb
poldc = ptryc
c.... See if maximum x value has been reached.
if (xtry .eq. xmax) go to 310
c.... Find next x value.
xtry = xtry + dx
xerr = xerr + tol * (abs (xtry) + abs (dx))
if (abs (xtry) .le. xerr) xtry = 0.0
if (abs (xtry - xmax) .le. xerr) xtry = xmax
if (xtry .gt. xmax) xtry = xmax
220 continue ! End of loop over trial values of x.
310 if (nx .eq. 0) go to 999
c.... Put all x values in ascending order.
do 330 n = 1, nx - 1
do 320 nn = n + 1, nx
if (xpoly(nn) .lt. xpoly(n)) then
xft = xpoly(n)
aft = atype(n)
xpoly(n) = xpoly(nn)
atype(n) = atype(nn)
xpoly(nn) = xft
atype(nn) = aft
endif
320 continue
330 continue
c.... Eliminate duplicate values of xpoly.
nn = 1
do 340 n = 2, nx
xdiff = xpoly(n) - xpoly(n-1)
xdiffe = tol * (abs (xpoly(n)) + abs (xpoly(n-1)))
if (abs (xdiff) .le. xdiffe) xdiff = 0.0
if (xdiff .ne. 0.0) then
nn = nn + 1
xpoly(nn) = xpoly(n)
atype(nn) = atype(n)
endif
340 continue
nx = nn
c.... Find values of polynomial, 1st derivative, 2nd derivative at x values.
do 420 n = 1, nx ! Loop over trial values of x.
pa(n) = a(0) ! Value of polynomial.
pb(n) = b(0) ! Value of 1st derivative of polynomial.
pc(n) = c(0) ! Value of 2nd derivative of polynomial.
perra = 0.0 ! Error in polynomial value.
perrb = 0.0 ! Error in 1st derivative value.
perrc = 0.0 ! Error in 2nd derivative value.
xperr = 0.0 ! Error in x**n, n = 0, nar.
xpow = 1.0 ! Value of x**n, n = 0, nar.
do 410 nn = 1, nar ! Loop over remaining coefficients.
xperr = abs (xpoly(n)) * (xperr + tol * abs (xpow))
xpow = xpow * xpoly(n)
pa(n) = pa(n) + a(nn) * xpow
pb(n) = pb(n) + b(nn) * xpow
pc(n) = pc(n) + c(nn) * xpow
perra = perra + abs (a(nn)) * (xperr + tol * abs (xpow))
perrb = perrb + abs (b(nn)) * (xperr + tol * abs (xpow))
perrc = perrc + abs (c(nn)) * (xperr + tol * abs (xpow))
410 continue ! End of loop over coefficients.
if (abs (pa(n)) .le. perra) pa(n) = 0.0
if (abs (pb(n)) .le. perrb) pb(n) = 0.0
if (abs (pc(n)) .le. perrc) pc(n) = 0.0
atype(n) = 'unknown'
if ((nar .ge. 3) .and. (pc(n) .eq. 0.0)) then ! Inflection point.
ntotc = ntotc + 1
atype(n) = 'inflectn'
endif
if (pb(n) .eq. 0.0) then ! Extremum.
ntotb = ntotb + 1
if (pc(n) .lt. 0.0) then
atype(n) = 'maximum'
elseif (pc(n) .gt. 0.0) then
atype(n) = 'minimum'
elseif (pc(n) .eq. 0.0) then
atype(n) = 'min/max'
endif
endif
if (pa(n) .eq. 0.0) then ! Root.
ntota = ntota + 1
if (atype(n) .eq. 'maximum') then
atype(n) = 'dbl max'
elseif (atype(n) .eq. 'minimum') then
atype(n) = 'dbl min'
elseif (atype(n) .eq. 'min/max') then
atype(n) = 'triple'
else
atype(n) = 'single'
endif
endif
420 continue ! End of loop over trial values of x.
999 continue
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptpolr results. nerr=',i3 )
cbug 9904 format (' # of roots =',i3,', extrema =',i3,
cbug & ', infl points =',i3,', all=',i4,'.')
cbug 9905 format (/ ' Event',15x,'x',19x,'P(x)',9x,'P''(x)',8x,'P''''(x)')
cbug 9906 format (' x (',a8,') =',1pe20.12,1p3e13.5)
cbug 9907 format ('Ran out of storage space for xpoly.')
cbug write ( 3, 9903) nerr
cbug if ((nerr .eq. 0) .or. (nerr .eq. 7)) then
cbug write ( 3, 9904) ntota, ntotb, ntotc, nx
cbug write ( 3, 9905)
cbug write ( 3, 9906) (atype(n), xpoly(n), pa(n), pb(n), pc(n),
cbug & n = 1, nx)
cbug endif
cbug if (nx .eq. nxmax) then
cbug write ( 3, 9907)
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine aptpolr. (+1 line)
end
UCRL-WEB-209832