module quadpack use const_and_precisions, only : wp_ implicit none contains subroutine dqags(f,a,b,epsabs,epsrel,result,abserr,neval,ier, & limit,lenw,last,iwork,work) !***begin prologue dqags !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a1a1 !***keywords automatic integrator, general-purpose, ! (end-point) singularities, extrapolation, ! globally adaptive !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & prog. div. - k.u.leuven !***purpose the routine calculates an approximation result to a given ! definite integral i = integral of f over (a,b), ! hopefully satisfying following claim for accuracy ! abs(i-result)<=max(epsabs,epsrel*abs(i)). !***description ! ! computation of a definite integral ! standard fortran subroutine ! real(8) version ! ! ! parameters ! on entry ! f - real(8) ! function subprogram defining the integrand ! function f(x). the actual name for f needs to be ! declared e x t e r n a l in the driver program. ! ! a - real(8) ! lower limit of integration ! ! b - real(8) ! upper limit of integration ! ! epsabs - real(8) ! absolute accuracy requested ! epsrel - real(8) ! relative accuracy requested ! if epsabs<=0 ! and epsrel0 abnormal termination of the routine ! the estimates for integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! error messages ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more sub- ! divisions by increasing the value of limit ! (and taking the according dimension ! adjustments into account. however, if ! this yields no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulties. if ! the position of a local difficulty can be ! determined (e.g. singularity, ! discontinuity within the interval) one ! will probably gain from splitting up the ! interval at this point and calling the ! integrator on the subranges. if possible, ! an appropriate special-purpose integrator ! should be used, which is designed for ! handling the type of difficulty involved. ! = 2 the occurrence of roundoff error is detec- ! ted, which prevents the requested ! tolerance from being achieved. ! the error may be under-estimated. ! = 3 extremely bad integrand behaviour ! occurs at some points of the integration ! interval. ! = 4 the algorithm does not converge. ! roundoff error is detected in the ! extrapolation table. it is presumed that ! the requested tolerance cannot be ! achieved, and that the returned result is ! the best which can be obtained. ! = 5 the integral is probably divergent, or ! slowly convergent. it must be noted that ! divergence can occur with any other value ! of ier. ! = 6 the input is invalid, because ! (epsabs<=0 and ! epsrel=1. ! if limit<1, the routine will end with ier = 6. ! ! lenw - integer ! dimensioning parameter for work ! lenw must be at least limit*4. ! if lenw=1.and.lenw>=limit*4) then ! ! prepare call for dqagse. ! l1 = limit+1 l2 = limit+l1 l3 = limit+l2 ! call dqagse(f,a,b,epsabs,epsrel,limit,result,abserr,neval, & ier,work(1),work(l1),work(l2),work(l3),iwork,last) ! ! call error handler if necessary. ! lvl = 0 end if if(ier==6) lvl = 1 if(ier/=0) print*,'habnormal return from dqags',ier,lvl end subroutine dqags subroutine dqagse(f,a,b,epsabs,epsrel,limit,result,abserr,neval, & ier,alist,blist,rlist,elist,iord,last) !***begin prologue dqagse !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a1a1 !***keywords automatic integrator, general-purpose, ! (end point) singularities, extrapolation, ! globally adaptive !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & progr. div. - k.u.leuven !***purpose the routine calculates an approximation result to a given ! definite integral i = integral of f over (a,b), ! hopefully satisfying following claim for accuracy ! abs(i-result)<=max(epsabs,epsrel*abs(i)). !***description ! ! computation of a definite integral ! standard fortran subroutine ! real(8) version ! ! parameters ! on entry ! f - real(8) ! function subprogram defining the integrand ! function f(x). the actual name for f needs to be ! declared e x t e r n a l in the driver program. ! ! a - real(8) ! lower limit of integration ! ! b - real(8) ! upper limit of integration ! ! epsabs - real(8) ! absolute accuracy requested ! epsrel - real(8) ! relative accuracy requested ! if epsabs<=0 ! and epsrel0 abnormal termination of the routine ! the estimates for integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! error messages ! = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more sub- ! divisions by increasing the value of limit ! (and taking the according dimension ! adjustments into account). however, if ! this yields no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulties. if ! the position of a local difficulty can be ! determined (e.g. singularity, ! discontinuity within the interval) one ! will probably gain from splitting up the ! interval at this point and calling the ! integrator on the subranges. if possible, ! an appropriate special-purpose integrator ! should be used, which is designed for ! handling the type of difficulty involved. ! = 2 the occurrence of roundoff error is detec- ! ted, which prevents the requested ! tolerance from being achieved. ! the error may be under-estimated. ! = 3 extremely bad integrand behaviour ! occurs at some points of the integration ! interval. ! = 4 the algorithm does not converge. ! roundoff error is detected in the ! extrapolation table. ! it is presumed that the requested ! tolerance cannot be achieved, and that the ! returned result is the best which can be ! obtained. ! = 5 the integral is probably divergent, or ! slowly convergent. it must be noted that ! divergence can occur with any other value ! of ier. ! = 6 the input is invalid, because ! epsabs<=0 and ! epsrelcomp_eps, uflow=>comp_tiny, & oflow=>comp_huge implicit none real(wp_), intent(in) :: a,b,epsabs,epsrel integer, intent(in) :: limit real(wp_), intent(out) :: result,abserr integer, intent(out) :: neval,ier,last real(wp_), dimension(limit), intent(inout) :: alist,blist,elist,rlist integer, dimension(limit), intent(inout) :: iord real(wp_), external :: f ! real(wp_) :: abseps,area,area1,area12,area2,a1,a2,b1,b2,correc,abs, & defabs,defab1,defab2,dmax1,dres,erlarg,erlast,errbnd,errmax, & error1,error2,erro12,errsum,ertest,resabs,reseps,small real(wp_) :: res3la(3),rlist2(52) integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, & ktmin,maxerr,nres,nrmax,numrl2 logical :: extrap,noext ! ! the dimension of rlist2 is determined by the value of ! limexp in subroutine dqelg (rlist2 should be of dimension ! (limexp+2) at least). ! ! list of major variables ! ----------------------- ! ! alist - list of left end points of all subintervals ! considered up to now ! blist - list of right end points of all subintervals ! considered up to now ! rlist(i) - approximation to the integral over ! (alist(i),blist(i)) ! rlist2 - array of dimension at least limexp+2 containing ! the part of the epsilon table which is still ! needed for further computations ! elist(i) - error estimate applying to rlist(i) ! maxerr - pointer to the interval with largest error ! estimate ! errmax - elist(maxerr) ! erlast - error on the interval currently subdivided ! (before that subdivision has taken place) ! area - sum of the integrals over the subintervals ! errsum - sum of the errors over the subintervals ! errbnd - requested accuracy max(epsabs,epsrel* ! abs(result)) ! *****1 - variable for the left interval ! *****2 - variable for the right interval ! last - index for subdivision ! nres - number of calls to the extrapolation routine ! numrl2 - number of elements currently in rlist2. if an ! appropriate approximation to the compounded ! integral has been obtained it is put in ! rlist2(numrl2) after numrl2 has been increased ! by one. ! small - length of the smallest interval considered up ! to now, multiplied by 1.5 ! erlarg - sum of the errors over the intervals larger ! than the smallest interval considered up to now ! extrap - logical variable denoting that the routine is ! attempting to perform extrapolation i.e. before ! subdividing the smallest interval we try to ! decrease the value of erlarg. ! noext - logical variable denoting that extrapolation ! is no longer allowed (true value) ! ! machine dependent constants ! --------------------------- ! ! epmach is the largest relative spacing. ! uflow is the smallest positive magnitude. ! oflow is the largest positive magnitude. ! !***first executable statement dqagse ! ! test on validity of parameters ! ------------------------------ ier = 0 neval = 0 last = 0 result = 0.0e+00_wp_ abserr = 0.0e+00_wp_ alist(1) = a blist(1) = b rlist(1) = 0.0e+00_wp_ elist(1) = 0.0e+00_wp_ if(epsabs<=0.0e+00_wp_.and.epsrelerrbnd) ier = 2 if(limit==1) ier = 1 if(ier/=0.or.(abserr<=errbnd.and.abserr/=resabs).or. & abserr==0.0e+00_wp_) go to 140 ! ! initialization ! -------------- ! rlist2(1) = result errmax = abserr maxerr = 1 area = result errsum = abserr abserr = oflow nrmax = 1 nres = 0 numrl2 = 2 ktmin = 0 extrap = .false. noext = .false. iroff1 = 0 iroff2 = 0 iroff3 = 0 ksgn = -1 if(dres>=(0.1e+01_wp_-0.5e+02_wp_*epmach)*defabs) ksgn = 1 ! ! main do-loop ! ------------ ! do last = 2,limit ! ! bisect the subinterval with the nrmax-th largest error ! estimate. ! a1 = alist(maxerr) b1 = 0.5e+00_wp_*(alist(maxerr)+blist(maxerr)) a2 = b1 b2 = blist(maxerr) erlast = errmax call dqk21(f,a1,b1,area1,error1,resabs,defab1) call dqk21(f,a2,b2,area2,error2,resabs,defab2) ! ! improve previous approximations to integral ! and error and test for accuracy. ! area12 = area1+area2 erro12 = error1+error2 errsum = errsum+erro12-errmax area = area+area12-rlist(maxerr) if(defab1/=error1.and.defab2/=error2) then if(abs(rlist(maxerr)-area12)<=0.1e-04_wp_*abs(area12) & .and.erro12>=0.99e+00_wp_*errmax) then if(extrap) iroff2 = iroff2+1 if(.not.extrap) iroff1 = iroff1+1 end if if(last>10.and.erro12>errmax) iroff3 = iroff3+1 end if rlist(maxerr) = area1 rlist(last) = area2 errbnd = dmax1(epsabs,epsrel*abs(area)) ! ! test for roundoff error and eventually set error flag. ! if(iroff1+iroff2>=10.or.iroff3>=20) ier = 2 if(iroff2>=5) ierro = 3 ! ! set error flag in the case that the number of subintervals ! equals limit. ! if(last==limit) ier = 1 ! ! set error flag in the case of bad integrand behaviour ! at a point of the integration range. ! if(dmax1(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* & (abs(a2)+0.1e+04_wp_*uflow)) ier = 4 ! ! append the newly-created intervals to the list. ! if(error2<=error1) then alist(last) = a2 blist(maxerr) = b1 blist(last) = b2 elist(maxerr) = error1 elist(last) = error2 else alist(maxerr) = a2 alist(last) = a1 blist(last) = b1 rlist(maxerr) = area2 rlist(last) = area1 elist(maxerr) = error2 elist(last) = error1 end if ! ! call subroutine dqpsrt to maintain the descending ordering ! in the list of error estimates and select the subinterval ! with nrmax-th largest error estimate (to be bisected next). ! call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) ! ***jump out of do-loop if(errsum<=errbnd) go to 115 ! ***jump out of do-loop if(ier/=0) exit if(last==2) then small = abs(b-a)*0.375e+00_wp_ erlarg = errsum ertest = errbnd rlist2(2) = area cycle end if if(noext) cycle erlarg = erlarg-erlast if(abs(b1-a1)>small) erlarg = erlarg+erro12 if(.not.extrap) then ! ! test whether the interval to be bisected next is the ! smallest interval. ! if(abs(blist(maxerr)-alist(maxerr))>small) cycle extrap = .true. nrmax = 2 end if if(ierro/=3.and.erlarg>ertest) then ! ! the smallest interval has the largest error. ! before bisecting decrease the sum of the errors over the ! larger intervals (erlarg) and perform extrapolation. ! id = nrmax jupbnd = last if(last>(2+limit/2)) jupbnd = limit+3-last do k = id,jupbnd maxerr = iord(nrmax) errmax = elist(maxerr) ! ***jump out of do-loop if(abs(blist(maxerr)-alist(maxerr))>small) go to 90 nrmax = nrmax+1 end do ! ! perform extrapolation. ! end if numrl2 = numrl2+1 rlist2(numrl2) = area call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres) ktmin = ktmin+1 if(ktmin>5.and.abserr<0.1e-02_wp_*errsum) ier = 5 if(absepserrsum) go to 115 if(area==0.0e+00_wp_) go to 130 go to 110 105 continue if(abserr/abs(result)>errsum/abs(area)) go to 115 ! ! test on divergence. ! 110 continue if(ksgn==(-1).and.dmax1(abs(result),abs(area))<= & defabs*0.1e-01_wp_) go to 130 if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ & .or.errsum>abs(area)) ier = 6 go to 130 ! ! compute global integral sum. ! 115 continue result = 0.0e+00_wp_ do k = 1,last result = result+rlist(k) end do abserr = errsum 130 continue if(ier>2) ier = ier-1 140 continue neval = 42*last-21 end subroutine dqagse subroutine dqelg(n,epstab,result,abserr,res3la,nres) !***begin prologue dqelg !***refer to dqagie,dqagoe,dqagpe,dqagse !***routines called (none) !***revision date 830518 (yymmdd) !***keywords epsilon algorithm, convergence acceleration, ! extrapolation !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math & progr. div. - k.u.leuven !***purpose the routine determines the limit of a given sequence of ! approximations, by means of the epsilon algorithm of ! p.wynn. an estimate of the absolute error is also given. ! the condensed epsilon table is computed. only those ! elements needed for the computation of the next diagonal ! are preserved. !***description ! ! epsilon algorithm ! standard fortran subroutine ! real(8) version ! ! parameters ! n - integer ! epstab(n) contains the new element in the ! first column of the epsilon table. ! ! epstab - real(8) ! vector of dimension 52 containing the elements ! of the two lower diagonals of the triangular ! epsilon table. the elements are numbered ! starting at the right-hand corner of the ! triangle. ! ! result - real(8) ! resulting approximation to the integral ! ! abserr - real(8) ! estimate of the absolute error computed from ! result and the 3 previous results ! ! res3la - real(8) ! vector of dimension 3 containing the last 3 ! results ! ! nres - integer ! number of calls to the routine ! (should be zero at first call) ! !***end prologue dqelg ! use const_and_precisions, only : epmach=>comp_eps, oflow=>comp_huge implicit none real(wp_), intent(out) :: abserr,result real(wp_), dimension(52), intent(inout) :: epstab real(wp_), dimension(3), intent(inout) :: res3la integer, intent(inout) :: n,nres real(wp_) :: abs,delta1,delta2,delta3,dmax1,epsinf,error, & err1,err2,err3,e0,e1,e1abs,e2,e3,res,ss,tol1,tol2,tol3 integer i,ib,ib2,ie,indx,k1,k2,k3,limexp,newelm,num ! ! list of major variables ! ----------------------- ! ! e0 - the 4 elements on which the computation of a new ! e1 element in the epsilon table is based ! e2 ! e3 e0 ! e3 e1 new ! e2 ! newelm - number of elements to be computed in the new ! diagonal ! error - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2) ! result - the element in the new diagonal with least value ! of error ! ! machine dependent constants ! --------------------------- ! ! epmach is the largest relative spacing. ! oflow is the largest positive magnitude. ! limexp is the maximum number of elements the epsilon ! table can contain. if this number is reached, the upper ! diagonal of the epsilon table is deleted. ! !***first executable statement dqelg nres = nres+1 abserr = oflow result = epstab(n) if(n<3) go to 100 limexp = 50 epstab(n+2) = epstab(n) newelm = (n-1)/2 epstab(n) = oflow num = n k1 = n do i = 1,newelm k2 = k1-1 k3 = k1-2 res = epstab(k1+2) e0 = epstab(k3) e1 = epstab(k2) e2 = res e1abs = abs(e1) delta2 = e2-e1 err2 = abs(delta2) tol2 = dmax1(abs(e2),e1abs)*epmach delta3 = e1-e0 err3 = abs(delta3) tol3 = dmax1(e1abs,abs(e0))*epmach if(err2<=tol2.and.err3<=tol3) then ! ! if e0, e1 and e2 are equal to within machine ! accuracy, convergence is assumed. ! result = e2 ! abserr = abs(e1-e0)+abs(e2-e1) ! result = res abserr = err2+err3 ! ***jump out of do-loop go to 100 end if e3 = epstab(k1) epstab(k1) = e1 delta1 = e1-e3 err1 = abs(delta1) tol1 = dmax1(e1abs,abs(e3))*epmach ! ! if two elements are very close to each other, omit ! a part of the table by adjusting the value of n ! if(err1<=tol1.or.err2<=tol2.or.err3<=tol3) go to 20 ss = 0.1e+01_wp_/delta1+0.1e+01_wp_/delta2-0.1e+01_wp_/delta3 epsinf = abs(ss*e1) ! ! test to detect irregular behaviour in the table, and ! eventually omit a part of the table adjusting the value ! of n. ! if(epsinf>0.1e-03_wp_) go to 30 ! ***jump out of do-loop 20 continue n = i+i-1 exit ! ! compute a new element and eventually adjust ! the value of result. ! 30 continue res = e1+0.1e+01_wp_/ss epstab(k1) = res k1 = k1-2 error = err2+abs(res-e2)+err3 if(error<=abserr) then abserr = error result = res end if end do ! ! shift the table. ! if(n==limexp) n = 2*(limexp/2)-1 ib = 1 if((num/2)*2==num) ib = 2 ie = newelm+1 do i=1,ie ib2 = ib+2 epstab(ib) = epstab(ib2) ib = ib2 end do if(num/=n) then indx = num-n+1 do i = 1,n epstab(i)= epstab(indx) indx = indx+1 end do end if if(nres<4) then res3la(nres) = result abserr = oflow else ! ! compute error estimate ! abserr = abs(result-res3la(3))+abs(result-res3la(2)) & +abs(result-res3la(1)) res3la(1) = res3la(2) res3la(2) = res3la(3) res3la(3) = result end if 100 continue abserr = dmax1(abserr,0.5e+01_wp_*epmach*abs(result)) end subroutine dqelg subroutine dqk21(f,a,b,result,abserr,resabs,resasc) !***begin prologue dqk21 !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a1a2 !***keywords 21-point gauss-kronrod rules !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & progr. div. - k.u.leuven !***purpose to compute i = integral of f over (a,b), with error ! estimate ! j = integral of abs(f) over (a,b) !***description ! ! integration rules ! standard fortran subroutine ! real(8) version ! ! parameters ! on entry ! f - real(8) ! function subprogram defining the integrand ! function f(x). the actual name for f needs to be ! declared e x t e r n a l in the driver program. ! ! a - real(8) ! lower limit of integration ! ! b - real(8) ! upper limit of integration ! ! on return ! result - real(8) ! approximation to the integral i ! result is computed by applying the 21-point ! kronrod rule (resk) obtained by optimal addition ! of abscissae to the 10-point gauss rule (resg). ! ! abserr - real(8) ! estimate of the modulus of the absolute error, ! which should not exceed abs(i-result) ! ! resabs - real(8) ! approximation to the integral j ! ! resasc - real(8) ! approximation to the integral of abs(f-i/(b-a)) ! over (a,b) ! !***references (none) !***routines called (none) !***end prologue dqk21 ! use const_and_precisions, only : epmach=>comp_eps, uflow=>comp_tiny implicit none real(wp_), intent(in) :: a,b real(wp_), intent(out) :: result,abserr,resabs,resasc real(wp_), external :: f real(wp_) :: absc,centr,abs,dhlgth,dmax1,dmin1,fc,fsum, & fval1,fval2,hlgth,resg,resk,reskh real(wp_), dimension(10) :: fv1,fv2 integer :: j,jtw,jtwm1 ! ! the abscissae and weights are given for the interval (-1,1). ! because of symmetry only the positive abscissae and their ! corresponding weights are given. ! ! xgk - abscissae of the 21-point kronrod rule ! xgk(2), xgk(4), ... abscissae of the 10-point ! gauss rule ! xgk(1), xgk(3), ... abscissae which are optimally ! added to the 10-point gauss rule ! ! wgk - weights of the 21-point kronrod rule ! ! wg - weights of the 10-point gauss rule ! ! ! gauss quadrature weights and kronron quadrature abscissae and weights ! as evaluated with 80 decimal digit arithmetic by l. w. fullerton, ! bell labs, nov. 1981. ! real(wp_), dimension(5), parameter :: & wg = (/ 0.066671344308688137593568809893332_wp_, & 0.149451349150580593145776339657697_wp_, & 0.219086362515982043995534934228163_wp_, & 0.269266719309996355091226921569469_wp_, & 0.295524224714752870173892994651338_wp_ /) ! real(wp_), dimension(11), parameter :: & xgk = (/ 0.995657163025808080735527280689003_wp_, & 0.973906528517171720077964012084452_wp_, & 0.930157491355708226001207180059508_wp_, & 0.865063366688984510732096688423493_wp_, & 0.780817726586416897063717578345042_wp_, & 0.679409568299024406234327365114874_wp_, & 0.562757134668604683339000099272694_wp_, & 0.433395394129247190799265943165784_wp_, & 0.294392862701460198131126603103866_wp_, & 0.148874338981631210884826001129720_wp_, & 0.000000000000000000000000000000000_wp_ /), & wgk = (/ 0.011694638867371874278064396062192_wp_, & 0.032558162307964727478818972459390_wp_, & 0.054755896574351996031381300244580_wp_, & 0.075039674810919952767043140916190_wp_, & 0.093125454583697605535065465083366_wp_, & 0.109387158802297641899210590325805_wp_, & 0.123491976262065851077958109831074_wp_, & 0.134709217311473325928054001771707_wp_, & 0.142775938577060080797094273138717_wp_, & 0.147739104901338491374841515972068_wp_, & 0.149445554002916905664936468389821_wp_ /) ! ! ! list of major variables ! ----------------------- ! ! centr - mid point of the interval ! hlgth - half-length of the interval ! absc - abscissa ! fval* - function value ! resg - result of the 10-point gauss formula ! resk - result of the 21-point kronrod formula ! reskh - approximation to the mean value of f over (a,b), ! i.e. to i/(b-a) ! ! ! machine dependent constants ! --------------------------- ! ! epmach is the largest relative spacing. ! uflow is the smallest positive magnitude. ! !***first executable statement dqk21 centr = 0.5e+00_wp_*(a+b) hlgth = 0.5e+00_wp_*(b-a) dhlgth = abs(hlgth) ! ! compute the 21-point kronrod approximation to ! the integral, and estimate the absolute error. ! resg = 0.0e+00_wp_ fc = f(centr) resk = wgk(11)*fc resabs = abs(resk) do j=1,5 jtw = 2*j absc = hlgth*xgk(jtw) fval1 = f(centr-absc) fval2 = f(centr+absc) fv1(jtw) = fval1 fv2(jtw) = fval2 fsum = fval1+fval2 resg = resg+wg(j)*fsum resk = resk+wgk(jtw)*fsum resabs = resabs+wgk(jtw)*(abs(fval1)+abs(fval2)) end do do j = 1,5 jtwm1 = 2*j-1 absc = hlgth*xgk(jtwm1) fval1 = f(centr-absc) fval2 = f(centr+absc) fv1(jtwm1) = fval1 fv2(jtwm1) = fval2 fsum = fval1+fval2 resk = resk+wgk(jtwm1)*fsum resabs = resabs+wgk(jtwm1)*(abs(fval1)+abs(fval2)) end do reskh = resk*0.5e+00_wp_ resasc = wgk(11)*abs(fc-reskh) do j=1,10 resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)) end do result = resk*hlgth resabs = resabs*dhlgth resasc = resasc*dhlgth abserr = abs((resk-resg)*hlgth) if(resasc/=0.0e+00_wp_.and.abserr/=0.0e+00_wp_) & abserr = resasc*dmin1(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_) if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 & ((epmach*0.5e+02_wp_)*resabs,abserr) end subroutine dqk21 subroutine dqpsrt(limit,last,maxerr,ermax,elist,iord,nrmax) !***begin prologue dqpsrt !***refer to dqage,dqagie,dqagpe,dqawse !***routines called (none) !***revision date 810101 (yymmdd) !***keywords sequential sorting !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & progr. div. - k.u.leuven !***purpose this routine maintains the descending ordering in the ! list of the local error estimated resulting from the ! interval subdivision process. at each call two error ! estimates are inserted using the sequential search ! method, top-down for the largest error estimate and ! bottom-up for the smallest error estimate. !***description ! ! ordering routine ! standard fortran subroutine ! real(8) version ! ! parameters (meaning at output) ! limit - integer ! maximum number of error estimates the list ! can contain ! ! last - integer ! number of error estimates currently in the list ! ! maxerr - integer ! maxerr points to the nrmax-th largest error ! estimate currently in the list ! ! ermax - real(8) ! nrmax-th largest error estimate ! ermax = elist(maxerr) ! ! elist - real(8) ! vector of dimension last containing ! the error estimates ! ! iord - integer ! vector of dimension last, the first k elements ! of which contain pointers to the error ! estimates, such that ! elist(iord(1)),..., elist(iord(k)) ! form a decreasing sequence, with ! k = last if last<=(limit/2+2), and ! k = limit+1-last otherwise ! ! nrmax - integer ! maxerr = iord(nrmax) ! !***end prologue dqpsrt ! implicit none integer, intent(in) :: last,limit real(wp_), intent(out) :: ermax integer, intent(inout) :: maxerr,nrmax real(wp_), dimension(last), intent(inout) :: elist integer, dimension(last), intent(inout) :: iord real(wp_) :: errmax,errmin integer :: i,ibeg,ido,isucc,j,jbnd,jupbn,k ! ! check whether the list contains more than ! two error estimates. ! !***first executable statement dqpsrt if(last<=2) then iord(1) = 1 iord(2) = 2 go to 90 end if ! ! this part of the routine is only executed if, due to a ! difficult integrand, subdivision increased the error ! estimate. in the normal case the insert procedure should ! start after the nrmax-th largest error estimate. ! errmax = elist(maxerr) if(nrmax/=1) then ido = nrmax-1 do i = 1,ido isucc = iord(nrmax-1) ! ***jump out of do-loop if(errmax<=elist(isucc)) exit iord(nrmax) = isucc nrmax = nrmax-1 end do end if ! ! compute the number of elements in the list to be maintained ! in descending order. this number depends on the number of ! subdivisions still allowed. ! jupbn = last if(last>(limit/2+2)) jupbn = limit+3-last errmin = elist(last) ! ! insert errmax by traversing the list top-down, ! starting comparison from the element elist(iord(nrmax+1)). ! jbnd = jupbn-1 ibeg = nrmax+1 do i=ibeg,jbnd isucc = iord(i) ! ***jump out of do-loop if(errmax>=elist(isucc)) then ! ! insert errmin by traversing the list bottom-up. ! iord(i-1) = maxerr k = jbnd do j=i,jbnd isucc = iord(k) ! ***jump out of do-loop if(errmin0 abnormal termination of the routine. the ! estimates for result and error are less ! reliable. it is assumed that the requested ! accuracy has not been achieved. ! error messages ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the value of ! limit (and taking the according dimension ! adjustments into account). however, if ! this yields no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulties. if ! the position of a local difficulty can be ! determined (e.g. singularity, ! discontinuity within the interval) one ! will probably gain from splitting up the ! interval at this point and calling the ! integrator on the subranges. if possible, ! an appropriate special-purpose integrator ! should be used, which is designed for ! handling the type of difficulty involved. ! = 2 the occurrence of roundoff error is ! detected, which prevents the requested ! tolerance from being achieved. ! the error may be under-estimated. ! = 3 extremely bad integrand behaviour occurs ! at some points of the integration ! interval. ! = 4 the algorithm does not converge. ! roundoff error is detected in the ! extrapolation table. ! it is assumed that the requested tolerance ! cannot be achieved, and that the returned ! result is the best which can be obtained. ! = 5 the integral is probably divergent, or ! slowly convergent. it must be noted that ! divergence can occur with any other value ! of ier. ! = 6 the input is invalid, because ! (epsabs<=0 and ! epsrel=1. ! if limit<1, the routine will end with ier = 6. ! ! lenw - integer ! dimensioning parameter for work ! lenw must be at least limit*4. ! if lenw=1.and.lenw>=limit*4) then ! ! prepare call for dqagie. ! l1 = limit+1 l2 = limit+l1 l3 = limit+l2 ! call dqagie(f,bound,inf,epsabs,epsrel,limit,result,abserr, & neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last) ! end if if(ier/=0) print*,'habnormal return from dqagi' end subroutine dqagi subroutine dqagie(f,bound,inf,epsabs,epsrel,limit,result,abserr, & neval,ier,alist,blist,rlist,elist,iord,last) !***begin prologue dqagie !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a3a1,h2a4a1 !***keywords automatic integrator, infinite intervals, ! general-purpose, transformation, extrapolation, ! globally adaptive !***author piessens,robert,appl. math & progr. div - k.u.leuven ! de doncker,elise,appl. math & progr. div - k.u.leuven !***purpose the routine calculates an approximation result to a given ! integral i = integral of f over (bound,+infinity) ! or i = integral of f over (-infinity,bound) ! or i = integral of f over (-infinity,+infinity), ! hopefully satisfying following claim for accuracy ! abs(i-result)<=max(epsabs,epsrel*abs(i)) !***description ! ! integration over infinite intervals ! standard fortran subroutine ! ! f - real(8) ! function subprogram defining the integrand ! function f(x). the actual name for f needs to be ! declared e x t e r n a l in the driver program. ! ! bound - real(8) ! finite bound of integration range ! (has no meaning if interval is doubly-infinite) ! ! inf - real(8) ! indicating the kind of integration range involved ! inf = 1 corresponds to (bound,+infinity), ! inf = -1 to (-infinity,bound), ! inf = 2 to (-infinity,+infinity). ! ! epsabs - real(8) ! absolute accuracy requested ! epsrel - real(8) ! relative accuracy requested ! if epsabs<=0 ! and epsrel=1 ! ! on return ! result - real(8) ! approximation to the integral ! ! abserr - real(8) ! estimate of the modulus of the absolute error, ! which should equal or exceed abs(i-result) ! ! neval - integer ! number of integrand evaluations ! ! ier - integer ! ier = 0 normal and reliable termination of the ! routine. it is assumed that the requested ! accuracy has been achieved. ! - ier>0 abnormal termination of the routine. the ! estimates for result and error are less ! reliable. it is assumed that the requested ! accuracy has not been achieved. ! error messages ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the value of ! limit (and taking the according dimension ! adjustments into account). however,if ! this yields no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulties. ! if the position of a local difficulty can ! be determined (e.g. singularity, ! discontinuity within the interval) one ! will probably gain from splitting up the ! interval at this point and calling the ! integrator on the subranges. if possible, ! an appropriate special-purpose integrator ! should be used, which is designed for ! handling the type of difficulty involved. ! = 2 the occurrence of roundoff error is ! detected, which prevents the requested ! tolerance from being achieved. ! the error may be under-estimated. ! = 3 extremely bad integrand behaviour occurs ! at some points of the integration ! interval. ! = 4 the algorithm does not converge. ! roundoff error is detected in the ! extrapolation table. ! it is assumed that the requested tolerance ! cannot be achieved, and that the returned ! result is the best which can be obtained. ! = 5 the integral is probably divergent, or ! slowly convergent. it must be noted that ! divergence can occur with any other value ! of ier. ! = 6 the input is invalid, because ! (epsabs<=0 and ! epsrelcomp_eps, uflow=>comp_tiny, & oflow=>comp_huge implicit none integer, intent(in) :: limit,inf real(wp_), intent(in) :: bound,epsabs,epsrel real(wp_), intent(out) :: result,abserr integer, intent(out) :: ier,neval,last real(wp_), dimension(limit), intent(inout) :: alist,blist,elist,rlist integer, dimension(limit), intent(inout) :: iord real(wp_), external :: f real(wp_) :: abseps,area,area1,area12,area2,a1,a2,boun,b1,b2,correc, & abs,defabs,defab1,defab2,dmax1,dres,erlarg,erlast,errbnd,errmax, & error1,error2,erro12,errsum,ertest,resabs,reseps,small real(wp_) :: res3la(3),rlist2(52) integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, & ktmin,maxerr,nres,nrmax,numrl2 logical :: extrap,noext ! ! the dimension of rlist2 is determined by the value of ! limexp in subroutine dqelg. ! ! ! list of major variables ! ----------------------- ! ! alist - list of left end points of all subintervals ! considered up to now ! blist - list of right end points of all subintervals ! considered up to now ! rlist(i) - approximation to the integral over ! (alist(i),blist(i)) ! rlist2 - array of dimension at least (limexp+2), ! containing the part of the epsilon table ! wich is still needed for further computations ! elist(i) - error estimate applying to rlist(i) ! maxerr - pointer to the interval with largest error ! estimate ! errmax - elist(maxerr) ! erlast - error on the interval currently subdivided ! (before that subdivision has taken place) ! area - sum of the integrals over the subintervals ! errsum - sum of the errors over the subintervals ! errbnd - requested accuracy max(epsabs,epsrel* ! abs(result)) ! *****1 - variable for the left subinterval ! *****2 - variable for the right subinterval ! last - index for subdivision ! nres - number of calls to the extrapolation routine ! numrl2 - number of elements currently in rlist2. if an ! appropriate approximation to the compounded ! integral has been obtained, it is put in ! rlist2(numrl2) after numrl2 has been increased ! by one. ! small - length of the smallest interval considered up ! to now, multiplied by 1.5 ! erlarg - sum of the errors over the intervals larger ! than the smallest interval considered up to now ! extrap - logical variable denoting that the routine ! is attempting to perform extrapolation. i.e. ! before subdividing the smallest interval we ! try to decrease the value of erlarg. ! noext - logical variable denoting that extrapolation ! is no longer allowed (true-value) ! ! machine dependent constants ! --------------------------- ! ! epmach is the largest relative spacing. ! uflow is the smallest positive magnitude. ! oflow is the largest positive magnitude. ! !***first executable statement dqagie ! ! test on validity of parameters ! ----------------------------- ! ier = 0 neval = 0 last = 0 result = 0.0e+00_wp_ abserr = 0.0e+00_wp_ alist(1) = 0.0e+00_wp_ blist(1) = 0.1e+01_wp_ rlist(1) = 0.0e+00_wp_ elist(1) = 0.0e+00_wp_ iord(1) = 0 if(epsabs<=0.0e+00_wp_.and.epsrelerrbnd) ier = 2 if(limit==1) ier = 1 if(ier/=0.or.(abserr<=errbnd.and.abserr/=resabs).or. & abserr==0.0e+00_wp_) go to 130 ! ! initialization ! -------------- ! rlist2(1) = result errmax = abserr maxerr = 1 area = result errsum = abserr abserr = oflow nrmax = 1 nres = 0 ktmin = 0 numrl2 = 2 extrap = .false. noext = .false. ierro = 0 iroff1 = 0 iroff2 = 0 iroff3 = 0 ksgn = -1 if(dres>=(0.1e+01_wp_-0.5e+02_wp_*epmach)*defabs) ksgn = 1 ! ! main do-loop ! ------------ ! do last = 2,limit ! ! bisect the subinterval with nrmax-th largest error estimate. ! a1 = alist(maxerr) b1 = 0.5e+00_wp_*(alist(maxerr)+blist(maxerr)) a2 = b1 b2 = blist(maxerr) erlast = errmax call dqk15i(f,boun,inf,a1,b1,area1,error1,resabs,defab1) call dqk15i(f,boun,inf,a2,b2,area2,error2,resabs,defab2) ! ! improve previous approximations to integral ! and error and test for accuracy. ! area12 = area1+area2 erro12 = error1+error2 errsum = errsum+erro12-errmax area = area+area12-rlist(maxerr) if(defab1/=error1.and.defab2/=error2) then if(abs(rlist(maxerr)-area12)<=0.1e-04_wp_*abs(area12) & .and.erro12>=0.99e+00_wp_*errmax) then if(extrap) iroff2 = iroff2+1 if(.not.extrap) iroff1 = iroff1+1 end if if(last>10.and.erro12>errmax) iroff3 = iroff3+1 end if rlist(maxerr) = area1 rlist(last) = area2 errbnd = dmax1(epsabs,epsrel*abs(area)) ! ! test for roundoff error and eventually set error flag. ! if(iroff1+iroff2>=10.or.iroff3>=20) ier = 2 if(iroff2>=5) ierro = 3 ! ! set error flag in the case that the number of ! subintervals equals limit. ! if(last==limit) ier = 1 ! ! set error flag in the case of bad integrand behaviour ! at some points of the integration range. ! if(dmax1(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* & (abs(a2)+0.1e+04_wp_*uflow)) ier = 4 ! ! append the newly-created intervals to the list. ! if(error2<=error1) then alist(last) = a2 blist(maxerr) = b1 blist(last) = b2 elist(maxerr) = error1 elist(last) = error2 else alist(maxerr) = a2 alist(last) = a1 blist(last) = b1 rlist(maxerr) = area2 rlist(last) = area1 elist(maxerr) = error2 elist(last) = error1 end if ! ! call subroutine dqpsrt to maintain the descending ordering ! in the list of error estimates and select the subinterval ! with nrmax-th largest error estimate (to be bisected next). ! call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) if(errsum<=errbnd) go to 115 if(ier/=0) exit if(last==2) then small = 0.375e+00_wp_ erlarg = errsum ertest = errbnd rlist2(2) = area cycle end if if(noext) cycle erlarg = erlarg-erlast if(abs(b1-a1)>small) erlarg = erlarg+erro12 if(.not.extrap) then ! ! test whether the interval to be bisected next is the ! smallest interval. ! if(abs(blist(maxerr)-alist(maxerr))>small) cycle extrap = .true. nrmax = 2 end if if(ierro/=3.and.erlarg>ertest) then ! ! the smallest interval has the largest error. ! before bisecting decrease the sum of the errors over the ! larger intervals (erlarg) and perform extrapolation. ! id = nrmax jupbnd = last if(last>(2+limit/2)) jupbnd = limit+3-last do k = id,jupbnd maxerr = iord(nrmax) errmax = elist(maxerr) if(abs(blist(maxerr)-alist(maxerr))>small) go to 90 nrmax = nrmax+1 end do end if ! ! perform extrapolation. ! numrl2 = numrl2+1 rlist2(numrl2) = area call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres) ktmin = ktmin+1 if(ktmin>5.and.abserr<0.1e-02_wp_*errsum) ier = 5 if(absepserrsum)go to 115 if(area==0.0e+00_wp_) go to 130 go to 110 105 continue if(abserr/abs(result)>errsum/abs(area)) go to 115 ! ! test on divergence ! 110 continue if(ksgn==(-1).and.dmax1(abs(result),abs(area))<= & defabs*0.1e-01_wp_) go to 130 if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ & .or.errsum>abs(area)) ier = 6 go to 130 ! ! compute global integral sum. ! 115 continue result = 0.0e+00_wp_ do k = 1,last result = result+rlist(k) end do abserr = errsum 130 continue neval = 30*last-15 if(inf==2) neval = 2*neval if(ier>2) ier=ier-1 end subroutine dqagie subroutine dqk15i(f,boun,inf,a,b,result,abserr,resabs,resasc) !***begin prologue dqk15i !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a3a2,h2a4a2 !***keywords 15-point transformed gauss-kronrod rules !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & progr. div. - k.u.leuven !***purpose the original (infinite integration range is mapped ! onto the interval (0,1) and (a,b) is a part of (0,1). ! it is the purpose to compute ! i = integral of transformed integrand over (a,b), ! j = integral of abs(transformed integrand) over (a,b). !***description ! ! integration rule ! standard fortran subroutine ! real(8) version ! ! parameters ! on entry ! f - real(8) ! fuction subprogram defining the integrand ! function f(x). the actual name for f needs to be ! declared e x t e r n a l in the calling program. ! ! boun - real(8) ! finite bound of original integration ! range (set to zero if inf = +2) ! ! inf - integer ! if inf = -1, the original interval is ! (-infinity,bound), ! if inf = +1, the original interval is ! (bound,+infinity), ! if inf = +2, the original interval is ! (-infinity,+infinity) and ! the integral is computed as the sum of two ! integrals, one over (-infinity,0) and one over ! (0,+infinity). ! ! a - real(8) ! lower limit for integration over subrange ! of (0,1) ! ! b - real(8) ! upper limit for integration over subrange ! of (0,1) ! ! on return ! result - real(8) ! approximation to the integral i ! result is computed by applying the 15-point ! kronrod rule(resk) obtained by optimal addition ! of abscissae to the 7-point gauss rule(resg). ! ! abserr - real(8) ! estimate of the modulus of the absolute error, ! which should equal or exceed abs(i-result) ! ! resabs - real(8) ! approximation to the integral j ! ! resasc - real(8) ! approximation to the integral of ! abs((transformed integrand)-i/(b-a)) over (a,b) ! !***references (none) !***routines called (none) !***end prologue dqk15i ! use const_and_precisions, only : epmach=>comp_eps, uflow=>comp_tiny implicit none real(wp_), intent(in) :: a,b,boun integer, intent(in) :: inf real(wp_), intent(out) :: result,abserr,resabs,resasc real(wp_), external :: f real(wp_) :: absc,absc1,absc2,centr,abs,dinf,dmax1,dmin1,fc,fsum, & fval1,fval2,hlgth,resg,resk,reskh,tabsc1,tabsc2 real(wp_), dimension(7) :: fv1,fv2 integer :: j ! ! the abscissae and weights are supplied for the interval ! (-1,1). because of symmetry only the positive abscissae and ! their corresponding weights are given. ! ! xgk - abscissae of the 15-point kronrod rule ! xgk(2), xgk(4), ... abscissae of the 7-point ! gauss rule ! xgk(1), xgk(3), ... abscissae which are optimally ! added to the 7-point gauss rule ! ! wgk - weights of the 15-point kronrod rule ! ! wg - weights of the 7-point gauss rule, corresponding ! to the abscissae xgk(2), xgk(4), ... ! wg(1), wg(3), ... are set to zero. ! real(wp_), dimension(8), parameter :: & wg = (/ 0.000000000000000000000000000000000_wp_, & 0.129484966168869693270611432679082_wp_, & 0.000000000000000000000000000000000_wp_, & 0.279705391489276667901467771423780_wp_, & 0.000000000000000000000000000000000_wp_, & 0.381830050505118944950369775488975_wp_, & 0.000000000000000000000000000000000_wp_, & 0.417959183673469387755102040816327_wp_ /), & xgk = (/ 0.991455371120812639206854697526329_wp_, & 0.949107912342758524526189684047851_wp_, & 0.864864423359769072789712788640926_wp_, & 0.741531185599394439863864773280788_wp_, & 0.586087235467691130294144838258730_wp_, & 0.405845151377397166906606412076961_wp_, & 0.207784955007898467600689403773245_wp_, & 0.000000000000000000000000000000000_wp_ /), & wgk = (/ 0.022935322010529224963732008058970_wp_, & 0.063092092629978553290700663189204_wp_, & 0.104790010322250183839876322541518_wp_, & 0.140653259715525918745189590510238_wp_, & 0.169004726639267902826583426598550_wp_, & 0.190350578064785409913256402421014_wp_, & 0.204432940075298892414161999234649_wp_, & 0.209482141084727828012999174891714_wp_ /) ! ! ! list of major variables ! ----------------------- ! ! centr - mid point of the interval ! hlgth - half-length of the interval ! absc* - abscissa ! tabsc* - transformed abscissa ! fval* - function value ! resg - result of the 7-point gauss formula ! resk - result of the 15-point kronrod formula ! reskh - approximation to the mean value of the transformed ! integrand over (a,b), i.e. to i/(b-a) ! ! machine dependent constants ! --------------------------- ! ! epmach is the largest relative spacing. ! uflow is the smallest positive magnitude. ! !***first executable statement dqk15i dinf = min0(1,inf) ! centr = 0.5e+00_wp_*(a+b) hlgth = 0.5e+00_wp_*(b-a) tabsc1 = boun+dinf*(0.1e+01_wp_-centr)/centr fval1 = f(tabsc1) if(inf==2) fval1 = fval1+f(-tabsc1) fc = (fval1/centr)/centr ! ! compute the 15-point kronrod approximation to ! the integral, and estimate the error. ! resg = wg(8)*fc resk = wgk(8)*fc resabs = abs(resk) do j=1,7 absc = hlgth*xgk(j) absc1 = centr-absc absc2 = centr+absc tabsc1 = boun+dinf*(0.1e+01_wp_-absc1)/absc1 tabsc2 = boun+dinf*(0.1e+01_wp_-absc2)/absc2 fval1 = f(tabsc1) fval2 = f(tabsc2) if(inf==2) fval1 = fval1+f(-tabsc1) if(inf==2) fval2 = fval2+f(-tabsc2) fval1 = (fval1/absc1)/absc1 fval2 = (fval2/absc2)/absc2 fv1(j) = fval1 fv2(j) = fval2 fsum = fval1+fval2 resg = resg+wg(j)*fsum resk = resk+wgk(j)*fsum resabs = resabs+wgk(j)*(abs(fval1)+abs(fval2)) end do reskh = resk*0.5e+00_wp_ resasc = wgk(8)*abs(fc-reskh) do j=1,7 resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)) end do result = resk*hlgth resasc = resasc*hlgth resabs = resabs*hlgth abserr = abs((resk-resg)*hlgth) if(resasc/=0.0e+00_wp_.and.abserr/=0._wp_) abserr = resasc* & dmin1(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_) if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 & ((epmach*0.5e+02_wp_)*resabs,abserr) end subroutine dqk15i subroutine dqagp(f,a,b,npts2,points,epsabs,epsrel,result,abserr, & neval,ier,leniw,lenw,last,iwork,work) !***begin prologue dqagp !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a2a1 !***keywords automatic integrator, general-purpose, ! singularities at user specified points, ! extrapolation, globally adaptive !***author piessens,robert,appl. math. & progr. div - k.u.leuven ! de doncker,elise,appl. math. & progr. div. - k.u.leuven !***purpose the routine calculates an approximation result to a given ! definite integral i = integral of f over (a,b), ! hopefully satisfying following claim for accuracy ! break points of the integration interval, where local ! difficulties of the integrand may occur (e.g. ! singularities, discontinuities), are provided by the user. !***description ! ! computation of a definite integral ! standard fortran subroutine ! double precision version ! ! parameters ! on entry ! f - double precision ! function subprogram defining the integrand ! function f(x). the actual name for f needs to be ! declared e x t e r n a l in the driver program. ! ! a - double precision ! lower limit of integration ! ! b - double precision ! upper limit of integration ! ! npts2 - integer ! number equal to two more than the number of ! user-supplied break points within the integration ! range, npts.ge.2. ! if npts2.lt.2, the routine will end with ier = 6. ! ! points - double precision ! vector of dimension npts2, the first (npts2-2) ! elements of which are the user provided break ! points. if these points do not constitute an ! ascending sequence there will be an automatic ! sorting. ! ! epsabs - double precision ! absolute accuracy requested ! epsrel - double precision ! relative accuracy requested ! if epsabs.le.0 ! and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), ! the routine will end with ier = 6. ! ! on return ! result - double precision ! approximation to the integral ! ! abserr - double precision ! estimate of the modulus of the absolute error, ! which should equal or exceed abs(i-result) ! ! neval - integer ! number of integrand evaluations ! ! ier - integer ! ier = 0 normal and reliable termination of the ! routine. it is assumed that the requested ! accuracy has been achieved. ! ier.gt.0 abnormal termination of the routine. ! the estimates for integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! error messages ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the value of ! limit (and taking the according dimension ! adjustments into account). however, if ! this yields no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulties. if ! the position of a local difficulty can be ! determined (i.e. singularity, ! discontinuity within the interval), it ! should be supplied to the routine as an ! element of the vector points. if necessary ! an appropriate special-purpose integrator ! must be used, which is designed for ! handling the type of difficulty involved. ! = 2 the occurrence of roundoff error is ! detected, which prevents the requested ! tolerance from being achieved. ! the error may be under-estimated. ! = 3 extremely bad integrand behaviour occurs ! at some points of the integration ! interval. ! = 4 the algorithm does not converge. ! roundoff error is detected in the ! extrapolation table. ! it is presumed that the requested ! tolerance cannot be achieved, and that ! the returned result is the best which ! can be obtained. ! = 5 the integral is probably divergent, or ! slowly convergent. it must be noted that ! divergence can occur with any other value ! of ier.gt.0. ! = 6 the input is invalid because ! npts2.lt.2 or ! break points are specified outside ! the integration range or ! (epsabs.le.0 and ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) ! result, abserr, neval, last are set to ! zero. exept when leniw or lenw or npts2 is ! invalid, iwork(1), iwork(limit+1), ! work(limit*2+1) and work(limit*3+1) ! are set to zero. ! work(1) is set to a and work(limit+1) ! to b (where limit = (leniw-npts2)/2). ! ! dimensioning parameters ! leniw - integer ! dimensioning parameter for iwork ! leniw determines limit = (leniw-npts2)/2, ! which is the maximum number of subintervals in the ! partition of the given integration interval (a,b), ! leniw.ge.(3*npts2-2). ! if leniw.lt.(3*npts2-2), the routine will end with ! ier = 6. ! ! lenw - integer ! dimensioning parameter for work ! lenw must be at least leniw*2-npts2. ! if lenw.lt.leniw*2-npts2, the routine will end ! with ier = 6. ! ! last - integer ! on return, last equals the number of subintervals ! produced in the subdivision process, which ! determines the number of significant elements ! actually in the work arrays. ! ! work arrays ! iwork - integer ! vector of dimension at least leniw. on return, ! the first k elements of which contain ! pointers to the error estimates over the ! subintervals, such that work(limit*3+iwork(1)),..., ! work(limit*3+iwork(k)) form a decreasing ! sequence, with k = last if last.le.(limit/2+2), and ! k = limit+1-last otherwise ! iwork(limit+1), ...,iwork(limit+last) contain the ! subdivision levels of the subintervals, i.e. ! if (aa,bb) is a subinterval of (p1,p2) ! where p1 as well as p2 is a user-provided ! break point or integration limit, then (aa,bb) has ! level l if abs(bb-aa) = abs(p2-p1)*2**(-l), ! iwork(limit*2+1), ..., iwork(limit*2+npts2) have ! no significance for the user, ! note that limit = (leniw-npts2)/2. ! ! work - double precision ! vector of dimension at least lenw ! on return ! work(1), ..., work(last) contain the left ! end points of the subintervals in the ! partition of (a,b), ! work(limit+1), ..., work(limit+last) contain ! the right end points, ! work(limit*2+1), ..., work(limit*2+last) contain ! the integral approximations over the subintervals, ! work(limit*3+1), ..., work(limit*3+last) ! contain the corresponding error estimates, ! work(limit*4+1), ..., work(limit*4+npts2) ! contain the integration limits and the ! break points sorted in an ascending sequence. ! note that limit = (leniw-npts2)/2. ! !***references (none) !***routines called dqagpe,xerror !***end prologue dqagp ! implicit none real(wp_), intent(in) :: a,b,epsabs,epsrel integer, intent(in) :: npts2,lenw,leniw real(wp_), intent(in), dimension(npts2) ::points real(wp_), intent(out) :: abserr,result integer, intent(out) :: neval,ier,last integer :: limit,lvl,l1,l2,l3,l4 ! real(wp_), dimension(lenw) :: work integer, dimension(leniw) :: iwork ! real(wp_), external :: f ! ! check validity of limit and lenw. ! !***first executable statement dqagp ier = 6 neval = 0 last = 0 result = 0.0_wp_ abserr = 0.0_wp_ if(leniw.ge.(3*npts2-2).and.lenw.ge.(leniw*2-npts2).and.npts2.ge.2) then ! ! prepare call for dqagpe. ! limit = (leniw-npts2)/2 l1 = limit+1 l2 = limit+l1 l3 = limit+l2 l4 = limit+l3 ! call dqagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,abserr, & neval,ier,work(1),work(l1),work(l2),work(l3),work(l4), & iwork(1),iwork(l1),iwork(l2),last) ! ! call error handler if necessary. ! lvl = 0 end if if(ier.eq.6) lvl = 1 if(ier.ne.0) print*,'habnormal return from dqaps',ier,lvl end subroutine dqagp subroutine dqagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result, & abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin,last) !***begin prologue dqagpe !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a2a1 !***keywords automatic integrator, general-purpose, ! singularities at user specified points, ! extrapolation, globally adaptive. !***author piessens,robert ,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & progr. div. - k.u.leuven !***purpose the routine calculates an approximation result to a given ! definite integral i = integral of f over (a,b), hopefully ! satisfying following claim for accuracy abs(i-result).le. ! max(epsabs,epsrel*abs(i)). break points of the integration ! interval, where local difficulties of the integrand may ! occur(e.g. singularities,discontinuities),provided by user. !***description ! ! computation of a definite integral ! standard fortran subroutine ! double precision version ! ! parameters ! on entry ! f - double precision ! function subprogram defining the integrand ! function f(x). the actual name for f needs to be ! declared e x t e r n a l in the driver program. ! ! a - double precision ! lower limit of integration ! ! b - double precision ! upper limit of integration ! ! npts2 - integer ! number equal to two more than the number of ! user-supplied break points within the integration ! range, npts2.ge.2. ! if npts2.lt.2, the routine will end with ier = 6. ! ! points - double precision ! vector of dimension npts2, the first (npts2-2) ! elements of which are the user provided break ! points. if these points do not constitute an ! ascending sequence there will be an automatic ! sorting. ! ! epsabs - double precision ! absolute accuracy requested ! epsrel - double precision ! relative accuracy requested ! if epsabs.le.0 ! and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), ! the routine will end with ier = 6. ! ! limit - integer ! gives an upper bound on the number of subintervals ! in the partition of (a,b), limit.ge.npts2 ! if limit.lt.npts2, the routine will end with ! ier = 6. ! ! on return ! result - double precision ! approximation to the integral ! ! abserr - double precision ! estimate of the modulus of the absolute error, ! which should equal or exceed abs(i-result) ! ! neval - integer ! number of integrand evaluations ! ! ier - integer ! ier = 0 normal and reliable termination of the ! routine. it is assumed that the requested ! accuracy has been achieved. ! ier.gt.0 abnormal termination of the routine. ! the estimates for integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! error messages ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the value of ! limit (and taking the according dimension ! adjustments into account). however, if ! this yields no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulties. if ! the position of a local difficulty can be ! determined (i.e. singularity, ! discontinuity within the interval), it ! should be supplied to the routine as an ! element of the vector points. if necessary ! an appropriate special-purpose integrator ! must be used, which is designed for ! handling the type of difficulty involved. ! = 2 the occurrence of roundoff error is ! detected, which prevents the requested ! tolerance from being achieved. ! the error may be under-estimated. ! = 3 extremely bad integrand behaviour occurs ! at some points of the integration ! interval. ! = 4 the algorithm does not converge. ! roundoff error is detected in the ! extrapolation table. it is presumed that ! the requested tolerance cannot be ! achieved, and that the returned result is ! the best which can be obtained. ! = 5 the integral is probably divergent, or ! slowly convergent. it must be noted that ! divergence can occur with any other value ! of ier.gt.0. ! = 6 the input is invalid because ! npts2.lt.2 or ! break points are specified outside ! the integration range or ! (epsabs.le.0 and ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) ! or limit.lt.npts2. ! result, abserr, neval, last, rlist(1), ! and elist(1) are set to zero. alist(1) and ! blist(1) are set to a and b respectively. ! ! alist - double precision ! vector of dimension at least limit, the first ! last elements of which are the left end points ! of the subintervals in the partition of the given ! integration range (a,b) ! ! blist - double precision ! vector of dimension at least limit, the first ! last elements of which are the right end points ! of the subintervals in the partition of the given ! integration range (a,b) ! ! rlist - double precision ! vector of dimension at least limit, the first ! last elements of which are the integral ! approximations on the subintervals ! ! elist - double precision ! vector of dimension at least limit, the first ! last elements of which are the moduli of the ! absolute error estimates on the subintervals ! ! pts - double precision ! vector of dimension at least npts2, containing the ! integration limits and the break points of the ! interval in ascending sequence. ! ! level - integer ! vector of dimension at least limit, containing the ! subdivision levels of the subinterval, i.e. if ! (aa,bb) is a subinterval of (p1,p2) where p1 as ! well as p2 is a user-provided break point or ! integration limit, then (aa,bb) has level l if ! abs(bb-aa) = abs(p2-p1)*2**(-l). ! ! ndin - integer ! vector of dimension at least npts2, after first ! integration over the intervals (pts(i)),pts(i+1), ! i = 0,1, ..., npts2-2, the error estimates over ! some of the intervals may have been increased ! artificially, in order to put their subdivision ! forward. if this happens for the subinterval ! numbered k, ndin(k) is put to 1, otherwise ! ndin(k) = 0. ! ! iord - integer ! vector of dimension at least limit, the first k ! elements of which are pointers to the ! error estimates over the subintervals, ! such that elist(iord(1)), ..., elist(iord(k)) ! form a decreasing sequence, with k = last ! if last.le.(limit/2+2), and k = limit+1-last ! otherwise ! ! last - integer ! number of subintervals actually produced in the ! subdivisions process ! !***references (none) !***routines called dqelg,dqk21,dqpsrt !***end prologue dqagpe use const_and_precisions, only : epmach=>comp_eps, uflow=>comp_tiny, & oflow=>comp_huge implicit none real(wp_) :: a,abseps,abserr,alist,area,area1,area12,area2,a1, & a2,b,blist,b1,b2,correc,dabs,defabs,defab1,defab2,dmax1,dmin1, & dres,elist,epsabs,epsrel,erlarg,erlast,errbnd, & errmax,error1,erro12,error2,errsum,ertest,points,pts, & resa,resabs,reseps,result,res3la,rlist,rlist2,sign,temp integer :: i,id,ier,ierro,ind1,ind2,iord,ip1,iroff1,iroff2,iroff3,j, & jlow,jupbnd,k,ksgn,ktmin,last,levcur,level,levmax,limit,maxerr, & ndin,neval,nint,nintp1,npts,npts2,nres,nrmax,numrl2 logical :: extrap,noext ! ! dimension alist(limit),blist(limit),elist(limit),iord(limit), & level(limit),ndin(npts2),points(npts2),pts(npts2),res3la(3), & rlist(limit),rlist2(52) ! real(wp_), external :: f ! ! the dimension of rlist2 is determined by the value of ! limexp in subroutine epsalg (rlist2 should be of dimension ! (limexp+2) at least). ! ! ! list of major variables ! ----------------------- ! ! alist - list of left end points of all subintervals ! considered up to now ! blist - list of right end points of all subintervals ! considered up to now ! rlist(i) - approximation to the integral over ! (alist(i),blist(i)) ! rlist2 - array of dimension at least limexp+2 ! containing the part of the epsilon table which ! is still needed for further computations ! elist(i) - error estimate applying to rlist(i) ! maxerr - pointer to the interval with largest error ! estimate ! errmax - elist(maxerr) ! erlast - error on the interval currently subdivided ! (before that subdivision has taken place) ! area - sum of the integrals over the subintervals ! errsum - sum of the errors over the subintervals ! errbnd - requested accuracy max(epsabs,epsrel* ! abs(result)) ! *****1 - variable for the left subinterval ! *****2 - variable for the right subinterval ! last - index for subdivision ! nres - number of calls to the extrapolation routine ! numrl2 - number of elements in rlist2. if an appropriate ! approximation to the compounded integral has ! been obtained, it is put in rlist2(numrl2) after ! numrl2 has been increased by one. ! erlarg - sum of the errors over the intervals larger ! than the smallest interval considered up to now ! extrap - logical variable denoting that the routine ! is attempting to perform extrapolation. i.e. ! before subdividing the smallest interval we ! try to decrease the value of erlarg. ! noext - logical variable denoting that extrapolation is ! no longer allowed (true-value) ! ! machine dependent constants ! --------------------------- ! ! epmach is the largest relative spacing. ! uflow is the smallest positive magnitude. ! oflow is the largest positive magnitude. ! !***first executable statement dqagpe ! ! test on validity of parameters ! ----------------------------- ! ier = 0 neval = 0 last = 0 result = 0.0_wp_ abserr = 0.0_wp_ alist(1) = a blist(1) = b rlist(1) = 0.0_wp_ elist(1) = 0.0_wp_ iord(1) = 0 level(1) = 0 npts = npts2-2 if(npts2.lt.2.or.limit.le.npts.or.(epsabs.le.0.0_wp_.and. & epsrel.lt.dmax1(0.5e+02_wp_*epmach,0.5e-28_wp_))) ier = 6 if(ier.eq.6) return ! ! if any break points are provided, sort them into an ! ascending sequence. ! sign = 1.0_wp_ if(a.gt.b) sign = -1.0_wp_ pts(1) = dmin1(a,b) do i = 1,npts pts(i+1) = points(i) end do pts(npts+2) = dmax1(a,b) nint = npts+1 a1 = pts(1) if(npts.ne.0) then nintp1 = nint+1 do i = 1,nint ip1 = i+1 do j = ip1,nintp1 if(pts(i).gt.pts(j)) then temp = pts(i) pts(i) = pts(j) pts(j) = temp end if end do end do if(pts(1).ne.dmin1(a,b).or.pts(nintp1).ne.dmax1(a,b)) ier = 6 if(ier.eq.6) return end if ! ! compute first integral and error approximations. ! ------------------------------------------------ ! resabs = 0.0_wp_ do i = 1,nint b1 = pts(i+1) call dqk21(f,a1,b1,area1,error1,defabs,resa) abserr = abserr+error1 result = result+area1 ndin(i) = 0 if(error1.eq.resa.and.error1.ne.0.0_wp_) ndin(i) = 1 resabs = resabs+defabs level(i) = 0 elist(i) = error1 alist(i) = a1 blist(i) = b1 rlist(i) = area1 iord(i) = i a1 = b1 end do errsum = 0.0_wp_ do i = 1,nint if(ndin(i).eq.1) elist(i) = abserr errsum = errsum+elist(i) end do ! ! test on accuracy. ! last = nint neval = 21*nint dres = dabs(result) errbnd = dmax1(epsabs,epsrel*dres) if(abserr.le.0.1d+03*epmach*resabs.and.abserr.gt.errbnd) ier = 2 if(nint.ne.1) then do i = 1,npts jlow = i+1 ind1 = iord(i) do j = jlow,nint ind2 = iord(j) if(elist(ind1).le.elist(ind2)) then ind1 = ind2 k = j end if end do if(ind1.ne.iord(i)) then iord(k) = iord(i) iord(i) = ind1 end if end do if(limit.lt.npts2) ier = 1 end if if(ier.eq.0.and.abserr.gt.errbnd) then ! ! initialization ! -------------- ! rlist2(1) = result maxerr = iord(1) errmax = elist(maxerr) area = result nrmax = 1 nres = 0 numrl2 = 1 ktmin = 0 extrap = .false. noext = .false. erlarg = errsum ertest = errbnd levmax = 1 iroff1 = 0 iroff2 = 0 iroff3 = 0 ierro = 0 abserr = oflow ksgn = -1 if(dres.ge.(0.1d+01-0.5d+02*epmach)*resabs) ksgn = 1 ! ! main do-loop ! ------------ ! do last = npts2,limit ! ! bisect the subinterval with the nrmax-th largest error ! estimate. ! levcur = level(maxerr)+1 a1 = alist(maxerr) b1 = 0.5d+00*(alist(maxerr)+blist(maxerr)) a2 = b1 b2 = blist(maxerr) erlast = errmax call dqk21(f,a1,b1,area1,error1,resa,defab1) call dqk21(f,a2,b2,area2,error2,resa,defab2) ! ! improve previous approximations to integral ! and error and test for accuracy. ! neval = neval+42 area12 = area1+area2 erro12 = error1+error2 errsum = errsum+erro12-errmax area = area+area12-rlist(maxerr) if(defab1.ne.error1.and.defab2.ne.error2) then if(dabs(rlist(maxerr)-area12).le.0.1d-04*dabs(area12) & .and.erro12.ge.0.99d+00*errmax) then if(extrap) iroff2 = iroff2+1 if(.not.extrap) iroff1 = iroff1+1 end if if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1 end if level(maxerr) = levcur level(last) = levcur rlist(maxerr) = area1 rlist(last) = area2 errbnd = dmax1(epsabs,epsrel*dabs(area)) ! ! test for roundoff error and eventually set error flag. ! if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2 if(iroff2.ge.5) ierro = 3 ! ! set error flag in the case that the number of ! subintervals equals limit. ! if(last.eq.limit) ier = 1 ! ! set error flag in the case of bad integrand behaviour ! at a point of the integration range ! if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach)* & (dabs(a2)+0.1d+04*uflow)) ier = 4 ! ! append the newly-created intervals to the list. ! if(error2.le.error1) then alist(last) = a2 blist(maxerr) = b1 blist(last) = b2 elist(maxerr) = error1 elist(last) = error2 else alist(maxerr) = a2 alist(last) = a1 blist(last) = b1 rlist(maxerr) = area2 rlist(last) = area1 elist(maxerr) = error2 elist(last) = error1 end if ! ! call subroutine dqpsrt to maintain the descending ordering ! in the list of error estimates and select the subinterval ! with nrmax-th largest error estimate (to be bisected next). ! call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) ! ***jump out of do-loop if(errsum.le.errbnd) go to 190 ! ***jump out of do-loop if(ier.ne.0) exit if(noext) cycle erlarg = erlarg-erlast if(levcur+1.le.levmax) erlarg = erlarg+erro12 if(.not.extrap) then ! ! test whether the interval to be bisected next is the ! smallest interval. ! if(level(maxerr)+1.le.levmax) cycle extrap = .true. nrmax = 2 end if if(ierro.ne.3.and.erlarg.gt.ertest) then ! ! the smallest interval has the largest error. ! before bisecting decrease the sum of the errors over ! the larger intervals (erlarg) and perform extrapolation. ! id = nrmax jupbnd = last if(last.gt.(2+limit/2)) jupbnd = limit+3-last do k = id,jupbnd maxerr = iord(nrmax) errmax = elist(maxerr) ! ***jump out of do-loop if(level(maxerr)+1.le.levmax) go to 160 nrmax = nrmax+1 end do end if ! ! perform extrapolation. ! numrl2 = numrl2+1 rlist2(numrl2) = area if(numrl2.gt.2) then call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres) ktmin = ktmin+1 if(ktmin.gt.5.and.abserr.lt.0.1d-02*errsum) ier = 5 if(abseps.lt.abserr) then ktmin = 0 abserr = abseps result = reseps correc = erlarg ertest = dmax1(epsabs,epsrel*dabs(reseps)) ! ***jump out of do-loop if(abserr.lt.ertest) exit end if ! ! prepare bisection of the smallest interval. ! if(numrl2.eq.1) noext = .true. if(ier.ge.5) exit end if maxerr = iord(1) errmax = elist(maxerr) nrmax = 1 extrap = .false. levmax = levmax+1 erlarg = errsum 160 continue end do ! ! set the final result. ! --------------------- ! ! if(abserr.eq.oflow) go to 190 if((ier+ierro).ne.0) then if(ierro.eq.3) abserr = abserr+correc if(ier.eq.0) ier = 3 if(result.ne.0.0d+00.and.area.ne.0.0d+00) then if(abserr/dabs(result).gt.errsum/dabs(area))go to 190 else if(abserr.gt.errsum)go to 190 if(area.eq.0.0d+00) go to 210 end if ! ! test on divergence. ! end if if(ksgn.ne.(-1).or.dmax1(dabs(result),dabs(area)).gt.resabs*0.1d-01) then if(0.1d-01.gt.(result/area).or.(result/area).gt.0.1d+03.or. & errsum.gt.dabs(area)) ier = 6 end if go to 210 ! ! compute global integral sum. ! 190 result = 0.0d+00 do k = 1,last result = result+rlist(k) end do abserr = errsum end if 210 if(ier.gt.2) ier = ier-1 result = result*sign end subroutine dqagpe ! ! ! Integration routine dqags.f from quadpack and dependencies: BEGIN ! Modified version for functions f(x,yi) with more than one variable ! ! subroutine dqagsmv(f,a,b,apar,np,epsabs,epsrel,result,abserr,neval,ier, & limit,lenw,last,iwork,work) !***begin prologue dqagsmv !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a1a1 !***keywords automatic integrator, general-purpose, ! (end-point) singularities, extrapolation, ! globally adaptive !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & prog. div. - k.u.leuven !***purpose the routine calculates an approximation result to a given ! definite integral i = integral of f over (a,b), ! hopefully satisfying following claim for accuracy ! abs(i-result)<=max(epsabs,epsrel*abs(i)). !***description ! ! computation of a definite integral ! standard fortran subroutine ! real(8) version ! ! ! parameters ! on entry ! f - real(8) ! function subprogram defining the integrand ! function f(x,apar,np). the actual name for f needs to be ! declared e x t e r n a l in the driver program. ! ! a - real(8) ! lower limit of integration ! ! b - real(8) ! upper limit of integration ! ! apar - array of parameters of the integrand function f ! ! np - number of parameters. size of apar ! ! epsabs - real(8) ! absolute accuracy requested ! epsrel - real(8) ! relative accuracy requested ! if epsabs<=0 ! and epsrel0 abnormal termination of the routine ! the estimates for integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! error messages ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more sub- ! divisions by increasing the value of limit ! (and taking the according dimension ! adjustments into account. however, if ! this yields no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulties. if ! the position of a local difficulty can be ! determined (e.g. singularity, ! discontinuity within the interval) one ! will probably gain from splitting up the ! interval at this point and calling the ! integrator on the subranges. if possible, ! an appropriate special-purpose integrator ! should be used, which is designed for ! handling the type of difficulty involved. ! = 2 the occurrence of roundoff error is detec- ! ted, which prevents the requested ! tolerance from being achieved. ! the error may be under-estimated. ! = 3 extremely bad integrand behaviour ! occurs at some points of the integration ! interval. ! = 4 the algorithm does not converge. ! roundoff error is detected in the ! extrapolation table. it is presumed that ! the requested tolerance cannot be ! achieved, and that the returned result is ! the best which can be obtained. ! = 5 the integral is probably divergent, or ! slowly convergent. it must be noted that ! divergence can occur with any other value ! of ier. ! = 6 the input is invalid, because ! (epsabs<=0 and ! epsrel=1. ! if limit<1, the routine will end with ier = 6. ! ! lenw - integer ! dimensioning parameter for work ! lenw must be at least limit*4. ! if lenw=1.and.lenw>=limit*4) then ! ! prepare call for dqagse. ! l1 = limit+1 l2 = limit+l1 l3 = limit+l2 ! call dqagsemv(f,a,b,apar,np,epsabs,epsrel,limit,result,abserr,neval, & ier,work(1),work(l1),work(l2),work(l3),iwork,last) ! ! call error handler if necessary. ! lvl = 0 end if if(ier==6) lvl = 1 if(ier/=0) print*,'habnormal return from dqags',ier,lvl end subroutine dqagsmv subroutine dqagsemv(f,a,b,apar,np,epsabs,epsrel,limit,result,abserr,neval, & ier,alist,blist,rlist,elist,iord,last) !***begin prologue dqagsemv !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a1a1 !***keywords automatic integrator, general-purpose, ! (end point) singularities, extrapolation, ! globally adaptive !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & progr. div. - k.u.leuven !***purpose the routine calculates an approximation result to a given ! definite integral i = integral of f over (a,b), ! hopefully satisfying following claim for accuracy ! abs(i-result)<=max(epsabs,epsrel*abs(i)). !***description ! ! computation of a definite integral ! standard fortran subroutine ! real(8) version ! ! parameters ! on entry ! f - real(8) ! function subprogram defining the integrand ! function f(x,apar,np). the actual name for f needs to be ! declared e x t e r n a l in the driver program. ! ! a - real(8) ! lower limit of integration ! ! b - real(8) ! upper limit of integration ! ! apar - array of parameters of the integrand function f ! ! np - number of parameters. size of apar ! ! epsabs - real(8) ! absolute accuracy requested ! epsrel - real(8) ! relative accuracy requested ! if epsabs<=0 ! and epsrel0 abnormal termination of the routine ! the estimates for integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! error messages ! = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more sub- ! divisions by increasing the value of limit ! (and taking the according dimension ! adjustments into account). however, if ! this yields no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulties. if ! the position of a local difficulty can be ! determined (e.g. singularity, ! discontinuity within the interval) one ! will probably gain from splitting up the ! interval at this point and calling the ! integrator on the subranges. if possible, ! an appropriate special-purpose integrator ! should be used, which is designed for ! handling the type of difficulty involved. ! = 2 the occurrence of roundoff error is detec- ! ted, which prevents the requested ! tolerance from being achieved. ! the error may be under-estimated. ! = 3 extremely bad integrand behaviour ! occurs at some points of the integration ! interval. ! = 4 the algorithm does not converge. ! roundoff error is detected in the ! extrapolation table. ! it is presumed that the requested ! tolerance cannot be achieved, and that the ! returned result is the best which can be ! obtained. ! = 5 the integral is probably divergent, or ! slowly convergent. it must be noted that ! divergence can occur with any other value ! of ier. ! = 6 the input is invalid, because ! epsabs<=0 and ! epsrelcomp_eps, uflow=>comp_tiny, & oflow=>comp_huge implicit none real(wp_), intent(in) :: a,b,epsabs,epsrel integer, intent(in) :: limit,np real(wp_), dimension(np), intent(in) :: apar real(wp_), intent(out) :: result,abserr integer, intent(out) :: neval,ier,last real(wp_), dimension(limit), intent(inout) :: alist,blist,elist,rlist integer, dimension(limit), intent(inout) :: iord real(wp_), external :: f ! real(wp_) :: abseps,area,area1,area12,area2,a1,a2,b1,b2,correc,abs, & defabs,defab1,defab2,dmax1,dres,erlarg,erlast,errbnd,errmax, & error1,error2,erro12,errsum,ertest,resabs,reseps,small real(wp_) :: res3la(3),rlist2(52) integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, & ktmin,maxerr,nres,nrmax,numrl2 logical :: extrap,noext ! ! the dimension of rlist2 is determined by the value of ! limexp in subroutine dqelg (rlist2 should be of dimension ! (limexp+2) at least). ! ! list of major variables ! ----------------------- ! ! alist - list of left end points of all subintervals ! considered up to now ! blist - list of right end points of all subintervals ! considered up to now ! rlist(i) - approximation to the integral over ! (alist(i),blist(i)) ! rlist2 - array of dimension at least limexp+2 containing ! the part of the epsilon table which is still ! needed for further computations ! elist(i) - error estimate applying to rlist(i) ! maxerr - pointer to the interval with largest error ! estimate ! errmax - elist(maxerr) ! erlast - error on the interval currently subdivided ! (before that subdivision has taken place) ! area - sum of the integrals over the subintervals ! errsum - sum of the errors over the subintervals ! errbnd - requested accuracy max(epsabs,epsrel* ! abs(result)) ! *****1 - variable for the left interval ! *****2 - variable for the right interval ! last - index for subdivision ! nres - number of calls to the extrapolation routine ! numrl2 - number of elements currently in rlist2. if an ! appropriate approximation to the compounded ! integral has been obtained it is put in ! rlist2(numrl2) after numrl2 has been increased ! by one. ! small - length of the smallest interval considered up ! to now, multiplied by 1.5 ! erlarg - sum of the errors over the intervals larger ! than the smallest interval considered up to now ! extrap - logical variable denoting that the routine is ! attempting to perform extrapolation i.e. before ! subdividing the smallest interval we try to ! decrease the value of erlarg. ! noext - logical variable denoting that extrapolation ! is no longer allowed (true value) ! ! machine dependent constants ! --------------------------- ! ! epmach is the largest relative spacing. ! uflow is the smallest positive magnitude. ! oflow is the largest positive magnitude. ! !***first executable statement dqagsemv ! ! test on validity of parameters ! ------------------------------ ier = 0 neval = 0 last = 0 result = 0.0e+00_wp_ abserr = 0.0e+00_wp_ alist(1) = a blist(1) = b rlist(1) = 0.0e+00_wp_ elist(1) = 0.0e+00_wp_ if(epsabs<=0.0e+00_wp_.and.epsrelerrbnd) ier = 2 if(limit==1) ier = 1 if(ier/=0.or.(abserr<=errbnd.and.abserr/=resabs).or. & abserr==0.0e+00_wp_) go to 140 ! ! initialization ! -------------- ! rlist2(1) = result errmax = abserr maxerr = 1 area = result errsum = abserr abserr = oflow nrmax = 1 nres = 0 numrl2 = 2 ktmin = 0 extrap = .false. noext = .false. iroff1 = 0 iroff2 = 0 iroff3 = 0 ksgn = -1 if(dres>=(0.1e+01_wp_-0.5e+02_wp_*epmach)*defabs) ksgn = 1 ! ! main do-loop ! ------------ ! do last = 2,limit ! ! bisect the subinterval with the nrmax-th largest error ! estimate. ! a1 = alist(maxerr) b1 = 0.5e+00_wp_*(alist(maxerr)+blist(maxerr)) a2 = b1 b2 = blist(maxerr) erlast = errmax call dqk21mv(f,a1,b1,apar,np,area1,error1,resabs,defab1) call dqk21mv(f,a2,b2,apar,np,area2,error2,resabs,defab2) ! ! improve previous approximations to integral ! and error and test for accuracy. ! area12 = area1+area2 erro12 = error1+error2 errsum = errsum+erro12-errmax area = area+area12-rlist(maxerr) if(defab1/=error1.and.defab2/=error2) then if(abs(rlist(maxerr)-area12)<=0.1e-04_wp_*abs(area12) & .and.erro12>=0.99e+00_wp_*errmax) then if(extrap) iroff2 = iroff2+1 if(.not.extrap) iroff1 = iroff1+1 end if if(last>10.and.erro12>errmax) iroff3 = iroff3+1 end if rlist(maxerr) = area1 rlist(last) = area2 errbnd = dmax1(epsabs,epsrel*abs(area)) ! ! test for roundoff error and eventually set error flag. ! if(iroff1+iroff2>=10.or.iroff3>=20) ier = 2 if(iroff2>=5) ierro = 3 ! ! set error flag in the case that the number of subintervals ! equals limit. ! if(last==limit) ier = 1 ! ! set error flag in the case of bad integrand behaviour ! at a point of the integration range. ! if(dmax1(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* & (abs(a2)+0.1e+04_wp_*uflow)) ier = 4 ! ! append the newly-created intervals to the list. ! if(error2<=error1) then alist(last) = a2 blist(maxerr) = b1 blist(last) = b2 elist(maxerr) = error1 elist(last) = error2 else alist(maxerr) = a2 alist(last) = a1 blist(last) = b1 rlist(maxerr) = area2 rlist(last) = area1 elist(maxerr) = error2 elist(last) = error1 end if ! ! call subroutine dqpsrt to maintain the descending ordering ! in the list of error estimates and select the subinterval ! with nrmax-th largest error estimate (to be bisected next). ! call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) ! ***jump out of do-loop if(errsum<=errbnd) go to 115 ! ***jump out of do-loop if(ier/=0) exit if(last==2) then small = abs(b-a)*0.375e+00_wp_ erlarg = errsum ertest = errbnd rlist2(2) = area cycle end if if(noext) cycle erlarg = erlarg-erlast if(abs(b1-a1)>small) erlarg = erlarg+erro12 if(.not.extrap) then ! ! test whether the interval to be bisected next is the ! smallest interval. ! if(abs(blist(maxerr)-alist(maxerr))>small) cycle extrap = .true. nrmax = 2 end if if(ierro/=3.and.erlarg>ertest) then ! ! the smallest interval has the largest error. ! before bisecting decrease the sum of the errors over the ! larger intervals (erlarg) and perform extrapolation. ! id = nrmax jupbnd = last if(last>(2+limit/2)) jupbnd = limit+3-last do k = id,jupbnd maxerr = iord(nrmax) errmax = elist(maxerr) ! ***jump out of do-loop if(abs(blist(maxerr)-alist(maxerr))>small) go to 90 nrmax = nrmax+1 end do ! ! perform extrapolation. ! end if numrl2 = numrl2+1 rlist2(numrl2) = area call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres) ktmin = ktmin+1 if(ktmin>5.and.abserr<0.1e-02_wp_*errsum) ier = 5 if(absepserrsum) go to 115 if(area==0.0e+00_wp_) go to 130 go to 110 105 continue if(abserr/abs(result)>errsum/abs(area)) go to 115 ! ! test on divergence. ! 110 continue if(ksgn==(-1).and.dmax1(abs(result),abs(area))<= & defabs*0.1e-01_wp_) go to 130 if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ & .or.errsum>abs(area)) ier = 6 go to 130 ! ! compute global integral sum. ! 115 continue result = 0.0e+00_wp_ do k = 1,last result = result+rlist(k) end do abserr = errsum 130 continue if(ier>2) ier = ier-1 140 continue neval = 42*last-21 end subroutine dqagsemv subroutine dqk21mv(f,a,b,apar,np,result,abserr,resabs,resasc) !***begin prologue dqk21mv !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a1a2 !***keywords 21-point gauss-kronrod rules !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & progr. div. - k.u.leuven !***purpose to compute i = integral of f over (a,b), with error ! estimate ! j = integral of abs(f) over (a,b) !***description ! ! integration rules ! standard fortran subroutine ! real(8) version ! ! parameters ! on entry ! f - real(8) ! function subprogram defining the integrand ! function f(x,apar,np). the actual name for f needs to be ! declared e x t e r n a l in the driver program. ! ! a - real(8) ! lower limit of integration ! ! b - real(8) ! upper limit of integration ! ! apar - array of parameters of the integrand function f ! ! np - number of parameters. size of apar ! ! on return ! result - real(8) ! approximation to the integral i ! result is computed by applying the 21-point ! kronrod rule (resk) obtained by optimal addition ! of abscissae to the 10-point gauss rule (resg). ! ! abserr - real(8) ! estimate of the modulus of the absolute error, ! which should not exceed abs(i-result) ! ! resabs - real(8) ! approximation to the integral j ! ! resasc - real(8) ! approximation to the integral of abs(f-i/(b-a)) ! over (a,b) ! !***references (none) !***routines called (none) !***end prologue dqk21mv ! use const_and_precisions, only : epmach=>comp_eps, uflow=>comp_tiny implicit none real(wp_), intent(in) :: a,b integer, intent(in) :: np real(wp_), dimension(np), intent(in) :: apar real(wp_), intent(out) :: result,abserr,resabs,resasc real(wp_), external :: f real(wp_) :: absc,centr,abs,dhlgth,dmax1,dmin1,fc,fsum, & fval1,fval2,hlgth,resg,resk,reskh real(wp_), dimension(10) :: fv1,fv2 integer :: j,jtw,jtwm1 ! ! the abscissae and weights are given for the interval (-1,1). ! because of symmetry only the positive abscissae and their ! corresponding weights are given. ! ! xgk - abscissae of the 21-point kronrod rule ! xgk(2), xgk(4), ... abscissae of the 10-point ! gauss rule ! xgk(1), xgk(3), ... abscissae which are optimally ! added to the 10-point gauss rule ! ! wgk - weights of the 21-point kronrod rule ! ! wg - weights of the 10-point gauss rule ! ! ! gauss quadrature weights and kronron quadrature abscissae and weights ! as evaluated with 80 decimal digit arithmetic by l. w. fullerton, ! bell labs, nov. 1981. ! real(wp_), dimension(5), parameter :: & wg = (/ 0.066671344308688137593568809893332_wp_, & 0.149451349150580593145776339657697_wp_, & 0.219086362515982043995534934228163_wp_, & 0.269266719309996355091226921569469_wp_, & 0.295524224714752870173892994651338_wp_ /) ! real(wp_), dimension(11), parameter :: & xgk = (/ 0.995657163025808080735527280689003_wp_, & 0.973906528517171720077964012084452_wp_, & 0.930157491355708226001207180059508_wp_, & 0.865063366688984510732096688423493_wp_, & 0.780817726586416897063717578345042_wp_, & 0.679409568299024406234327365114874_wp_, & 0.562757134668604683339000099272694_wp_, & 0.433395394129247190799265943165784_wp_, & 0.294392862701460198131126603103866_wp_, & 0.148874338981631210884826001129720_wp_, & 0.000000000000000000000000000000000_wp_ /), & wgk = (/ 0.011694638867371874278064396062192_wp_, & 0.032558162307964727478818972459390_wp_, & 0.054755896574351996031381300244580_wp_, & 0.075039674810919952767043140916190_wp_, & 0.093125454583697605535065465083366_wp_, & 0.109387158802297641899210590325805_wp_, & 0.123491976262065851077958109831074_wp_, & 0.134709217311473325928054001771707_wp_, & 0.142775938577060080797094273138717_wp_, & 0.147739104901338491374841515972068_wp_, & 0.149445554002916905664936468389821_wp_ /) ! ! ! list of major variables ! ----------------------- ! ! centr - mid point of the interval ! hlgth - half-length of the interval ! absc - abscissa ! fval* - function value ! resg - result of the 10-point gauss formula ! resk - result of the 21-point kronrod formula ! reskh - approximation to the mean value of f over (a,b), ! i.e. to i/(b-a) ! ! ! machine dependent constants ! --------------------------- ! ! epmach is the largest relative spacing. ! uflow is the smallest positive magnitude. ! !***first executable statement dqk21mv centr = 0.5e+00_wp_*(a+b) hlgth = 0.5e+00_wp_*(b-a) dhlgth = abs(hlgth) ! ! compute the 21-point kronrod approximation to ! the integral, and estimate the absolute error. ! resg = 0.0e+00_wp_ fc = f(centr,apar,np) resk = wgk(11)*fc resabs = abs(resk) do j=1,5 jtw = 2*j absc = hlgth*xgk(jtw) fval1 = f(centr-absc,apar,np) fval2 = f(centr+absc,apar,np) fv1(jtw) = fval1 fv2(jtw) = fval2 fsum = fval1+fval2 resg = resg+wg(j)*fsum resk = resk+wgk(jtw)*fsum resabs = resabs+wgk(jtw)*(abs(fval1)+abs(fval2)) end do do j = 1,5 jtwm1 = 2*j-1 absc = hlgth*xgk(jtwm1) fval1 = f(centr-absc,apar,np) fval2 = f(centr+absc,apar,np) fv1(jtwm1) = fval1 fv2(jtwm1) = fval2 fsum = fval1+fval2 resk = resk+wgk(jtwm1)*fsum resabs = resabs+wgk(jtwm1)*(abs(fval1)+abs(fval2)) end do reskh = resk*0.5e+00_wp_ resasc = wgk(11)*abs(fc-reskh) do j=1,10 resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)) end do result = resk*hlgth resabs = resabs*dhlgth resasc = resasc*dhlgth abserr = abs((resk-resg)*hlgth) if(resasc/=0.0e+00_wp_.and.abserr/=0.0e+00_wp_) & abserr = resasc*dmin1(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_) if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 & ((epmach*0.5e+02_wp_)*resabs,abserr) end subroutine dqk21mv subroutine dqagimv(f,bound,inf,apar,np,epsabs,epsrel,result,abserr,neval, & ier,limit,lenw,last,iwork,work) !***begin prologue dqagimv !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a3a1,h2a4a1 !***keywords automatic integrator, infinite intervals, ! general-purpose, transformation, extrapolation, ! globally adaptive !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & progr. div. -k.u.leuven !***purpose the routine calculates an approximation result to a given ! integral i = integral of f over (bound,+infinity) ! or i = integral of f over (-infinity,bound) ! or i = integral of f over (-infinity,+infinity) ! hopefully satisfying following claim for accuracy ! abs(i-result)<=max(epsabs,epsrel*abs(i)). !***description ! ! integration over infinite intervals ! standard fortran subroutine ! ! parameters ! on entry ! f - real(8) ! function subprogram defining the integrand ! function f(x,apar,np). the actual name for f needs to be ! declared e x t e r n a l in the driver program. ! ! bound - real(8) ! finite bound of integration range ! (has no meaning if interval is doubly-infinite) ! ! inf - integer ! indicating the kind of integration range involved ! inf = 1 corresponds to (bound,+infinity), ! inf = -1 to (-infinity,bound), ! inf = 2 to (-infinity,+infinity). ! ! apar - array of parameters of the integrand function f ! ! np - number of parameters. size of apar ! ! epsabs - real(8) ! absolute accuracy requested ! epsrel - real(8) ! relative accuracy requested ! if epsabs<=0 ! and epsrel0 abnormal termination of the routine. the ! estimates for result and error are less ! reliable. it is assumed that the requested ! accuracy has not been achieved. ! error messages ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the value of ! limit (and taking the according dimension ! adjustments into account). however, if ! this yields no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulties. if ! the position of a local difficulty can be ! determined (e.g. singularity, ! discontinuity within the interval) one ! will probably gain from splitting up the ! interval at this point and calling the ! integrator on the subranges. if possible, ! an appropriate special-purpose integrator ! should be used, which is designed for ! handling the type of difficulty involved. ! = 2 the occurrence of roundoff error is ! detected, which prevents the requested ! tolerance from being achieved. ! the error may be under-estimated. ! = 3 extremely bad integrand behaviour occurs ! at some points of the integration ! interval. ! = 4 the algorithm does not converge. ! roundoff error is detected in the ! extrapolation table. ! it is assumed that the requested tolerance ! cannot be achieved, and that the returned ! result is the best which can be obtained. ! = 5 the integral is probably divergent, or ! slowly convergent. it must be noted that ! divergence can occur with any other value ! of ier. ! = 6 the input is invalid, because ! (epsabs<=0 and ! epsrel=1. ! if limit<1, the routine will end with ier = 6. ! ! lenw - integer ! dimensioning parameter for work ! lenw must be at least limit*4. ! if lenw=1.and.lenw>=limit*4) then ! ! prepare call for dqagie. ! l1 = limit+1 l2 = limit+l1 l3 = limit+l2 ! call dqagiemv(f,bound,inf,apar,np,epsabs,epsrel,limit,result,abserr, & neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last) ! end if if(ier/=0) print*,'habnormal return from dqagi' end subroutine dqagimv subroutine dqagiemv(f,bound,inf,apar,np,epsabs,epsrel,limit,result,abserr, & neval,ier,alist,blist,rlist,elist,iord,last) !***begin prologue dqagiemv !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a3a1,h2a4a1 !***keywords automatic integrator, infinite intervals, ! general-purpose, transformation, extrapolation, ! globally adaptive !***author piessens,robert,appl. math & progr. div - k.u.leuven ! de doncker,elise,appl. math & progr. div - k.u.leuven !***purpose the routine calculates an approximation result to a given ! integral i = integral of f over (bound,+infinity) ! or i = integral of f over (-infinity,bound) ! or i = integral of f over (-infinity,+infinity), ! hopefully satisfying following claim for accuracy ! abs(i-result)<=max(epsabs,epsrel*abs(i)) !***description ! ! integration over infinite intervals ! standard fortran subroutine ! ! f - real(8) ! function subprogram defining the integrand ! function f(x,apar,np). the actual name for f needs to be ! declared e x t e r n a l in the driver program. ! ! bound - real(8) ! finite bound of integration range ! (has no meaning if interval is doubly-infinite) ! ! inf - real(8) ! indicating the kind of integration range involved ! inf = 1 corresponds to (bound,+infinity), ! inf = -1 to (-infinity,bound), ! inf = 2 to (-infinity,+infinity). ! ! apar - array of parameters of the integrand function f ! ! np - number of parameters. size of apar ! ! epsabs - real(8) ! absolute accuracy requested ! epsrel - real(8) ! relative accuracy requested ! if epsabs<=0 ! and epsrel=1 ! ! on return ! result - real(8) ! approximation to the integral ! ! abserr - real(8) ! estimate of the modulus of the absolute error, ! which should equal or exceed abs(i-result) ! ! neval - integer ! number of integrand evaluations ! ! ier - integer ! ier = 0 normal and reliable termination of the ! routine. it is assumed that the requested ! accuracy has been achieved. ! - ier>0 abnormal termination of the routine. the ! estimates for result and error are less ! reliable. it is assumed that the requested ! accuracy has not been achieved. ! error messages ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the value of ! limit (and taking the according dimension ! adjustments into account). however,if ! this yields no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulties. ! if the position of a local difficulty can ! be determined (e.g. singularity, ! discontinuity within the interval) one ! will probably gain from splitting up the ! interval at this point and calling the ! integrator on the subranges. if possible, ! an appropriate special-purpose integrator ! should be used, which is designed for ! handling the type of difficulty involved. ! = 2 the occurrence of roundoff error is ! detected, which prevents the requested ! tolerance from being achieved. ! the error may be under-estimated. ! = 3 extremely bad integrand behaviour occurs ! at some points of the integration ! interval. ! = 4 the algorithm does not converge. ! roundoff error is detected in the ! extrapolation table. ! it is assumed that the requested tolerance ! cannot be achieved, and that the returned ! result is the best which can be obtained. ! = 5 the integral is probably divergent, or ! slowly convergent. it must be noted that ! divergence can occur with any other value ! of ier. ! = 6 the input is invalid, because ! (epsabs<=0 and ! epsrelcomp_eps, uflow=>comp_tiny, & oflow=>comp_huge implicit none integer, intent(in) :: limit,inf,np real(wp_), intent(in) :: bound,epsabs,epsrel real(wp_), dimension(np), intent(in) :: apar real(wp_), intent(out) :: result,abserr integer, intent(out) :: ier,neval,last real(wp_), dimension(limit), intent(inout) :: alist,blist,elist,rlist integer, dimension(limit), intent(inout) :: iord real(wp_), external :: f real(wp_) :: abseps,area,area1,area12,area2,a1,a2,boun,b1,b2,correc, & abs,defabs,defab1,defab2,dmax1,dres,erlarg,erlast,errbnd,errmax, & error1,error2,erro12,errsum,ertest,resabs,reseps,small real(wp_) :: res3la(3),rlist2(52) integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, & ktmin,maxerr,nres,nrmax,numrl2 logical :: extrap,noext ! ! the dimension of rlist2 is determined by the value of ! limexp in subroutine dqelg. ! ! ! list of major variables ! ----------------------- ! ! alist - list of left end points of all subintervals ! considered up to now ! blist - list of right end points of all subintervals ! considered up to now ! rlist(i) - approximation to the integral over ! (alist(i),blist(i)) ! rlist2 - array of dimension at least (limexp+2), ! containing the part of the epsilon table ! wich is still needed for further computations ! elist(i) - error estimate applying to rlist(i) ! maxerr - pointer to the interval with largest error ! estimate ! errmax - elist(maxerr) ! erlast - error on the interval currently subdivided ! (before that subdivision has taken place) ! area - sum of the integrals over the subintervals ! errsum - sum of the errors over the subintervals ! errbnd - requested accuracy max(epsabs,epsrel* ! abs(result)) ! *****1 - variable for the left subinterval ! *****2 - variable for the right subinterval ! last - index for subdivision ! nres - number of calls to the extrapolation routine ! numrl2 - number of elements currently in rlist2. if an ! appropriate approximation to the compounded ! integral has been obtained, it is put in ! rlist2(numrl2) after numrl2 has been increased ! by one. ! small - length of the smallest interval considered up ! to now, multiplied by 1.5 ! erlarg - sum of the errors over the intervals larger ! than the smallest interval considered up to now ! extrap - logical variable denoting that the routine ! is attempting to perform extrapolation. i.e. ! before subdividing the smallest interval we ! try to decrease the value of erlarg. ! noext - logical variable denoting that extrapolation ! is no longer allowed (true-value) ! ! machine dependent constants ! --------------------------- ! ! epmach is the largest relative spacing. ! uflow is the smallest positive magnitude. ! oflow is the largest positive magnitude. ! !***first executable statement dqagie ! ! test on validity of parameters ! ----------------------------- ! ier = 0 neval = 0 last = 0 result = 0.0e+00_wp_ abserr = 0.0e+00_wp_ alist(1) = 0.0e+00_wp_ blist(1) = 0.1e+01_wp_ rlist(1) = 0.0e+00_wp_ elist(1) = 0.0e+00_wp_ iord(1) = 0 if(epsabs<=0.0e+00_wp_.and.epsrelerrbnd) ier = 2 if(limit==1) ier = 1 if(ier/=0.or.(abserr<=errbnd.and.abserr/=resabs).or. & abserr==0.0e+00_wp_) go to 130 ! ! initialization ! -------------- ! rlist2(1) = result errmax = abserr maxerr = 1 area = result errsum = abserr abserr = oflow nrmax = 1 nres = 0 ktmin = 0 numrl2 = 2 extrap = .false. noext = .false. ierro = 0 iroff1 = 0 iroff2 = 0 iroff3 = 0 ksgn = -1 if(dres>=(0.1e+01_wp_-0.5e+02_wp_*epmach)*defabs) ksgn = 1 ! ! main do-loop ! ------------ ! do last = 2,limit ! ! bisect the subinterval with nrmax-th largest error estimate. ! a1 = alist(maxerr) b1 = 0.5e+00_wp_*(alist(maxerr)+blist(maxerr)) a2 = b1 b2 = blist(maxerr) erlast = errmax call dqk15imv(f,boun,inf,a1,b1,apar,np,area1,error1,resabs,defab1) call dqk15imv(f,boun,inf,a2,b2,apar,np,area2,error2,resabs,defab2) ! ! improve previous approximations to integral ! and error and test for accuracy. ! area12 = area1+area2 erro12 = error1+error2 errsum = errsum+erro12-errmax area = area+area12-rlist(maxerr) if(defab1/=error1.and.defab2/=error2) then if(abs(rlist(maxerr)-area12)<=0.1e-04_wp_*abs(area12) & .and.erro12>=0.99e+00_wp_*errmax) then if(extrap) iroff2 = iroff2+1 if(.not.extrap) iroff1 = iroff1+1 end if if(last>10.and.erro12>errmax) iroff3 = iroff3+1 end if rlist(maxerr) = area1 rlist(last) = area2 errbnd = dmax1(epsabs,epsrel*abs(area)) ! ! test for roundoff error and eventually set error flag. ! if(iroff1+iroff2>=10.or.iroff3>=20) ier = 2 if(iroff2>=5) ierro = 3 ! ! set error flag in the case that the number of ! subintervals equals limit. ! if(last==limit) ier = 1 ! ! set error flag in the case of bad integrand behaviour ! at some points of the integration range. ! if(dmax1(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* & (abs(a2)+0.1e+04_wp_*uflow)) ier = 4 ! ! append the newly-created intervals to the list. ! if(error2<=error1) then alist(last) = a2 blist(maxerr) = b1 blist(last) = b2 elist(maxerr) = error1 elist(last) = error2 else alist(maxerr) = a2 alist(last) = a1 blist(last) = b1 rlist(maxerr) = area2 rlist(last) = area1 elist(maxerr) = error2 elist(last) = error1 end if ! ! call subroutine dqpsrt to maintain the descending ordering ! in the list of error estimates and select the subinterval ! with nrmax-th largest error estimate (to be bisected next). ! call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) if(errsum<=errbnd) go to 115 if(ier/=0) exit if(last==2) then small = 0.375e+00_wp_ erlarg = errsum ertest = errbnd rlist2(2) = area cycle end if if(noext) cycle erlarg = erlarg-erlast if(abs(b1-a1)>small) erlarg = erlarg+erro12 if(.not.extrap) then ! ! test whether the interval to be bisected next is the ! smallest interval. ! if(abs(blist(maxerr)-alist(maxerr))>small) cycle extrap = .true. nrmax = 2 end if if(ierro/=3.and.erlarg>ertest) then ! ! the smallest interval has the largest error. ! before bisecting decrease the sum of the errors over the ! larger intervals (erlarg) and perform extrapolation. ! id = nrmax jupbnd = last if(last>(2+limit/2)) jupbnd = limit+3-last do k = id,jupbnd maxerr = iord(nrmax) errmax = elist(maxerr) if(abs(blist(maxerr)-alist(maxerr))>small) go to 90 nrmax = nrmax+1 end do end if ! ! perform extrapolation. ! numrl2 = numrl2+1 rlist2(numrl2) = area call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres) ktmin = ktmin+1 if(ktmin>5.and.abserr<0.1e-02_wp_*errsum) ier = 5 if(absepserrsum)go to 115 if(area==0.0e+00_wp_) go to 130 go to 110 105 continue if(abserr/abs(result)>errsum/abs(area)) go to 115 ! ! test on divergence ! 110 continue if(ksgn==(-1).and.dmax1(abs(result),abs(area))<= & defabs*0.1e-01_wp_) go to 130 if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ & .or.errsum>abs(area)) ier = 6 go to 130 ! ! compute global integral sum. ! 115 continue result = 0.0e+00_wp_ do k = 1,last result = result+rlist(k) end do abserr = errsum 130 continue neval = 30*last-15 if(inf==2) neval = 2*neval if(ier>2) ier=ier-1 end subroutine dqagiemv subroutine dqk15imv(f,boun,inf,a,b,apar,np,result,abserr,resabs,resasc) !***begin prologue dqk15imv !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a3a2,h2a4a2 !***keywords 15-point transformed gauss-kronrod rules !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & progr. div. - k.u.leuven !***purpose the original (infinite integration range is mapped ! onto the interval (0,1) and (a,b) is a part of (0,1). ! it is the purpose to compute ! i = integral of transformed integrand over (a,b), ! j = integral of abs(transformed integrand) over (a,b). !***description ! ! integration rule ! standard fortran subroutine ! real(8) version ! ! parameters ! on entry ! f - real(8) ! fuction subprogram defining the integrand ! function f(x,apar,np). the actual name for f needs to be ! declared e x t e r n a l in the calling program. ! ! boun - real(8) ! finite bound of original integration ! range (set to zero if inf = +2) ! ! inf - integer ! if inf = -1, the original interval is ! (-infinity,bound), ! if inf = +1, the original interval is ! (bound,+infinity), ! if inf = +2, the original interval is ! (-infinity,+infinity) and ! the integral is computed as the sum of two ! integrals, one over (-infinity,0) and one over ! (0,+infinity). ! ! a - real(8) ! lower limit for integration over subrange ! of (0,1) ! ! b - real(8) ! upper limit for integration over subrange ! of (0,1) ! ! apar - array of parameters of the integrand function f ! ! np - number of parameters. size of apar ! ! on return ! result - real(8) ! approximation to the integral i ! result is computed by applying the 15-point ! kronrod rule(resk) obtained by optimal addition ! of abscissae to the 7-point gauss rule(resg). ! ! abserr - real(8) ! estimate of the modulus of the absolute error, ! which should equal or exceed abs(i-result) ! ! resabs - real(8) ! approximation to the integral j ! ! resasc - real(8) ! approximation to the integral of ! abs((transformed integrand)-i/(b-a)) over (a,b) ! !***references (none) !***routines called (none) !***end prologue dqk15imv ! use const_and_precisions, only : epmach=>comp_eps, uflow=>comp_tiny implicit none real(wp_), intent(in) :: a,b,boun integer, intent(in) :: inf,np real(wp_), dimension(np), intent(in) :: apar real(wp_), intent(out) :: result,abserr,resabs,resasc real(wp_), external :: f real(wp_) :: absc,absc1,absc2,centr,abs,dinf,dmax1,dmin1,fc,fsum, & fval1,fval2,hlgth,resg,resk,reskh,tabsc1,tabsc2 real(wp_), dimension(7) :: fv1,fv2 integer :: j ! ! the abscissae and weights are supplied for the interval ! (-1,1). because of symmetry only the positive abscissae and ! their corresponding weights are given. ! ! xgk - abscissae of the 15-point kronrod rule ! xgk(2), xgk(4), ... abscissae of the 7-point ! gauss rule ! xgk(1), xgk(3), ... abscissae which are optimally ! added to the 7-point gauss rule ! ! wgk - weights of the 15-point kronrod rule ! ! wg - weights of the 7-point gauss rule, corresponding ! to the abscissae xgk(2), xgk(4), ... ! wg(1), wg(3), ... are set to zero. ! real(wp_), dimension(8), parameter :: & wg = (/ 0.000000000000000000000000000000000_wp_, & 0.129484966168869693270611432679082_wp_, & 0.000000000000000000000000000000000_wp_, & 0.279705391489276667901467771423780_wp_, & 0.000000000000000000000000000000000_wp_, & 0.381830050505118944950369775488975_wp_, & 0.000000000000000000000000000000000_wp_, & 0.417959183673469387755102040816327_wp_ /), & xgk = (/ 0.991455371120812639206854697526329_wp_, & 0.949107912342758524526189684047851_wp_, & 0.864864423359769072789712788640926_wp_, & 0.741531185599394439863864773280788_wp_, & 0.586087235467691130294144838258730_wp_, & 0.405845151377397166906606412076961_wp_, & 0.207784955007898467600689403773245_wp_, & 0.000000000000000000000000000000000_wp_ /), & wgk = (/ 0.022935322010529224963732008058970_wp_, & 0.063092092629978553290700663189204_wp_, & 0.104790010322250183839876322541518_wp_, & 0.140653259715525918745189590510238_wp_, & 0.169004726639267902826583426598550_wp_, & 0.190350578064785409913256402421014_wp_, & 0.204432940075298892414161999234649_wp_, & 0.209482141084727828012999174891714_wp_ /) ! ! ! list of major variables ! ----------------------- ! ! centr - mid point of the interval ! hlgth - half-length of the interval ! absc* - abscissa ! tabsc* - transformed abscissa ! fval* - function value ! resg - result of the 7-point gauss formula ! resk - result of the 15-point kronrod formula ! reskh - approximation to the mean value of the transformed ! integrand over (a,b), i.e. to i/(b-a) ! ! machine dependent constants ! --------------------------- ! ! epmach is the largest relative spacing. ! uflow is the smallest positive magnitude. ! !***first executable statement dqk15imv dinf = min0(1,inf) ! centr = 0.5e+00_wp_*(a+b) hlgth = 0.5e+00_wp_*(b-a) tabsc1 = boun+dinf*(0.1e+01_wp_-centr)/centr fval1 = f(tabsc1,apar,np) if(inf==2) fval1 = fval1+f(-tabsc1,apar,np) fc = (fval1/centr)/centr ! ! compute the 15-point kronrod approximation to ! the integral, and estimate the error. ! resg = wg(8)*fc resk = wgk(8)*fc resabs = abs(resk) do j=1,7 absc = hlgth*xgk(j) absc1 = centr-absc absc2 = centr+absc tabsc1 = boun+dinf*(0.1e+01_wp_-absc1)/absc1 tabsc2 = boun+dinf*(0.1e+01_wp_-absc2)/absc2 fval1 = f(tabsc1,apar,np) fval2 = f(tabsc2,apar,np) if(inf==2) fval1 = fval1+f(-tabsc1,apar,np) if(inf==2) fval2 = fval2+f(-tabsc2,apar,np) fval1 = (fval1/absc1)/absc1 fval2 = (fval2/absc2)/absc2 fv1(j) = fval1 fv2(j) = fval2 fsum = fval1+fval2 resg = resg+wg(j)*fsum resk = resk+wgk(j)*fsum resabs = resabs+wgk(j)*(abs(fval1)+abs(fval2)) end do reskh = resk*0.5e+00_wp_ resasc = wgk(8)*abs(fc-reskh) do j=1,7 resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)) end do result = resk*hlgth resasc = resasc*hlgth resabs = resabs*hlgth abserr = abs((resk-resg)*hlgth) if(resasc/=0.0e+00_wp_.and.abserr/=0._wp_) abserr = resasc* & dmin1(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_) if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 & ((epmach*0.5e+02_wp_)*resabs,abserr) end subroutine dqk15imv end module quadpack