subroutine aptsliv (xa, pa, xb, pb, np, xran, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTSLIV c c call aptsliv (xa, pa, xb, pb, np, xran, nerr) c c Version: aptsliv Updated 1990 December 3 14:20. c aptsliv Originated 1990 February 6 11:40. 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 np linear distribution c functions having probabilities pa at xa, and pb at xb. c Flag nerr indicates any input error. c c If both pa and pb are non-negative, the expected value of x is c= fa * xa + fb * xb, c where fa = (2.0 * pa + pb) / (3.0 * (pa + pb)) c and fb = (pa + 2.0 * pb) / (3.0 * (pa + pb)). c c If pa is negative, and pb is positive, then c no values of x will be sampled between xa and c xa' = (xa * pb - xb * pa) / (pb - pa), and c = fa' * xa + fb' * xb, c where fa' = pb / (3.0 * (pb - pa)) c and fb' = (2.0 * pb - 3.0 * pa) / (3.0 * (pb - pa)). c c If pa is positive, and pb is negative, then c no values of x will be sampled between xb and c xb' = (xb * pa - xa * pb) / (pa - pb), and c = fa' * xa + fb' * xb, c where fa' = (2.0 * pa - 3.0 * pb) / (3.0 * (pa - pb)) c and fb' = pa / (3.0 * (pa - pb)). c c If both pa and pb are non-positive, then the returned value c of xran will be -1.e+99. c c Input: xa, pa, xb, pb, np. c c Output: xran, nerr. c c Glossary: c c nerr Output Indicates an input error, if not 0. c 1 if np is not positive. c c np Input Size of arrays xa, pa, xb, pb, xran. c c pa, pb Input Relative probabilities of random variables xa and xb, c respectively. Probability p(x) is linear in x. c The values of pa and pb need not be normalized. c A value of x will be sampled only from the positive c range of each probability distribution function. c If neither pa nor pb is positive, the value of c xran will be -1.e+99. Size np. c To find pa and pb, when is given: c xma = (2.0 * xa + xb) / 3.0, c xmb = (xa + 2.0 * xb) / 3.0. c For between xa and xma, pb is negative, and: c pb = -((xma - ) / ( - xa)) * pa. c For between xma and xmb: c ( - xma) * pa = (xmb - ) * pb. c For between xmb and xb, pa is negative, and: c pa = -(( - xmb) / (xb - )) * pb. c c xa, xb Input Values of random variable x with relative probabilities c pa and pb, respectively. Size np. c c xran Output Randomly sampled value of x in range from xa to xb. c Value returned is -1.e+99 if neither pa nor pb c is positive. Size np. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Probability at xa. dimension pa (1) c---- Probability at xb. dimension pb (1) c---- Value of x with probability pa. dimension xa (1) c---- Value of x with probability pb. dimension xb (1) c---- Randomly sampled value of x. dimension xran (1) c.... Local variables. c---- A very big number. common /laptsliv/ big c---- A very small number. common /laptsliv/ fuz c---- Index in local array. common /laptsliv/ n c---- First index of subset of data. common /laptsliv/ n1 c---- Last index of subset of data. common /laptsliv/ n2 c---- Size of current subset of data. common /laptsliv/ ns c---- Index in external array. common /laptsliv/ nn c---- Probability of x1. common /laptsliv/ p1 c---- Value of p(x) at xxa. common /laptsliv/ ppa (64) c---- Value of p(x) at xxb. common /laptsliv/ ppb (64) c---- Sum of ppa and ppb. common /laptsliv/ psum c---- First random number. common /laptsliv/ ranf1 (64) c---- Second random number. common /laptsliv/ ranf2 (64) c---- First random number ranf1. common /laptsliv/ u1 c---- Value 1.0 - u1. common /laptsliv/ um c---- First random value of x. common /laptsliv/ x1 c---- Second random value of x. common /laptsliv/ x2 c---- Maximum of xa, xb. common /laptsliv/ xmax c---- Minimum of xa, xb. common /laptsliv/ xmin c---- Adjusted value of xa. common /laptsliv/ xxa (64) c---- Adjusted value of xb. common /laptsliv/ xxb (64) c---- Value of x where p(x) = 0.0. common /laptsliv/ xzero cbugc***DEBUG begins. cbugc---- Mean value of x. cbug common /laptsliv/ xmean cbugc---- Standard deviation from mean x. cbug common /laptsliv/ xdev cbug 9901 format (/ 'aptsliv sampling from a linear distribution. np=',i3 / cbug & (i3,' xa= ',1pe22.14,' pa= ',1pe22.14 / cbug & ' xb= ',1pe22.14,' pb= ',1pe22.14)) cbug write (3, 9901) np, (n, xa(n), pa(n), xb(n), pb(n), n = 1, np) cbugc***DEBUG ends. c.... initialize. c---- A very big number. big = 1.e+99 c---- A very small number. fuz = 1.e-99 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 data. n1 = 1 n2 = min (np, 64) 110 ns = n2 - n1 + 1 c.... Adjust the ranges of x to include only positive probabilities. c.... If not possible, will return x = -big. c---- Loop over local arrays. do 120 n = 1, ns nn = n + n1 - 1 xzero = (xa(nn) * pb(nn) - xb(nn) * pa(nn)) / & (pb(nn) - pa(nn) + fuz) xmax = amax1 (xa(nn), xb(nn)) xmin = amin1 (xa(nn), xb(nn)) xzero = amax1 (xmin, amin1 (xzero, xmax)) if (pa(nn) .ge. 0.0) then xxa(n) = xa(nn) ppa(n) = pa(nn) else xxa(n) = xzero ppa(n) = 0.0 endif if (pb(nn) .ge. 0.0) then xxb(n) = xb(nn) ppb(n) = pb(nn) else xxb(n) = xzero ppb(n) = 0.0 endif c---- End of loop over local arrays. 120 continue cbugc***DEBUG begins. cbug 9902 format (/ ' Eliminated negative probabilities:' / cbug & i3,' xxa=',1pe22.14,' ppa=',1pe22.14 / cbug & ' xxb=',1pe22.14,' ppb=',1pe22.14) cbugc---- Loop over local arrays. cbug do 125 n = 1, ns cbug nn = n + n1 - 1 cbug if ((xxa(n) .ne. xa(nn)) .or. (xxb(n) .ne. xb(nn))) then cbug write ( 3, 9902) nn, xxa(n), ppa(n), xxb(n), ppb(n) cbug endif cbugc---- End of loop over local arrays. cbug 125 continue cbugc***DEBUG ends. c.... Generate ns pairs of random numbers. c---- Loop over local arrays. do 130 n = 1, ns ranf1(n) = ranf( ) c---- End of loop over local arrays. 130 continue c---- Loop over local arrays. do 140 n = 1, ns ranf2(n) = ranf( ) c---- End of loop over local arrays. 140 continue c.... Sample values of x from the distributions p(x). c---- Loop over local arrays. do 150 n = 1, ns nn = n + n1 - 1 u1 = ranf1(n) um = 1.0 - u1 x1 = um * xxa(n) + u1 * xxb(n) x2 = u1 * xxa(n) + um * xxb(n) p1 = um * ppa(n) + u1 * ppb(n) psum = ppa(n) + ppb(n) if (p1 .ge. ranf2(n) * psum) then xran(nn) = x1 else xran(nn) = x2 endif if (psum .le. 0.0) then xran(nn) = -big endif c---- End of loop over local arrays. 150 continue c.... See if all data subsets are done. c---- Do another subset of data. if (n2 .lt. np) then n1 = n2 + 1 n2 = min (np, n1 + 63) go to 110 endif cbugc***DEBUG begins. cbug 9903 format (/ 'aptsliv results: x') cbug write ( 3, 9903) cbug call aptmean (xran, np, 1.e-11, xmean, xdev, nerr) cbugc***DEBUG ends. 210 return c.... End of subroutine aptsliv. (+1 line.) end UCRL-WEB-209832