V3FIT
dqag.f
1  subroutine dqag(f,a,b,epsabs,epsrel,key,result,abserr,neval,ier,
2  * limit,lenw,last,iwork,work)
3 c***begin prologue dqag
4 c***date written 800101 (yymmdd)
5 c***revision date 830518 (yymmdd)
6 c***category no. h2a1a1
7 c***keywords automatic integrator, general-purpose,
8 c integrand examinator, globally adaptive,
9 c gauss-kronrod
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 abs(i-result)le.max(epsabs,epsrel*abs(i)).
16 c***description
17 c
18 c computation of a definite integral
19 c standard fortran subroutine
20 c double precision version
21 c
22 c f - double precision
23 c function subprogam 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 a - double precision
28 c lower limit of integration
29 c
30 c b - double precision
31 c upper limit of integration
32 c
33 c epsabs - double precision
34 c absolute accoracy requested
35 c epsrel - double precision
36 c relative accuracy requested
37 c if epsabs.le.0
38 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
39 c the routine will end with ier = 6.
40 c
41 c key - integer
42 c key for choice of local integration rule
43 c a gauss-kronrod pair is used with
44 c 7 - 15 points if key.lt.2,
45 c 10 - 21 points if key = 2,
46 c 15 - 31 points if key = 3,
47 c 20 - 41 points if key = 4,
48 c 25 - 51 points if key = 5,
49 c 30 - 61 points if key.gt.5.
50 c
51 c on return
52 c result - double precision
53 c approximation to the integral
54 c
55 c abserr - double precision
56 c estimate of the modulus of the absolute error,
57 c which should equal or exceed abs(i-result)
58 c
59 c neval - integer
60 c number of integrand evaluations
61 c
62 c ier - integer
63 c ier = 0 normal and reliable termination of the
64 c routine. it is assumed that the requested
65 c accuracy has been achieved.
66 c ier.gt.0 abnormal termination of the routine
67 c the estimates for result and error are
68 c less reliable. it is assumed that the
69 c requested accuracy has not been achieved.
70 c error messages
71 c ier = 1 maximum number of subdivisions allowed
72 c has been achieved. one can allow more
73 c subdivisions by increasing the value of
74 c limit (and taking the according dimension
75 c adjustments into account). however, if
76 c this yield no improvement it is advised
77 c to analyze the integrand in order to
78 c determine the integration difficulaties.
79 c if the position of a local difficulty can
80 c be determined (i.e.singularity,
81 c discontinuity within the interval) one
82 c will probably gain from splitting up the
83 c interval at this point and calling the
84 c integrator on the subranges. if possible,
85 c an appropriate special-purpose integrator
86 c should be used which is designed for
87 c handling the type of difficulty involved.
88 c = 2 the occurrence of roundoff error is
89 c detected, which prevents the requested
90 c tolerance from being achieved.
91 c = 3 extremely bad integrand behaviour occurs
92 c at some points of the integration
93 c interval.
94 c = 6 the input is invalid, because
95 c (epsabs.le.0 and
96 c epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
97 c or limit.lt.1 or lenw.lt.limit*4.
98 c result, abserr, neval, last are set
99 c to zero.
100 c except when lenw is invalid, iwork(1),
101 c work(limit*2+1) and work(limit*3+1) are
102 c set to zero, work(1) is set to a and
103 c work(limit+1) to b.
104 c
105 c dimensioning parameters
106 c limit - integer
107 c dimensioning parameter for iwork
108 c limit determines the maximum number of subintervals
109 c in the partition of the given integration interval
110 c (a,b), limit.ge.1.
111 c if limit.lt.1, the routine will end with ier = 6.
112 c
113 c lenw - integer
114 c dimensioning parameter for work
115 c lenw must be at least limit*4.
116 c if lenw.lt.limit*4, the routine will end with
117 c ier = 6.
118 c
119 c last - integer
120 c on return, last equals the number of subintervals
121 c produced in the subdiviosion process, which
122 c determines the number of significant elements
123 c actually in the work arrays.
124 c
125 c work arrays
126 c iwork - integer
127 c vector of dimension at least limit, the first k
128 c elements of which contain pointers to the error
129 c estimates over the subintervals, such that
130 c work(limit*3+iwork(1)),... , work(limit*3+iwork(k))
131 c form a decreasing sequence with k = last if
132 c last.le.(limit/2+2), and k = limit+1-last otherwise
133 c
134 c work - double precision
135 c vector of dimension at least lenw
136 c on return
137 c work(1), ..., work(last) contain the left end
138 c points of the subintervals in the partition of
139 c (a,b),
140 c work(limit+1), ..., work(limit+last) contain the
141 c right end points,
142 c work(limit*2+1), ..., work(limit*2+last) contain
143 c the integral approximations over the subintervals,
144 c work(limit*3+1), ..., work(limit*3+last) contain
145 c the error estimates.
146 c
147 c***references (none)
148 c***routines called dqage,xerror
149 c***end prologue dqag
150  double precision a,abserr,b,epsabs,epsrel,f,result,work
151  integer ier,iwork,key,last,lenw,limit,lvl,l1,l2,l3,neval
152 c
153  dimension iwork(limit),work(lenw)
154 c
155  external f
156 c
157 c check validity of lenw.
158 c
159 c***first executable statement dqag
160  ier = 6
161  neval = 0
162  last = 0
163  result = 0.0d+00
164  abserr = 0.0d+00
165  if(limit.lt.1.or.lenw.lt.limit*4) go to 10
166 c
167 c prepare call for dqage.
168 c
169  l1 = limit+1
170  l2 = limit+l1
171  l3 = limit+l2
172 c
173  call dqage(f,a,b,epsabs,epsrel,key,limit,result,abserr,neval,
174  * ier,work(1),work(l1),work(l2),work(l3),iwork,last)
175 c
176 c call error handler if necessary.
177 c
178  lvl = 0
179 10 if(ier.eq.6) lvl = 1
180 c SAL 04/30/14 xerror not inlucdes/couldn't be found at linking
181 c if(ier.ne.0) call xerror(26habnormal return from dqag ,26,ier,lvl)
182  return
183  end
184  subroutine dqage(f,a,b,epsabs,epsrel,key,limit,result,abserr,
185  * neval,ier,alist,blist,rlist,elist,iord,last)
186 c***begin prologue dqage
187 c***date written 800101 (yymmdd)
188 c***revision date 830518 (yymmdd)
189 c***category no. h2a1a1
190 c***keywords automatic integrator, general-purpose,
191 c integrand examinator, globally adaptive,
192 c gauss-kronrod
193 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
194 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
195 c***purpose the routine calculates an approximation result to a given
196 c definite integral i = integral of f over (a,b),
197 c hopefully satisfying following claim for accuracy
198 c abs(i-reslt).le.max(epsabs,epsrel*abs(i)).
199 c***description
200 c
201 c computation of a definite integral
202 c standard fortran subroutine
203 c double precision version
204 c
205 c parameters
206 c on entry
207 c f - double precision
208 c function subprogram defining the integrand
209 c function f(x). the actual name for f needs to be
210 c declared e x t e r n a l in the driver program.
211 c
212 c a - double precision
213 c lower limit of integration
214 c
215 c b - double precision
216 c upper limit of integration
217 c
218 c epsabs - double precision
219 c absolute accuracy requested
220 c epsrel - double precision
221 c relative accuracy requested
222 c if epsabs.le.0
223 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
224 c the routine will end with ier = 6.
225 c
226 c key - integer
227 c key for choice of local integration rule
228 c a gauss-kronrod pair is used with
229 c 7 - 15 points if key.lt.2,
230 c 10 - 21 points if key = 2,
231 c 15 - 31 points if key = 3,
232 c 20 - 41 points if key = 4,
233 c 25 - 51 points if key = 5,
234 c 30 - 61 points if key.gt.5.
235 c
236 c limit - integer
237 c gives an upperbound on the number of subintervals
238 c in the partition of (a,b), limit.ge.1.
239 c
240 c on return
241 c result - double precision
242 c approximation to the integral
243 c
244 c abserr - double precision
245 c estimate of the modulus of the absolute error,
246 c which should equal or exceed abs(i-result)
247 c
248 c neval - integer
249 c number of integrand evaluations
250 c
251 c ier - integer
252 c ier = 0 normal and reliable termination of the
253 c routine. it is assumed that the requested
254 c accuracy has been achieved.
255 c ier.gt.0 abnormal termination of the routine
256 c the estimates for result and error are
257 c less reliable. it is assumed that the
258 c requested accuracy has not been achieved.
259 c error messages
260 c ier = 1 maximum number of subdivisions allowed
261 c has been achieved. one can allow more
262 c subdivisions by increasing the value
263 c of limit.
264 c however, if this yields no improvement it
265 c is rather advised to analyze the integrand
266 c in order to determine the integration
267 c difficulties. if the position of a local
268 c difficulty can be determined(e.g.
269 c singularity, discontinuity within the
270 c interval) one will probably gain from
271 c splitting up the interval at this point
272 c and calling the integrator on the
273 c subranges. if possible, an appropriate
274 c special-purpose integrator should be used
275 c which is designed for handling the type of
276 c difficulty involved.
277 c = 2 the occurrence of roundoff error is
278 c detected, which prevents the requested
279 c tolerance from being achieved.
280 c = 3 extremely bad integrand behaviour occurs
281 c at some points of the integration
282 c interval.
283 c = 6 the input is invalid, because
284 c (epsabs.le.0 and
285 c epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
286 c result, abserr, neval, last, rlist(1) ,
287 c elist(1) and iord(1) are set to zero.
288 c alist(1) and blist(1) are set to a and b
289 c respectively.
290 c
291 c alist - double precision
292 c vector of dimension at least limit, the first
293 c last elements of which are the left
294 c end points of the subintervals in the partition
295 c of the given integration range (a,b)
296 c
297 c blist - double precision
298 c vector of dimension at least limit, the first
299 c last elements of which are the right
300 c end points of the subintervals in the partition
301 c of the given integration range (a,b)
302 c
303 c rlist - double precision
304 c vector of dimension at least limit, the first
305 c last elements of which are the
306 c integral approximations on the subintervals
307 c
308 c elist - double precision
309 c vector of dimension at least limit, the first
310 c last elements of which are the moduli of the
311 c absolute error estimates on the subintervals
312 c
313 c iord - integer
314 c vector of dimension at least limit, the first k
315 c elements of which are pointers to the
316 c error estimates over the subintervals,
317 c such that elist(iord(1)), ...,
318 c elist(iord(k)) form a decreasing sequence,
319 c with k = last if last.le.(limit/2+2), and
320 c k = limit+1-last otherwise
321 c
322 c last - integer
323 c number of subintervals actually produced in the
324 c subdivision process
325 c
326 c***references (none)
327 c***routines called d1mach,dqk15,dqk21,dqk31,
328 c dqk41,dqk51,dqk61,dqpsrt
329 c***end prologue dqage
330 c
331  double precision a,abserr,alist,area,area1,area12,area2,a1,a2,b,
332  * blist,b1,b2,dabs,defabs,defab1,defab2,dmax1,d1mach,elist,epmach,
333  * epsabs,epsrel,errbnd,errmax,error1,error2,erro12,errsum,f,
334  * resabs,result,rlist,uflow
335  integer ier,iord,iroff1,iroff2,k,key,keyf,last,limit,maxerr,neval,
336  * nrmax
337 c
338  dimension alist(limit),blist(limit),elist(limit),iord(limit),
339  * rlist(limit)
340 c
341  external f
342 c
343 c list of major variables
344 c -----------------------
345 c
346 c alist - list of left end points of all subintervals
347 c considered up to now
348 c blist - list of right end points of all subintervals
349 c considered up to now
350 c rlist(i) - approximation to the integral over
351 c (alist(i),blist(i))
352 c elist(i) - error estimate applying to rlist(i)
353 c maxerr - pointer to the interval with largest
354 c error estimate
355 c errmax - elist(maxerr)
356 c area - sum of the integrals over the subintervals
357 c errsum - sum of the errors over the subintervals
358 c errbnd - requested accuracy max(epsabs,epsrel*
359 c abs(result))
360 c *****1 - variable for the left subinterval
361 c *****2 - variable for the right subinterval
362 c last - index for subdivision
363 c
364 c
365 c machine dependent constants
366 c ---------------------------
367 c
368 c epmach is the largest relative spacing.
369 c uflow is the smallest positive magnitude.
370 c
371 c***first executable statement dqage
372  epmach = d1mach(4)
373  uflow = d1mach(1)
374 c
375 c test on validity of parameters
376 c ------------------------------
377 c
378  ier = 0
379  neval = 0
380  last = 0
381  result = 0.0d+00
382  abserr = 0.0d+00
383  alist(1) = a
384  blist(1) = b
385  rlist(1) = 0.0d+00
386  elist(1) = 0.0d+00
387  iord(1) = 0
388  if(epsabs.le.0.0d+00.and.
389  * epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28)) ier = 6
390  if(ier.eq.6) go to 999
391 c
392 c first approximation to the integral
393 c -----------------------------------
394 c
395  keyf = key
396  if(key.le.0) keyf = 1
397  if(key.ge.7) keyf = 6
398  neval = 0
399  if(keyf.eq.1) call dqk15(f,a,b,result,abserr,defabs,resabs)
400  if(keyf.eq.2) call dqk21(f,a,b,result,abserr,defabs,resabs)
401  if(keyf.eq.3) call dqk31(f,a,b,result,abserr,defabs,resabs)
402  if(keyf.eq.4) call dqk41(f,a,b,result,abserr,defabs,resabs)
403  if(keyf.eq.5) call dqk51(f,a,b,result,abserr,defabs,resabs)
404  if(keyf.eq.6) call dqk61(f,a,b,result,abserr,defabs,resabs)
405  last = 1
406  rlist(1) = result
407  elist(1) = abserr
408  iord(1) = 1
409 c
410 c test on accuracy.
411 c
412  errbnd = dmax1(epsabs,epsrel*dabs(result))
413  if(abserr.le.0.5d+02*epmach*defabs.and.abserr.gt.errbnd) ier = 2
414  if(limit.eq.1) ier = 1
415  if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs)
416  * .or.abserr.eq.0.0d+00) go to 60
417 c
418 c initialization
419 c --------------
420 c
421 c
422  errmax = abserr
423  maxerr = 1
424  area = result
425  errsum = abserr
426  nrmax = 1
427  iroff1 = 0
428  iroff2 = 0
429 c
430 c main do-loop
431 c ------------
432 c
433  do 30 last = 2,limit
434 c
435 c bisect the subinterval with the largest error estimate.
436 c
437  a1 = alist(maxerr)
438  b1 = 0.5d+00*(alist(maxerr)+blist(maxerr))
439  a2 = b1
440  b2 = blist(maxerr)
441  if(keyf.eq.1) call dqk15(f,a1,b1,area1,error1,resabs,defab1)
442  if(keyf.eq.2) call dqk21(f,a1,b1,area1,error1,resabs,defab1)
443  if(keyf.eq.3) call dqk31(f,a1,b1,area1,error1,resabs,defab1)
444  if(keyf.eq.4) call dqk41(f,a1,b1,area1,error1,resabs,defab1)
445  if(keyf.eq.5) call dqk51(f,a1,b1,area1,error1,resabs,defab1)
446  if(keyf.eq.6) call dqk61(f,a1,b1,area1,error1,resabs,defab1)
447  if(keyf.eq.1) call dqk15(f,a2,b2,area2,error2,resabs,defab2)
448  if(keyf.eq.2) call dqk21(f,a2,b2,area2,error2,resabs,defab2)
449  if(keyf.eq.3) call dqk31(f,a2,b2,area2,error2,resabs,defab2)
450  if(keyf.eq.4) call dqk41(f,a2,b2,area2,error2,resabs,defab2)
451  if(keyf.eq.5) call dqk51(f,a2,b2,area2,error2,resabs,defab2)
452  if(keyf.eq.6) call dqk61(f,a2,b2,area2,error2,resabs,defab2)
453 c
454 c improve previous approximations to integral
455 c and error and test for accuracy.
456 c
457  neval = neval+1
458  area12 = area1+area2
459  erro12 = error1+error2
460  errsum = errsum+erro12-errmax
461  area = area+area12-rlist(maxerr)
462  if(defab1.eq.error1.or.defab2.eq.error2) go to 5
463  if(dabs(rlist(maxerr)-area12).le.0.1d-04*dabs(area12)
464  * .and.erro12.ge.0.99d+00*errmax) iroff1 = iroff1+1
465  if(last.gt.10.and.erro12.gt.errmax) iroff2 = iroff2+1
466  5 rlist(maxerr) = area1
467  rlist(last) = area2
468  errbnd = dmax1(epsabs,epsrel*dabs(area))
469  if(errsum.le.errbnd) go to 8
470 c
471 c test for roundoff error and eventually set error flag.
472 c
473  if(iroff1.ge.6.or.iroff2.ge.20) ier = 2
474 c
475 c set error flag in the case that the number of subintervals
476 c equals limit.
477 c
478  if(last.eq.limit) ier = 1
479 c
480 c set error flag in the case of bad integrand behaviour
481 c at a point of the integration range.
482 c
483  if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*
484  * epmach)*(dabs(a2)+0.1d+04*uflow)) ier = 3
485 c
486 c append the newly-created intervals to the list.
487 c
488  8 if(error2.gt.error1) go to 10
489  alist(last) = a2
490  blist(maxerr) = b1
491  blist(last) = b2
492  elist(maxerr) = error1
493  elist(last) = error2
494  go to 20
495  10 alist(maxerr) = a2
496  alist(last) = a1
497  blist(last) = b1
498  rlist(maxerr) = area2
499  rlist(last) = area1
500  elist(maxerr) = error2
501  elist(last) = error1
502 c
503 c call subroutine dqpsrt to maintain the descending ordering
504 c in the list of error estimates and select the subinterval
505 c with the largest error estimate (to be bisected next).
506 c
507  20 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
508 c ***jump out of do-loop
509  if(ier.ne.0.or.errsum.le.errbnd) go to 40
510  30 continue
511 c
512 c compute final result.
513 c ---------------------
514 c
515  40 result = 0.0d+00
516  do 50 k=1,last
517  result = result+rlist(k)
518  50 continue
519  abserr = errsum
520  60 if(keyf.ne.1) neval = (10*keyf+1)*(2*neval+1)
521  if(keyf.eq.1) neval = 30*neval+15
522  999 return
523  end
524  subroutine dqk15(f,a,b,result,abserr,resabs,resasc)
525 c***begin prologue dqk15
526 c***date written 800101 (yymmdd)
527 c***revision date 830518 (yymmdd)
528 c***category no. h2a1a2
529 c***keywords 15-point gauss-kronrod rules
530 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
531 c de doncker,elise,appl. math. & progr. div - k.u.leuven
532 c***purpose to compute i = integral of f over (a,b), with error
533 c estimate
534 c j = integral of abs(f) over (a,b)
535 c***description
536 c
537 c integration rules
538 c standard fortran subroutine
539 c double precision version
540 c
541 c parameters
542 c on entry
543 c f - double precision
544 c function subprogram defining the integrand
545 c function f(x). the actual name for f needs to be
546 c declared e x t e r n a l in the calling program.
547 c
548 c a - double precision
549 c lower limit of integration
550 c
551 c b - double precision
552 c upper limit of integration
553 c
554 c on return
555 c result - double precision
556 c approximation to the integral i
557 c result is computed by applying the 15-point
558 c kronrod rule (resk) obtained by optimal addition
559 c of abscissae to the7-point gauss rule(resg).
560 c
561 c abserr - double precision
562 c estimate of the modulus of the absolute error,
563 c which should not exceed abs(i-result)
564 c
565 c resabs - double precision
566 c approximation to the integral j
567 c
568 c resasc - double precision
569 c approximation to the integral of abs(f-i/(b-a))
570 c over (a,b)
571 c
572 c***references (none)
573 c***routines called d1mach
574 c***end prologue dqk15
575 c
576  double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
577  * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
578  * resg,resk,reskh,result,uflow,wg,wgk,xgk
579  integer j,jtw,jtwm1
580  external f
581 c
582  dimension fv1(7),fv2(7),wg(4),wgk(8),xgk(8)
583 c
584 c the abscissae and weights are given for the interval (-1,1).
585 c because of symmetry only the positive abscissae and their
586 c corresponding weights are given.
587 c
588 c xgk - abscissae of the 15-point kronrod rule
589 c xgk(2), xgk(4), ... abscissae of the 7-point
590 c gauss rule
591 c xgk(1), xgk(3), ... abscissae which are optimally
592 c added to the 7-point gauss rule
593 c
594 c wgk - weights of the 15-point kronrod rule
595 c
596 c wg - weights of the 7-point gauss rule
597 c
598 c
599 c gauss quadrature weights and kronron quadrature abscissae and weights
600 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
601 c bell labs, nov. 1981.
602 c
603  data wg( 1) / 0.1294849661 6886969327 0611432679 082 d0 /
604  data wg( 2) / 0.2797053914 8927666790 1467771423 780 d0 /
605  data wg( 3) / 0.3818300505 0511894495 0369775488 975 d0 /
606  data wg( 4) / 0.4179591836 7346938775 5102040816 327 d0 /
607 c
608  data xgk( 1) / 0.9914553711 2081263920 6854697526 329 d0 /
609  data xgk( 2) / 0.9491079123 4275852452 6189684047 851 d0 /
610  data xgk( 3) / 0.8648644233 5976907278 9712788640 926 d0 /
611  data xgk( 4) / 0.7415311855 9939443986 3864773280 788 d0 /
612  data xgk( 5) / 0.5860872354 6769113029 4144838258 730 d0 /
613  data xgk( 6) / 0.4058451513 7739716690 6606412076 961 d0 /
614  data xgk( 7) / 0.2077849550 0789846760 0689403773 245 d0 /
615  data xgk( 8) / 0.0000000000 0000000000 0000000000 000 d0 /
616 c
617  data wgk( 1) / 0.0229353220 1052922496 3732008058 970 d0 /
618  data wgk( 2) / 0.0630920926 2997855329 0700663189 204 d0 /
619  data wgk( 3) / 0.1047900103 2225018383 9876322541 518 d0 /
620  data wgk( 4) / 0.1406532597 1552591874 5189590510 238 d0 /
621  data wgk( 5) / 0.1690047266 3926790282 6583426598 550 d0 /
622  data wgk( 6) / 0.1903505780 6478540991 3256402421 014 d0 /
623  data wgk( 7) / 0.2044329400 7529889241 4161999234 649 d0 /
624  data wgk( 8) / 0.2094821410 8472782801 2999174891 714 d0 /
625 c
626 c
627 c list of major variables
628 c -----------------------
629 c
630 c centr - mid point of the interval
631 c hlgth - half-length of the interval
632 c absc - abscissa
633 c fval* - function value
634 c resg - result of the 7-point gauss formula
635 c resk - result of the 15-point kronrod formula
636 c reskh - approximation to the mean value of f over (a,b),
637 c i.e. to i/(b-a)
638 c
639 c machine dependent constants
640 c ---------------------------
641 c
642 c epmach is the largest relative spacing.
643 c uflow is the smallest positive magnitude.
644 c
645 c***first executable statement dqk15
646  epmach = d1mach(4)
647  uflow = d1mach(1)
648 c
649  centr = 0.5d+00*(a+b)
650  hlgth = 0.5d+00*(b-a)
651  dhlgth = dabs(hlgth)
652 c
653 c compute the 15-point kronrod approximation to
654 c the integral, and estimate the absolute error.
655 c
656  fc = f(centr)
657  resg = fc*wg(4)
658  resk = fc*wgk(8)
659  resabs = dabs(resk)
660  do 10 j=1,3
661  jtw = j*2
662  absc = hlgth*xgk(jtw)
663  fval1 = f(centr-absc)
664  fval2 = f(centr+absc)
665  fv1(jtw) = fval1
666  fv2(jtw) = fval2
667  fsum = fval1+fval2
668  resg = resg+wg(j)*fsum
669  resk = resk+wgk(jtw)*fsum
670  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
671  10 continue
672  do 15 j = 1,4
673  jtwm1 = j*2-1
674  absc = hlgth*xgk(jtwm1)
675  fval1 = f(centr-absc)
676  fval2 = f(centr+absc)
677  fv1(jtwm1) = fval1
678  fv2(jtwm1) = fval2
679  fsum = fval1+fval2
680  resk = resk+wgk(jtwm1)*fsum
681  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
682  15 continue
683  reskh = resk*0.5d+00
684  resasc = wgk(8)*dabs(fc-reskh)
685  do 20 j=1,7
686  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
687  20 continue
688  result = resk*hlgth
689  resabs = resabs*dhlgth
690  resasc = resasc*dhlgth
691  abserr = dabs((resk-resg)*hlgth)
692  if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
693  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
694  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
695  * ((epmach*0.5d+02)*resabs,abserr)
696  return
697  end
698  subroutine dqk21(f,a,b,result,abserr,resabs,resasc)
699 c***begin prologue dqk21
700 c***date written 800101 (yymmdd)
701 c***revision date 830518 (yymmdd)
702 c***category no. h2a1a2
703 c***keywords 21-point gauss-kronrod rules
704 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
705 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
706 c***purpose to compute i = integral of f over (a,b), with error
707 c estimate
708 c j = integral of abs(f) over (a,b)
709 c***description
710 c
711 c integration rules
712 c standard fortran subroutine
713 c double precision version
714 c
715 c parameters
716 c on entry
717 c f - double precision
718 c function subprogram defining the integrand
719 c function f(x). the actual name for f needs to be
720 c declared e x t e r n a l in the driver program.
721 c
722 c a - double precision
723 c lower limit of integration
724 c
725 c b - double precision
726 c upper limit of integration
727 c
728 c on return
729 c result - double precision
730 c approximation to the integral i
731 c result is computed by applying the 21-point
732 c kronrod rule (resk) obtained by optimal addition
733 c of abscissae to the 10-point gauss rule (resg).
734 c
735 c abserr - double precision
736 c estimate of the modulus of the absolute error,
737 c which should not exceed abs(i-result)
738 c
739 c resabs - double precision
740 c approximation to the integral j
741 c
742 c resasc - double precision
743 c approximation to the integral of abs(f-i/(b-a))
744 c over (a,b)
745 c
746 c***references (none)
747 c***routines called d1mach
748 c***end prologue dqk21
749 c
750  double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
751  * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
752  * resg,resk,reskh,result,uflow,wg,wgk,xgk
753  integer j,jtw,jtwm1
754  external f
755 c
756  dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
757 c
758 c the abscissae and weights are given for the interval (-1,1).
759 c because of symmetry only the positive abscissae and their
760 c corresponding weights are given.
761 c
762 c xgk - abscissae of the 21-point kronrod rule
763 c xgk(2), xgk(4), ... abscissae of the 10-point
764 c gauss rule
765 c xgk(1), xgk(3), ... abscissae which are optimally
766 c added to the 10-point gauss rule
767 c
768 c wgk - weights of the 21-point kronrod rule
769 c
770 c wg - weights of the 10-point gauss rule
771 c
772 c
773 c gauss quadrature weights and kronron quadrature abscissae and weights
774 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
775 c bell labs, nov. 1981.
776 c
777  data wg( 1) / 0.0666713443 0868813759 3568809893 332 d0 /
778  data wg( 2) / 0.1494513491 5058059314 5776339657 697 d0 /
779  data wg( 3) / 0.2190863625 1598204399 5534934228 163 d0 /
780  data wg( 4) / 0.2692667193 0999635509 1226921569 469 d0 /
781  data wg( 5) / 0.2955242247 1475287017 3892994651 338 d0 /
782 c
783  data xgk( 1) / 0.9956571630 2580808073 5527280689 003 d0 /
784  data xgk( 2) / 0.9739065285 1717172007 7964012084 452 d0 /
785  data xgk( 3) / 0.9301574913 5570822600 1207180059 508 d0 /
786  data xgk( 4) / 0.8650633666 8898451073 2096688423 493 d0 /
787  data xgk( 5) / 0.7808177265 8641689706 3717578345 042 d0 /
788  data xgk( 6) / 0.6794095682 9902440623 4327365114 874 d0 /
789  data xgk( 7) / 0.5627571346 6860468333 9000099272 694 d0 /
790  data xgk( 8) / 0.4333953941 2924719079 9265943165 784 d0 /
791  data xgk( 9) / 0.2943928627 0146019813 1126603103 866 d0 /
792  data xgk( 10) / 0.1488743389 8163121088 4826001129 720 d0 /
793  data xgk( 11) / 0.0000000000 0000000000 0000000000 000 d0 /
794 c
795  data wgk( 1) / 0.0116946388 6737187427 8064396062 192 d0 /
796  data wgk( 2) / 0.0325581623 0796472747 8818972459 390 d0 /
797  data wgk( 3) / 0.0547558965 7435199603 1381300244 580 d0 /
798  data wgk( 4) / 0.0750396748 1091995276 7043140916 190 d0 /
799  data wgk( 5) / 0.0931254545 8369760553 5065465083 366 d0 /
800  data wgk( 6) / 0.1093871588 0229764189 9210590325 805 d0 /
801  data wgk( 7) / 0.1234919762 6206585107 7958109831 074 d0 /
802  data wgk( 8) / 0.1347092173 1147332592 8054001771 707 d0 /
803  data wgk( 9) / 0.1427759385 7706008079 7094273138 717 d0 /
804  data wgk( 10) / 0.1477391049 0133849137 4841515972 068 d0 /
805  data wgk( 11) / 0.1494455540 0291690566 4936468389 821 d0 /
806 c
807 c
808 c list of major variables
809 c -----------------------
810 c
811 c centr - mid point of the interval
812 c hlgth - half-length of the interval
813 c absc - abscissa
814 c fval* - function value
815 c resg - result of the 10-point gauss formula
816 c resk - result of the 21-point kronrod formula
817 c reskh - approximation to the mean value of f over (a,b),
818 c i.e. to i/(b-a)
819 c
820 c
821 c machine dependent constants
822 c ---------------------------
823 c
824 c epmach is the largest relative spacing.
825 c uflow is the smallest positive magnitude.
826 c
827 c***first executable statement dqk21
828  epmach = d1mach(4)
829  uflow = d1mach(1)
830 c
831  centr = 0.5d+00*(a+b)
832  hlgth = 0.5d+00*(b-a)
833  dhlgth = dabs(hlgth)
834 c
835 c compute the 21-point kronrod approximation to
836 c the integral, and estimate the absolute error.
837 c
838  resg = 0.0d+00
839  fc = f(centr)
840  resk = wgk(11)*fc
841  resabs = dabs(resk)
842  do 10 j=1,5
843  jtw = 2*j
844  absc = hlgth*xgk(jtw)
845  fval1 = f(centr-absc)
846  fval2 = f(centr+absc)
847  fv1(jtw) = fval1
848  fv2(jtw) = fval2
849  fsum = fval1+fval2
850  resg = resg+wg(j)*fsum
851  resk = resk+wgk(jtw)*fsum
852  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
853  10 continue
854  do 15 j = 1,5
855  jtwm1 = 2*j-1
856  absc = hlgth*xgk(jtwm1)
857  fval1 = f(centr-absc)
858  fval2 = f(centr+absc)
859  fv1(jtwm1) = fval1
860  fv2(jtwm1) = fval2
861  fsum = fval1+fval2
862  resk = resk+wgk(jtwm1)*fsum
863  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
864  15 continue
865  reskh = resk*0.5d+00
866  resasc = wgk(11)*dabs(fc-reskh)
867  do 20 j=1,10
868  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
869  20 continue
870  result = resk*hlgth
871  resabs = resabs*dhlgth
872  resasc = resasc*dhlgth
873  abserr = dabs((resk-resg)*hlgth)
874  if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
875  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
876  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
877  * ((epmach*0.5d+02)*resabs,abserr)
878  return
879  end
880  subroutine dqk31(f,a,b,result,abserr,resabs,resasc)
881 c***begin prologue dqk31
882 c***date written 800101 (yymmdd)
883 c***revision date 830518 (yymmdd)
884 c***category no. h2a1a2
885 c***keywords 31-point gauss-kronrod rules
886 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
887 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
888 c***purpose to compute i = integral of f over (a,b) with error
889 c estimate
890 c j = integral of abs(f) over (a,b)
891 c***description
892 c
893 c integration rules
894 c standard fortran subroutine
895 c double precision version
896 c
897 c parameters
898 c on entry
899 c f - double precision
900 c function subprogram defining the integrand
901 c function f(x). the actual name for f needs to be
902 c declared e x t e r n a l in the calling program.
903 c
904 c a - double precision
905 c lower limit of integration
906 c
907 c b - double precision
908 c upper limit of integration
909 c
910 c on return
911 c result - double precision
912 c approximation to the integral i
913 c result is computed by applying the 31-point
914 c gauss-kronrod rule (resk), obtained by optimal
915 c addition of abscissae to the 15-point gauss
916 c rule (resg).
917 c
918 c abserr - double precison
919 c estimate of the modulus of the modulus,
920 c which should not exceed abs(i-result)
921 c
922 c resabs - double precision
923 c approximation to the integral j
924 c
925 c resasc - double precision
926 c approximation to the integral of abs(f-i/(b-a))
927 c over (a,b)
928 c
929 c***references (none)
930 c***routines called d1mach
931 c***end prologue dqk31
932  double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
933  * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
934  * resg,resk,reskh,result,uflow,wg,wgk,xgk
935  integer j,jtw,jtwm1
936  external f
937 c
938  dimension fv1(15),fv2(15),xgk(16),wgk(16),wg(8)
939 c
940 c the abscissae and weights are given for the interval (-1,1).
941 c because of symmetry only the positive abscissae and their
942 c corresponding weights are given.
943 c
944 c xgk - abscissae of the 31-point kronrod rule
945 c xgk(2), xgk(4), ... abscissae of the 15-point
946 c gauss rule
947 c xgk(1), xgk(3), ... abscissae which are optimally
948 c added to the 15-point gauss rule
949 c
950 c wgk - weights of the 31-point kronrod rule
951 c
952 c wg - weights of the 15-point gauss rule
953 c
954 c
955 c gauss quadrature weights and kronron quadrature abscissae and weights
956 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
957 c bell labs, nov. 1981.
958 c
959  data wg( 1) / 0.0307532419 9611726835 4628393577 204 d0 /
960  data wg( 2) / 0.0703660474 8810812470 9267416450 667 d0 /
961  data wg( 3) / 0.1071592204 6717193501 1869546685 869 d0 /
962  data wg( 4) / 0.1395706779 2615431444 7804794511 028 d0 /
963  data wg( 5) / 0.1662692058 1699393355 3200860481 209 d0 /
964  data wg( 6) / 0.1861610000 1556221102 6800561866 423 d0 /
965  data wg( 7) / 0.1984314853 2711157645 6118326443 839 d0 /
966  data wg( 8) / 0.2025782419 2556127288 0620199967 519 d0 /
967 c
968  data xgk( 1) / 0.9980022986 9339706028 5172840152 271 d0 /
969  data xgk( 2) / 0.9879925180 2048542848 9565718586 613 d0 /
970  data xgk( 3) / 0.9677390756 7913913425 7347978784 337 d0 /
971  data xgk( 4) / 0.9372733924 0070590430 7758947710 209 d0 /
972  data xgk( 5) / 0.8972645323 4408190088 2509656454 496 d0 /
973  data xgk( 6) / 0.8482065834 1042721620 0648320774 217 d0 /
974  data xgk( 7) / 0.7904185014 4246593296 7649294817 947 d0 /
975  data xgk( 8) / 0.7244177313 6017004741 6186054613 938 d0 /
976  data xgk( 9) / 0.6509967412 9741697053 3735895313 275 d0 /
977  data xgk( 10) / 0.5709721726 0853884753 7226737253 911 d0 /
978  data xgk( 11) / 0.4850818636 4023968069 3655740232 351 d0 /
979  data xgk( 12) / 0.3941513470 7756336989 7207370981 045 d0 /
980  data xgk( 13) / 0.2991800071 5316881216 6780024266 389 d0 /
981  data xgk( 14) / 0.2011940939 9743452230 0628303394 596 d0 /
982  data xgk( 15) / 0.1011420669 1871749902 7074231447 392 d0 /
983  data xgk( 16) / 0.0000000000 0000000000 0000000000 000 d0 /
984 c
985  data wgk( 1) / 0.0053774798 7292334898 7792051430 128 d0 /
986  data wgk( 2) / 0.0150079473 2931612253 8374763075 807 d0 /
987  data wgk( 3) / 0.0254608473 2671532018 6874001019 653 d0 /
988  data wgk( 4) / 0.0353463607 9137584622 2037948478 360 d0 /
989  data wgk( 5) / 0.0445897513 2476487660 8227299373 280 d0 /
990  data wgk( 6) / 0.0534815246 9092808726 5343147239 430 d0 /
991  data wgk( 7) / 0.0620095678 0067064028 5139230960 803 d0 /
992  data wgk( 8) / 0.0698541213 1872825870 9520077099 147 d0 /
993  data wgk( 9) / 0.0768496807 5772037889 4432777482 659 d0 /
994  data wgk( 10) / 0.0830805028 2313302103 8289247286 104 d0 /
995  data wgk( 11) / 0.0885644430 5621177064 7275443693 774 d0 /
996  data wgk( 12) / 0.0931265981 7082532122 5486872747 346 d0 /
997  data wgk( 13) / 0.0966427269 8362367850 5179907627 589 d0 /
998  data wgk( 14) / 0.0991735987 2179195933 2393173484 603 d0 /
999  data wgk( 15) / 0.1007698455 2387559504 4946662617 570 d0 /
1000  data wgk( 16) / 0.1013300070 1479154901 7374792767 493 d0 /
1001 c
1002 c
1003 c list of major variables
1004 c -----------------------
1005 c centr - mid point of the interval
1006 c hlgth - half-length of the interval
1007 c absc - abscissa
1008 c fval* - function value
1009 c resg - result of the 15-point gauss formula
1010 c resk - result of the 31-point kronrod formula
1011 c reskh - approximation to the mean value of f over (a,b),
1012 c i.e. to i/(b-a)
1013 c
1014 c machine dependent constants
1015 c ---------------------------
1016 c epmach is the largest relative spacing.
1017 c uflow is the smallest positive magnitude.
1018 c***first executable statement dqk31
1019  epmach = d1mach(4)
1020  uflow = d1mach(1)
1021 c
1022  centr = 0.5d+00*(a+b)
1023  hlgth = 0.5d+00*(b-a)
1024  dhlgth = dabs(hlgth)
1025 c
1026 c compute the 31-point kronrod approximation to
1027 c the integral, and estimate the absolute error.
1028 c
1029  fc = f(centr)
1030  resg = wg(8)*fc
1031  resk = wgk(16)*fc
1032  resabs = dabs(resk)
1033  do 10 j=1,7
1034  jtw = j*2
1035  absc = hlgth*xgk(jtw)
1036  fval1 = f(centr-absc)
1037  fval2 = f(centr+absc)
1038  fv1(jtw) = fval1
1039  fv2(jtw) = fval2
1040  fsum = fval1+fval2
1041  resg = resg+wg(j)*fsum
1042  resk = resk+wgk(jtw)*fsum
1043  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1044  10 continue
1045  do 15 j = 1,8
1046  jtwm1 = j*2-1
1047  absc = hlgth*xgk(jtwm1)
1048  fval1 = f(centr-absc)
1049  fval2 = f(centr+absc)
1050  fv1(jtwm1) = fval1
1051  fv2(jtwm1) = fval2
1052  fsum = fval1+fval2
1053  resk = resk+wgk(jtwm1)*fsum
1054  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1055  15 continue
1056  reskh = resk*0.5d+00
1057  resasc = wgk(16)*dabs(fc-reskh)
1058  do 20 j=1,15
1059  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1060  20 continue
1061  result = resk*hlgth
1062  resabs = resabs*dhlgth
1063  resasc = resasc*dhlgth
1064  abserr = dabs((resk-resg)*hlgth)
1065  if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
1066  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1067  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1068  * ((epmach*0.5d+02)*resabs,abserr)
1069  return
1070  end
1071  subroutine dqk41(f,a,b,result,abserr,resabs,resasc)
1072 c***begin prologue dqk41
1073 c***date written 800101 (yymmdd)
1074 c***revision date 830518 (yymmdd)
1075 c***category no. h2a1a2
1076 c***keywords 41-point gauss-kronrod rules
1077 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1078 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
1079 c***purpose to compute i = integral of f over (a,b), with error
1080 c estimate
1081 c j = integral of abs(f) over (a,b)
1082 c***description
1083 c
1084 c integration rules
1085 c standard fortran subroutine
1086 c double precision version
1087 c
1088 c parameters
1089 c on entry
1090 c f - double precision
1091 c function subprogram defining the integrand
1092 c function f(x). the actual name for f needs to be
1093 c declared e x t e r n a l in the calling program.
1094 c
1095 c a - double precision
1096 c lower limit of integration
1097 c
1098 c b - double precision
1099 c upper limit of integration
1100 c
1101 c on return
1102 c result - double precision
1103 c approximation to the integral i
1104 c result is computed by applying the 41-point
1105 c gauss-kronrod rule (resk) obtained by optimal
1106 c addition of abscissae to the 20-point gauss
1107 c rule (resg).
1108 c
1109 c abserr - double precision
1110 c estimate of the modulus of the absolute error,
1111 c which should not exceed abs(i-result)
1112 c
1113 c resabs - double precision
1114 c approximation to the integral j
1115 c
1116 c resasc - double precision
1117 c approximation to the integal of abs(f-i/(b-a))
1118 c over (a,b)
1119 c
1120 c***references (none)
1121 c***routines called d1mach
1122 c***end prologue dqk41
1123 c
1124  double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
1125  * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
1126  * resg,resk,reskh,result,uflow,wg,wgk,xgk
1127  integer j,jtw,jtwm1
1128  external f
1129 c
1130  dimension fv1(20),fv2(20),xgk(21),wgk(21),wg(10)
1131 c
1132 c the abscissae and weights are given for the interval (-1,1).
1133 c because of symmetry only the positive abscissae and their
1134 c corresponding weights are given.
1135 c
1136 c xgk - abscissae of the 41-point gauss-kronrod rule
1137 c xgk(2), xgk(4), ... abscissae of the 20-point
1138 c gauss rule
1139 c xgk(1), xgk(3), ... abscissae which are optimally
1140 c added to the 20-point gauss rule
1141 c
1142 c wgk - weights of the 41-point gauss-kronrod rule
1143 c
1144 c wg - weights of the 20-point gauss rule
1145 c
1146 c
1147 c gauss quadrature weights and kronron quadrature abscissae and weights
1148 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
1149 c bell labs, nov. 1981.
1150 c
1151  data wg( 1) / 0.0176140071 3915211831 1861962351 853 d0 /
1152  data wg( 2) / 0.0406014298 0038694133 1039952274 932 d0 /
1153  data wg( 3) / 0.0626720483 3410906356 9506535187 042 d0 /
1154  data wg( 4) / 0.0832767415 7670474872 4758143222 046 d0 /
1155  data wg( 5) / 0.1019301198 1724043503 6750135480 350 d0 /
1156  data wg( 6) / 0.1181945319 6151841731 2377377711 382 d0 /
1157  data wg( 7) / 0.1316886384 4917662689 8494499748 163 d0 /
1158  data wg( 8) / 0.1420961093 1838205132 9298325067 165 d0 /
1159  data wg( 9) / 0.1491729864 7260374678 7828737001 969 d0 /
1160  data wg( 10) / 0.1527533871 3072585069 8084331955 098 d0 /
1161 c
1162  data xgk( 1) / 0.9988590315 8827766383 8315576545 863 d0 /
1163  data xgk( 2) / 0.9931285991 8509492478 6122388471 320 d0 /
1164  data xgk( 3) / 0.9815078774 5025025919 3342994720 217 d0 /
1165  data xgk( 4) / 0.9639719272 7791379126 7666131197 277 d0 /
1166  data xgk( 5) / 0.9408226338 3175475351 9982722212 443 d0 /
1167  data xgk( 6) / 0.9122344282 5132590586 7752441203 298 d0 /
1168  data xgk( 7) / 0.8782768112 5228197607 7442995113 078 d0 /
1169  data xgk( 8) / 0.8391169718 2221882339 4529061701 521 d0 /
1170  data xgk( 9) / 0.7950414288 3755119835 0638833272 788 d0 /
1171  data xgk( 10) / 0.7463319064 6015079261 4305070355 642 d0 /
1172  data xgk( 11) / 0.6932376563 3475138480 5490711845 932 d0 /
1173  data xgk( 12) / 0.6360536807 2651502545 2836696226 286 d0 /
1174  data xgk( 13) / 0.5751404468 1971031534 2946036586 425 d0 /
1175  data xgk( 14) / 0.5108670019 5082709800 4364050955 251 d0 /
1176  data xgk( 15) / 0.4435931752 3872510319 9992213492 640 d0 /
1177  data xgk( 16) / 0.3737060887 1541956067 2548177024 927 d0 /
1178  data xgk( 17) / 0.3016278681 1491300432 0555356858 592 d0 /
1179  data xgk( 18) / 0.2277858511 4164507808 0496195368 575 d0 /
1180  data xgk( 19) / 0.1526054652 4092267550 5220241022 678 d0 /
1181  data xgk( 20) / 0.0765265211 3349733375 4640409398 838 d0 /
1182  data xgk( 21) / 0.0000000000 0000000000 0000000000 000 d0 /
1183 c
1184  data wgk( 1) / 0.0030735837 1852053150 1218293246 031 d0 /
1185  data wgk( 2) / 0.0086002698 5564294219 8661787950 102 d0 /
1186  data wgk( 3) / 0.0146261692 5697125298 3787960308 868 d0 /
1187  data wgk( 4) / 0.0203883734 6126652359 8010231432 755 d0 /
1188  data wgk( 5) / 0.0258821336 0495115883 4505067096 153 d0 /
1189  data wgk( 6) / 0.0312873067 7703279895 8543119323 801 d0 /
1190  data wgk( 7) / 0.0366001697 5820079803 0557240707 211 d0 /
1191  data wgk( 8) / 0.0416688733 2797368626 3788305936 895 d0 /
1192  data wgk( 9) / 0.0464348218 6749767472 0231880926 108 d0 /
1193  data wgk( 10) / 0.0509445739 2372869193 2707670050 345 d0 /
1194  data wgk( 11) / 0.0551951053 4828599474 4832372419 777 d0 /
1195  data wgk( 12) / 0.0591114008 8063957237 4967220648 594 d0 /
1196  data wgk( 13) / 0.0626532375 5478116802 5870122174 255 d0 /
1197  data wgk( 14) / 0.0658345971 3361842211 1563556969 398 d0 /
1198  data wgk( 15) / 0.0686486729 2852161934 5623411885 368 d0 /
1199  data wgk( 16) / 0.0710544235 5344406830 5790361723 210 d0 /
1200  data wgk( 17) / 0.0730306903 3278666749 5189417658 913 d0 /
1201  data wgk( 18) / 0.0745828754 0049918898 6581418362 488 d0 /
1202  data wgk( 19) / 0.0757044976 8455667465 9542775376 617 d0 /
1203  data wgk( 20) / 0.0763778676 7208073670 5502835038 061 d0 /
1204  data wgk( 21) / 0.0766007119 1799965644 5049901530 102 d0 /
1205 c
1206 c
1207 c list of major variables
1208 c -----------------------
1209 c
1210 c centr - mid point of the interval
1211 c hlgth - half-length of the interval
1212 c absc - abscissa
1213 c fval* - function value
1214 c resg - result of the 20-point gauss formula
1215 c resk - result of the 41-point kronrod formula
1216 c reskh - approximation to mean value of f over (a,b), i.e.
1217 c to i/(b-a)
1218 c
1219 c machine dependent constants
1220 c ---------------------------
1221 c
1222 c epmach is the largest relative spacing.
1223 c uflow is the smallest positive magnitude.
1224 c
1225 c***first executable statement dqk41
1226  epmach = d1mach(4)
1227  uflow = d1mach(1)
1228 c
1229  centr = 0.5d+00*(a+b)
1230  hlgth = 0.5d+00*(b-a)
1231  dhlgth = dabs(hlgth)
1232 c
1233 c compute the 41-point gauss-kronrod approximation to
1234 c the integral, and estimate the absolute error.
1235 c
1236  resg = 0.0d+00
1237  fc = f(centr)
1238  resk = wgk(21)*fc
1239  resabs = dabs(resk)
1240  do 10 j=1,10
1241  jtw = j*2
1242  absc = hlgth*xgk(jtw)
1243  fval1 = f(centr-absc)
1244  fval2 = f(centr+absc)
1245  fv1(jtw) = fval1
1246  fv2(jtw) = fval2
1247  fsum = fval1+fval2
1248  resg = resg+wg(j)*fsum
1249  resk = resk+wgk(jtw)*fsum
1250  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1251  10 continue
1252  do 15 j = 1,10
1253  jtwm1 = j*2-1
1254  absc = hlgth*xgk(jtwm1)
1255  fval1 = f(centr-absc)
1256  fval2 = f(centr+absc)
1257  fv1(jtwm1) = fval1
1258  fv2(jtwm1) = fval2
1259  fsum = fval1+fval2
1260  resk = resk+wgk(jtwm1)*fsum
1261  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1262  15 continue
1263  reskh = resk*0.5d+00
1264  resasc = wgk(21)*dabs(fc-reskh)
1265  do 20 j=1,20
1266  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1267  20 continue
1268  result = resk*hlgth
1269  resabs = resabs*dhlgth
1270  resasc = resasc*dhlgth
1271  abserr = dabs((resk-resg)*hlgth)
1272  if(resasc.ne.0.0d+00.and.abserr.ne.0.d+00)
1273  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1274  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1275  * ((epmach*0.5d+02)*resabs,abserr)
1276  return
1277  end
1278  subroutine dqk51(f,a,b,result,abserr,resabs,resasc)
1279 c***begin prologue dqk51
1280 c***date written 800101 (yymmdd)
1281 c***revision date 830518 (yymmdd)
1282 c***category no. h2a1a2
1283 c***keywords 51-point gauss-kronrod rules
1284 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1285 c de doncker,elise,appl. math & progr. div. - k.u.leuven
1286 c***purpose to compute i = integral of f over (a,b) with error
1287 c estimate
1288 c j = integral of abs(f) over (a,b)
1289 c***description
1290 c
1291 c integration rules
1292 c standard fortran subroutine
1293 c double precision version
1294 c
1295 c parameters
1296 c on entry
1297 c f - double precision
1298 c function subroutine defining the integrand
1299 c function f(x). the actual name for f needs to be
1300 c declared e x t e r n a l in the calling program.
1301 c
1302 c a - double precision
1303 c lower limit of integration
1304 c
1305 c b - double precision
1306 c upper limit of integration
1307 c
1308 c on return
1309 c result - double precision
1310 c approximation to the integral i
1311 c result is computed by applying the 51-point
1312 c kronrod rule (resk) obtained by optimal addition
1313 c of abscissae to the 25-point gauss rule (resg).
1314 c
1315 c abserr - double precision
1316 c estimate of the modulus of the absolute error,
1317 c which should not exceed abs(i-result)
1318 c
1319 c resabs - double precision
1320 c approximation to the integral j
1321 c
1322 c resasc - double precision
1323 c approximation to the integral of abs(f-i/(b-a))
1324 c over (a,b)
1325 c
1326 c***references (none)
1327 c***routines called d1mach
1328 c***end prologue dqk51
1329 c
1330  double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
1331  * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
1332  * resg,resk,reskh,result,uflow,wg,wgk,xgk
1333  integer j,jtw,jtwm1
1334  external f
1335 c
1336  dimension fv1(25),fv2(25),xgk(26),wgk(26),wg(13)
1337 c
1338 c the abscissae and weights are given for the interval (-1,1).
1339 c because of symmetry only the positive abscissae and their
1340 c corresponding weights are given.
1341 c
1342 c xgk - abscissae of the 51-point kronrod rule
1343 c xgk(2), xgk(4), ... abscissae of the 25-point
1344 c gauss rule
1345 c xgk(1), xgk(3), ... abscissae which are optimally
1346 c added to the 25-point gauss rule
1347 c
1348 c wgk - weights of the 51-point kronrod rule
1349 c
1350 c wg - weights of the 25-point gauss rule
1351 c
1352 c
1353 c gauss quadrature weights and kronron quadrature abscissae and weights
1354 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
1355 c bell labs, nov. 1981.
1356 c
1357  data wg( 1) / 0.0113937985 0102628794 7902964113 235 d0 /
1358  data wg( 2) / 0.0263549866 1503213726 1901815295 299 d0 /
1359  data wg( 3) / 0.0409391567 0130631265 5623487711 646 d0 /
1360  data wg( 4) / 0.0549046959 7583519192 5936891540 473 d0 /
1361  data wg( 5) / 0.0680383338 1235691720 7187185656 708 d0 /
1362  data wg( 6) / 0.0801407003 3500101801 3234959669 111 d0 /
1363  data wg( 7) / 0.0910282619 8296364981 1497220702 892 d0 /
1364  data wg( 8) / 0.1005359490 6705064420 2206890392 686 d0 /
1365  data wg( 9) / 0.1085196244 7426365311 6093957050 117 d0 /
1366  data wg( 10) / 0.1148582591 4571164833 9325545869 556 d0 /
1367  data wg( 11) / 0.1194557635 3578477222 8178126512 901 d0 /
1368  data wg( 12) / 0.1222424429 9031004168 8959518945 852 d0 /
1369  data wg( 13) / 0.1231760537 2671545120 3902873079 050 d0 /
1370 c
1371  data xgk( 1) / 0.9992621049 9260983419 3457486540 341 d0 /
1372  data xgk( 2) / 0.9955569697 9049809790 8784946893 902 d0 /
1373  data xgk( 3) / 0.9880357945 3407724763 7331014577 406 d0 /
1374  data xgk( 4) / 0.9766639214 5951751149 8315386479 594 d0 /
1375  data xgk( 5) / 0.9616149864 2584251241 8130033660 167 d0 /
1376  data xgk( 6) / 0.9429745712 2897433941 4011169658 471 d0 /
1377  data xgk( 7) / 0.9207471152 8170156174 6346084546 331 d0 /
1378  data xgk( 8) / 0.8949919978 7827536885 1042006782 805 d0 /
1379  data xgk( 9) / 0.8658470652 9327559544 8996969588 340 d0 /
1380  data xgk( 10) / 0.8334426287 6083400142 1021108693 570 d0 /
1381  data xgk( 11) / 0.7978737979 9850005941 0410904994 307 d0 /
1382  data xgk( 12) / 0.7592592630 3735763057 7282865204 361 d0 /
1383  data xgk( 13) / 0.7177664068 1308438818 6654079773 298 d0 /
1384  data xgk( 14) / 0.6735663684 7346836448 5120633247 622 d0 /
1385  data xgk( 15) / 0.6268100990 1031741278 8122681624 518 d0 /
1386  data xgk( 16) / 0.5776629302 4122296772 3689841612 654 d0 /
1387  data xgk( 17) / 0.5263252843 3471918259 9623778158 010 d0 /
1388  data xgk( 18) / 0.4730027314 4571496052 2182115009 192 d0 /
1389  data xgk( 19) / 0.4178853821 9303774885 1814394594 572 d0 /
1390  data xgk( 20) / 0.3611723058 0938783773 5821730127 641 d0 /
1391  data xgk( 21) / 0.3030895389 3110783016 7478909980 339 d0 /
1392  data xgk( 22) / 0.2438668837 2098843204 5190362797 452 d0 /
1393  data xgk( 23) / 0.1837189394 2104889201 5969888759 528 d0 /
1394  data xgk( 24) / 0.1228646926 1071039638 7359818808 037 d0 /
1395  data xgk( 25) / 0.0615444830 0568507888 6546392366 797 d0 /
1396  data xgk( 26) / 0.0000000000 0000000000 0000000000 000 d0 /
1397 c
1398  data wgk( 1) / 0.0019873838 9233031592 6507851882 843 d0 /
1399  data wgk( 2) / 0.0055619321 3535671375 8040236901 066 d0 /
1400  data wgk( 3) / 0.0094739733 8617415160 7207710523 655 d0 /
1401  data wgk( 4) / 0.0132362291 9557167481 3656405846 976 d0 /
1402  data wgk( 5) / 0.0168478177 0912829823 1516667536 336 d0 /
1403  data wgk( 6) / 0.0204353711 4588283545 6568292235 939 d0 /
1404  data wgk( 7) / 0.0240099456 0695321622 0092489164 881 d0 /
1405  data wgk( 8) / 0.0274753175 8785173780 2948455517 811 d0 /
1406  data wgk( 9) / 0.0307923001 6738748889 1109020215 229 d0 /
1407  data wgk( 10) / 0.0340021302 7432933783 6748795229 551 d0 /
1408  data wgk( 11) / 0.0371162714 8341554356 0330625367 620 d0 /
1409  data wgk( 12) / 0.0400838255 0403238207 4839284467 076 d0 /
1410  data wgk( 13) / 0.0428728450 2017004947 6895792439 495 d0 /
1411  data wgk( 14) / 0.0455029130 4992178890 9870584752 660 d0 /
1412  data wgk( 15) / 0.0479825371 3883671390 6392255756 915 d0 /
1413  data wgk( 16) / 0.0502776790 8071567196 3325259433 440 d0 /
1414  data wgk( 17) / 0.0523628858 0640747586 4366712137 873 d0 /
1415  data wgk( 18) / 0.0542511298 8854549014 4543370459 876 d0 /
1416  data wgk( 19) / 0.0559508112 2041231730 8240686382 747 d0 /
1417  data wgk( 20) / 0.0574371163 6156783285 3582693939 506 d0 /
1418  data wgk( 21) / 0.0586896800 2239420796 1974175856 788 d0 /
1419  data wgk( 22) / 0.0597203403 2417405997 9099291932 562 d0 /
1420  data wgk( 23) / 0.0605394553 7604586294 5360267517 565 d0 /
1421  data wgk( 24) / 0.0611285097 1705304830 5859030416 293 d0 /
1422  data wgk( 25) / 0.0614711898 7142531666 1544131965 264 d0 /
1423 c note: wgk (26) was calculated from the values of wgk(1..25)
1424  data wgk( 26) / 0.0615808180 6783293507 8759824240 066 d0 /
1425 c
1426 c
1427 c list of major variables
1428 c -----------------------
1429 c
1430 c centr - mid point of the interval
1431 c hlgth - half-length of the interval
1432 c absc - abscissa
1433 c fval* - function value
1434 c resg - result of the 25-point gauss formula
1435 c resk - result of the 51-point kronrod formula
1436 c reskh - approximation to the mean value of f over (a,b),
1437 c i.e. to i/(b-a)
1438 c
1439 c machine dependent constants
1440 c ---------------------------
1441 c
1442 c epmach is the largest relative spacing.
1443 c uflow is the smallest positive magnitude.
1444 c
1445 c***first executable statement dqk51
1446  epmach = d1mach(4)
1447  uflow = d1mach(1)
1448 c
1449  centr = 0.5d+00*(a+b)
1450  hlgth = 0.5d+00*(b-a)
1451  dhlgth = dabs(hlgth)
1452 c
1453 c compute the 51-point kronrod approximation to
1454 c the integral, and estimate the absolute error.
1455 c
1456  fc = f(centr)
1457  resg = wg(13)*fc
1458  resk = wgk(26)*fc
1459  resabs = dabs(resk)
1460  do 10 j=1,12
1461  jtw = j*2
1462  absc = hlgth*xgk(jtw)
1463  fval1 = f(centr-absc)
1464  fval2 = f(centr+absc)
1465  fv1(jtw) = fval1
1466  fv2(jtw) = fval2
1467  fsum = fval1+fval2
1468  resg = resg+wg(j)*fsum
1469  resk = resk+wgk(jtw)*fsum
1470  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1471  10 continue
1472  do 15 j = 1,13
1473  jtwm1 = j*2-1
1474  absc = hlgth*xgk(jtwm1)
1475  fval1 = f(centr-absc)
1476  fval2 = f(centr+absc)
1477  fv1(jtwm1) = fval1
1478  fv2(jtwm1) = fval2
1479  fsum = fval1+fval2
1480  resk = resk+wgk(jtwm1)*fsum
1481  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1482  15 continue
1483  reskh = resk*0.5d+00
1484  resasc = wgk(26)*dabs(fc-reskh)
1485  do 20 j=1,25
1486  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1487  20 continue
1488  result = resk*hlgth
1489  resabs = resabs*dhlgth
1490  resasc = resasc*dhlgth
1491  abserr = dabs((resk-resg)*hlgth)
1492  if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
1493  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1494  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1495  * ((epmach*0.5d+02)*resabs,abserr)
1496  return
1497  end
1498  subroutine dqk61(f,a,b,result,abserr,resabs,resasc)
1499 c***begin prologue dqk61
1500 c***date written 800101 (yymmdd)
1501 c***revision date 830518 (yymmdd)
1502 c***category no. h2a1a2
1503 c***keywords 61-point gauss-kronrod rules
1504 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1505 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
1506 c***purpose to compute i = integral of f over (a,b) with error
1507 c estimate
1508 c j = integral of dabs(f) over (a,b)
1509 c***description
1510 c
1511 c integration rule
1512 c standard fortran subroutine
1513 c double precision version
1514 c
1515 c
1516 c parameters
1517 c on entry
1518 c f - double precision
1519 c function subprogram defining the integrand
1520 c function f(x). the actual name for f needs to be
1521 c declared e x t e r n a l in the calling program.
1522 c
1523 c a - double precision
1524 c lower limit of integration
1525 c
1526 c b - double precision
1527 c upper limit of integration
1528 c
1529 c on return
1530 c result - double precision
1531 c approximation to the integral i
1532 c result is computed by applying the 61-point
1533 c kronrod rule (resk) obtained by optimal addition of
1534 c abscissae to the 30-point gauss rule (resg).
1535 c
1536 c abserr - double precision
1537 c estimate of the modulus of the absolute error,
1538 c which should equal or exceed dabs(i-result)
1539 c
1540 c resabs - double precision
1541 c approximation to the integral j
1542 c
1543 c resasc - double precision
1544 c approximation to the integral of dabs(f-i/(b-a))
1545 c
1546 c
1547 c***references (none)
1548 c***routines called d1mach
1549 c***end prologue dqk61
1550 c
1551  double precision a,dabsc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
1552  * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
1553  * resg,resk,reskh,result,uflow,wg,wgk,xgk
1554  integer j,jtw,jtwm1
1555  external f
1556 c
1557  dimension fv1(30),fv2(30),xgk(31),wgk(31),wg(15)
1558 c
1559 c the abscissae and weights are given for the
1560 c interval (-1,1). because of symmetry only the positive
1561 c abscissae and their corresponding weights are given.
1562 c
1563 c xgk - abscissae of the 61-point kronrod rule
1564 c xgk(2), xgk(4) ... abscissae of the 30-point
1565 c gauss rule
1566 c xgk(1), xgk(3) ... optimally added abscissae
1567 c to the 30-point gauss rule
1568 c
1569 c wgk - weights of the 61-point kronrod rule
1570 c
1571 c wg - weigths of the 30-point gauss rule
1572 c
1573 c
1574 c gauss quadrature weights and kronron quadrature abscissae and weights
1575 c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
1576 c bell labs, nov. 1981.
1577 c
1578  data wg( 1) / 0.0079681924 9616660561 5465883474 674 d0 /
1579  data wg( 2) / 0.0184664683 1109095914 2302131912 047 d0 /
1580  data wg( 3) / 0.0287847078 8332336934 9719179611 292 d0 /
1581  data wg( 4) / 0.0387991925 6962704959 6801936446 348 d0 /
1582  data wg( 5) / 0.0484026728 3059405290 2938140422 808 d0 /
1583  data wg( 6) / 0.0574931562 1761906648 1721689402 056 d0 /
1584  data wg( 7) / 0.0659742298 8218049512 8128515115 962 d0 /
1585  data wg( 8) / 0.0737559747 3770520626 8243850022 191 d0 /
1586  data wg( 9) / 0.0807558952 2942021535 4694938460 530 d0 /
1587  data wg( 10) / 0.0868997872 0108297980 2387530715 126 d0 /
1588  data wg( 11) / 0.0921225222 3778612871 7632707087 619 d0 /
1589  data wg( 12) / 0.0963687371 7464425963 9468626351 810 d0 /
1590  data wg( 13) / 0.0995934205 8679526706 2780282103 569 d0 /
1591  data wg( 14) / 0.1017623897 4840550459 6428952168 554 d0 /
1592  data wg( 15) / 0.1028526528 9355884034 1285636705 415 d0 /
1593 c
1594  data xgk( 1) / 0.9994844100 5049063757 1325895705 811 d0 /
1595  data xgk( 2) / 0.9968934840 7464954027 1630050918 695 d0 /
1596  data xgk( 3) / 0.9916309968 7040459485 8628366109 486 d0 /
1597  data xgk( 4) / 0.9836681232 7974720997 0032581605 663 d0 /
1598  data xgk( 5) / 0.9731163225 0112626837 4693868423 707 d0 /
1599  data xgk( 6) / 0.9600218649 6830751221 6871025581 798 d0 /
1600  data xgk( 7) / 0.9443744447 4855997941 5831324037 439 d0 /
1601  data xgk( 8) / 0.9262000474 2927432587 9324277080 474 d0 /
1602  data xgk( 9) / 0.9055733076 9990779854 6522558925 958 d0 /
1603  data xgk( 10) / 0.8825605357 9205268154 3116462530 226 d0 /
1604  data xgk( 11) / 0.8572052335 4606109895 8658510658 944 d0 /
1605  data xgk( 12) / 0.8295657623 8276839744 2898119732 502 d0 /
1606  data xgk( 13) / 0.7997278358 2183908301 3668942322 683 d0 /
1607  data xgk( 14) / 0.7677774321 0482619491 7977340974 503 d0 /
1608  data xgk( 15) / 0.7337900624 5322680472 6171131369 528 d0 /
1609  data xgk( 16) / 0.6978504947 9331579693 2292388026 640 d0 /
1610  data xgk( 17) / 0.6600610641 2662696137 0053668149 271 d0 /
1611  data xgk( 18) / 0.6205261829 8924286114 0477556431 189 d0 /
1612  data xgk( 19) / 0.5793452358 2636169175 6024932172 540 d0 /
1613  data xgk( 20) / 0.5366241481 4201989926 4169793311 073 d0 /
1614  data xgk( 21) / 0.4924804678 6177857499 3693061207 709 d0 /
1615  data xgk( 22) / 0.4470337695 3808917678 0609900322 854 d0 /
1616  data xgk( 23) / 0.4004012548 3039439253 5476211542 661 d0 /
1617  data xgk( 24) / 0.3527047255 3087811347 1037207089 374 d0 /
1618  data xgk( 25) / 0.3040732022 7362507737 2677107199 257 d0 /
1619  data xgk( 26) / 0.2546369261 6788984643 9805129817 805 d0 /
1620  data xgk( 27) / 0.2045251166 8230989143 8957671002 025 d0 /
1621  data xgk( 28) / 0.1538699136 0858354696 3794672743 256 d0 /
1622  data xgk( 29) / 0.1028069379 6673703014 7096751318 001 d0 /
1623  data xgk( 30) / 0.0514718425 5531769583 3025213166 723 d0 /
1624  data xgk( 31) / 0.0000000000 0000000000 0000000000 000 d0 /
1625 c
1626  data wgk( 1) / 0.0013890136 9867700762 4551591226 760 d0 /
1627  data wgk( 2) / 0.0038904611 2709988405 1267201844 516 d0 /
1628  data wgk( 3) / 0.0066307039 1593129217 3319826369 750 d0 /
1629  data wgk( 4) / 0.0092732796 5951776342 8441146892 024 d0 /
1630  data wgk( 5) / 0.0118230152 5349634174 2232898853 251 d0 /
1631  data wgk( 6) / 0.0143697295 0704580481 2451432443 580 d0 /
1632  data wgk( 7) / 0.0169208891 8905327262 7572289420 322 d0 /
1633  data wgk( 8) / 0.0194141411 9394238117 3408951050 128 d0 /
1634  data wgk( 9) / 0.0218280358 2160919229 7167485738 339 d0 /
1635  data wgk( 10) / 0.0241911620 7808060136 5686370725 232 d0 /
1636  data wgk( 11) / 0.0265099548 8233310161 0601709335 075 d0 /
1637  data wgk( 12) / 0.0287540487 6504129284 3978785354 334 d0 /
1638  data wgk( 13) / 0.0309072575 6238776247 2884252943 092 d0 /
1639  data wgk( 14) / 0.0329814470 5748372603 1814191016 854 d0 /
1640  data wgk( 15) / 0.0349793380 2806002413 7499670731 468 d0 /
1641  data wgk( 16) / 0.0368823646 5182122922 3911065617 136 d0 /
1642  data wgk( 17) / 0.0386789456 2472759295 0348651532 281 d0 /
1643  data wgk( 18) / 0.0403745389 5153595911 1995279752 468 d0 /
1644  data wgk( 19) / 0.0419698102 1516424614 7147541285 970 d0 /
1645  data wgk( 20) / 0.0434525397 0135606931 6831728117 073 d0 /
1646  data wgk( 21) / 0.0448148001 3316266319 2355551616 723 d0 /
1647  data wgk( 22) / 0.0460592382 7100698811 6271735559 374 d0 /
1648  data wgk( 23) / 0.0471855465 6929915394 5261478181 099 d0 /
1649  data wgk( 24) / 0.0481858617 5708712914 0779492298 305 d0 /
1650  data wgk( 25) / 0.0490554345 5502977888 7528165367 238 d0 /
1651  data wgk( 26) / 0.0497956834 2707420635 7811569379 942 d0 /
1652  data wgk( 27) / 0.0504059214 0278234684 0893085653 585 d0 /
1653  data wgk( 28) / 0.0508817958 9874960649 2297473049 805 d0 /
1654  data wgk( 29) / 0.0512215478 4925877217 0656282604 944 d0 /
1655  data wgk( 30) / 0.0514261285 3745902593 3862879215 781 d0 /
1656  data wgk( 31) / 0.0514947294 2945156755 8340433647 099 d0 /
1657 c
1658 c list of major variables
1659 c -----------------------
1660 c
1661 c centr - mid point of the interval
1662 c hlgth - half-length of the interval
1663 c dabsc - abscissa
1664 c fval* - function value
1665 c resg - result of the 30-point gauss rule
1666 c resk - result of the 61-point kronrod rule
1667 c reskh - approximation to the mean value of f
1668 c over (a,b), i.e. to i/(b-a)
1669 c
1670 c machine dependent constants
1671 c ---------------------------
1672 c
1673 c epmach is the largest relative spacing.
1674 c uflow is the smallest positive magnitude.
1675 c
1676  epmach = d1mach(4)
1677  uflow = d1mach(1)
1678 c
1679  centr = 0.5d+00*(b+a)
1680  hlgth = 0.5d+00*(b-a)
1681  dhlgth = dabs(hlgth)
1682 c
1683 c compute the 61-point kronrod approximation to the
1684 c integral, and estimate the absolute error.
1685 c
1686 c***first executable statement dqk61
1687  resg = 0.0d+00
1688  fc = f(centr)
1689  resk = wgk(31)*fc
1690  resabs = dabs(resk)
1691  do 10 j=1,15
1692  jtw = j*2
1693  dabsc = hlgth*xgk(jtw)
1694  fval1 = f(centr-dabsc)
1695  fval2 = f(centr+dabsc)
1696  fv1(jtw) = fval1
1697  fv2(jtw) = fval2
1698  fsum = fval1+fval2
1699  resg = resg+wg(j)*fsum
1700  resk = resk+wgk(jtw)*fsum
1701  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1702  10 continue
1703  do 15 j=1,15
1704  jtwm1 = j*2-1
1705  dabsc = hlgth*xgk(jtwm1)
1706  fval1 = f(centr-dabsc)
1707  fval2 = f(centr+dabsc)
1708  fv1(jtwm1) = fval1
1709  fv2(jtwm1) = fval2
1710  fsum = fval1+fval2
1711  resk = resk+wgk(jtwm1)*fsum
1712  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1713  15 continue
1714  reskh = resk*0.5d+00
1715  resasc = wgk(31)*dabs(fc-reskh)
1716  do 20 j=1,30
1717  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1718  20 continue
1719  result = resk*hlgth
1720  resabs = resabs*dhlgth
1721  resasc = resasc*dhlgth
1722  abserr = dabs((resk-resg)*hlgth)
1723  if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
1724  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1725  if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1726  * ((epmach*0.5d+02)*resabs,abserr)
1727  return
1728  end
1729  subroutine dqpsrt(limit,last,maxerr,ermax,elist,iord,nrmax)
1730 c***begin prologue dqpsrt
1731 c***refer to dqage,dqagie,dqagpe,dqawse
1732 c***routines called (none)
1733 c***revision date 810101 (yymmdd)
1734 c***keywords sequential sorting
1735 c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1736 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
1737 c***purpose this routine maintains the descending ordering in the
1738 c list of the local error estimated resulting from the
1739 c interval subdivision process. at each call two error
1740 c estimates are inserted using the sequential search
1741 c method, top-down for the largest error estimate and
1742 c bottom-up for the smallest error estimate.
1743 c***description
1744 c
1745 c ordering routine
1746 c standard fortran subroutine
1747 c double precision version
1748 c
1749 c parameters (meaning at output)
1750 c limit - integer
1751 c maximum number of error estimates the list
1752 c can contain
1753 c
1754 c last - integer
1755 c number of error estimates currently in the list
1756 c
1757 c maxerr - integer
1758 c maxerr points to the nrmax-th largest error
1759 c estimate currently in the list
1760 c
1761 c ermax - double precision
1762 c nrmax-th largest error estimate
1763 c ermax = elist(maxerr)
1764 c
1765 c elist - double precision
1766 c vector of dimension last containing
1767 c the error estimates
1768 c
1769 c iord - integer
1770 c vector of dimension last, the first k elements
1771 c of which contain pointers to the error
1772 c estimates, such that
1773 c elist(iord(1)),..., elist(iord(k))
1774 c form a decreasing sequence, with
1775 c k = last if last.le.(limit/2+2), and
1776 c k = limit+1-last otherwise
1777 c
1778 c nrmax - integer
1779 c maxerr = iord(nrmax)
1780 c
1781 c***end prologue dqpsrt
1782 c
1783  double precision elist,ermax,errmax,errmin
1784  integer i,ibeg,ido,iord,isucc,j,jbnd,jupbn,k,last,limit,maxerr,
1785  * nrmax
1786  dimension elist(last),iord(last)
1787 c
1788 c check whether the list contains more than
1789 c two error estimates.
1790 c
1791 c***first executable statement dqpsrt
1792  if(last.gt.2) go to 10
1793  iord(1) = 1
1794  iord(2) = 2
1795  go to 90
1796 c
1797 c this part of the routine is only executed if, due to a
1798 c difficult integrand, subdivision increased the error
1799 c estimate. in the normal case the insert procedure should
1800 c start after the nrmax-th largest error estimate.
1801 c
1802  10 errmax = elist(maxerr)
1803  if(nrmax.eq.1) go to 30
1804  ido = nrmax-1
1805  do 20 i = 1,ido
1806  isucc = iord(nrmax-1)
1807 c ***jump out of do-loop
1808  if(errmax.le.elist(isucc)) go to 30
1809  iord(nrmax) = isucc
1810  nrmax = nrmax-1
1811  20 continue
1812 c
1813 c compute the number of elements in the list to be maintained
1814 c in descending order. this number depends on the number of
1815 c subdivisions still allowed.
1816 c
1817  30 jupbn = last
1818  if(last.gt.(limit/2+2)) jupbn = limit+3-last
1819  errmin = elist(last)
1820 c
1821 c insert errmax by traversing the list top-down,
1822 c starting comparison from the element elist(iord(nrmax+1)).
1823 c
1824  jbnd = jupbn-1
1825  ibeg = nrmax+1
1826  if(ibeg.gt.jbnd) go to 50
1827  do 40 i=ibeg,jbnd
1828  isucc = iord(i)
1829 c ***jump out of do-loop
1830  if(errmax.ge.elist(isucc)) go to 60
1831  iord(i-1) = isucc
1832  40 continue
1833  50 iord(jbnd) = maxerr
1834  iord(jupbn) = last
1835  go to 90
1836 c
1837 c insert errmin by traversing the list bottom-up.
1838 c
1839  60 iord(i-1) = maxerr
1840  k = jbnd
1841  do 70 j=i,jbnd
1842  isucc = iord(k)
1843 c ***jump out of do-loop
1844  if(errmin.lt.elist(isucc)) go to 80
1845  iord(k+1) = isucc
1846  k = k-1
1847  70 continue
1848  iord(i) = last
1849  go to 90
1850  80 iord(k+1) = last
1851 c
1852 c set maxerr and ermax.
1853 c
1854  90 maxerr = iord(nrmax)
1855  ermax = elist(maxerr)
1856  return
1857  end
1858 
1859  DOUBLE PRECISION FUNCTION d1mach(I)
1860  INTEGER I
1861 C
1862 C DOUBLE-PRECISION MACHINE CONSTANTS
1863 C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
1864 C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
1865 C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
1866 C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
1867 C D1MACH( 5) = LOG10(B)
1868 C
1869  INTEGER SMALL(2)
1870  INTEGER LARGE(2)
1871  INTEGER RIGHT(2)
1872  INTEGER DIVER(2)
1873  INTEGER LOG10(2)
1874  INTEGER SC, CRAY1(38), J
1875  COMMON /d9mach/ cray1
1876  SAVE small, large, right, diver, log10, sc
1877  DOUBLE PRECISION DMACH(5)
1878  equivalence(dmach(1),small(1))
1879  equivalence(dmach(2),large(1))
1880  equivalence(dmach(3),right(1))
1881  equivalence(dmach(4),diver(1))
1882  equivalence(dmach(5),log10(1))
1883 C THIS VERSION ADAPTS AUTOMATICALLY TO MOST CURRENT MACHINES.
1884 C R1MACH CAN HANDLE AUTO-DOUBLE COMPILING, BUT THIS VERSION OF
1885 C D1MACH DOES NOT, BECAUSE WE DO NOT HAVE QUAD CONSTANTS FOR
1886 C MANY MACHINES YET.
1887 C TO COMPILE ON OLDER MACHINES, ADD A C IN COLUMN 1
1888 C ON THE NEXT LINE
1889  DATA sc/0/
1890 C AND REMOVE THE C FROM COLUMN 1 IN ONE OF THE SECTIONS BELOW.
1891 C CONSTANTS FOR EVEN OLDER MACHINES CAN BE OBTAINED BY
1892 C mail netlib@research.bell-labs.com
1893 C send old1mach from blas
1894 C PLEASE SEND CORRECTIONS TO dmg OR ehg@bell-labs.com.
1895 C
1896 C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
1897 C DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 /
1898 C DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 /
1899 C DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 /
1900 C DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 /
1901 C DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /, SC/987/
1902 C
1903 C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
1904 C 32-BIT INTEGERS.
1905 C DATA SMALL(1),SMALL(2) / 8388608, 0 /
1906 C DATA LARGE(1),LARGE(2) / 2147483647, -1 /
1907 C DATA RIGHT(1),RIGHT(2) / 612368384, 0 /
1908 C DATA DIVER(1),DIVER(2) / 620756992, 0 /
1909 C DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /, SC/987/
1910 C
1911 C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
1912 C DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 /
1913 C DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 /
1914 C DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 /
1915 C DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 /
1916 C DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /, SC/987/
1917 C
1918 C ON FIRST CALL, IF NO DATA UNCOMMENTED, TEST MACHINE TYPES.
1919  IF (sc .NE. 987) THEN
1920  dmach(1) = 1.d13
1921  IF ( small(1) .EQ. 1117925532
1922  * .AND. small(2) .EQ. -448790528) THEN
1923 * *** IEEE BIG ENDIAN ***
1924  small(1) = 1048576
1925  small(2) = 0
1926  large(1) = 2146435071
1927  large(2) = -1
1928  right(1) = 1017118720
1929  right(2) = 0
1930  diver(1) = 1018167296
1931  diver(2) = 0
1932  log10(1) = 1070810131
1933  log10(2) = 1352628735
1934  ELSE IF ( small(2) .EQ. 1117925532
1935  * .AND. small(1) .EQ. -448790528) THEN
1936 * *** IEEE LITTLE ENDIAN ***
1937  small(2) = 1048576
1938  small(1) = 0
1939  large(2) = 2146435071
1940  large(1) = -1
1941  right(2) = 1017118720
1942  right(1) = 0
1943  diver(2) = 1018167296
1944  diver(1) = 0
1945  log10(2) = 1070810131
1946  log10(1) = 1352628735
1947  ELSE IF ( small(1) .EQ. -2065213935
1948  * .AND. small(2) .EQ. 10752) THEN
1949 * *** VAX WITH D_FLOATING ***
1950  small(1) = 128
1951  small(2) = 0
1952  large(1) = -32769
1953  large(2) = -1
1954  right(1) = 9344
1955  right(2) = 0
1956  diver(1) = 9472
1957  diver(2) = 0
1958  log10(1) = 546979738
1959  log10(2) = -805796613
1960  ELSE IF ( small(1) .EQ. 1267827943
1961  * .AND. small(2) .EQ. 704643072) THEN
1962 * *** IBM MAINFRAME ***
1963  small(1) = 1048576
1964  small(2) = 0
1965  large(1) = 2147483647
1966  large(2) = -1
1967  right(1) = 856686592
1968  right(2) = 0
1969  diver(1) = 873463808
1970  diver(2) = 0
1971  log10(1) = 1091781651
1972  log10(2) = 1352628735
1973  ELSE IF ( small(1) .EQ. 1120022684
1974  * .AND. small(2) .EQ. -448790528) THEN
1975 * *** CONVEX C-1 ***
1976  small(1) = 1048576
1977  small(2) = 0
1978  large(1) = 2147483647
1979  large(2) = -1
1980  right(1) = 1019215872
1981  right(2) = 0
1982  diver(1) = 1020264448
1983  diver(2) = 0
1984  log10(1) = 1072907283
1985  log10(2) = 1352628735
1986  ELSE IF ( small(1) .EQ. 815547074
1987  * .AND. small(2) .EQ. 58688) THEN
1988 * *** VAX G-FLOATING ***
1989  small(1) = 16
1990  small(2) = 0
1991  large(1) = -32769
1992  large(2) = -1
1993  right(1) = 15552
1994  right(2) = 0
1995  diver(1) = 15568
1996  diver(2) = 0
1997  log10(1) = 1142112243
1998  log10(2) = 2046775455
1999  ELSE
2000  dmach(2) = 1.d27 + 1
2001  dmach(3) = 1.d27
2002  large(2) = large(2) - right(2)
2003  IF (large(2) .EQ. 64 .AND. small(2) .EQ. 0) THEN
2004  cray1(1) = 67291416
2005  DO 10 j = 1, 20
2006  cray1(j+1) = cray1(j) + cray1(j)
2007  10 CONTINUE
2008  cray1(22) = cray1(21) + 321322
2009  DO 20 j = 22, 37
2010  cray1(j+1) = cray1(j) + cray1(j)
2011  20 CONTINUE
2012  IF (cray1(38) .EQ. small(1)) THEN
2013 * *** CRAY ***
2014  CALL i1mcry(small(1), j, 8285, 8388608, 0)
2015  small(2) = 0
2016  CALL i1mcry(large(1), j, 24574, 16777215, 16777215)
2017  CALL i1mcry(large(2), j, 0, 16777215, 16777214)
2018  CALL i1mcry(right(1), j, 16291, 8388608, 0)
2019  right(2) = 0
2020  CALL i1mcry(diver(1), j, 16292, 8388608, 0)
2021  diver(2) = 0
2022  CALL i1mcry(log10(1), j, 16383, 10100890, 8715215)
2023  CALL i1mcry(log10(2), j, 0, 16226447, 9001388)
2024  ELSE
2025  WRITE(*,9000)
2026  stop 779
2027  END IF
2028  ELSE
2029  WRITE(*,9000)
2030  stop 779
2031  END IF
2032  END IF
2033  sc = 987
2034  END IF
2035 * SANITY CHECK
2036  IF (dmach(4) .GE. 1.0d0) stop 778
2037  IF (i .LT. 1 .OR. i .GT. 5) THEN
2038  WRITE(*,*) 'D1MACH(I): I =',i,' is out of bounds.'
2039  stop
2040  END IF
2041  d1mach = dmach(i)
2042  RETURN
2043  9000 FORMAT(/' Adjust D1MACH by uncommenting data statements'/
2044  *' appropriate for your machine.')
2045  END FUNCTION d1mach
2046 
2047  SUBROUTINE i1mcry(A, A1, B, C, D)
2048 **** SPECIAL COMPUTATION FOR OLD CRAY MACHINES ****
2049  INTEGER A, A1, B, C, D
2050  A1 = 16777216*b + c
2051  a = 16777216*a1 + d
2052  END SUBROUTINE
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17