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