subroutine aptalsb (fbin, nalb, nbins, np, frstr, nb, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTALSB c c call aptalsb (fbin, nalb, nbins, np, frstr, nb, nerr) c c Version: aptalsb Updated 1990 November 27 14:00. c aptalsb Originated 1990 July 9 10:40. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To randomly sample np bin numbers nb from nbins aliased bins, c each with probability fbin(n) of selecting bin n, and c probability 1.0 - fbin(n) of selecting bin nalb(n). c Any fraction frstr of sampling may be striated, from 0.0 to 1.0. c Flag nerr indicates any input error. c c WARNING: STRIATED SAMPLING INTRODUCES CORRELATIONS BETWEEN SAMPLE INDICES c AND BIN INDICES. USE ONLY WHEN SUCH CORRELATIONS ACCEPTABLE. c c Note: Use subroutine aptalin to convert a piecewise linear probability c distribution to the aliased form required by aptalsb. c c Timing: For a test problem with 8 bins, 1000 samples, the cpu time was c 970 microseconds unstriated, 664 microseconds striated. c c Input: fbin, nalb, nbins, np, frstr. c c Output: nb, nerr. c c Glossary: c c fbin Input For an initially selected bin n, fbin(n) is the c probability of keeping that selection, and c 1.0 - fbin(n) is the probability of selecting bin c nalb(n) instead. Size nbins. c c frstr Input The fraction (from 0.0 to 1.0) of samples to striate c over the bins. Remaining samples will be selected c randomly over the bins. c c nalb Input For bin n, nalb(n) is the index of the alternate, c or aliased bin. Size nbins. c c nb Output The randomly sampled bin numbers. Size np. c c nbins Input The number of bins. c c nerr Output Indicates an input error, if not 0. c 1 if nbins is not positive. c 2 if np is not positive. c 3 if frstr is less than 0.0 or more than 1.0. c c np Input Size of array nb. Number of samples. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Probability of keeping initial bin. dimension fbin (1) c---- Index of alias bin. dimension nalb (1) c---- Index of sampled bin. dimension nb (1) c.... Local variables. c---- Index in internal array. common /laptalsb/ n c---- First index of subset of data. common /laptalsb/ n1 c---- Last index of subset of data. common /laptalsb/ n2 c---- Index in external array. common /laptalsb/ nn c---- Size of current subset of data. common /laptalsb/ ns c---- Number of striated samples. common /laptalsb/ nstrd c---- Number of unstriated samples. common /laptalsb/ nunst c---- Trial value of bin index. common /laptalsb/ nx c---- Random numbers in range 0.0 to 1.0. common /laptalsb/ ranfp (64) c---- Striated random number, 0.0 to 1.0. common /laptalsb/ ranst c---- Unrounded, floating point nx. common /laptalsb/ xn cbugc***DEBUG begins. cbugc---- Fraction of samples in bin. cbug common /laptalsb/ fb cbugc---- Index in loop over bins. cbug common /laptalsb/ nbin cbug 9901 format (/ 'aptalsb sampling from aliased bins.' / cbug & 'np=',i5,' frstr=',1pe22.14 / cbug & (i3,' fbin=',1pe22.14,' nalb=',i3)) cbug write (3, 9901) np, frstr, (n, fbin(n), nalb(n), n = 1, nbins) cbugc***DEBUG ends. c.... initialize. nerr = 0 c.... Test for input errors. if (nbins .le. 0) then nerr = 1 go to 210 endif if (np .le. 0) then nerr = 2 go to 210 endif if ((frstr .lt. 0.0) .or. (frstr .gt. 1.0)) then nerr = 3 go to 210 endif c.... Find the number of striated and unstriated samples. nstrd = frstr * np + 0.5 nunst = np - nstrd c.... See if any unstriated samples are needed. c---- Need unstriated samples. if (nunst .gt. 0) then c.... Set up the indices of the first subset of samples. n1 = 1 n2 = min (nunst, 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 numbers. c---- Loop over samples. do 130 n = 1, ns nn = n + n1 - 1 xn = 1.0 + ranfp(n) * nbins nx = xn if ((xn - nx) .lt. fbin(nx)) then nb(nn) = nx else nb(nn) = nalb(nx) 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. nunst) then n1 = n2 + 1 n2 = min (nunst, n1 + 63) c---- End of loop over subset of samples. go to 110 endif c---- Tested nunst. endif c.... See if any striated samples are needed. c---- Need striated samples. if (nstrd .gt. 0) then c.... Set up the indices of the first subset of samples. n1 = 1 + nunst n2 = min (np, n1 + 63) c.... Loop over the subset of samples. c---- Loop over subset of samples. 140 ns = n2 - n1 + 1 c.... Generate the needed random numbers. c---- Loop over samples. do 150 n = 1, ns ranfp(n) = ranf( ) c---- End of loop over samples. 150 continue c.... Randomly sample bin numbers, using striated sampling. c---- Loop over samples. do 160 n = 1, ns nn = n + n1 - 1 ranst = (nn - nunst - ranfp(n)) / nstrd xn = 1.0 + ranst * nbins nx = xn if ((xn - nx) .lt. fbin(nx)) then nb(nn) = nx else nb(nn) = nalb(nx) endif c---- End of loop over samples. 160 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 140 endif c---- Tested nstrd. endif cbugc***DEBUG begins. cbug 9902 format (/ 'aptalsb results:' / cbug & ' nunst=',i5,' nstrd=',i5) cbug 9903 format (5(' n=',i5,' nb=',i3)) cbug write ( 3, 9902) nunst, nstrd cbug write ( 3, 9903) (n, nb(n), n = 1, np) cbug cbug 9904 format (/ ' summary by bins:') cbug 9905 format (i5,' samples=',i5,' fraction=',1pe22.14) cbug write ( 3, 9904) cbugc---- Loop over bins. cbug do 190 nbin = 1, nbins cbug ns = 0 cbugc---- Loop over samples. cbug do 185 n = 1, np cbug if (nb(n) .eq. nbin) then cbug ns = ns + 1 cbug endif cbugc---- End of loop over samples. cbug 185 continue cbug fb = ns / (np + 1.e-99) cbug write ( 3, 9905) nbin, ns, fb cbugc---- End of loop over bins. cbug 190 continue cbugc***DEBUG ends. 210 return c.... End of subroutine aptalsb. (+1 line.) end UCRL-WEB-209832