subroutine aptsliv (xa, pa, xb, pb, np, xran, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSLIV
c
c call aptsliv (xa, pa, xb, pb, np, xran, nerr)
c
c Version: aptsliv Updated 1990 December 3 14:20.
c aptsliv Originated 1990 February 6 11:40.
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 np linear distribution
c functions having probabilities pa at xa, and 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
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
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 the returned value
c of xran will be -1.e+99.
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
c np Input Size of arrays xa, pa, xb, pb, 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 A value of x will be sampled only from the positive
c range of each probability distribution function.
c If neither pa nor pb is positive, the value of
c xran will be -1.e+99. Size np.
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. Size np.
c
c xran Output Randomly sampled value of x in range from xa to xb.
c Value returned is -1.e+99 if neither pa nor pb
c is positive. Size np.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Probability at xa.
dimension pa (1)
c---- Probability at xb.
dimension pb (1)
c---- Value of x with probability pa.
dimension xa (1)
c---- Value of x with probability pb.
dimension xb (1)
c---- Randomly sampled value of x.
dimension xran (1)
c.... Local variables.
c---- A very big number.
common /laptsliv/ big
c---- A very small number.
common /laptsliv/ fuz
c---- Index in local array.
common /laptsliv/ n
c---- First index of subset of data.
common /laptsliv/ n1
c---- Last index of subset of data.
common /laptsliv/ n2
c---- Size of current subset of data.
common /laptsliv/ ns
c---- Index in external array.
common /laptsliv/ nn
c---- Probability of x1.
common /laptsliv/ p1
c---- Value of p(x) at xxa.
common /laptsliv/ ppa (64)
c---- Value of p(x) at xxb.
common /laptsliv/ ppb (64)
c---- Sum of ppa and ppb.
common /laptsliv/ psum
c---- First random number.
common /laptsliv/ ranf1 (64)
c---- Second random number.
common /laptsliv/ ranf2 (64)
c---- First random number ranf1.
common /laptsliv/ u1
c---- Value 1.0 - u1.
common /laptsliv/ um
c---- First random value of x.
common /laptsliv/ x1
c---- Second random value of x.
common /laptsliv/ x2
c---- Maximum of xa, xb.
common /laptsliv/ xmax
c---- Minimum of xa, xb.
common /laptsliv/ xmin
c---- Adjusted value of xa.
common /laptsliv/ xxa (64)
c---- Adjusted value of xb.
common /laptsliv/ xxb (64)
c---- Value of x where p(x) = 0.0.
common /laptsliv/ xzero
cbugc***DEBUG begins.
cbugc---- Mean value of x.
cbug common /laptsliv/ xmean
cbugc---- Standard deviation from mean x.
cbug common /laptsliv/ xdev
cbug 9901 format (/ 'aptsliv sampling from a linear distribution. np=',i3 /
cbug & (i3,' xa= ',1pe22.14,' pa= ',1pe22.14 /
cbug & ' xb= ',1pe22.14,' pb= ',1pe22.14))
cbug write (3, 9901) np, (n, xa(n), pa(n), xb(n), pb(n), n = 1, np)
cbugc***DEBUG ends.
c.... initialize.
c---- A very big number.
big = 1.e+99
c---- A very small number.
fuz = 1.e-99
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Adjust the ranges of x to include only positive probabilities.
c.... If not possible, will return x = -big.
c---- Loop over local arrays.
do 120 n = 1, ns
nn = n + n1 - 1
xzero = (xa(nn) * pb(nn) - xb(nn) * pa(nn)) /
& (pb(nn) - pa(nn) + fuz)
xmax = amax1 (xa(nn), xb(nn))
xmin = amin1 (xa(nn), xb(nn))
xzero = amax1 (xmin, amin1 (xzero, xmax))
if (pa(nn) .ge. 0.0) then
xxa(n) = xa(nn)
ppa(n) = pa(nn)
else
xxa(n) = xzero
ppa(n) = 0.0
endif
if (pb(nn) .ge. 0.0) then
xxb(n) = xb(nn)
ppb(n) = pb(nn)
else
xxb(n) = xzero
ppb(n) = 0.0
endif
c---- End of loop over local arrays.
120 continue
cbugc***DEBUG begins.
cbug 9902 format (/ ' Eliminated negative probabilities:' /
cbug & i3,' xxa=',1pe22.14,' ppa=',1pe22.14 /
cbug & ' xxb=',1pe22.14,' ppb=',1pe22.14)
cbugc---- Loop over local arrays.
cbug do 125 n = 1, ns
cbug nn = n + n1 - 1
cbug if ((xxa(n) .ne. xa(nn)) .or. (xxb(n) .ne. xb(nn))) then
cbug write ( 3, 9902) nn, xxa(n), ppa(n), xxb(n), ppb(n)
cbug endif
cbugc---- End of loop over local arrays.
cbug 125 continue
cbugc***DEBUG ends.
c.... Generate ns pairs of random numbers.
c---- Loop over local arrays.
do 130 n = 1, ns
ranf1(n) = ranf( )
c---- End of loop over local arrays.
130 continue
c---- Loop over local arrays.
do 140 n = 1, ns
ranf2(n) = ranf( )
c---- End of loop over local arrays.
140 continue
c.... Sample values of x from the distributions p(x).
c---- Loop over local arrays.
do 150 n = 1, ns
nn = n + n1 - 1
u1 = ranf1(n)
um = 1.0 - u1
x1 = um * xxa(n) + u1 * xxb(n)
x2 = u1 * xxa(n) + um * xxb(n)
p1 = um * ppa(n) + u1 * ppb(n)
psum = ppa(n) + ppb(n)
if (p1 .ge. ranf2(n) * psum) then
xran(nn) = x1
else
xran(nn) = x2
endif
if (psum .le. 0.0) then
xran(nn) = -big
endif
c---- End of loop over local arrays.
150 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 (/ 'aptsliv results: x')
cbug write ( 3, 9903)
cbug call aptmean (xran, np, 1.e-11, xmean, xdev, nerr)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptsliv. (+1 line.)
end
UCRL-WEB-209832