subroutine aptalst (fbin, nalb, istart, nbins, incr, np, nb, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTALST
c
c call aptalst (fbin, nalb, istart, nbins, incr, np, nb, nerr)
c
c Version: aptalst Updated 1990 November 27 14:00.
c aptalst Originated 1990 August 22 13:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To randomly sample np probability distribution function (PDF)
c bin indices nb. The indices are in packed PDF bin arrays,
c which include the aliased PDFs for all of the samples, and
c each sample is selected from its own specified PDF.
c For the n'th sample, the selection is from nbins(n) bins,
c for which the index in the packed arrays fbin and nalb
c starts at index istart(n), with an increment of incr.
c Each of the bins is initially sampled with equal probability,
c but the bin with packed index nx has a probability fbin(nx)
c of using nb(n) = nx, and a probability 1.0 - fbin(nx) of using
c nb(n) = istart(n) + incr * (nalb(nx) - 1), instead.
c Flag nerr indicates any input error.
c
c Note: Use subroutine aptalin to convert a piecewise linear PDF
c to the aliased form required by aptalst.
c Use subroutine aptalsb to randomly sample bins from a single
c aliased PDF.
c
c Timing: For a test problem with 9 bins, 1000 samples, the cpu time was
c 1130 microseconds (incr = 1), 1320 microseconds (incr = 2).
c
c Input: fbin, nalb, istart, nbins, incr, np.
c
c Output: nb, nerr.
c
c Glossary:
c
c fbin Input The probability of keeping the initially selected bin.
c The value fbin(1) is the first fbin for the first bin
c of the first PDF in the packed aliased PDF arrays.
c For the n'th sample, and the initially selected index
c nx in the packed PDF bin arrays, the probability of
c using nb(n) = nx is fbin(nx), and the probability of
c using nb(n) = istart(n) + incr * (nalb(nx) - 1) is
c 1.0 - fbin(nx).
c The fbin values for the n'th sample have indices from
c istart(n) to istart(n) + incr * (nbins(n) - 1), in
c increments of incr.
c
c incr Input The increment between memory locations for the values
c of fbin, or of nalb, for two adjacent bins in a
c single aliased PDF. Must be the same for all of the
c aliased PDFs needed for all np samples.
c
c istart Input For the n'th sample, istart(n) is the initial index, in
c the packed aliased PDF bin arrays, of the values of
c fbin and nalb for the first bin of the PDF from which
c the n'th sample is to be selected.
c
c nalb Input The bin number (1 to nbins(n)) of the aliased bin.
c The value nalb(1) is the first nalb for the first bin
c of the first PDF in the packed aliased PDF arrays.
c For the n'th sample, and the initially selected index
c nx in the packed PDF bin arrays, the probability of
c using nb(n) = nx is fbin(nx), and the probability of
c using nb(n) = istart(n) + incr * (nalb(nx) - 1) is
c 1.0 - fbin(nx).
c The nalb values for the n'th sample have indices from
c istart(n) to istart(n) + incr * (nbins(n) - 1), in
c increments of incr.
c
c nb Output The randomly sampled bin indices in the packed PDF bin
c arrays. For the n'th sample, nb(n) will be between
c istart(n) and istart(n) + incr * (nbins(n) - 1).
c The equivalent bin number within the sampled PDF is
c 1 + (nb(n) - istart(n)) / incr.
c The sequential bin number in the packed PDF tables is
c 1 + (nb(n) - 1) / incr.
c
c nbins Input For the n'th sample, nbins(n) is the number of bins
c in the aliased PDF from which the sample is to be
c selected.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c
c np Input The number of values of nb to be randomly sampled,
c each from the aliased PDF specified for that
c sample. The size of arrays istart, nbins, nb.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Probability of keeping initial bin.
dimension fbin (1)
c---- Initial PDF index in fbin, nalb.
dimension istart (1)
c---- Bin number of aliased bin.
dimension nalb (1)
c---- Packed array index of sampled bin.
dimension nb (1)
c---- Number of bins in an aliased PDF.
dimension nbins (1)
c.... Local variables.
c---- Index in internal array.
common /laptalst/ n
c---- First index of subset of data.
common /laptalst/ n1
c---- Last index of subset of data.
common /laptalst/ n2
c---- Packed array index of aliased bin.
common /laptalst/ na
c---- Randomly sampled local bin number.
common /laptalst/ nbin
c---- Index in external array.
common /laptalst/ nn
c---- Size of current subset of data.
common /laptalst/ ns
c---- Trial value of packed array bin index.
common /laptalst/ nx
c---- Random numbers in range 0.0 to 1.0.
common /laptalst/ ranfp (64)
c---- A second random number.
common /laptalst/ ranfp2
c---- Unrounded, floating point nbin.
common /laptalst/ xbin
c---- Unrounded, floating point nx.
common /laptalst/ xn
cbugc***DEBUG begins.
cbugc---- Highest starting index.
cbug common /laptalst/ istartmx
cbugc---- Highest packed bin index.
cbug common /laptalst/ nbmx
cbugc---- Sequential bin index in packed arrays.
cbug common /laptalst/ nseq
cbugc---- Fraction of samples in bin.
cbug common /laptalst/ fb
cbug 9901 format (/ 'aptalst sampling from aliased bins.',' incr=',i5)
cbug 9902 format (/ ' packed bins:')
cbug 9903 format (i5,' nseq=',i5,' fbin=',1pe22.14,' nalb=',i5)
cbug 9904 format (/ ' n, istart, nbins:' /
cbug & (i6,i5,i5,i6,i5,i5,i6,i5,i5,i6,i5,i5,i6,i5,i5))
cbug write ( 3, 9901) incr
cbug istartmx = istart(1)
cbug nbmx = istart(1) + incr * (nbins(1) - 1)
cbug do 140 n = 2, np
cbug if (istart(n) .gt. istartmx) then
cbug istartmx = istart(n)
cbug nbmx = istart(n) + incr * (nbins(n) - 1)
cbug endif
cbug 140 continue
cbug write ( 3, 9902)
cbug do 150 nx = 1, nbmx, incr
cbug nseq = 1 + (nx - 1) / incr
cbug write ( 3, 9903) nx, nseq, fbin(nx), nalb(nx)
cbug 150 continue
cbug write ( 3, 9904) (n, istart(n), nbins(n), n = 1, np)
cbugc***DEBUG ends.
c.... initialize.
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 samples.
n1 = 1
n2 = min (np, 64)
c.... Loop over the subset of samples.
c---- Loop over subset of samples.
110 ns = n2 - n1 + 1
c.... Generate the needed random numbers.
c---- Loop over samples.
do 120 n = 1, ns
ranfp(n) = ranf( )
c---- End of loop over samples.
120 continue
c.... Randomly sample bin indices in the packed arrays.
c---- Loop over samples.
do 130 n = 1, ns
nn = n + n1 - 1
xbin = 1 + ranfp(n) * nbins(nn)
nbin = xbin
ranfp2 = xbin - nbin
nx = istart(nn) + incr * (nbin - 1)
na = istart(nn) + incr * (nalb(nx) - 1)
if (ranfp2 .lt. fbin(nx)) then
nb(nn) = nx
else
nb(nn) = na
endif
c---- End of loop over samples.
130 continue
c.... See if all subsets of samples are done.
c---- Do another subset of data.
if (n2 .lt. np) then
n1 = n2 + 1
n2 = min (np, n1 + 63)
c---- End of loop over subset of samples.
go to 110
endif
cbugc***DEBUG begins.
cbug 9905 format (/ 'aptalst results:')
cbug 9906 format (/ ' summary by packed bins:')
cbug 9907 format (i5,' nseq=',i5,' samples=',i5,' fraction=',1pe22.14)
cbug 9908 format (5(' n=',i5,' nb=',i3))
cbug write ( 3, 9905)
cbug write ( 3, 9906)
cbugc---- Loop over packed bins.
cbug do 170 nx = 1, nbmx, incr
cbug nseq = 1 + (nx - 1) / incr
cbug ns = 0
cbugc---- Loop over samples.
cbug do 160 nn = 1, np
cbug if (nb(nn) .eq. nx) then
cbug ns = ns + 1
cbug endif
cbugc---- End of loop over samples.
cbug 160 continue
cbug fb = ns / (np + 1.e-99)
cbug write ( 3, 9907) nx, nseq, ns, fb
cbugc---- End of loop over packed bins.
cbug 170 continue
cbug write ( 3, 9908) (n, nb(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptalst. (+1 line.)
end
UCRL-WEB-209832