]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:57 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | SUBROUTINE FUNPRE (FUNC,XFCUM,X2LOW,X2HIGH) | |
11 | C F. JAMES, MAY, 1976 | |
12 | C MODIFIED OCT, 1980 TO ADD PRINTOUT OF INTEGRAL | |
13 | C MODIFIED DEC., 1980 TO DELETE LEADING AND TRAILING | |
14 | C RANGES OF X WHERE FUNCTION IS ZERO. | |
15 | C MODIFIED JUNE,1982 TO FIX POSSIBLE INFINITE LOOP. | |
16 | C | |
17 | C PREPARES THE USER FUNCTION "FUNC" FOR FUNRAN | |
18 | C BY FINDING THE PERCENTILES | |
19 | C (IN EFFECT, INVERTING THE CUMULATIVE DISTRIBUTION) | |
20 | EXTERNAL FUNC | |
21 | COMMON/FUNINT/TFTOT | |
22 | DIMENSION XFCUM(100) | |
23 | DATA NBINS / 99/ | |
24 | DATA NZ / 10/ | |
25 | DATA MAXZ / 20/ | |
26 | DATA NITMAX / 6 / | |
27 | DATA IFUNC/0/ | |
28 | IFUNC = IFUNC + 1 | |
29 | C FIND MACHINE ACCURACY | |
30 | COMP1 = 1.0 | |
31 | DO 200 I= 1, 100 | |
32 | COMP1 = COMP1*0.5 | |
33 | COMP2 = 1.0 - COMP1 | |
34 | IF(COMP2 .EQ. 1.0) GOTO 210 | |
35 | 200 CONTINUE | |
36 | COMP1 = 1.0E-10 | |
37 | 210 PRECIS = COMP1 | |
38 | C FIND RANGE WHERE FUNCTION IS NON-ZERO. | |
39 | CALL FUNZER(FUNC,X2LOW,X2HIGH,XLOW,XHIGH) | |
40 | XRANGE = XHIGH-XLOW | |
41 | IF(XRANGE .LE. 0.) GOTO 900 | |
42 | RTEPS = MAX(0.0001,PRECIS*10.) | |
43 | TFTOT = GAUSS(FUNC,XLOW,XHIGH,RTEPS) | |
44 | C PRINT OUT VALUE OF NORMALIZATION INTEGRAL | |
45 | WRITE(6,1003) IFUNC,XLOW,XHIGH,TFTOT | |
46 | RTEPS = 0.001 | |
47 | IF(TFTOT .LE. 0.) GOTO 900 | |
48 | TPCTIL = TFTOT/NBINS | |
49 | TZ = TPCTIL/NZ | |
50 | TZMAX = TZ * 2. | |
51 | XFCUM(1) = XLOW | |
52 | XFCUM(NBINS+1) = XHIGH | |
53 | X = XLOW | |
54 | F = FUNC(X) | |
55 | IF(F .LT. 0.) GOTO 900 | |
56 | NBINM1 = NBINS - 1 | |
57 | C LOOP OVER BINS (HUNDREDTH PERCENTILES) | |
58 | DO 600 IBIN = 1, NBINM1 | |
59 | TCUM = 0. | |
60 | X1 = X | |
61 | F1 = F | |
62 | DXMAX = (XHIGH -X) / NZ | |
63 | FMIN = TZ/DXMAX | |
64 | FMINZ = FMIN | |
65 | C LOOP OVER TRAPEZOIDS WITHIN A SUPPOSED PERCENTILE | |
66 | DO 500 IZ= 1, MAXZ | |
67 | XINCR = TZ/MAX(F1,FMIN,FMINZ) | |
68 | 350 X = X1 + XINCR | |
69 | F = FUNC(X) | |
70 | IF(F .LT. 0.) GOTO 900 | |
71 | TINCR = (X-X1) * 0.5 * (F+F1) | |
72 | IF(TINCR .LT. TZMAX) GOTO 370 | |
73 | XINCR = XINCR * 0.5 | |
74 | GOTO 350 | |
75 | 370 CONTINUE | |
76 | TCUM = TCUM + TINCR | |
77 | IF(TCUM .GE. TPCTIL*0.99) GOTO 520 | |
78 | FMINZ = TZ*F/ (TPCTIL-TCUM) | |
79 | F1 = F | |
80 | X1 = X | |
81 | 500 CONTINUE | |
82 | WRITE(6,2000) | |
83 | 2000 FORMAT('0 FAILURE TO FIND TRAPEZOID HELP') | |
84 | C END OF TRAPEZOID LOOP | |
85 | C ADJUST TRAPEZOID INTEGRAL BY GAUSS WITH NEWTON CORRECTION | |
86 | 520 CONTINUE | |
87 | X1 = XFCUM(IBIN) | |
88 | XBEST = X | |
89 | DTBEST = TPCTIL | |
90 | TPART = TPCTIL | |
91 | C ALLOW FOR MAXIMUM NITMAX MORE ITERATIONS ON GAUSS | |
92 | DO 550 IHOME= 1, NITMAX | |
93 | XINCR = (TPCTIL-TPART) / MAX(F,FMIN) | |
94 | 535 X = XBEST + XINCR | |
95 | X2 = X | |
96 | TPART2 = GAUSS(FUNC,X1,X2,RTEPS) | |
97 | DTPAR2 = TPART2-TPCTIL | |
98 | DTABS = ABS(DTPAR2) | |
99 | IF(ABS(XINCR) .LT. PRECIS) GOTO 545 | |
100 | IF(DTABS .LT. DTBEST) GOTO 545 | |
101 | XINCR = XINCR * 0.5 | |
102 | GOTO 535 | |
103 | 545 DTBEST = DTABS | |
104 | XBEST = X | |
105 | IF(DTABS .LT. RTEPS*TPCTIL) GOTO 580 | |
106 | TPART = TPART2 | |
107 | F = FUNC(X) | |
108 | 550 CONTINUE | |
109 | IHOME = NITMAX | |
110 | C | |
111 | 580 CONTINUE | |
112 | XFCUM(IBIN+1) = X | |
113 | F = FUNC(X) | |
114 | IF(F .LT. 0.) GOTO 900 | |
115 | 600 CONTINUE | |
116 | C END OF LOOP OVER BINS | |
117 | X1 = XFCUM(NBINS) | |
118 | X2 = XHIGH | |
119 | TPART = GAUSS(FUNC,X1,X2,RTEPS) | |
120 | ABERR = ABS(TPART-TPCTIL)/TFTOT | |
121 | WRITE(6,1001) IFUNC,ABERR | |
122 | IF(ABERR .GT. RTEPS) WRITE(6,1002) | |
123 | RETURN | |
124 | 900 WRITE(6,1000) X,F,XLOW,XHIGH | |
125 | RETURN | |
126 | 1000 FORMAT('0FUNPRE FINDS NEGATIVE FUNCTION VALUE OR RANGE OF X'/ | |
127 | + ,' X=',E15.6,', F=',E15.6,20X,'XLOW=',E15.6,' XHIGH=',E15.6/) | |
128 | 1001 FORMAT(' SUBROUTINE FUNPRE HAS PREPARED USER FUNCTION NUMBER',I4, | |
129 | + ' FOR FUNRAN'/' MAXIMUM RELATIVE ERROR IN CUMULATIVE', | |
130 | + 'DISTRIBUTION WILL BE',E15.5) | |
131 | 1002 FORMAT('+',80X,'WARNING,THIS MAY BE TOO BIG'//) | |
132 | 1003 FORMAT('0SUBROUTINE FUNPRE FINDS THE INTEGRAL OF USER FUNCTION', | |
133 | + I2,' FROM ',E12.5,' TO ',E12.5,' IS ',E14.6) | |
134 | END |