MODULE toms519 ! Code converted using TO_F90 by Alan Miller ! Date: 2000-06-15 Time: 00:33:15 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), PARAMETER, PRIVATE :: one = 1.0_dp, zero = 0.0_dp CONTAINS FUNCTION pln(n, s2, b2, s1, b1) RESULT(fn_val) ! COMPUTES PROB FOR BOUNDARIES OF FORM S*Y + B. N IS SAMPLE SIZE, ! S2,B2 ARE PARAMETERS FOR UPPER BOUNDARY. ! CALLS SUBROUTINE RAKK. ! TO CALL DURB, CHANGE AS INDICATED BELOW. IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: s2 REAL (dp), INTENT(IN) :: b2 REAL (dp), INTENT(IN) :: s1 REAL (dp), INTENT(IN) :: b1 REAL (dp) :: fn_val REAL (dp) :: yy(2,202), s(2), b(2), an, ans, c INTEGER :: i, k, n1 ! MAX SIZE FOR N IS 200 s(1) = s1 s(2) = s2 b(1) = b1 b(2) = b2 IF (n <= 0) THEN fn_val = zero RETURN END IF IF (s(1) < zero) THEN s(1) = zero END IF IF (s(2) < zero) THEN b(2) = s(2) + b(2) s(2) = zero END IF n1 = 1 + n an = n DO i = 1, 2 DO k = 1, n1 c = REAL(k-1) / an - b(i) yy(i,k) = -0.5 IF (c > zero) THEN yy(i,k) = 1.1 IF (s(i) > c) yy(i,k) = c / s(i) END IF END DO END DO CALL rakk(n, yy, ans) ! TO CALL DURB CHANGE PRECEDING STATEMENT TO ! CALL DURB (N, YY, ANS) fn_val = ans RETURN END FUNCTION pln SUBROUTINE rakk(n, yy, ans) ! N IS SAMPLE SIZE. YY(2,.),YY(1,.) CONTAIN, RESPECTIVELY, ! UPPER AND LOWER BOUNDARY LEVEL POINTS CORRESPONDING TO 0,1/N,2/N,...,1. ! MUST BE N+1 OF EACH. ROUTINE ASSUMES THEY ARE NONDECREASING. IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: yy(2,202) REAL (dp), INTENT(OUT) :: ans REAL (dp) :: convi, frng(41,2), p, rrngs, sss(41), s, sqrng, sqrrng, t, tf, & tmlt, w(201,2), y0, y1 INTEGER :: i, i1, iext, j, j1, j11, jbot0, jbotp, jtop0, jtop1, jtopp, & k, k0, k01, k0p1, k1, k1p1, kbli, kbloc(41), l, l0, l1, lb, & lb1, ldf, lfg, ll, lt, m1, m2, mbav, mt, n1, nbloc, nnn REAL (dp), PARAMETER :: rng = HUGE(one), rrng = TINY(one) INTEGER, PARAMETER :: maxnb = 40 ! DATA rng, rrng, sqrng, sqrrng, maxnb /1.e292, 1.e-292, 1.e146, 1.e-146, 40/ ! RNG AND RRNG ARE MACHINE DEPENDENT CONSTANTS SELECTED ! SO RRNG=1./RNG >= SMALLEST POSITIVE FLOATING ! NUMBER AND RNG <= LARGEST POSITIVE FLOATING ! NUMBER OF MACHINE. SQRNG AND SQRRNG ARE SQUARE ROOTS OF RNG AND RRNG. ! MAXNB IS THE MAX NUMBER OF BLOCKS INTO WHICH W(.,.) WILL BE PARTITIONED. ! IF MAXNB IS INCREASED REVISE DIMENSION STATEMENT ACCORDINGLY. sqrng = SQRT(rng) sqrrng = SQRT(rrng) n1 = n + 1 tf = one nnn = 1 iext = 1 m2 = 1 m1 = 2 w(1,m2) = sqrng nbloc = 1 kbloc(1) = 1 kbloc(2) = 2 frng(1,m2) = one ans = zero IF (n1 <= 1) RETURN yy(2,n1+1) = one y1 = zero IF (yy(1,1) <= zero .OR. yy(2,n1) >= one .OR. yy(1,n1) <= one .OR. & yy(2,1) >= zero) RETURN jbotp = 1 jtopp = 1 jtop1 = 1 ! ENTRY POINT FOR SUBSEQUENT ITERATIONS. 10 jbot0 = jbotp jtop0 = jtop1 mt = m1 m1 = m2 m2 = mt y0 = y1 ! COMPUTE SUMMATION INDICES, JBOTP,JTOPP FOR NEXT ITERATION. 20 IF (yy(2,jtopp+1) <= y0) THEN jtopp = jtopp + 1 GO TO 20 END IF y1 = MIN(yy(2,jtopp+1), yy(1,jbotp)) IF (y1 < one) THEN 30 IF (yy(1,jbotp) > y1) GO TO 40 jbotp = jbotp + 1 GO TO 30 END IF y1 = one ! SET EXIT FLAG. NEXT ITERATION IS FINAL ONE. iext = 2 jbotp = n1 40 IF (jbotp > jtopp) RETURN ! RETURN WITH ANS=0 SINCE NO PATHS BETWEEN BOUNDARIES. jtop1 = jtopp p = y1 - y0 kbloc(nbloc+1) = jtop1 + 1 DO l = jbot0, jtop1 w(l,m2) = zero END DO l0 = kbloc(nbloc) l1 = jtop0 DO l = l0, l1 IF (w(l,m1) <= sqrrng) THEN jtop0 = l - 1 EXIT END IF END DO 70 IF (jbot0 >= kbloc(2)) THEN ! THE NEXT FEW STATEMENTS ELIMINATE BLOCKS NO LONGER USED. DO i = 1, nbloc frng(i,m1) = frng(i+1,m1) kbloc(i) = kbloc(i+1) END DO nbloc = nbloc - 1 CALL accum(rrng/frng(1,m1),tf,nnn) GO TO 70 END IF kbloc(1) = jbot0 mbav = (maxnb-nbloc+1) / 2 IF (p <= (one/REAL(4*n1))) mbav = 2 mbav = MIN(maxnb, nbloc+mbav) ! MBAV IS UPPER LIMIT ON NUMBER OF BLOCKS DURING NEXT ITERATION. ! MBAV IS SET TO LIMIT NUMBER OF NEW BLOCKS WHEN P IS SMALL. convi = one i1 = 1 j1 = 1 90 i = i + 1 ! Dummy instruction loop320: & DO i = i1, nbloc kbli = kbloc(i) IF (i /= i1) THEN frng(i,m2) = MIN(frng(i,m1),one/(sqrng*w(kbli-1,m2))) convi = (convi/frng(i-j1+1,m1)) * frng(i,m2) IF (j1 /= 1) THEN k1 = kbloc(i) - kbloc(i-j1+2) + 1 k0 = kbloc(i-1) - kbloc(i-j1+1) + 1 IF (k1 < k0) THEN k1p1 = k1 + 1 DO k = k1p1, k0 convi = convi * REAL(k) / p END DO ELSE IF (k1 > k0) THEN k0p1 = k0 + 1 DO k = k0p1, k1 convi = convi * p / REAL(k) END DO END IF END IF END IF IF (convi <= rrng) THEN DO j = j1, j11 ! CHANGE UPPER LIMIT IF EARLY EXIT FROM PREVIOUS LOOP convi = frng(i,m2) * rng k0 = MAX(1, kbloc(i-1) - MIN(jtop0+1, kbloc(i-j+1))+2) k1 = kbloc(i) - MIN(jtop0+1, kbloc(i-j+1)) + 1 DO k = k0, k1 convi = convi * p / REAL(k) END DO convi = sss(j) * convi IF (convi > rrng) GO TO 170 END DO SELECT CASE ( iext ) CASE ( 1) GO TO 10 CASE ( 2) GO TO 360 END SELECT 170 j1 = j + 1 END IF DO j = j1, i lfg = 1 ! LFG=2 IF NEXT ITERATION ON J BECOMES NECESSARY. lb1 = kbloc(i) k01 = kbloc(i) - kbloc(i-j+1) + 1 s = convi IF (j /= j1) THEN IF (tmlt == zero) THEN CYCLE loop320 ELSE IF (tmlt < zero) THEN s = -tmlt * frng(i-j+2,m1) IF (s > rrng) GO TO 200 END IF s = tmlt * (frng(i-j+2,m1)*rng) END IF 200 k0 = MAX(kbloc(i)-MIN(kbloc(i-j+2)-1,jtop0)+1, 1) k1 = kbloc(i+1) - kbloc(i-j+1) tmlt = zero rrngs = one IF (s < zero) rrngs = -rrng sss(j) = s DO k = k0, k1 ldf = 1 - k lb = MAX(kbloc(i),kbloc(i-j+1)-ldf,lb1) LT = MIN(kbloc(i+1)-1, MIN(jtop0,kbloc(i-j+2)-1)-ldf) IF (lb > lb1) lfg = 2 IF (lb <= LT) THEN DO l = lb, LT IF (w(l,m2) /= (w(l+ldf,m1)*rrngs)*s + w(l,m2)) THEN lb1 = l ! WHEN THE SUMMANDS BECOME SO SMALL THAT THE VALUE OF W(L,M2) IS UNCHANGED, ! LB1 IS INCREASED. THIS PREVENTS FURTHER ADDITIONS TO SUCH W(L,M2). GO TO 230 END IF IF (w(l,m2) <= rrng) THEN IF (k-k01 > 0) THEN GO TO 220 ELSE GO TO 280 END IF END IF END DO END IF lb1 = LT + 1 IF (lb1 /= kbloc(i+1) .OR. k <= k01) GO TO 280 220 j11 = j IF (w(kbli,m2) <= rrng) lfg = 2 SELECT CASE ( lfg ) CASE ( 1) CYCLE loop320 CASE ( 2) GO TO 300 END SELECT 230 IF (s < zero) THEN DO l = lb1, LT w(l,m2) = (w(l+ldf,m1)*rrngs) * s + w(l,m2) END DO ELSE IF (s > zero) THEN DO l = lb1, LT w(l,m2) = w(l+ldf,m1) * s + w(l,m2) END DO ELSE GO TO 300 END IF 280 t = s s = s * p / REAL(k) IF (s == 0) THEN IF (t <= zero) EXIT s = -(t*rng) * p / REAL(k) rrngs = -rrng END IF ! THE NEGATIVE BIT IN S,TMLT IS USED TO ! INDICATE ITS VALUE IS RNG TIMES THE USUAL VALUE. IF (k == k01) tmlt = s ! THE ITERATION ON K MUST CONTINUE UNTIL TMLT IS SET. ! TMLT USED IN NEXT ITERATION ON J. END DO 300 j11 = j END DO END DO loop320 ! CREATE NEW BLOCK IF NECESSARY. IF (nbloc >= mbav) THEN SELECT CASE ( iext ) CASE ( 1) GO TO 10 CASE ( 2) GO TO 360 END SELECT END IF l1 = kbloc(nbloc+1) - kbloc(nbloc) DO l = 1, l1 ll = kbloc(nbloc+1) - l IF (w(ll,m2) >= sqrrng) EXIT END DO IF (ll == jtop1) THEN SELECT CASE ( iext ) CASE ( 1) GO TO 10 CASE ( 2) GO TO 360 END SELECT END IF kbloc(nbloc+2) = kbloc(nbloc+1) kbloc(nbloc+1) = ll + 1 nbloc = nbloc + 1 l0 = kbloc(nbloc) DO l = l0, jtop1 w(l,m2) = zero END DO i1 = nbloc frng(nbloc,m2) = one convi = zero GO TO 90 360 CALL accum(w(n1,m2), tf, nnn) CALL accum(sqrrng,tf,nnn) IF (nbloc /= 1) THEN DO i = 2, nbloc CALL accum(rrng/frng(i,m2), tf, nnn) END DO END IF ans = tf IF (n < nnn) THEN DO j = n1, nnn ans = ans / REAL(j) END DO ELSE IF (n > nnn) THEN nnn = nnn + 1 DO j = nnn, n ans = ans * REAL(j) END DO END IF RETURN END SUBROUTINE rakk SUBROUTINE accum(f, tf, nnn) ! THIS SUBROUTINE ANCILLARY TO RAKK AND ACCUMULATES THE ! MULTIPLIERS USED TO KEEP W(.,.) IN RANGE. IMPLICIT NONE REAL (dp), INTENT(IN) :: f REAL (dp), INTENT(IN OUT) :: tf INTEGER, INTENT(IN OUT) :: nnn tf = tf * f DO IF (tf > one .OR. tf == zero) EXIT nnn = nnn + 1 tf = tf * REAL(nnn) END DO RETURN END SUBROUTINE accum SUBROUTINE durb(n,yy,ans) ! DURBINS METHOD. J.APPL.PROB., 8(1971), PP431-453, ! FORMULAS (36),(37),(40). ! N IS SAMPLE SIZE. YY(2,.),YY(1,.) ! CONTAIN, RESPECTIVELY, UPPER AND LOWER BOUNDARY ! LEVEL POINTS CORRESPONDING TO 0,1/N,2/N,...,1. MUST ! BE N+1 OF EACH. ROUTINE ASSUMES THEY ARE NONDECREASING. IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: yy(2,202) REAL (dp), INTENT(OUT) :: ans REAL (dp) :: alpha(201), beta(201) ! DIMENSION STATEMENT PERMITS MAX SAMPLE SIZE OF 200 INTEGER :: i, i1, ip, j, k, km1, n1 REAL (dp) :: comb n1 = n + 1 ans = zero IF (yy(1,1) <= zero .OR. yy(2,n1) >= one .OR. yy(1,n1) <= one .OR. & yy(2,1) >= zero) RETURN ip = 1 DO i = 1, n1 IF (yy(2,i) >= zero) THEN ip = i EXIT END IF yy(2,i) = zero END DO DO i = 1, n1 IF (yy(1,i) > one) THEN yy(1,i) = one END IF END DO alpha(1) = zero beta(1) = one loop90: DO k = 2, n1 alpha(k) = yy(2,k)**(k-1) IF (ip <= (k-1)) THEN comb = one i1 = k - ip DO i = 1, i1 comb = comb * REAL(k-i) / REAL(i) ! COMB IS BINOMIAL COEFFICIENT K-1 OVER I-1. alpha(k) = alpha(k) - alpha(k-i) * (yy(2,k) - yy(2,k-i))**i * comb END DO END IF comb = one DO j = 1, n1 IF (yy(2,k) <= yy(1,j)) EXIT alpha(k) = alpha(k) - beta(j) * (yy(2,k) - yy(1,j))**(k-j) * comb ! COMB IS BINOMIAL COEFFICIENT K-1 OVER J-1. comb = comb * REAL(k-j) / REAL(j) END DO beta(k) = yy(1,k) ** (k-1) IF (ip <= k) THEN comb = one i1 = k - ip + 1 DO i = 1, i1 beta(k) = beta(k) - alpha(k-i+1) * (yy(1,k) - yy(2,k-i+1))**(i-1) * comb ! COMB IS BINOMIAL COEFFICIENT K-1 OVER I-1. comb = comb * REAL(k-i) / REAL(i) END DO comb = one END IF km1 = k - 1 DO j = 1, km1 IF (yy(1,j) == one) CYCLE loop90 ! COMB IS BINOMIAL COEFFICIENT K-1 OVER J-1. beta(k) = beta(k) - beta(j) * (yy(1,k) - yy(1,j))**(k-j) * comb comb = comb * REAL(k-j) / REAL(j) END DO END DO loop90 comb = one i1 = n1 - ip + 1 DO i = 1, i1 ans = ans + (one - yy(2,n1-i+1))**(i-1) * alpha(n1-i+1) * comb comb = comb * REAL(n1-i) / REAL(i) END DO comb = one DO j = 1, n1 IF (yy(1,j) == one) EXIT ans = ans + (one - yy(1,j))**(n1-j) * beta(j) * comb ! COMB IS BINOMIAL COEFFICIENT N OVER J-1. comb = comb * REAL(n1-j) / REAL(j) END DO ! FOR PROB OF COMPLEMENTARY EVENT CHANGE NEXT STATEMENT TO ! CONTINUE ans = one - ans RETURN END SUBROUTINE durb END MODULE toms519