GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
dqagie.f
Go to the documentation of this file.
1  SUBROUTINE dqagie(F,BOUND,INF,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
2  * neval,ier,alist,blist,rlist,elist,iord,last)
3 C***BEGIN PROLOGUE DQAGIE
4 C***DATE WRITTEN 800101 (YYMMDD)
5 C***REVISION DATE 830518 (YYMMDD)
6 C***CATEGORY NO. H2A3A1,H2A4A1
7 C***KEYWORDS AUTOMATIC INTEGRATOR, INFINITE INTERVALS,
8 C GENERAL-PURPOSE, TRANSFORMATION, EXTRAPOLATION,
9 C 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 INTEGRAL I = INTEGRAL OF F OVER (BOUND,+INFINITY)
14 C OR I = INTEGRAL OF F OVER (-INFINITY,BOUND)
15 C OR I = INTEGRAL OF F OVER (-INFINITY,+INFINITY),
16 C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
17 C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I))
18 C***DESCRIPTION
19 C
20 C INTEGRATION OVER INFINITE INTERVALS
21 C STANDARD FORTRAN SUBROUTINE
22 C
23 C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND
24 C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
25 C DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
26 C
27 C BOUND - DOUBLE PRECISION
28 C FINITE BOUND OF INTEGRATION RANGE
29 C (HAS NO MEANING IF INTERVAL IS DOUBLY-INFINITE)
30 C
31 C INF - DOUBLE PRECISION
32 C INDICATING THE KIND OF INTEGRATION RANGE INVOLVED
33 C INF = 1 CORRESPONDS TO (BOUND,+INFINITY),
34 C INF = -1 TO (-INFINITY,BOUND),
35 C INF = 2 TO (-INFINITY,+INFINITY).
36 C
37 C EPSABS - DOUBLE PRECISION
38 C ABSOLUTE ACCURACY REQUESTED
39 C EPSREL - DOUBLE PRECISION
40 C RELATIVE ACCURACY REQUESTED
41 C IF EPSABS.LE.0
42 C AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
43 C THE ROUTINE WILL END WITH IER = 6.
44 C
45 C LIMIT - INTEGER
46 C GIVES AN UPPER BOUND ON THE NUMBER OF SUBINTERVALS
47 C IN THE PARTITION OF (A,B), LIMIT.GE.1
48 C
49 C ON RETURN
50 C RESULT - DOUBLE PRECISION
51 C APPROXIMATION TO THE INTEGRAL
52 C
53 C ABSERR - DOUBLE PRECISION
54 C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
55 C WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
56 C
57 C NEVAL - INTEGER
58 C NUMBER OF INTEGRAND EVALUATIONS
59 C
60 C IER - INTEGER
61 C IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
62 C ROUTINE. IT IS ASSUMED THAT THE REQUESTED
63 C ACCURACY HAS BEEN ACHIEVED.
64 C IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE. THE
65 C ESTIMATES FOR RESULT AND ERROR ARE LESS
66 C RELIABLE. IT IS ASSUMED THAT THE REQUESTED
67 C ACCURACY HAS NOT BEEN ACHIEVED.
68 C IER.LT.0 EXIT REQUESTED FROM USER-SUPPLIED
69 C FUNCTION.
70 C
71 C ERROR MESSAGES
72 C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
73 C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
74 C SUBDIVISIONS BY INCREASING THE VALUE OF
75 C LIMIT (AND TAKING THE ACCORDING DIMENSION
76 C ADJUSTMENTS INTO ACCOUNT). HOWEVER,IF
77 C THIS YIELDS NO IMPROVEMENT IT IS ADVISED
78 C TO ANALYZE THE INTEGRAND IN ORDER TO
79 C DETERMINE THE INTEGRATION DIFFICULTIES.
80 C IF THE POSITION OF A LOCAL DIFFICULTY CAN
81 C BE DETERMINED (E.G. SINGULARITY,
82 C DISCONTINUITY WITHIN THE INTERVAL) ONE
83 C WILL PROBABLY GAIN FROM SPLITTING UP THE
84 C INTERVAL AT THIS POINT AND CALLING THE
85 C INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
86 C AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
87 C SHOULD BE USED, WHICH IS DESIGNED FOR
88 C HANDLING THE TYPE OF DIFFICULTY INVOLVED.
89 C = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
90 C DETECTED, WHICH PREVENTS THE REQUESTED
91 C TOLERANCE FROM BEING ACHIEVED.
92 C THE ERROR MAY BE UNDER-ESTIMATED.
93 C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
94 C AT SOME POINTS OF THE INTEGRATION
95 C INTERVAL.
96 C = 4 THE ALGORITHM DOES NOT CONVERGE.
97 C ROUNDOFF ERROR IS DETECTED IN THE
98 C EXTRAPOLATION TABLE.
99 C IT IS ASSUMED THAT THE REQUESTED TOLERANCE
100 C CANNOT BE ACHIEVED, AND THAT THE RETURNED
101 C RESULT IS THE BEST WHICH CAN BE OBTAINED.
102 C = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
103 C SLOWLY CONVERGENT. IT MUST BE NOTED THAT
104 C DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
105 C OF IER.
106 C = 6 THE INPUT IS INVALID, BECAUSE
107 C (EPSABS.LE.0 AND
108 C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
109 C RESULT, ABSERR, NEVAL, LAST, RLIST(1),
110 C ELIST(1) AND IORD(1) ARE SET TO ZERO.
111 C ALIST(1) AND BLIST(1) ARE SET TO 0
112 C AND 1 RESPECTIVELY.
113 C
114 C ALIST - DOUBLE PRECISION
115 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
116 C LAST ELEMENTS OF WHICH ARE THE LEFT
117 C END POINTS OF THE SUBINTERVALS IN THE PARTITION
118 C OF THE TRANSFORMED INTEGRATION RANGE (0,1).
119 C
120 C BLIST - DOUBLE PRECISION
121 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
122 C LAST ELEMENTS OF WHICH ARE THE RIGHT
123 C END POINTS OF THE SUBINTERVALS IN THE PARTITION
124 C OF THE TRANSFORMED INTEGRATION RANGE (0,1).
125 C
126 C RLIST - DOUBLE PRECISION
127 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
128 C LAST ELEMENTS OF WHICH ARE THE INTEGRAL
129 C APPROXIMATIONS ON THE SUBINTERVALS
130 C
131 C ELIST - DOUBLE PRECISION
132 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
133 C LAST ELEMENTS OF WHICH ARE THE MODULI OF THE
134 C ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
135 C
136 C IORD - INTEGER
137 C VECTOR OF DIMENSION LIMIT, THE FIRST K
138 C ELEMENTS OF WHICH ARE POINTERS TO THE
139 C ERROR ESTIMATES OVER THE SUBINTERVALS,
140 C SUCH THAT ELIST(IORD(1)), ..., ELIST(IORD(K))
141 C FORM A DECREASING SEQUENCE, WITH K = LAST
142 C IF LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST
143 C OTHERWISE
144 C
145 C LAST - INTEGER
146 C NUMBER OF SUBINTERVALS ACTUALLY PRODUCED
147 C IN THE SUBDIVISION PROCESS
148 C
149 C***REFERENCES (NONE)
150 C***ROUTINES CALLED D1MACH,DQELG,DQK15I,DQPSRT
151 C***END PROLOGUE DQAGIE
152  DOUBLE PRECISION abseps,abserr,alist,area,area1,area12,area2,a1,
153  * a2,blist,boun,bound,b1,b2,correc,dabs,defabs,defab1,defab2,
154  * dmax1,dres,d1mach,elist,epmach,epsabs,epsrel,erlarg,erlast,
155  * errbnd,errmax,error1,error2,erro12,errsum,ertest,oflow,resabs,
156  * reseps,result,res3la,rlist,rlist2,small,uflow
157  INTEGER id,ier,ierro,inf,iord,iroff1,iroff2,iroff3,jupbnd,k,ksgn,
158  * ktmin,last,limit,maxerr,neval,nres,nrmax,numrl2
159  LOGICAL extrap,noext
160 C
161  dimension alist(limit),blist(limit),elist(limit),iord(limit),
162  * res3la(3),rlist(limit),rlist2(52)
163 C
164  EXTERNAL f
165 C
166 C THE DIMENSION OF RLIST2 IS DETERMINED BY THE VALUE OF
167 C LIMEXP IN SUBROUTINE DQELG.
168 C
169 C
170 C LIST OF MAJOR VARIABLES
171 C -----------------------
172 C
173 C ALIST - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
174 C CONSIDERED UP TO NOW
175 C BLIST - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
176 C CONSIDERED UP TO NOW
177 C RLIST(I) - APPROXIMATION TO THE INTEGRAL OVER
178 C (ALIST(I),BLIST(I))
179 C RLIST2 - ARRAY OF DIMENSION AT LEAST (LIMEXP+2),
180 C CONTAINING THE PART OF THE EPSILON TABLE
181 C WICH IS STILL NEEDED FOR FURTHER COMPUTATIONS
182 C ELIST(I) - ERROR ESTIMATE APPLYING TO RLIST(I)
183 C MAXERR - POINTER TO THE INTERVAL WITH LARGEST ERROR
184 C ESTIMATE
185 C ERRMAX - ELIST(MAXERR)
186 C ERLAST - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED
187 C (BEFORE THAT SUBDIVISION HAS TAKEN PLACE)
188 C AREA - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
189 C ERRSUM - SUM OF THE ERRORS OVER THE SUBINTERVALS
190 C ERRBND - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
191 C ABS(RESULT))
192 C *****1 - VARIABLE FOR THE LEFT SUBINTERVAL
193 C *****2 - VARIABLE FOR THE RIGHT SUBINTERVAL
194 C LAST - INDEX FOR SUBDIVISION
195 C NRES - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE
196 C NUMRL2 - NUMBER OF ELEMENTS CURRENTLY IN RLIST2. IF AN
197 C APPROPRIATE APPROXIMATION TO THE COMPOUNDED
198 C INTEGRAL HAS BEEN OBTAINED, IT IS PUT IN
199 C RLIST2(NUMRL2) AFTER NUMRL2 HAS BEEN INCREASED
200 C BY ONE.
201 C SMALL - LENGTH OF THE SMALLEST INTERVAL CONSIDERED UP
202 C TO NOW, MULTIPLIED BY 1.5
203 C ERLARG - SUM OF THE ERRORS OVER THE INTERVALS LARGER
204 C THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW
205 C EXTRAP - LOGICAL VARIABLE DENOTING THAT THE ROUTINE
206 C IS ATTEMPTING TO PERFORM EXTRAPOLATION. I.E.
207 C BEFORE SUBDIVIDING THE SMALLEST INTERVAL WE
208 C TRY TO DECREASE THE VALUE OF ERLARG.
209 C NOEXT - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION
210 C IS NO LONGER ALLOWED (TRUE-VALUE)
211 C
212 C MACHINE DEPENDENT CONSTANTS
213 C ---------------------------
214 C
215 C EPMACH IS THE LARGEST RELATIVE SPACING.
216 C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
217 C OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
218 C
219 C***FIRST EXECUTABLE STATEMENT DQAGIE
220  epmach = d1mach(4)
221 C
222 C TEST ON VALIDITY OF PARAMETERS
223 C -----------------------------
224 C
225  ier = 0
226  neval = 0
227  last = 0
228  result = 0.0d+00
229  abserr = 0.0d+00
230  alist(1) = 0.0d+00
231  blist(1) = 0.1d+01
232  rlist(1) = 0.0d+00
233  elist(1) = 0.0d+00
234  iord(1) = 0
235  IF(epsabs.LE.0.0d+00.AND.epsrel.LT.dmax1(0.5d+02*epmach,0.5d-28))
236  * ier = 6
237  IF(ier.EQ.6) go to 999
238 C
239 C
240 C FIRST APPROXIMATION TO THE INTEGRAL
241 C -----------------------------------
242 C
243 C DETERMINE THE INTERVAL TO BE MAPPED ONTO (0,1).
244 C IF INF = 2 THE INTEGRAL IS COMPUTED AS I = I1+I2, WHERE
245 C I1 = INTEGRAL OF F OVER (-INFINITY,0),
246 C I2 = INTEGRAL OF F OVER (0,+INFINITY).
247 C
248  boun = bound
249  IF(inf.EQ.2) boun = 0.0d+00
250  CALL dqk15i(f,boun,inf,0.0d+00,0.1d+01,result,abserr,
251  * defabs,resabs,ier)
252  IF (ier .LT. 0) RETURN
253 C
254 C TEST ON ACCURACY
255 C
256  last = 1
257  rlist(1) = result
258  elist(1) = abserr
259  iord(1) = 1
260  dres = dabs(result)
261  errbnd = dmax1(epsabs,epsrel*dres)
262  IF(abserr.LE.1.0d+02*epmach*defabs.AND.abserr.GT.errbnd) ier = 2
263  IF(limit.EQ.1) ier = 1
264  IF(ier.NE.0.OR.(abserr.LE.errbnd.AND.abserr.NE.resabs).OR.
265  * abserr.EQ.0.0d+00) go to 130
266 C
267 C INITIALIZATION
268 C --------------
269 C
270  uflow = d1mach(1)
271  oflow = d1mach(2)
272  rlist2(1) = result
273  errmax = abserr
274  maxerr = 1
275  area = result
276  errsum = abserr
277  abserr = oflow
278  nrmax = 1
279  nres = 0
280  ktmin = 0
281  numrl2 = 2
282  extrap = .false.
283  noext = .false.
284  ierro = 0
285  iroff1 = 0
286  iroff2 = 0
287  iroff3 = 0
288  ksgn = -1
289  IF(dres.GE.(0.1d+01-0.5d+02*epmach)*defabs) ksgn = 1
290 C
291 C MAIN DO-LOOP
292 C ------------
293 C
294  DO 90 last = 2,limit
295 C
296 C BISECT THE SUBINTERVAL WITH NRMAX-TH LARGEST ERROR ESTIMATE.
297 C
298  a1 = alist(maxerr)
299  b1 = 0.5d+00*(alist(maxerr)+blist(maxerr))
300  a2 = b1
301  b2 = blist(maxerr)
302  erlast = errmax
303  CALL dqk15i(f,boun,inf,a1,b1,area1,error1,resabs,defab1,ier)
304  IF (ier .LT. 0) RETURN
305  CALL dqk15i(f,boun,inf,a2,b2,area2,error2,resabs,defab2,ier)
306  IF (ier .LT. 0) RETURN
307 C
308 C IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
309 C AND ERROR AND TEST FOR ACCURACY.
310 C
311  area12 = area1+area2
312  erro12 = error1+error2
313  errsum = errsum+erro12-errmax
314  area = area+area12-rlist(maxerr)
315  IF(defab1.EQ.error1.OR.defab2.EQ.error2)go to 15
316  IF(dabs(rlist(maxerr)-area12).GT.0.1d-04*dabs(area12)
317  * .OR.erro12.LT.0.99d+00*errmax) go to 10
318  IF(extrap) iroff2 = iroff2+1
319  IF(.NOT.extrap) iroff1 = iroff1+1
320  10 IF(last.GT.10.AND.erro12.GT.errmax) iroff3 = iroff3+1
321  15 rlist(maxerr) = area1
322  rlist(last) = area2
323  errbnd = dmax1(epsabs,epsrel*dabs(area))
324 C
325 C TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
326 C
327  IF(iroff1+iroff2.GE.10.OR.iroff3.GE.20) ier = 2
328  IF(iroff2.GE.5) ierro = 3
329 C
330 C SET ERROR FLAG IN THE CASE THAT THE NUMBER OF
331 C SUBINTERVALS EQUALS LIMIT.
332 C
333  IF(last.EQ.limit) ier = 1
334 C
335 C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
336 C AT SOME POINTS OF THE INTEGRATION RANGE.
337 C
338  IF(dmax1(dabs(a1),dabs(b2)).LE.(0.1d+01+0.1d+03*epmach)*
339  * (dabs(a2)+0.1d+04*uflow)) ier = 4
340 C
341 C APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
342 C
343  IF(error2.GT.error1) go to 20
344  alist(last) = a2
345  blist(maxerr) = b1
346  blist(last) = b2
347  elist(maxerr) = error1
348  elist(last) = error2
349  go to 30
350  20 alist(maxerr) = a2
351  alist(last) = a1
352  blist(last) = b1
353  rlist(maxerr) = area2
354  rlist(last) = area1
355  elist(maxerr) = error2
356  elist(last) = error1
357 C
358 C CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
359 C IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
360 C WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
361 C
362  30 CALL dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
363  IF(errsum.LE.errbnd) go to 115
364  IF(ier.NE.0) go to 100
365  IF(last.EQ.2) go to 80
366  IF(noext) go to 90
367  erlarg = erlarg-erlast
368  IF(dabs(b1-a1).GT.small) erlarg = erlarg+erro12
369  IF(extrap) go to 40
370 C
371 C TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
372 C SMALLEST INTERVAL.
373 C
374  IF(dabs(blist(maxerr)-alist(maxerr)).GT.small) go to 90
375  extrap = .true.
376  nrmax = 2
377  40 IF(ierro.EQ.3.OR.erlarg.LE.ertest) go to 60
378 C
379 C THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
380 C BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER THE
381 C LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
382 C
383  id = nrmax
384  jupbnd = last
385  IF(last.GT.(2+limit/2)) jupbnd = limit+3-last
386  DO 50 k = id,jupbnd
387  maxerr = iord(nrmax)
388  errmax = elist(maxerr)
389  IF(dabs(blist(maxerr)-alist(maxerr)).GT.small) go to 90
390  nrmax = nrmax+1
391  50 CONTINUE
392 C
393 C PERFORM EXTRAPOLATION.
394 C
395  60 numrl2 = numrl2+1
396  rlist2(numrl2) = area
397  CALL dqelg(numrl2,rlist2,reseps,abseps,res3la,nres)
398  ktmin = ktmin+1
399  IF(ktmin.GT.5.AND.abserr.LT.0.1d-02*errsum) ier = 5
400  IF(abseps.GE.abserr) go to 70
401  ktmin = 0
402  abserr = abseps
403  result = reseps
404  correc = erlarg
405  ertest = dmax1(epsabs,epsrel*dabs(reseps))
406  IF(abserr.LE.ertest) go to 100
407 C
408 C PREPARE BISECTION OF THE SMALLEST INTERVAL.
409 C
410  70 IF(numrl2.EQ.1) noext = .true.
411  IF(ier.EQ.5) go to 100
412  maxerr = iord(1)
413  errmax = elist(maxerr)
414  nrmax = 1
415  extrap = .false.
416  small = small*0.5d+00
417  erlarg = errsum
418  go to 90
419  80 small = 0.375d+00
420  erlarg = errsum
421  ertest = errbnd
422  rlist2(2) = area
423  90 CONTINUE
424 C
425 C SET FINAL RESULT AND ERROR ESTIMATE.
426 C ------------------------------------
427 C
428  100 IF(abserr.EQ.oflow) go to 115
429  IF((ier+ierro).EQ.0) go to 110
430  IF(ierro.EQ.3) abserr = abserr+correc
431  IF(ier.EQ.0) ier = 3
432  IF(result.NE.0.0d+00.AND.area.NE.0.0d+00)go to 105
433  IF(abserr.GT.errsum)go to 115
434  IF(area.EQ.0.0d+00) go to 130
435  go to 110
436  105 IF(abserr/dabs(result).GT.errsum/dabs(area))go to 115
437 C
438 C TEST ON DIVERGENCE
439 C
440  110 IF(ksgn.EQ.(-1).AND.dmax1(dabs(result),dabs(area)).LE.
441  * defabs*0.1d-01) go to 130
442  IF(0.1d-01.GT.(result/area).OR.(result/area).GT.0.1d+03.
443  *or.errsum.GT.dabs(area)) ier = 6
444  go to 130
445 C
446 C COMPUTE GLOBAL INTEGRAL SUM.
447 C
448  115 result = 0.0d+00
449  DO 120 k = 1,last
450  result = result+rlist(k)
451  120 CONTINUE
452  abserr = errsum
453  130 neval = 30*last-15
454  IF(inf.EQ.2) neval = 2*neval
455  IF(ier.GT.2) ier=ier-1
456  999 RETURN
457  END