GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
double precision function d1mach(i)
Definition: d1mach.f:23
subroutine dqagie(F, BOUND, INF, EPSABS, EPSREL, LIMIT, RESULT, ABSERR, NEVAL, IER, ALIST, BLIST, RLIST, ELIST, IORD, LAST)
Definition: dqagie.f:3
subroutine dqelg(N, EPSTAB, RESULT, ABSERR, RES3LA, NRES)
Definition: dqelg.f:2
subroutine dqk15i(F, BOUN, INF, A, B, RESULT, ABSERR, RESABS, RESASC, IERR)
Definition: dqk15i.f:3
subroutine dqpsrt(LIMIT, LAST, MAXERR, ERMAX, ELIST, IORD, NRMAX)
Definition: dqpsrt.f:2