subroutine aptpole (x, nx, a, na, tol, pa, pb, pc, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPOLE
c
c call aptpole (x, nx, a, na, tol, pa, pb, pc, nerr)
c
c Version: aptpole Updated 1997 August 12 14:20.
c aptpole Originated 1997 August 12 14:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find, for the arguments x(m), m = 1, nx, the values of the
c polynomial equation in x, P(x) = SUM {a(n)*x**n}, (n = 0, na),
c and its first and second derivatives P'(x) and P''(x),
c returned as the arrays pa, pb and pc, respectively.
c Flag nerr, if non-zero, indicates an input error.
c
c Method: The first derivative P'(x) = sum (b(n)*x**n, (n = 0, nb),
c where b(n) = (n + 1) * a(n+1), n = 0, nb, where nb = na -1.
c The second derivative P''(x) = sum (c(n)*x**n, (n = 0, nc),
c where c(n) = (n + 1) * b(n+1), n = 0, nc, where nc = nb -1.
c
c Input: x, nx, a, na, tol.
c
c Output: 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
c nerr Output Indicates no input errors were found, if zero.
c Otherwise:
c 1 nx was not positive.
c 2 na was not positive.
c
c na Input The largest exponent of x in the polynomial equation
c P(x), for which the coefficient is non-zero.
c Must be positive.
c
c nx Input The size of array x.
c
c pa Output The value of P(x). Size nx or larger.
c
c pb Output The value of P'(x). Size nx or larger.
c
c pc Output The value of P''(x). Size nx 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 pa, pb and pc less than the relative
c error in their calculation, based on tol, are
c truncated to zero.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
dimension a(0:1) ! The coefficients of x**n in P(x).
dimension x(1) ! Arguments of P(x), P'(x) and P''(x).
dimension pa(1) ! Value of P(x).
dimension pb(1) ! Value of P'(x).
dimension pc(1) ! Value of P''(x).
c.... Local variables.
dimension b(0:64) ! The coefficients of x**n in P'(x).
dimension c(0:64) ! The coefficients of x**n in P''(x).
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptpole finding values of P(x), P''(x), P''''(x)')
cbug 9902 format ('tol=',1pe22.14,' nx=',i5,' na=',i5)
cbug 9903 format ('n=',i5,' x=',1pe22.14)
cbug 9904 format ('n=',i5,' a=',1pe22.14)
cbug write ( 3, 9901)
cbug write ( 3, 9902) tol, nx, na
cbug if (nx .gt. 0) then
cbug write ( 3, 9903) (n, x(n), n = 1, nx)
cbug endif
cbug if (na .gt. 0) then
cbug write ( 3, 9904) (n, a(n), n = 0, na)
cbug endif
cbugc***DEBUG ends.
c.... Test for input errors.
nerr = 0
if (nx .le. 0) then
nerr = 1
go to 210
endif
if (na .le. 0) then
nerr = 2
go to 210
endif
c.... Find the coefficients of P'(x).
nb = na - 1
if (nb .ge. 0) then
do 120 n = 0, nb
b(n) = (n + 1) * a(n+1)
120 continue
endif
b(na) = 0.0
c.... Find the coefficients of P''(x).
nc = nb - 1
if (nc .ge. 0) then
do 130 n = 0, nc
c(n) = (n + 1) * b(n+1)
130 continue
endif
c(nb) = 0.0
c(na) = 0.0
c.... Find the values of P(x), P'(x) and P''(x) at each value of x.
do 150 n = 1, nx ! Loop over x array.
pa(n) = a(0) ! Value of P(x).
pb(n) = b(0) ! Value of P'(x).
pc(n) = c(0) ! Value of P''(x).
perra = 0.0 ! Error in P(x).
perrb = 0.0 ! Error in P'(x).
perrc = 0.0 ! Error in P''(x).
xperr = 0.0 ! Error in x**n.
xpow = 1.0 ! Initial value of x**n.
do 140 nn = 1, na ! Loop over remaining coefficients.
xperr = abs (x(n)) * (xperr + tol * abs (xpow))
xpow = xpow * x(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))
140 continue ! End of loop over coefficients.
c.... Truncate small values to zero.
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
150 continue ! End of loop over x array.
210 continue
cbugc***DEBUG begins.
cbug 9905 format ('aptpole results. nerr=',i3)
cbug 9906 format ('n,x,pa,pb,pc=',i3,1p4e16.8)
cbug write ( 3, 9905) nerr
cbug if (nerr .eq. 0) then
cbug write ( 3, 9906) (n, x(n), pa(n), pb(n), pc(n), n = 1, nx)
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine aptpole. (+1 line)
end
UCRL-WEB-209832