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