subroutine aptslid (xa, pa, xb, pb, np, xran, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSLID
c
c call aptslid (xa, pa, xb, pb, np, xran, nerr)
c
c Version: aptslid Updated 1990 December 3 14:20.
c aptslid Originated 1990 February 2 14:10.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To randomly sample np values of x from a linear distribution
c function having probability pa at xa, and probability pb at xb.
c Flag nerr indicates any input error.
c
c If both pa and pb are non-negative, the expected value of x is
c = fa * xa + fb * xb,
c where fa = (2.0 * pa + pb) / (3.0 * (pa + pb))
c and fb = (pa + 2.0 * pb) / (3.0 * (pa + pb)).
c
c If pa is negative, and pb is positive, then nerr = 3,
c no values of x will be sampled between xa and
c xa' = (xa * pb - xb * pa) / (pb - pa), and
c = fa' * xa + fb' * xb,
c where fa' = pb / (3.0 * (pb - pa))
c and fb' = (2.0 * pb - 3.0 * pa) / (3.0 * (pb - pa)).
c
c If pa is positive, and pb is negative, then nerr = 3,
c no values of x will be sampled between xb and
c xb' = (xb * pa - xa * pb) / (pa - pb), and
c = fa' * xa + fb' * xb,
c where fa' = (2.0 * pa - 3.0 * pb) / (3.0 * (pa - pb))
c and fb' = pa / (3.0 * (pa - pb)).
c
c If both pa and pb are non-positive, then nerr = 2, and
c no values of x will be sampled.
c
c History: 1990 February 5 13:20. Modified to eliminate sampling from
c any range of x with negative probability.
c
c Input: xa, pa, xb, pb, np.
c
c Output: xran, nerr.
c
c Glossary:
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if both pa and pb are non-positive.
c No x values are sampled.
c 3 if either pa or pb is negative, and the other is
c positive. Values of x are sampled only from the
c positive part of the probability distribution.
c
c np Input Size of array xran.
c
c pa, pb Input Relative probabilities of random variables xa and xb,
c respectively. Probability p(x) is linear in x.
c The values of pa and pb need not be normalized.
c If neither pa nor pb is positive (nerr = 2), no
c values of x will be sampled. If either pa or pb is
c negative, and the other is positive, (nerr = 3),
c x values will only be sampled from the positive
c part of the probability distribution.
c To find pa and pb, when is given:
c xma = (2.0 * xa + xb) / 3.0,
c xmb = (xa + 2.0 * xb) / 3.0.
c For between xa and xma, pb is negative, and:
c pb = -((xma - ) / ( - xa)) * pa.
c For between xma and xmb:
c ( - xma) * pa = (xmb - ) * pb.
c For between xmb and xb, pa is negative, and:
c pa = -(( - xmb) / (xb - )) * pb.
c
c xa, xb Input Values of random variable x with relative probabilities
c pa and pb, respectively.
c
c xran Output Randomly sampled value of x in range from xa to xb.
c Size np.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Randomly sampled value of x.
dimension xran (1)
c.... Local variables.
c---- Index in local array.
common /laptslid/ n
c---- First index of subset of data.
common /laptslid/ n1
c---- Last index of subset of data.
common /laptslid/ n2
c---- Index in external array.
common /laptslid/ nn
c---- Size of current subset of data.
common /laptslid/ ns
c---- Probability of x1.
common /laptslid/ p1
c---- Value of p(x) at xxa.
common /laptslid/ ppa
c---- Value of p(x) at xxb.
common /laptslid/ ppb
c---- Sum of ppa and ppb.
common /laptslid/ psum
c---- First random number.
common /laptslid/ ranf1 (64)
c---- Second random number.
common /laptslid/ ranf2 (64)
c---- First random number ranf1.
common /laptslid/ u1
c---- Value 1.0 - u1.
common /laptslid/ um
c---- First random value of x.
common /laptslid/ x1
c---- Second random value of x.
common /laptslid/ x2
c---- Adjusted value of xa.
common /laptslid/ xxa
c---- Adjusted value of xb.
common /laptslid/ xxb
cbugc***DEBUG begins.
cbugc---- Standard deviation from mean x.
cbug common /laptslid/ xdev
cbugc---- Mean value of x.
cbug common /laptslid/ xmean
cbug 9901 format (/ 'aptslid sampling from a linear distribution. np=',i3 /
cbug & ' xa= ',1pe22.14,' pa= ',1pe22.14 /
cbug & ' xb= ',1pe22.14,' pb= ',1pe22.14)
cbug write (3, 9901) np, xa, pa, xb, pb
cbugc***DEBUG ends.
c.... initialize.
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
if ((pa .le. 0.0) .and. (pb .le. 0.0)) then
nerr = 2
go to 210
endif
c.... Adjust the range of x to include only positive probabilities.
xxa = xa
xxb = xb
ppa = pa
ppb = pb
if (pa .lt. 0.0) then
nerr = 3
xxa = (xa * pb - xb * pa) / (pb - pa)
ppa = 0.0
elseif (pb .lt. 0.0) then
nerr = 3
xxb = (xa * pb - xb * pa) / (pb - pa)
ppb = 0.0
endif
cbugc***DEBUG begins.
cbug 9902 format (/ ' Eliminated negative probabilities:' /
cbug & ' xxa=',1pe22.14,' ppa=',1pe22.14 /
cbug & ' xxb=',1pe22.14,' ppb=',1pe22.14)
cbug if ((xxa .ne. xa) .or. (xxb .ne. xb)) then
cbug write ( 3, 9902) xxa, ppa, xxb, ppb
cbug endif
cbugc***DEBUG ends.
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Generate ns pairs of random numbers.
c---- Loop over local arrays.
do 120 n = 1, ns
ranf1(n) = ranf( )
c---- End of loop over local arrays.
120 continue
c---- Loop over local arrays.
do 130 n = 1, ns
ranf2(n) = ranf( )
c---- End of loop over local arrays.
130 continue
c.... Sample values of x from the distribution p(x).
c---- Loop over local arrays.
do 140 n = 1, ns
nn = n + n1 - 1
u1 = ranf1(n)
um = 1.0 - u1
x1 = um * xxa + u1 * xxb
x2 = u1 * xxa + um * xxb
p1 = um * ppa + u1 * ppb
psum = ppa + ppb
if (p1 .ge. ranf2(n) * psum) then
xran(nn) = x1
else
xran(nn) = x2
endif
c---- End of loop over local arrays.
140 continue
c.... See if all data subsets are done.
c---- Do another subset of data.
if (n2 .lt. np) then
n1 = n2 + 1
n2 = min (np, n1 + 63)
go to 110
endif
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptslid results: xran')
cbug write ( 3, 9903)
cbug call aptmean (xran, np, 1.e-11, xmean, xdev, nerr)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptslid. (+1 line.)
end
UCRL-WEB-209832