|
2329
|
1 SUBROUTINE DQAGP(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,RESULT,ABSERR, |
|
|
2 * NEVAL,IER,LENIW,LENW,LAST,IWORK,WORK) |
|
|
3 C***BEGIN PROLOGUE DQAGP |
|
|
4 C***DATE WRITTEN 800101 (YYMMDD) |
|
|
5 C***REVISION DATE 830518 (YYMMDD) |
|
|
6 C***CATEGORY NO. H2A2A1 |
|
|
7 C***KEYWORDS AUTOMATIC INTEGRATOR, GENERAL-PURPOSE, |
|
|
8 C SINGULARITIES AT USER SPECIFIED POINTS, |
|
|
9 C EXTRAPOLATION, GLOBALLY ADAPTIVE |
|
|
10 C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV - K.U.LEUVEN |
|
|
11 C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN |
|
|
12 C***PURPOSE THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN |
|
|
13 C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), |
|
|
14 C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY |
|
|
15 C BREAK POINTS OF THE INTEGRATION INTERVAL, WHERE LOCAL |
|
|
16 C DIFFICULTIES OF THE INTEGRAND MAY OCCUR (E.G. |
|
|
17 C SINGULARITIES, DISCONTINUITIES), ARE PROVIDED BY THE USER. |
|
|
18 C***DESCRIPTION |
|
|
19 C |
|
|
20 C COMPUTATION OF A DEFINITE INTEGRAL |
|
|
21 C STANDARD FORTRAN SUBROUTINE |
|
|
22 C DOUBLE PRECISION VERSION |
|
|
23 C |
|
|
24 C PARAMETERS |
|
|
25 C ON ENTRY |
|
3135
|
26 C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND |
|
2329
|
27 C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE |
|
|
28 C DECLARED E X T E R N A L IN THE DRIVER PROGRAM. |
|
|
29 C |
|
|
30 C A - DOUBLE PRECISION |
|
|
31 C LOWER LIMIT OF INTEGRATION |
|
|
32 C |
|
|
33 C B - DOUBLE PRECISION |
|
|
34 C UPPER LIMIT OF INTEGRATION |
|
|
35 C |
|
|
36 C NPTS2 - INTEGER |
|
|
37 C NUMBER EQUAL TO TWO MORE THAN THE NUMBER OF |
|
|
38 C USER-SUPPLIED BREAK POINTS WITHIN THE INTEGRATION |
|
|
39 C RANGE, NPTS.GE.2. |
|
|
40 C IF NPTS2.LT.2, THE ROUTINE WILL END WITH IER = 6. |
|
|
41 C |
|
|
42 C POINTS - DOUBLE PRECISION |
|
|
43 C VECTOR OF DIMENSION NPTS2, THE FIRST (NPTS2-2) |
|
|
44 C ELEMENTS OF WHICH ARE THE USER PROVIDED BREAK |
|
|
45 C POINTS. IF THESE POINTS DO NOT CONSTITUTE AN |
|
|
46 C ASCENDING SEQUENCE THERE WILL BE AN AUTOMATIC |
|
|
47 C SORTING. |
|
|
48 C |
|
|
49 C EPSABS - DOUBLE PRECISION |
|
|
50 C ABSOLUTE ACCURACY REQUESTED |
|
|
51 C EPSREL - DOUBLE PRECISION |
|
|
52 C RELATIVE ACCURACY REQUESTED |
|
|
53 C IF EPSABS.LE.0 |
|
|
54 C AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), |
|
|
55 C THE ROUTINE WILL END WITH IER = 6. |
|
|
56 C |
|
|
57 C ON RETURN |
|
|
58 C RESULT - DOUBLE PRECISION |
|
|
59 C APPROXIMATION TO THE INTEGRAL |
|
|
60 C |
|
|
61 C ABSERR - DOUBLE PRECISION |
|
|
62 C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, |
|
|
63 C WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT) |
|
|
64 C |
|
|
65 C NEVAL - INTEGER |
|
|
66 C NUMBER OF INTEGRAND EVALUATIONS |
|
|
67 C |
|
|
68 C IER - INTEGER |
|
|
69 C IER = 0 NORMAL AND RELIABLE TERMINATION OF THE |
|
|
70 C ROUTINE. IT IS ASSUMED THAT THE REQUESTED |
|
|
71 C ACCURACY HAS BEEN ACHIEVED. |
|
|
72 C IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE. |
|
|
73 C THE ESTIMATES FOR INTEGRAL AND ERROR ARE |
|
|
74 C LESS RELIABLE. IT IS ASSUMED THAT THE |
|
|
75 C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED. |
|
|
76 C ERROR MESSAGES |
|
|
77 C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED |
|
|
78 C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE |
|
|
79 C SUBDIVISIONS BY INCREASING THE VALUE OF |
|
|
80 C LIMIT (AND TAKING THE ACCORDING DIMENSION |
|
|
81 C ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF |
|
|
82 C THIS YIELDS NO IMPROVEMENT IT IS ADVISED |
|
|
83 C TO ANALYZE THE INTEGRAND IN ORDER TO |
|
|
84 C DETERMINE THE INTEGRATION DIFFICULTIES. IF |
|
|
85 C THE POSITION OF A LOCAL DIFFICULTY CAN BE |
|
|
86 C DETERMINED (I.E. SINGULARITY, |
|
|
87 C DISCONTINUITY WITHIN THE INTERVAL), IT |
|
|
88 C SHOULD BE SUPPLIED TO THE ROUTINE AS AN |
|
|
89 C ELEMENT OF THE VECTOR POINTS. IF NECESSARY |
|
|
90 C AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR |
|
|
91 C MUST BE USED, WHICH IS DESIGNED FOR |
|
|
92 C HANDLING THE TYPE OF DIFFICULTY INVOLVED. |
|
|
93 C = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS |
|
|
94 C DETECTED, WHICH PREVENTS THE REQUESTED |
|
|
95 C TOLERANCE FROM BEING ACHIEVED. |
|
|
96 C THE ERROR MAY BE UNDER-ESTIMATED. |
|
|
97 C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS |
|
|
98 C AT SOME POINTS OF THE INTEGRATION |
|
|
99 C INTERVAL. |
|
|
100 C = 4 THE ALGORITHM DOES NOT CONVERGE. |
|
|
101 C ROUNDOFF ERROR IS DETECTED IN THE |
|
|
102 C EXTRAPOLATION TABLE. |
|
|
103 C IT IS PRESUMED THAT THE REQUESTED |
|
|
104 C TOLERANCE CANNOT BE ACHIEVED, AND THAT |
|
|
105 C THE RETURNED RESULT IS THE BEST WHICH |
|
|
106 C CAN BE OBTAINED. |
|
|
107 C = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR |
|
|
108 C SLOWLY CONVERGENT. IT MUST BE NOTED THAT |
|
|
109 C DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE |
|
|
110 C OF IER.GT.0. |
|
|
111 C = 6 THE INPUT IS INVALID BECAUSE |
|
|
112 C NPTS2.LT.2 OR |
|
|
113 C BREAK POINTS ARE SPECIFIED OUTSIDE |
|
|
114 C THE INTEGRATION RANGE OR |
|
|
115 C (EPSABS.LE.0 AND |
|
|
116 C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28)) |
|
|
117 C RESULT, ABSERR, NEVAL, LAST ARE SET TO |
|
|
118 C ZERO. EXEPT WHEN LENIW OR LENW OR NPTS2 IS |
|
|
119 C INVALID, IWORK(1), IWORK(LIMIT+1), |
|
|
120 C WORK(LIMIT*2+1) AND WORK(LIMIT*3+1) |
|
|
121 C ARE SET TO ZERO. |
|
|
122 C WORK(1) IS SET TO A AND WORK(LIMIT+1) |
|
|
123 C TO B (WHERE LIMIT = (LENIW-NPTS2)/2). |
|
|
124 C |
|
|
125 C DIMENSIONING PARAMETERS |
|
|
126 C LENIW - INTEGER |
|
|
127 C DIMENSIONING PARAMETER FOR IWORK |
|
|
128 C LENIW DETERMINES LIMIT = (LENIW-NPTS2)/2, |
|
|
129 C WHICH IS THE MAXIMUM NUMBER OF SUBINTERVALS IN THE |
|
|
130 C PARTITION OF THE GIVEN INTEGRATION INTERVAL (A,B), |
|
|
131 C LENIW.GE.(3*NPTS2-2). |
|
|
132 C IF LENIW.LT.(3*NPTS2-2), THE ROUTINE WILL END WITH |
|
|
133 C IER = 6. |
|
|
134 C |
|
|
135 C LENW - INTEGER |
|
|
136 C DIMENSIONING PARAMETER FOR WORK |
|
|
137 C LENW MUST BE AT LEAST LENIW*2-NPTS2. |
|
|
138 C IF LENW.LT.LENIW*2-NPTS2, THE ROUTINE WILL END |
|
|
139 C WITH IER = 6. |
|
|
140 C |
|
|
141 C LAST - INTEGER |
|
|
142 C ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS |
|
|
143 C PRODUCED IN THE SUBDIVISION PROCESS, WHICH |
|
|
144 C DETERMINES THE NUMBER OF SIGNIFICANT ELEMENTS |
|
|
145 C ACTUALLY IN THE WORK ARRAYS. |
|
|
146 C |
|
|
147 C WORK ARRAYS |
|
|
148 C IWORK - INTEGER |
|
|
149 C VECTOR OF DIMENSION AT LEAST LENIW. ON RETURN, |
|
|
150 C THE FIRST K ELEMENTS OF WHICH CONTAIN |
|
|
151 C POINTERS TO THE ERROR ESTIMATES OVER THE |
|
|
152 C SUBINTERVALS, SUCH THAT WORK(LIMIT*3+IWORK(1)),..., |
|
|
153 C WORK(LIMIT*3+IWORK(K)) FORM A DECREASING |
|
|
154 C SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND |
|
|
155 C K = LIMIT+1-LAST OTHERWISE |
|
|
156 C IWORK(LIMIT+1), ...,IWORK(LIMIT+LAST) CONTAIN THE |
|
|
157 C SUBDIVISION LEVELS OF THE SUBINTERVALS, I.E. |
|
|
158 C IF (AA,BB) IS A SUBINTERVAL OF (P1,P2) |
|
|
159 C WHERE P1 AS WELL AS P2 IS A USER-PROVIDED |
|
|
160 C BREAK POINT OR INTEGRATION LIMIT, THEN (AA,BB) HAS |
|
|
161 C LEVEL L IF ABS(BB-AA) = ABS(P2-P1)*2**(-L), |
|
|
162 C IWORK(LIMIT*2+1), ..., IWORK(LIMIT*2+NPTS2) HAVE |
|
|
163 C NO SIGNIFICANCE FOR THE USER, |
|
|
164 C NOTE THAT LIMIT = (LENIW-NPTS2)/2. |
|
|
165 C |
|
|
166 C WORK - DOUBLE PRECISION |
|
|
167 C VECTOR OF DIMENSION AT LEAST LENW |
|
|
168 C ON RETURN |
|
|
169 C WORK(1), ..., WORK(LAST) CONTAIN THE LEFT |
|
|
170 C END POINTS OF THE SUBINTERVALS IN THE |
|
|
171 C PARTITION OF (A,B), |
|
|
172 C WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN |
|
|
173 C THE RIGHT END POINTS, |
|
|
174 C WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN |
|
|
175 C THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS, |
|
|
176 C WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST) |
|
|
177 C CONTAIN THE CORRESPONDING ERROR ESTIMATES, |
|
|
178 C WORK(LIMIT*4+1), ..., WORK(LIMIT*4+NPTS2) |
|
|
179 C CONTAIN THE INTEGRATION LIMITS AND THE |
|
|
180 C BREAK POINTS SORTED IN AN ASCENDING SEQUENCE. |
|
|
181 C NOTE THAT LIMIT = (LENIW-NPTS2)/2. |
|
|
182 C |
|
|
183 C***REFERENCES (NONE) |
|
|
184 C***ROUTINES CALLED DQAGPE,XERROR |
|
|
185 C***END PROLOGUE DQAGP |
|
|
186 C |
|
3136
|
187 DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,POINTS,RESULT,WORK |
|
2329
|
188 INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,NEVAL, |
|
|
189 * NPTS2 |
|
|
190 C |
|
|
191 DIMENSION IWORK(LENIW),POINTS(NPTS2),WORK(LENW) |
|
|
192 C |
|
|
193 EXTERNAL F |
|
|
194 C |
|
|
195 C CHECK VALIDITY OF LIMIT AND LENW. |
|
|
196 C |
|
|
197 C***FIRST EXECUTABLE STATEMENT DQAGP |
|
|
198 IER = 6 |
|
|
199 NEVAL = 0 |
|
|
200 LAST = 0 |
|
|
201 RESULT = 0.0D+00 |
|
|
202 ABSERR = 0.0D+00 |
|
|
203 IF(LENIW.LT.(3*NPTS2-2).OR.LENW.LT.(LENIW*2-NPTS2).OR.NPTS2.LT.2) |
|
|
204 * GO TO 10 |
|
|
205 C |
|
|
206 C PREPARE CALL FOR DQAGPE. |
|
|
207 C |
|
|
208 LIMIT = (LENIW-NPTS2)/2 |
|
|
209 L1 = LIMIT+1 |
|
|
210 L2 = LIMIT+L1 |
|
|
211 L3 = LIMIT+L2 |
|
|
212 L4 = LIMIT+L3 |
|
|
213 C |
|
|
214 CALL DQAGPE(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,LIMIT,RESULT,ABSERR, |
|
|
215 * NEVAL,IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),WORK(L4), |
|
|
216 * IWORK(1),IWORK(L1),IWORK(L2),LAST) |
|
|
217 C |
|
|
218 C CALL ERROR HANDLER IF NECESSARY. |
|
|
219 C |
|
|
220 LVL = 0 |
|
|
221 10 IF(IER.EQ.6) LVL = 1 |
|
4040
|
222 IF(IER.GT.0) CALL XERROR('ABNORMAL RETURN FROM DQAGP',26,IER,LVL) |
|
2329
|
223 RETURN |
|
|
224 END |
