subroutine apteqxs (xep, nep, np, frstr, x, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTEQXS c c call apteqxs (xep, nep, np, frstr, x, nerr) c c Version: apteqxs Updated 1990 August 14 10:00. c apteqxs Originated 1990 August 14 10:00. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To randomly sample np values of x from the nep equal-probability c bins bounded by the nep + 1 values of xep. Within each bin, c values of x are sampled uniformly. 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 X VALUES. USE ONLY WHEN SUCH CORRELATIONS ARE ACCEPTABLE. c c Note: Use subroutine apteqxn to convert a piecewise linear probability c distribution to the equal-probability table required by apteqxs. c Use subroutine apteqin to convert a piecewise linear probability c distribution to an equal-probability table of bin indices, to c allow sampling from the original linear distribution within each c bin, using subroutines apteqsb, aptalsh, aptalsl. c c Timing: For a test problem with 128 bins, 1000 samples, the cpu time c was 972 microseconds unstriated, 1109 microseconds striated. c c Input: xep, nep, np, frstr. c c Output: x, 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 nep Input The number of equal-probability bins. 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 x. Number of samples. c c x Output The randomly sampled values of the random variable. c Size np. Sampled uniformly in each bin. c c xep Input A bin boundary. The n'th bin is bounded by xep(n) and c xep(n+1). Size nep + 1. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Randomly sampled values. dimension x (1) c---- Boundaries of equal-probability bins. dimension xep (1) c.... Local variables. c---- Index in internal array. common /lapteqxs/ n c---- First index of subset of data. common /lapteqxs/ n1 c---- Last index of subset of data. common /lapteqxs/ n2 c---- Index in external array. common /lapteqxs/ nn c---- Size of current subset of data. common /lapteqxs/ ns c---- Number of striated samples. common /lapteqxs/ nstrd c---- Number of unstriated samples. common /lapteqxs/ nunst c---- Sampled bin index. common /lapteqxs/ nx c---- Sampled bin index + 1. common /lapteqxs/ nx1 c---- Random numbers in range 0.0 to 1.0. common /lapteqxs/ ranfp (64) c---- Striated random number, 0.0 to 1.0. common /lapteqxs/ ranst c---- Floating point sampled bin index. common /lapteqxs/ xn cbugc***DEBUG begins. cbugc---- Standard deviation of mean. cbug common /lapteqxs/ xdev cbugc---- Mean value of x. cbug common /lapteqxs/ xmean cbug 9901 format (/ 'apteqxs sampling from equal-probability bins.' / cbug & 'np=',i5,' frstr=',1pe22.14) cbug 9902 format (3(i4,1pe22.14)) cbug write ( 3, 9901) np, frstr cbug write ( 3, 9902) (n, xep(n), n = 1, nep + 1) 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 bins and x values (both uniformly). c---- Loop over samples. do 130 n = 1, ns nn = n + n1 - 1 xn = 1.0 + ranfp(n) * nep nx = xn x(nn) = xep(nx) + (xn - nx) * (xep(nx+1) - xep(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 bins and x values, 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 * nep nx = xn c---- For vectorization. ??? nx1 = xn + 1.0 x(nn) = xep(nx) + (xn - nx) * (xep(nx1) - xep(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 (/ 'apteqxs results:' / cbug & ' nunst=',i5,' nstrd=',i5) cbug write ( 3, 9904) nunst, nstrd cbug call aptmean (x, np, 1.e-11, xmean, xdev, nerr) cbugc***DEBUG ends. 210 return c.... End of subroutine apteqxs. (+1 line.) end UCRL-WEB-209832