subroutine aptgrat (t1, nt, sumt, tol, ratt, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTGRAT
c
c call aptgrat (t1, nt, sumt, tol, ratt, nerr)
c
c Version: aptgrat Updated 1994 August 8 10:00.
c aptgrat Originated 1994 August 8 10:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the geometric series with first term t1, a total of nt terms
c with sum sumt, to find the ratio ratt between successive terms,
c to within a relative error tol.
c Flag nerr indicates any input error.
c Method: solve F = (ratt**nt - 1) / (ratt - 1) = sumt / t1
c by Newtonian iteration, with first guess from approximation
c F ~= 1 + (nt - 2) * ratt**((nt - 1)/2) + ratt**(nt -1).
c
c Input: t1, nt, sumt
c
c Output: ratt, nerr
c
c Glossary:
c
c nerr Output If not 0, indicates an input error.
c 1 if t1 is not positive.
c 2 if nt is not 2 or more.
c 3 if sumt is not greater than t1.
c 4 if more than 100 iterations failed to converge.
c
c nt Input Number of terms in the geometric series.
c
c tol Input Maximum relative error of solution.
c Method will not converge if zero.
c On Cray computers, recommend 1.e-11.
c
c t1 Input The first term of the geometric series.
c If zero, replaced by sumt / nt, and a value of
c ratt = 1 is returned, with nerr = -1.
c
c ratt Output The ratio between successive terms of the series.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptgrat finding ratio of terms. nt=',i5,
cbug & ' tol=',1pe22.14 /
cbug & ' t1 =' 1pe22.14,' sumt=',1pe22.14 )
cbug write ( 3, 9901) nt, tol, t1, sumt
cbug write ( *, 9901) nt, tol, t1, sumt
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
ratt = -1.e99
c.... Test for input errors.
if (t1 .le. 0.0) then
nerr = 1
go to 210
endif
if (nt .le. 1) then
nerr = 2
go to 210
endif
if (sumt .le. t1) then
nerr = 3
go to 210
endif
c.... Find the ratio between terms, by Newtonian iteration.
fn = nt
q = sumt / t1
ntry = 0
if (q .eq. fn) then
ratt = 1.0
go to 210
endif
c.... First quess: solve approximate quadratic for ratt**((nt-1)/2).
ratt = 1.0 - 0.5 * fn + sqrt (q - fn * (1.0 - 0.25 * fn))
ratt = ratt**(2.0 / (fn - 1.0))
nm1 = nt - 1
cbugc***DEBUG begins.
cbug 9801 format (' aptgrat DEBUG. Initial guess = ',1pe22.14)
cbug write ( 3, 9801) ratt
cbug write ( *, 9801) ratt
cbugc***DEBUG ends.
110 ntry = ntry + 1
rnm1 = ratt**nm1
rm1 = ratt - 1.0
qq = fn
if (abs (rm1) .gt. tol) then
qq = (ratt * rnm1 - 1.0) / rm1
endif
cbugc***DEBUG begins.
cbug 9802 format (' aptgrat DEBUG. q =',1pe22.14 /
cbug & ' qq =',1pe22.14 /
cbug & ' diff=',1pe22.14 )
cbug diff = qq - q
cbug write ( 3, 9802) q, qq, diff
cbug write ( *, 9802) q, qq, diff
cbugc***DEBUG ends.
if (abs (qq - q) .le. tol * q) go to 210
if (ntry .gt. 100) then
nerr = 4
go to 210
endif
ratt = ratt + rm1 * (qq - q) / (qq - fn * rnm1)
cbugc***DEBUG begins.
cbug 9803 format (' aptgrat DEBUG. Next guess = ',1pe22.14)
cbug write ( 3, 9803) ratt
cbug write ( *, 9803) ratt
cbugc***DEBUG ends.
go to 110
210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptgrat results: ',' nerr=',i3 /
cbug & ' ratt=',1pe22.14 )
cbug write ( 3, 9902) nerr, ratt
cbug write ( *, 9902) nerr, ratt
cbugc***DEBUG ends.
return
c.... End of subroutine aptgrat. (+1 line.)
end
UCRL-WEB-209832