subroutine aptcumt (pcum, istart, nbins, incr, np, nb, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCUMT
c
c call aptcumt (pcum, istart, nbins, incr, np, nb, nerr)
c
c Version: aptcumt Updated 1990 August 24 11:00.
c aptcumt Originated 1990 August 24 11:00.
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 cumulative PDFs for all of the samples, and
c each sample is selected from its own specified cumulative PDF.
c For the n'th sample, the selection is from nbins(n) bins,
c for which the index in the packed array pcum starts at
c istart(n), with an increment of incr (the same for all n).
c For the packed index nx, there is a probability pcum(nx) of
c sampling any nb(n) from istart(n) to incr * (nbins(n) - 1),
c with an increment of incr.
c Flag nerr indicates any input error.
c
c Note: Use subroutine aptcump to find the cumulative PDF for a
c histogram or piecewise linear PDF.
c Use subroutine aptcums to randomly sample bins from a single
c cumulative PDF.
c
c Timing: For a test problem with 8 bins, 1000 samples, the cpu time was
c 3845 microseconds (incr = 1), 6816 microseconds (incr = 2).
c
c Input: pcum, istart, nbins, incr, np.
c
c Output: nb, nerr.
c
c Glossary:
c
c incr Input The increment between memory locations for the values
c of fbin, or of nalb, for two adjacent bins in the
c same cumulative PDF. Must be the same for all of the
c cumulative PDFs needed for all np samples.
c
c istart Input For the n'th sample, istart(n) is the initial index,
c in the packed cumulative PDF bin arrays, of the pcum
c for the first bin of the PDF from which the n'th
c sample is to be selected.
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 cumulative 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 cumulative PDF specified for that
c sample. The size of arrays istart, nbins, nb.
c
c pcum Input The relative probability of selecting any of the bins
c in the current PDF, up to and including the bin with
c value pcum. Absolute probabilities may be obtained
c by dividing each pcum by the pcum for the last bin
c in the same PDF.
c The value pcum(1) is the first pcum for the first bin
c of the first PDF in the packed cumulative PDF arrays.
c The pcum values for the n'th sample have indices from
c istart(n) to istart(n) + incr * (nbins(n) - 1), in
c increments of incr. For any such index nx, pcum(nx)
c is the relative probability of selecting nb(n)
c between istart(n) and nx.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Initial PDF index in pcum.
dimension istart (1)
c---- Packed array index of sampled bin.
dimension nb (1)
c---- Number of bins in an cumulative PDF.
dimension nbins (1)
c---- Cumulative probability, bins 1 to n.
dimension pcum (1)
c.... Local variables.
c---- Index in internal array.
common /laptcumt/ n
c---- First index of subset of data.
common /laptcumt/ n1
c---- Last index of subset of data.
common /laptcumt/ n2
c---- Index in external array.
common /laptcumt/ nn
c---- Size of current subset of data.
common /laptcumt/ ns
c---- Trial value of packed array bin index.
common /laptcumt/ nx
c---- Final packed index in one CPDF.
common /laptcumt/ nxmax
c---- Sampled relative cum probability.
common /laptcumt/ plook
c---- Random numbers in range 0.0 to 1.0.
common /laptcumt/ ranfp (64)
cbugc***DEBUG begins.
cbugc---- Highest starting index.
cbug common /laptcumt/ istartmx
cbugc---- Highest packed bin index.
cbug common /laptcumt/ nbmx
cbugc---- Sequential bin index in packed arrays.
cbug common /laptcumt/ nseq
cbugc---- Fraction of samples in bin.
cbug common /laptcumt/ fb
cbug 9901 format (/ 'aptcumt sampling from cumulative bins.',' incr=',i5)
cbug 9902 format (/ ' packed bins:')
cbug 9903 format (i5,' nseq=',i5,' pcum=',1pe22.14)
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 160 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 160 continue
cbug write ( 3, 9902)
cbug do 170 nx = 1, nbmx, incr
cbug nseq = 1 + (nx - 1) / incr
cbug write ( 3, 9903) nx, nseq, pcum(nx)
cbug 170 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 cumulative PDFs.
c---- Can use luf.
if (incr .eq. 1) then
c---- Loop over samples.
do 130 n = 1, ns
nn = n + n1 - 1
nxmax = istart(nn) + incr * (nbins(nn) - 1)
plook = ranfp(n) * pcum(nxmax)
nb(nn) = luf (plook, pcum(istart(nn)), nbins(nn) - 1)
c---- End of loop over samples.
130 continue
c---- Replace luf with table look-up loop.
else
c---- Loop over samples.
do 150 n = 1, ns
nn = n + n1 - 1
nxmax = istart(nn) + incr * (nbins(nn) - 1)
plook = ranfp(n) * pcum(nxmax)
do 140 nx = istart(n), nxmax, incr
c++++ Found table interval.
if (pcum(nx) .ge. plook) then
go to 145
endif
140 continue
145 nb(nn) = nx
c---- End of loop over samples.
150 continue
endif
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 (/ 'aptcumt 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 190 nx = 1, nbmx, incr
cbug nseq = 1 + (nx - 1) / incr
cbug ns = 0
cbugc---- Loop over samples.
cbug do 180 nn = 1, np
cbug if (nb(nn) .eq. nx) then
cbug ns = ns + 1
cbug endif
cbugc---- End of loop over samples.
cbug 180 continue
cbug fb = ns / (np + 1.e-99)
cbug write ( 3, 9907) nx, nseq, ns, fb
cbugc---- End of loop over packed bins.
cbug 190 continue
cbug write ( 3, 9908) (n, nb(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptcumt. (+1 line.)
end
UCRL-WEB-209832