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