subroutine aptpolf (xroot, a, na, tol, nz, b, nb, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPOLF
c
c call aptpolf (xroot, a, na, tol, nz, b, nb, nerr)
c
c Version: aptpolf Updated 1997 September 25 13:40.
c aptpolf Originated 1997 September 25 13:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To factor a real root xroot, and any roots at x = 0, out of a
c polynomial equation P(x) of order na, with real coefficients
c a(n), n = 0, na. P(x) = sum (n = 1, na) {a(n) * x**n} = 0.
c and find the coefficients b(n), n = 1, nb of the reduced
c polynomial equation.
c
c Input: xroot, a, na, tol.
c
c Output: nz, b, nb, nerr.
c
c Glossary:
c
c a Input The coefficients of the polynomial equation P(x) of
c order na: 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, and need not be specified.
c
c b Output The coefficients of the reduced polynomial equation of
c order nb: R(x) = sum {b(n) * x**n} (n = 0, nb).
c Must be dimensioned a(0:na-1) or larger.
c R(x) = P(x) / (x - xroot), normalized to b(nb) = 1.
c If not calculated -1.e99 will be returned.
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 2 or more, or nerr = 1 will be returned
c
c nb Input The largest exponent of x in the reduced polynomial
c equation.
c
c nerr Output Indicates no input errors were found, if zero.
c Otherwise:
c 1 if na less than 2.
c 2 if xroot is zero.
c 3 if xroot is not a root.
c
c nz Output The number of leading coefficients of P(x) that are
c zero. There are nz roots at zero.
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-15 is recommended.
c
c xroot Input A root of the equation P(x) = 0. If not a root,
c nerr = 3 will be returned. If zero, nerr = 2 will
c be returned.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
dimension a(0:64) ! Coefficients of polynomial equation.
dimension b(0:64) ! Coefficients of factored equation.
c.... Local variables.
dimension berr(0:64) ! Error in calculating b (first method).
dimension c(0:64) ! Coefficients of factored equation.
dimension cerr(0:64) ! Error in calculating c (second method).
cbugc***DEBUG begins.
cbug 9901 format ('aptpolf factoring out root of polynomial.' /
cbug & ' xroot =',1pe22.14,' order=',i3,' tol=',1pe22.14)
cbug 9902 format (' a(',i2,') =',1pe22.14)
cbug write ( 3, 9901) xroot, na, tol
cbug if (na .ge. 0) then
cbug write ( 3, 9902) (n, a(n), n = 0, na)
cbug endif
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
nz = 0
do 110 n = 0, na - 1
b(n) = -1.e99
110 continue
c.... Eliminate trailing zero coefficients.
nax = na
120 if (a(nax) .eq. 0.0) then
nax = nax - 1
cbugcbugc***DEBUG begins.
cbugcbug 9930 format (' Reduced order to',i3,', trailing zero.')
cbugcbug write ( 3, 9930) nax
cbugcbugc***DEBUG ends.
if (nax .gt. 0) go to 120
endif
c.... Eliminate leading zero coefficients.
130 if ((nax .gt. 0) .and. (a(0) .eq. 0.0)) then ! A root is zero.
nax = nax - 1
nz = nz + 1
cbugcbugc***DEBUG begins.
cbugcbug 9931 format (' Reduced order to',i3,', roots at zero=',i3)
cbugcbug write ( 3, 9931) nax, nz
cbugcbugc***DEBUG ends.
do 140 n = 0, nax
a(n) = a(n+1)
140 continue
go to 130
endif ! Tested for leading zero coefficient.
nb = nax - 1
c.... Test for input errors.
if (nax .lt. 2) then
nerr = 1
go to 999
endif
if (xroot .eq. 0.0) then
nerr = 2
go to 999
endif
c.... See if xroot is really a root of the polynomial.
call aptpole (xroot, 1, a, nax, tol, pa, pb, pc, nerr)
if (pa .ne. 0.0) then
nerr = 3
go to 999
endif
c.... Find the coefficients b(n) (n = 0, nax - 1) of the reduced equation.
c.... b(n) = a(n+1) + a(n+2)*r + ... + a(N)*r**(N-n-1), or
c.... b(n) = -[a(0) + a(1)*r + ... + a(n)*r**n] / r**(n+1).
polya = a(0) ! Value of polynomial.
perra = 0.0 ! Error in polynomial value.
xperr = 0.0 ! Error in x**n, n = 0, nax.
xpow = 1.0 ! Value of x**n, n = 0, nax.
b(0) = -a(0) / xroot ! Coefficient of reduced equation.
berr(0) = 0.0 ! Error in calculating b.
do 220 n = 1, nax ! Loop over remaining coefficients.
xperr = abs (xroot) * (xperr + tol * abs (xpow))
xpow = xpow * xroot
polya = polya + a(n) * xpow
perra = perra + abs (a(n)) * (xperr + tol * abs (xpow))
b(n) = -polya / (xpow * xroot)
berrx = perra + abs (polya) * (tol + xperr / abs (xpow))
berr(n) = berrx / abs (xpow * xroot)
if (abs (b(n)) .le. berr(n)) b(n) = 0.0
220 continue ! End of loop over coefficients.
if (abs (polya) .le. perra) polya = 0.0
cbugcbugc***DEBUG begins.
cbugcbug 8220 format ('Polynomial, error=',1p2e20.12)
cbugcbug
cbugcbug write ( 6, 8220) polya, perra
cbugcbug write (iout, 8220) polya, perra
cbugcbugc***DEBUG ends.
c.... Find the coefficients of the reduced equation by the second method.
c.... Keep the coefficients with the lesser error.
do 320 n = 0, nb
c(n) = 0.0
cerr(n) = 0.0
xperr = 0.0
xpow = 1.0
do 310 nn = n + 1, nax
c(n) = c(n) + a(nn) * xpow
cerr(n) = cerr(n) + abs (a(nn)) * (xperr + tol * abs (xpow))
xperr = abs (xroot) * (xperr + tol * abs (xpow))
xpow = xpow * xroot
310 continue
if (abs (c(n)) .le. cerr(n)) c(n) = 0.0
if (cerr(n) .lt. berr(n)) then
b(n) = c(n)
berr(n) = cerr(n)
endif
320 continue
cbugcbugc***DEBUG begins.
cbugcbug 8320 format (' b(',i2,'), error =',1p2e20.12)
cbugcbug write ( 6, 8320) (n, b(n), berr(n), n = 0, nax - 1)
cbugcbug write (iout, 8320) (n, b(n), berr(n), n = 0, nax - 1)
cbugcbugc***DEBUG ends.
999 continue
cbugc***DEBUG begins.
cbug 9904 format ('aptpolf results. Roots at zero =',i2,
cbug & '. Reduced order = ',i3,' nerr =',i3)
cbug 9905 format (' b(',i2,') =',1pe22.14)
cbug write ( 3, 9904) nz, nb, nerr
cbug if (nerr .eq. 0) then
cbug write ( 3, 9905) (n, b(n), n = 0, nb)
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine aptpolf. (+1 line)
end
UCRL-WEB-209832