subroutine apteqsb (nbep, nep, np, frstr, nb, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTEQSB c c call apteqsb (nbep, nep, np, frstr, nb, nerr) c c Version: apteqsb Updated 1990 August 7 15:50. c apteqsb Originated 1990 August 7 15:50. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To randomly sample np bin indices nb from the array nbep(n), c n = 1, nep, where each value of n has equal probability. 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 apteqin to convert a piecewise linear probability c distribution to the equal-probability table required by apteqsb. c c Timing: For a test problem with 1000 bins, 1000 samples, the cpu time c was 740 microseconds unstriated, 829 microseconds striated. c c Input: nbep, nep, np, frstr. c c Output: nb, nerr. c c Glossary: 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 nb Output The randomly sampled bin numbers. Size np. c c nbep Input A bin index. Size nep. c c nep Input The number of values of nbep. c c nerr Output Indicates an input error, if not 0. c 1 if nep 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---- Randomly sampled bin index. dimension nb (1) c---- Equally probable bin indices. dimension nbep (1) c.... Local variables. c---- Index in internal array. common /lapteqsb/ n c---- First index of subset of data. common /lapteqsb/ n1 c---- Last index of subset of data. common /lapteqsb/ n2 c---- Index in external array. common /lapteqsb/ nn c---- Size of current subset of data. common /lapteqsb/ ns c---- Number of striated samples. common /lapteqsb/ nstrd c---- Number of unstriated samples. common /lapteqsb/ nunst c---- Sampled value of bin index. common /lapteqsb/ nx c---- Random numbers in range 0.0 to 1.0. common /lapteqsb/ ranfp (64) c---- Striated random number, 0.0 to 1.0. common /lapteqsb/ ranst cbugc***DEBUG begins. cbugc---- First index with same nbep. cbug common /lapteqsb/ nfirst cbugc---- Last index with same nbep. cbug common /lapteqsb/ nlast cbugc---- Fraction of samples in bin. cbug common /lapteqsb/ fb cbugc---- Index in loop over bins. cbug common /lapteqsb/ nbin cbugc---- Max value of nbeq. cbug common /lapteqsb/ nbinmx cbug 9901 format (/ 'apteqsb sampling from equal-probability bins.' / cbug & 'np=',i5,' frstr=',1pe22.14) cbug 9902 format (i5,' to',i5,' nbep=',i5) cbug write ( 3, 9901) np, frstr cbug nfirst = 1 cbug do 170 n = 1, nep cbug if (nbep(n) .ne. nbep(nfirst)) then cbug nlast = n - 1 cbug write ( 3, 9902) nfirst, nlast, nbep(nfirst) cbug nfirst = n cbug endif cbug 170 continue cbug nlast = nep cbug write ( 3, 9902) nfirst, nlast, nbep(nfirst) cbugc***DEBUG ends. c.... initialize. nerr = 0 c.... Test for input errors. if (nep .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 nx = 1.0 + ranfp(n) * nep nb(nn) = nbep(nx) 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 nx = 1.0 + ranst * nep nb(nn) = nbep(nx) 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 9904 format (/ 'apteqsb results:' / cbug & ' nunst=',i5,' nstrd=',i5) cbug 9905 format (4(' n=',i5,' nb=',i5)) cbug write ( 3, 9904) nunst, nstrd cbug write ( 3, 9905) (n, nb(n), n = 1, np) cbug cbug 9906 format (/ ' summary by bins:') cbug 9907 format (i5,' samples=',i5,' fraction=',1pe22.14) cbug write ( 3, 9906) cbug nbinmx = 0 cbug do 180 nbin = 1, nep cbug nbinmx = max0 (nbinmx, nbep(nbin)) cbug 180 continue cbugc---- Loop over bins. cbug do 190 nbin = 1, nbinmx 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 if (ns .ge. 1) then cbug write ( 3, 9907) nbin, ns, fb cbug endif cbugc---- End of loop over bins. cbug 190 continue cbugc***DEBUG ends. 210 return c.... End of subroutine apteqsb. (+1 line.) end UCRL-WEB-209832