subroutine aaaa !*****************************************************************************80 ! !! AAAA is a dummy subroutine with QUADPACK documentation in its comments. ! ! 1. introduction ! ! quadpack is a fortran subroutine package for the numerical ! computation of definite one-dimensional integrals. it originated ! from a joint project of r. piessens and e. de doncker (appl. ! math. and progr. div.- k.u.leuven, belgium), c. ueberhuber (inst. ! fuer math.- techn.u.wien, austria), and d. kahaner (nation. bur. ! of standards- washington d.c., u.s.a.). ! ! 2. survey ! ! - qags : is an integrator based on globally adaptive interval ! subdivision in connection with extrapolation (de doncker, ! 1978) by the epsilon algorithm (wynn, 1956). ! ! - qagp : serves the same purposes as qags, but also allows ! for eventual user-supplied information, i.e. the ! abscissae of internal singularities, discontinuities ! and other difficulties of the integrand function. ! the algorithm is a modification of that in qags. ! ! - qagi : handles integration over infinite intervals. the ! infinite range is mapped onto a finite interval and ! then the same strategy as in qags is applied. ! ! - qawo : is a routine for the integration of cos(omega*x)*f(x) ! or sin(omega*x)*f(x) over a finite interval (a,b). ! omega is is specified by the user ! the rule evaluation component is based on the ! modified clenshaw-curtis technique. ! an adaptive subdivision scheme is used connected with ! an extrapolation procedure, which is a modification ! of that in qags and provides the possibility to deal ! even with singularities in f. ! ! - qawf : calculates the fourier cosine or fourier sine ! transform of f(x), for user-supplied interval (a, ! infinity), omega, and f. the procedure of qawo is ! used on successive finite intervals, and convergence ! acceleration by means of the epsilon algorithm (wynn, ! 1956) is applied to the series of the integral ! contributions. ! ! - qaws : integrates w(x)*f(x) over (a,b) with a < b finite, ! and w(x) = ((x-a)**alfa)*((b-x)**beta)*v(x) ! where v(x) = 1 or log(x-a) or log(b-x) ! or log(x-a)*log(b-x) ! and -1 < alfa, -1 < beta. ! the user specifies a, b, alfa, beta and the type of ! the function v. ! a globally adaptive subdivision strategy is applied, ! with modified clenshaw-curtis integration on the ! subintervals which contain a or b. ! ! - qawc : computes the cauchy principal value of f(x)/(x-c) ! over a finite interval (a,b) and for ! user-determined c. ! the strategy is globally adaptive, and modified ! clenshaw-curtis integration is used on the subranges ! which contain the point x = c. ! ! each of the routines above also has a "more detailed" version ! with a name ending in e, as qage. these provide more ! information and control than the easier versions. ! ! ! the preceeding routines are all automatic. that is, the user ! inputs his problem and an error tolerance. the routine ! attempts to perform the integration to within the requested ! absolute or relative error. ! there are, in addition, a number of non-automatic integrators. ! these are most useful when the problem is such that the ! user knows that a fixed rule will provide the accuracy ! required. typically they return an error estimate but make ! no attempt to satisfy any particular input error request. ! ! qk15 ! qk21 ! qk31 ! qk41 ! qk51 ! qk61 ! estimate the integral on [a,b] using 15, 21,..., 61 ! point rule and return an error estimate. ! qk15i 15 point rule for (semi)infinite interval. ! qk15w 15 point rule for special singular weight functions. ! qc25c 25 point rule for cauchy principal values ! qc25o 25 point rule for sin/cos integrand. ! qmomo integrates k-th degree chebychev polynomial times ! function with various explicit singularities. ! ! 3. guidelines for the use of quadpack ! ! here it is not our purpose to investigate the question when ! automatic quadrature should be used. we shall rather attempt ! to help the user who already made the decision to use quadpack, ! with selecting an appropriate routine or a combination of ! several routines for handling his problem. ! ! for both quadrature over finite and over infinite intervals, ! one of the first questions to be answered by the user is ! related to the amount of computer time he wants to spend, ! versus his -own- time which would be needed, for example, for ! manual subdivision of the interval or other analytic ! manipulations. ! ! (1) the user may not care about computer time, or not be ! willing to do any analysis of the problem. especially when ! only one or a few integrals must be calculated, this attitude ! can be perfectly reasonable. in this case it is clear that ! either the most sophisticated of the routines for finite ! intervals, qags, must be used, or its analogue for infinite ! intervals, qagi. these routines are able to cope with ! rather difficult, even with improper integrals. ! this way of proceeding may be expensive. but the integrator ! is supposed to give you an answer in return, with additional ! information in the case of a failure, through its error ! estimate and flag. yet it must be stressed that the programs ! cannot be totally reliable. ! ! (2) the user may want to examine the integrand function. ! if bad local difficulties occur, such as a discontinuity, a ! singularity, derivative singularity or high peak at one or ! more points within the interval, the first advice is to ! split up the interval at these points. the integrand must ! then be examinated over each of the subintervals separately, ! so that a suitable integrator can be selected for each of ! them. if this yields problems involving relative accuracies ! to be imposed on -finite- subintervals, one can make use of ! qagp, which must be provided with the positions of the local ! difficulties. however, if strong singularities are present ! and a high accuracy is requested, application of qags on the ! subintervals may yield a better result. ! ! for quadrature over finite intervals we thus dispose of qags ! and ! - qng for well-behaved integrands, ! - qag for functions with an oscillating behavior of a non ! specific type, ! - qawo for functions, eventually singular, containing a ! factor cos(omega*x) or sin(omega*x) where omega is known, ! - qaws for integrands with algebraico-logarithmic end point ! singularities of known type, ! - qawc for cauchy principal values. ! ! remark ! ! on return, the work arrays in the argument lists of the ! adaptive integrators contain information about the interval ! subdivision process and hence about the integrand behavior: ! the end points of the subintervals, the local integral ! contributions and error estimates, and eventually other ! characteristics. for this reason, and because of its simple ! globally adaptive nature, the routine qag in particular is ! well-suited for integrand examination. difficult spots can ! be located by investigating the error estimates on the ! subintervals. ! ! for infinite intervals we provide only one general-purpose ! routine, qagi. it is based on the qags algorithm applied ! after a transformation of the original interval into (0,1). ! yet it may eventuate that another type of transformation is ! more appropriate, or one might prefer to break up the ! original interval and use qagi only on the infinite part ! and so on. these kinds of actions suggest a combined use of ! different quadpack integrators. note that, when the only ! difficulty is an integrand singularity at the finite ! integration limit, it will in general not be necessary to ! break up the interval, as qagi deals with several types of ! singularity at the boundary point of the integration range. ! it also handles slowly convergent improper integrals, on ! the condition that the integrand does not oscillate over ! the entire infinite interval. if it does we would advise ! to sum succeeding positive and negative contributions to ! the integral -e.g. integrate between the zeros- with one ! or more of the finite-range integrators, and apply ! convergence acceleration eventually by means of quadpack ! subroutine qelg which implements the epsilon algorithm. ! such quadrature problems include the fourier transform as ! a special case. yet for the latter we have an automatic ! integrator available, qawf. ! return end subroutine qag ( f, a, b, epsabs, epsrel, key, result, abserr, neval, ier ) !*****************************************************************************80 ! !! QAG approximates an integral over a finite interval. ! ! Discussion: ! ! The routine calculates an approximation RESULT to a definite integral ! I = integral of F over (A,B), ! hopefully satisfying ! || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ). ! ! QAG is a simple globally adaptive integrator using the strategy of ! Aind (Piessens, 1973). It is possible to choose between 6 pairs of ! Gauss-Kronrod quadrature formulae for the rule evaluation component. ! The pairs of high degree of precision are suitable for handling ! integration difficulties due to a strongly oscillating integrand. ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, B, the limits of integration. ! ! Input, real EPSABS, EPSREL, the absolute and relative accuracy requested. ! ! Input, integer KEY, chooses the order of the local integration rule: ! 1, 7 Gauss points, 15 Gauss-Kronrod points, ! 2, 10 Gauss points, 21 Gauss-Kronrod points, ! 3, 15 Gauss points, 31 Gauss-Kronrod points, ! 4, 20 Gauss points, 41 Gauss-Kronrod points, ! 5, 25 Gauss points, 51 Gauss-Kronrod points, ! 6, 30 Gauss points, 61 Gauss-Kronrod points. ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! Output, integer IER, return code. ! 0, normal and reliable termination of the routine. It is assumed that the ! requested accuracy has been achieved. ! 1, maximum number of subdivisions allowed has been achieved. One can ! allow more subdivisions by increasing the value of LIMIT in QAG. ! However, if this yields no improvement it is advised to analyze the ! integrand to determine the integration difficulties. If the position ! of a local difficulty can be determined, such as a singularity or ! 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. ! 3, extremely bad integrand behavior occurs at some points of the ! integration interval. ! 6, the input is invalid, because EPSABS < 0 and EPSREL < 0. ! ! Local parameters: ! ! LIMIT is the maximum number of subintervals allowed in ! the subdivision process of QAGE. ! implicit none integer, parameter :: limit = 500 real a real abserr real alist(limit) real b real blist(limit) real elist(limit) real epsabs real epsrel real, external :: f integer ier integer iord(limit) integer key integer last integer neval real result real rlist(limit) call qage ( f, a, b, epsabs, epsrel, key, limit, result, abserr, neval, & ier, alist, blist, rlist, elist, iord, last ) return end subroutine qage ( f, a, b, epsabs, epsrel, key, limit, result, abserr, neval, & ier, alist, blist, rlist, elist, iord, last ) !*****************************************************************************80 ! !! QAGE estimates a definite integral. ! ! Discussion: ! ! The routine calculates an approximation RESULT to a definite integral ! I = integral of F over (A,B), ! hopefully satisfying ! || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ). ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, B, the limits of integration. ! ! Input, real EPSABS, EPSREL, the absolute and relative accuracy requested. ! ! Input, integer KEY, chooses the order of the local integration rule: ! 1, 7 Gauss points, 15 Gauss-Kronrod points, ! 2, 10 Gauss points, 21 Gauss-Kronrod points, ! 3, 15 Gauss points, 31 Gauss-Kronrod points, ! 4, 20 Gauss points, 41 Gauss-Kronrod points, ! 5, 25 Gauss points, 51 Gauss-Kronrod points, ! 6, 30 Gauss points, 61 Gauss-Kronrod points. ! ! Input, integer LIMIT, the maximum number of subintervals that ! can be used. ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! Output, integer IER, return code. ! 0, normal and reliable termination of the routine. It is assumed that the ! requested accuracy has been achieved. ! 1, maximum number of subdivisions allowed has been achieved. One can ! allow more subdivisions by increasing the value of LIMIT in QAG. ! However, if this yields no improvement it is advised to analyze the ! integrand to determine the integration difficulties. If the position ! of a local difficulty can be determined, such as a singularity or ! 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. ! 3, extremely bad integrand behavior occurs at some points of the ! integration interval. ! 6, the input is invalid, because EPSABS < 0 and EPSREL < 0. ! ! Workspace, real ALIST(LIMIT), BLIST(LIMIT), contains in entries 1 ! through LAST the left and right ends of the partition subintervals. ! ! Workspace, real RLIST(LIMIT), contains in entries 1 through LAST ! the integral approximations on the subintervals. ! ! Workspace, real ELIST(LIMIT), contains in entries 1 through LAST ! the absolute error estimates on the subintervals. ! ! Output, integer IORD(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 <= (limit/2+2), and k = limit+1-last otherwise. ! ! Output, integer LAST, the number of subintervals actually produced ! in the subdivision process. ! ! Local parameters: ! ! 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 ! elist(i) - error estimate applying to rlist(i) ! maxerr - pointer to the interval with largest error estimate ! errmax - elist(maxerr) ! 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 ! implicit none integer limit real a real abserr real alist(limit) real area real area1 real area12 real area2 real a1 real a2 real b real blist(limit) real b1 real b2 real c real defabs real defab1 real defab2 real elist(limit) real epsabs real epsrel real errbnd real errmax real error1 real error2 real erro12 real errsum real, external :: f integer ier integer iord(limit) integer iroff1 integer iroff2 integer key integer keyf integer last integer maxerr integer neval integer nrmax real resabs real result real rlist(limit) ! ! Test on validity of parameters. ! ier = 0 neval = 0 last = 0 result = 0.0e+00 abserr = 0.0e+00 alist(1) = a blist(1) = b rlist(1) = 0.0e+00 elist(1) = 0.0e+00 iord(1) = 0 if ( epsabs < 0.0e+00 .and. epsrel < 0.0e+00 ) then ier = 6 return end if ! ! First approximation to the integral. ! keyf = key keyf = max ( keyf, 1 ) keyf = min ( keyf, 6 ) c = keyf neval = 0 if ( keyf == 1 ) then call qk15 ( f, a, b, result, abserr, defabs, resabs ) else if ( keyf == 2 ) then call qk21 ( f, a, b, result, abserr, defabs, resabs ) else if ( keyf == 3 ) then call qk31 ( f, a, b, result, abserr, defabs, resabs ) else if ( keyf == 4 ) then call qk41 ( f, a, b, result, abserr, defabs, resabs ) else if ( keyf == 5 ) then call qk51 ( f, a, b, result, abserr, defabs, resabs ) else if ( keyf == 6 ) then call qk61 ( f, a, b, result, abserr, defabs, resabs ) end if last = 1 rlist(1) = result elist(1) = abserr iord(1) = 1 ! ! Test on accuracy. ! errbnd = max ( epsabs, epsrel * abs ( result ) ) if ( abserr <= 5.0e+01 * epsilon ( defabs ) * defabs .and. & errbnd < abserr ) then ier = 2 end if if ( limit == 1 ) then ier = 1 end if if ( ier /= 0 .or. & ( abserr <= errbnd .and. abserr /= resabs ) .or. & abserr == 0.0e+00 ) then if ( keyf /= 1 ) then neval = (10*keyf+1) * (2*neval+1) else neval = 30 * neval + 15 end if return end if ! ! Initialization. ! errmax = abserr maxerr = 1 area = result errsum = abserr nrmax = 1 iroff1 = 0 iroff2 = 0 do last = 2, limit ! ! Bisect the subinterval with the largest error estimate. ! a1 = alist(maxerr) b1 = 0.5E+00 * ( alist(maxerr) + blist(maxerr) ) a2 = b1 b2 = blist(maxerr) if ( keyf == 1 ) then call qk15 ( f, a1, b1, area1, error1, resabs, defab1 ) else if ( keyf == 2 ) then call qk21 ( f, a1, b1, area1, error1, resabs, defab1 ) else if ( keyf == 3 ) then call qk31 ( f, a1, b1, area1, error1, resabs, defab1 ) else if ( keyf == 4 ) then call qk41 ( f, a1, b1, area1, error1, resabs, defab1) else if ( keyf == 5 ) then call qk51 ( f, a1, b1, area1, error1, resabs, defab1 ) else if ( keyf == 6 ) then call qk61 ( f, a1, b1, area1, error1, resabs, defab1 ) end if if ( keyf == 1 ) then call qk15 ( f, a2, b2, area2, error2, resabs, defab2 ) else if ( keyf == 2 ) then call qk21 ( f, a2, b2, area2, error2, resabs, defab2 ) else if ( keyf == 3 ) then call qk31 ( f, a2, b2, area2, error2, resabs, defab2 ) else if ( keyf == 4 ) then call qk41 ( f, a2, b2, area2, error2, resabs, defab2 ) else if ( keyf == 5 ) then call qk51 ( f, a2, b2, area2, error2, resabs, defab2 ) else if ( keyf == 6 ) then call qk61 ( f, a2, b2, area2, error2, resabs, defab2 ) end if ! ! Improve previous approximations to integral and error and ! test for accuracy. ! neval = neval + 1 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 ) <= 1.0e-05 * abs ( area12 ) & .and. 9.9e-01 * errmax <= erro12 ) then iroff1 = iroff1 + 1 end if if ( 10 < last .and. errmax < erro12 ) then iroff2 = iroff2 + 1 end if end if rlist(maxerr) = area1 rlist(last) = area2 errbnd = max ( epsabs, epsrel * abs ( area ) ) ! ! Test for roundoff error and eventually set error flag. ! if ( errbnd < errsum ) then if ( 6 <= iroff1 .or. 20 <= iroff2 ) then ier = 2 end if ! ! Set error flag in the case that the number of subintervals ! equals limit. ! if ( last == limit ) then ier = 1 end if ! ! Set error flag in the case of bad integrand behavior ! at a point of the integration range. ! if ( max ( abs ( a1 ), abs ( b2 ) ) <= ( 1.0e+00 + c * 1.0e+03 * & epsilon ( a1 ) ) * ( abs ( a2 ) + 1.0e+04 * tiny ( a2 ) ) ) then ier = 3 end if end if ! ! 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 QSORT to maintain the descending ordering ! in the list of error estimates and select the subinterval ! with the largest error estimate (to be bisected next). ! call qsort ( limit, last, maxerr, errmax, elist, iord, nrmax ) if ( ier /= 0 .or. errsum <= errbnd ) then exit end if end do ! ! Compute final result. ! result = sum ( rlist(1:last) ) abserr = errsum if ( keyf /= 1 ) then neval = ( 10 * keyf + 1 ) * ( 2 * neval + 1 ) else neval = 30 * neval + 15 end if return end subroutine qagi ( f, bound, inf, epsabs, epsrel, result, abserr, neval, ier ) !*****************************************************************************80 ! !! QAGI estimates an integral over a semi-infinite or infinite interval. ! ! Discussion: ! ! The routine calculates an approximation RESULT to a definite integral ! I = integral of F over (A, +Infinity), ! or ! I = integral of F over (-Infinity,A) ! or ! I = integral of F over (-Infinity,+Infinity), ! hopefully satisfying ! || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ). ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real BOUND, the value of the finite endpoint of the integration ! range, if any, that is, if INF is 1 or -1. ! ! Input, integer INF, indicates the type of integration range. ! 1: ( BOUND, +Infinity), ! -1: ( -Infinity, BOUND), ! 2: ( -Infinity, +Infinity). ! ! Input, real EPSABS, EPSREL, the absolute and relative accuracy requested. ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! Output, integer IER, error indicator. ! 0, normal and reliable termination of the routine. It is assumed that ! the requested accuracy has been achieved. ! > 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. ! 1, maximum number of subdivisions allowed has been achieved. One can ! allow more subdivisions by increasing the data value of LIMIT in QAGI ! (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 behavior 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 INF /= 1 and INF /= -1 and INF /= 2, or ! epsabs < 0 and epsrel < 0. result, abserr, neval are set to zero. ! ! Local parameters: ! ! the dimension of rlist2 is determined by the value of ! limexp in QEXTR. ! ! 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 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) ! implicit none integer, parameter :: limit = 500 real abseps real abserr real alist(limit) real area real area1 real area12 real area2 real a1 real a2 real blist(limit) real boun real bound real b1 real b2 real correc real defabs real defab1 real defab2 real dres real elist(limit) real epsabs real epsrel real erlarg real erlast real errbnd real errmax real error1 real error2 real erro12 real errsum real ertest logical extrap real, external :: f integer id integer ier integer ierro integer inf integer iord(limit) integer iroff1 integer iroff2 integer iroff3 integer jupbnd integer k integer ksgn integer ktmin integer last integer maxerr integer neval logical noext integer nres integer nrmax integer numrl2 real resabs real reseps real result real res3la(3) real rlist(limit) real rlist2(52) real small ! ! Test on validity of parameters. ! ier = 0 neval = 0 last = 0 result = 0.0e+00 abserr = 0.0e+00 alist(1) = 0.0e+00 blist(1) = 1.0e+00 rlist(1) = 0.0e+00 elist(1) = 0.0e+00 iord(1) = 0 if ( epsabs < 0.0e+00 .and. epsrel < 0.0e+00 ) then ier = 6 return end if ! ! First approximation to the integral. ! ! Determine the interval to be mapped onto (0,1). ! If INF = 2 the integral is computed as i = i1+i2, where ! i1 = integral of f over (-infinity,0), ! i2 = integral of f over (0,+infinity). ! if ( inf == 2 ) then boun = 0.0e+00 else boun = bound end if call qk15i ( f, boun, inf, 0.0e+00, 1.0e+00, result, abserr, defabs, resabs ) ! ! Test on accuracy. ! last = 1 rlist(1) = result elist(1) = abserr iord(1) = 1 dres = abs ( result ) errbnd = max ( epsabs, epsrel * dres ) if ( abserr <= 100.0E+00 * epsilon ( defabs ) * defabs .and. & errbnd < abserr ) then ier = 2 end if if ( limit == 1 ) then ier = 1 end if if ( ier /= 0 .or. (abserr <= errbnd .and. abserr /= resabs ) .or. & abserr == 0.0e+00 ) go to 130 ! ! Initialization. ! rlist2(1) = result errmax = abserr maxerr = 1 area = result errsum = abserr abserr = huge ( abserr ) nrmax = 1 nres = 0 ktmin = 0 numrl2 = 2 extrap = .false. noext = .false. ierro = 0 iroff1 = 0 iroff2 = 0 iroff3 = 0 if ( ( 1.0e+00 - 5.0e+01 * epsilon ( defabs ) ) * defabs <= dres ) then ksgn = 1 else ksgn = -1 end if do last = 2, limit ! ! Bisect the subinterval with nrmax-th largest error estimate. ! a1 = alist(maxerr) b1 = 5.0e-01 * ( alist(maxerr) + blist(maxerr) ) a2 = b1 b2 = blist(maxerr) erlast = errmax call qk15i ( f, boun, inf, a1, b1, area1, error1, resabs, defab1 ) call qk15i ( 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 ) <= 1.0e-05 * abs ( area12 ) & .and. 9.9e-01 * errmax <= erro12 ) then if ( extrap ) then iroff2 = iroff2 + 1 end if if ( .not. extrap ) then iroff1 = iroff1 + 1 end if end if if ( 10 < last .and. errmax < erro12 ) then iroff3 = iroff3 + 1 end if end if rlist(maxerr) = area1 rlist(last) = area2 errbnd = max ( epsabs, epsrel * abs ( area ) ) ! ! Test for roundoff error and eventually set error flag. ! if ( 10 <= iroff1 + iroff2 .or. 20 <= iroff3 ) then ier = 2 end if if ( 5 <= iroff2 ) then ierro = 3 end if ! ! Set error flag in the case that the number of subintervals equals LIMIT. ! if ( last == limit ) then ier = 1 end if ! ! Set error flag in the case of bad integrand behavior ! at some points of the integration range. ! if ( max ( abs(a1), abs(b2) ) <= (1.0e+00 + 1.0e+03 * epsilon ( a1 ) ) * & ( abs(a2) + 1.0e+03 * tiny ( a2 ) )) then ier = 4 end if ! ! 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 QSORT 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 qsort ( limit, last, maxerr, errmax, elist, iord, nrmax ) if ( errsum <= errbnd ) go to 115 if ( ier /= 0 ) then exit end if if ( last == 2 ) then small = 3.75e-01 erlarg = errsum ertest = errbnd rlist2(2) = area cycle end if if ( noext ) then cycle end if erlarg = erlarg - erlast if ( small < abs ( b1 - a1 ) ) then erlarg = erlarg + erro12 end if ! ! Test whether the interval to be bisected next is the ! smallest interval. ! if ( .not. extrap ) then if ( small < abs ( blist(maxerr) - alist(maxerr) ) ) then cycle end if extrap = .true. nrmax = 2 end if if ( ierro == 3 .or. erlarg <= ertest ) then go to 60 end if ! ! 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 ( (2+limit/2) < last ) then jupbnd = limit + 3 - last end if do k = id, jupbnd maxerr = iord(nrmax) errmax = elist(maxerr) if ( small < abs ( blist(maxerr) - alist(maxerr) ) ) then go to 90 end if nrmax = nrmax + 1 end do ! ! Extrapolate. ! 60 continue numrl2 = numrl2 + 1 rlist2(numrl2) = area call qextr ( numrl2, rlist2, reseps, abseps, res3la, nres ) ktmin = ktmin+1 if ( 5 < ktmin .and. abserr < 1.0e-03 * errsum ) then ier = 5 end if if ( abseps < abserr ) then ktmin = 0 abserr = abseps result = reseps correc = erlarg ertest = max ( epsabs, epsrel * abs(reseps) ) if ( abserr <= ertest ) then exit end if end if ! ! Prepare bisection of the smallest interval. ! if ( numrl2 == 1 ) then noext = .true. end if if ( ier == 5 ) then exit end if maxerr = iord(1) errmax = elist(maxerr) nrmax = 1 extrap = .false. small = small * 5.0e-01 erlarg = errsum 90 continue end do ! ! Set final result and error estimate. ! if ( abserr == huge ( abserr ) ) then go to 115 end if if ( ( ier + ierro ) == 0 ) then go to 110 end if if ( ierro == 3 ) then abserr = abserr + correc end if if ( ier == 0 ) then ier = 3 end if if ( result /= 0.0e+00 .and. area /= 0.0e+00) then go to 105 end if if ( errsum < abserr ) then go to 115 end if if ( area == 0.0e+00 ) then go to 130 end if go to 110 105 continue if ( errsum / abs ( area ) < abserr / abs ( result ) ) then go to 115 end if ! ! Test on divergence ! 110 continue if ( ksgn == (-1) .and. & max ( abs(result), abs(area) ) <= defabs * 1.0e-02) go to 130 if ( 1.0e-02 > (result/area) .or. & (result/area) > 1.0e+02 .or. & errsum > abs(area)) then ier = 6 end if go to 130 ! ! Compute global integral sum. ! 115 continue result = sum ( rlist(1:last) ) abserr = errsum 130 continue neval = 30 * last - 15 if ( inf == 2 ) then neval = 2 * neval end if if ( 2 < ier ) then ier = ier - 1 end if return end subroutine qagp ( f, a, b, npts2, points, epsabs, epsrel, result, abserr, & neval, ier ) !*****************************************************************************80 ! !! QAGP computes a definite integral. ! ! Discussion: ! ! The routine calculates an approximation RESULT to a definite integral ! I = integral of F over (A,B), ! hopefully satisfying ! || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ). ! ! Interior break points of the integration interval, ! where local difficulties of the integrand may occur, such as ! singularities or discontinuities, are provided by the user. ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, B, the limits of integration. ! ! Input, integer NPTS2, the number of user-supplied break points within ! the integration range, plus 2. NPTS2 must be at least 2. ! ! Input/output, real POINTS(NPTS2), contains the user provided interior ! breakpoints in entries 1 through NPTS2-2. If these points are not ! in ascending order on input, they will be sorted. ! ! Input, real EPSABS, EPSREL, the absolute and relative accuracy requested. ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! Output, integer IER, return flag. ! 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 integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the data value ! of limit in qagp(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 behavior 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 > 0. ! = 6 the input is invalid because ! npts2 < 2 or ! break points are specified outside ! the integration range or ! epsabs < 0 and epsrel < 0, ! or limit < npts2. ! result, abserr, neval are set to zero. ! ! Local parameters: ! ! the dimension of rlist2 is determined by the value of ! limexp in QEXTR (rlist2 should be of dimension ! (limexp+2) at least). ! ! 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 ! 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) ! implicit none integer, parameter :: limit = 500 real a real abseps real abserr real alist(limit) real area real area1 real area12 real area2 real a1 real a2 real b real blist(limit) real b1 real b2 real correc real defabs real defab1 real defab2 real dres real elist(limit) real epsabs real epsrel real erlarg real erlast real errbnd real errmax real error1 real erro12 real error2 real errsum real ertest logical extrap real, external :: f integer i integer id integer ier integer ierro integer ind1 integer ind2 integer iord(limit) integer iroff1 integer iroff2 integer iroff3 integer j integer jlow integer jupbnd integer k integer ksgn integer ktmin integer last integer levcur integer level(limit) integer levmax integer maxerr integer ndin(40) integer neval integer nint logical noext integer npts integer npts2 integer nres integer nrmax integer numrl2 real points(40) real pts(40) real resa real resabs real reseps real result real res3la(3) real rlist(limit) real rlist2(52) real sign real temp ! ! Test on validity of parameters. ! ier = 0 neval = 0 last = 0 result = 0.0e+00 abserr = 0.0e+00 alist(1) = a blist(1) = b rlist(1) = 0.0e+00 elist(1) = 0.0e+00 iord(1) = 0 level(1) = 0 npts = npts2 - 2 if ( npts2 < 2 ) then ier = 6 return else if ( limit <= npts .or. ( epsabs < 0.0e+00 .and. & epsrel < 0.0e+00) ) then ier = 6 return end if ! ! If any break points are provided, sort them into an ! ascending sequence. ! if ( b < a ) then sign = -1.0e+00 else sign = +1.0E+00 end if pts(1) = min ( a, b ) do i = 1, npts pts(i+1) = points(i) end do pts(npts+2) = max ( a, b ) nint = npts+1 a1 = pts(1) if ( npts /= 0 ) then do i = 1, nint do j = i+1, nint+1 if ( pts(j) < pts(i) ) then temp = pts(i) pts(i) = pts(j) pts(j) = temp end if end do end do if ( pts(1) /= min ( a, b ) .or. pts(nint+1) /= max ( a, b ) ) then ier = 6 return end if end if ! ! Compute first integral and error approximations. ! resabs = 0.0e+00 do i = 1, nint b1 = pts(i+1) call qk21 ( f, a1, b1, area1, error1, defabs, resa ) abserr = abserr + error1 result = result + area1 ndin(i) = 0 if ( error1 == resa .and. error1 /= 0.0e+00 ) then ndin(i) = 1 end if 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.0e+00 do i = 1, nint if ( ndin(i) == 1 ) then elist(i) = abserr end if errsum = errsum + elist(i) end do ! ! Test on accuracy. ! last = nint neval = 21 * nint dres = abs ( result ) errbnd = max ( epsabs, epsrel * dres ) if ( abserr <= 1.0e+02 * epsilon ( resabs ) * resabs .and. & abserr > errbnd ) then ier = 2 end if if ( nint /= 1 ) then do i = 1, npts jlow = i+1 ind1 = iord(i) do j = jlow, nint ind2 = iord(j) if ( elist(ind1) <= elist(ind2) ) then ind1 = ind2 k = j end if end do if ( ind1 /= iord(i) ) then iord(k) = iord(i) iord(i) = ind1 end if end do if ( limit < npts2 ) then ier = 1 end if end if if ( ier /= 0 .or. abserr <= errbnd ) then return end if ! ! 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 = huge ( abserr ) if ( dres >= ( 1.0e+00 - 0.5E+00 * epsilon ( resabs ) ) * resabs ) then ksgn = 1 else ksgn = -1 end if do last = npts2, limit ! ! Bisect the subinterval with the NRMAX-th largest error estimate. ! levcur = level(maxerr) + 1 a1 = alist(maxerr) b1 = 0.5E+00 * ( alist(maxerr) + blist(maxerr) ) a2 = b1 b2 = blist(maxerr) erlast = errmax call qk21 ( f, a1, b1, area1, error1, resa, defab1 ) call qk21 ( 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 /= error1 .and. defab2 /= error2 ) then if ( abs ( rlist ( maxerr ) - area12 ) <= 1.0e-05 * abs(area12) .and. & erro12 >= 9.9e-01 * errmax ) then if ( extrap ) then iroff2 = iroff2+1 else iroff1 = iroff1+1 end if end if if ( last > 10 .and. erro12 > errmax ) then iroff3 = iroff3 + 1 end if end if level(maxerr) = levcur level(last) = levcur rlist(maxerr) = area1 rlist(last) = area2 errbnd = max ( epsabs, epsrel * abs ( area ) ) ! ! Test for roundoff error and eventually set error flag. ! if ( 10 <= iroff1 + iroff2 .or. 20 <= iroff3 ) then ier = 2 end if if ( 5 <= iroff2 ) then ierro = 3 end if ! ! Set error flag in the case that the number of subintervals ! equals limit. ! if ( last == limit ) then ier = 1 end if ! ! Set error flag in the case of bad integrand behavior ! at a point of the integration range ! if ( max ( abs(a1), abs(b2)) <= ( 1.0e+00 + 1.0e+03 * epsilon ( a1 ) )* & ( abs ( a2 ) + 1.0e+03 * tiny ( a2 ) ) ) then ier = 4 end if ! ! 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 QSORT 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 qsort ( limit, last, maxerr, errmax, elist, iord, nrmax ) if ( errsum <= errbnd ) then go to 190 end if if ( ier /= 0 ) then exit end if if ( noext ) then cycle end if erlarg = erlarg - erlast if ( levcur+1 <= levmax ) then erlarg = erlarg + erro12 end if ! ! Test whether the interval to be bisected next is the ! smallest interval. ! if ( .not. extrap ) then if ( level(maxerr)+1 <= levmax ) then cycle end if extrap = .true. nrmax = 2 end if ! ! The smallest interval has the largest error. ! Before bisecting decrease the sum of the errors over the ! larger intervals (erlarg) and perform extrapolation. ! if ( ierro /= 3 .and. erlarg > ertest ) then id = nrmax jupbnd = last if ( last > (2+limit/2) ) then jupbnd = limit+3-last end if do k = id, jupbnd maxerr = iord(nrmax) errmax = elist(maxerr) if ( level(maxerr)+1 <= levmax ) then go to 160 end if nrmax = nrmax + 1 end do end if ! ! Perform extrapolation. ! numrl2 = numrl2 + 1 rlist2(numrl2) = area if ( numrl2 <= 2 ) then go to 155 end if call qextr ( numrl2, rlist2, reseps, abseps, res3la, nres ) ktmin = ktmin+1 if ( 5 < ktmin .and. abserr < 1.0e-03 * errsum ) then ier = 5 end if if ( abseps < abserr ) then ktmin = 0 abserr = abseps result = reseps correc = erlarg ertest = max ( epsabs, epsrel * abs(reseps) ) if ( abserr < ertest ) then exit end if end if ! ! Prepare bisection of the smallest interval. ! if ( numrl2 == 1 ) then noext = .true. end if if ( 5 <= ier ) then exit end if 155 continue 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 == huge ( abserr ) ) then go to 190 end if if ( ( ier + ierro ) == 0 ) then go to 180 end if if ( ierro == 3 ) then abserr = abserr + correc end if if ( ier == 0 ) then ier = 3 end if if ( result /= 0.0e+00 .and. area /= 0.0e+00 ) then go to 175 end if if ( errsum < abserr ) then go to 190 end if if ( area == 0.0e+00 ) then go to 210 end if go to 180 175 continue if ( abserr / abs(result) > errsum / abs(area) ) then go to 190 end if ! ! Test on divergence. ! 180 continue if ( ksgn == (-1) .and. max ( abs(result),abs(area)) <= & resabs*1.0e-02 ) go to 210 if ( 1.0e-02 > (result/area) .or. (result/area) > 1.0e+02 .or. & errsum > abs(area) ) then ier = 6 end if go to 210 ! ! Compute global integral sum. ! 190 continue result = sum ( rlist(1:last) ) abserr = errsum 210 continue if ( 2 < ier ) then ier = ier - 1 end if result = result * sign return end subroutine qags ( f, a, b, epsabs, epsrel, result, abserr, neval, ier ) !*****************************************************************************80 ! !! QAGS estimates the integral of a function. ! ! Discussion: ! ! The routine calculates an approximation RESULT to a definite integral ! I = integral of F over (A,B), ! hopefully satisfying ! || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ). ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, B, the limits of integration. ! ! Input, real EPSABS, EPSREL, the absolute and relative accuracy requested. ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! Output, integer IER, error flag. ! 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 integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more sub- ! divisions by increasing the data value of ! limit in qags (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 sub- ! ranges. 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 behavior 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 < 0, ! result, abserr and neval are set to zero. ! ! Local Parameters: ! ! 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) ! implicit none integer, parameter :: limit = 500 real a real abseps real abserr real alist(limit) real area real area1 real area12 real area2 real a1 real a2 real b real blist(limit) real b1 real b2 real correc real defabs real defab1 real defab2 real dres real elist(limit) real epsabs real epsrel real erlarg real erlast real errbnd real errmax real error1 real error2 real erro12 real errsum real ertest logical extrap real, external :: f integer id integer ier integer ierro integer iord(limit) integer iroff1 integer iroff2 integer iroff3 integer jupbnd integer k integer ksgn integer ktmin integer last logical noext integer maxerr integer neval integer nres integer nrmax integer numrl2 real resabs real reseps real result real res3la(3) real rlist(limit) real rlist2(52) real small ! ! The dimension of rlist2 is determined by the value of ! limexp in QEXTR (rlist2 should be of dimension ! (limexp+2) at least). ! ! Test on validity of parameters. ! ier = 0 neval = 0 last = 0 result = 0.0e+00 abserr = 0.0e+00 alist(1) = a blist(1) = b rlist(1) = 0.0e+00 elist(1) = 0.0e+00 if ( epsabs < 0.0e+00 .and. epsrel < 0.0e+00 ) then ier = 6 return end if ! ! First approximation to the integral. ! ierro = 0 call qk21 ( f, a, b, result, abserr, defabs, resabs ) ! ! Test on accuracy. ! dres = abs ( result ) errbnd = max ( epsabs, epsrel * dres ) last = 1 rlist(1) = result elist(1) = abserr iord(1) = 1 if ( abserr <= 1.0e+02 * epsilon ( defabs ) * defabs .and. & abserr > errbnd ) then ier = 2 end if if ( limit == 1 ) then ier = 1 end if if ( ier /= 0 .or. (abserr <= errbnd .and. abserr /= resabs ) .or. & abserr == 0.0e+00 ) go to 140 ! ! Initialization. ! rlist2(1) = result errmax = abserr maxerr = 1 area = result errsum = abserr abserr = huge ( abserr ) nrmax = 1 nres = 0 numrl2 = 2 ktmin = 0 extrap = .false. noext = .false. iroff1 = 0 iroff2 = 0 iroff3 = 0 if ( dres >= (1.0e+00-5.0e+01* epsilon ( defabs ) ) * defabs ) then ksgn = 1 else ksgn = -1 end if do last = 2, limit ! ! Bisect the subinterval with the nrmax-th largest error estimate. ! a1 = alist(maxerr) b1 = 5.0e-01 * ( alist(maxerr) + blist(maxerr) ) a2 = b1 b2 = blist(maxerr) erlast = errmax call qk21 ( f, a1, b1, area1, error1, resabs, defab1 ) call qk21 ( 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 .or. defab2 == error2 ) go to 15 if ( abs ( rlist(maxerr) - area12) > 1.0e-05 * abs(area12) & .or. erro12 < 9.9e-01 * errmax ) go to 10 if ( extrap ) then iroff2 = iroff2+1 else iroff1 = iroff1+1 end if 10 continue if ( last > 10 .and. erro12 > errmax ) then iroff3 = iroff3+1 end if 15 continue rlist(maxerr) = area1 rlist(last) = area2 errbnd = max ( epsabs, epsrel*abs(area) ) ! ! Test for roundoff error and eventually set error flag. ! if ( iroff1+iroff2 >= 10 .or. iroff3 >= 20 ) then ier = 2 end if if ( iroff2 >= 5 ) then ierro = 3 end if ! ! Set error flag in the case that the number of subintervals ! equals limit. ! if ( last == limit ) then ier = 1 end if ! ! Set error flag in the case of bad integrand behavior ! at a point of the integration range. ! if ( max ( abs(a1),abs(b2)) <= (1.0e+00+1.0e+03* epsilon ( a1 ) )* & (abs(a2)+1.0e+03* tiny ( a2 ) ) ) then ier = 4 end if ! ! 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 QSORT 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 qsort ( limit, last, maxerr, errmax, elist, iord, nrmax ) if ( errsum <= errbnd ) go to 115 if ( ier /= 0 ) then exit end if if ( last == 2 ) go to 80 if ( noext ) go to 90 erlarg = erlarg-erlast if ( abs(b1-a1) > small ) then erlarg = erlarg+erro12 end if ! ! Test whether the interval to be bisected next is the ! smallest interval. ! if ( .not. extrap ) then if ( abs(blist(maxerr)-alist(maxerr)) > small ) go to 90 extrap = .true. nrmax = 2 end if !40 continue ! ! The smallest interval has the largest error. ! Before bisecting decrease the sum of the errors over the ! larger intervals (erlarg) and perform extrapolation. ! if ( ierro /= 3 .and. erlarg > ertest ) then id = nrmax jupbnd = last if ( last > (2+limit/2) ) then jupbnd = limit+3-last end if do k = id, jupbnd maxerr = iord(nrmax) errmax = elist(maxerr) if ( abs(blist(maxerr)-alist(maxerr)) > small ) then go to 90 end if nrmax = nrmax+1 end do end if ! ! Perform extrapolation. ! !60 continue numrl2 = numrl2+1 rlist2(numrl2) = area call qextr ( numrl2, rlist2, reseps, abseps, res3la, nres ) ktmin = ktmin+1 if ( ktmin > 5 .and. abserr < 1.0e-03 * errsum ) then ier = 5 end if if ( abseps < abserr ) then ktmin = 0 abserr = abseps result = reseps correc = erlarg ertest = max ( epsabs,epsrel*abs(reseps)) if ( abserr <= ertest ) then exit end if end if ! ! Prepare bisection of the smallest interval. ! if ( numrl2 == 1 ) then noext = .true. end if if ( ier == 5 ) then exit end if maxerr = iord(1) errmax = elist(maxerr) nrmax = 1 extrap = .false. small = small * 5.0e-01 erlarg = errsum go to 90 80 continue small = abs ( b - a ) * 3.75e-01 erlarg = errsum ertest = errbnd rlist2(2) = area 90 continue end do ! ! Set final result and error estimate. ! if ( abserr == huge ( abserr ) ) then go to 115 end if if ( ier + ierro == 0 ) then go to 110 end if if ( ierro == 3 ) then abserr = abserr + correc end if if ( ier == 0 ) then ier = 3 end if if ( result /= 0.0e+00.and.area /= 0.0e+00 ) then go to 105 end if if ( abserr > errsum ) go to 115 if ( area == 0.0e+00 ) 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.max ( abs(result),abs(area)) <= & defabs*1.0e-02 ) go to 130 if ( 1.0e-02 > (result/area) .or. (result/area) > 1.0e+02 & .or. errsum > abs(area) ) then ier = 6 end if go to 130 ! ! Compute global integral sum. ! 115 continue result = sum ( rlist(1:last) ) abserr = errsum 130 continue if ( 2 < ier ) then ier = ier - 1 end if 140 continue neval = 42*last-21 return end subroutine qawc ( f, a, b, c, epsabs, epsrel, result, abserr, neval, ier ) !*****************************************************************************80 ! !! QAWC computes a Cauchy principal value. ! ! Discussion: ! ! The routine calculates an approximation RESULT to a Cauchy principal ! value ! I = integral of F*W over (A,B), ! with ! W(X) = 1 / (X-C), ! with C distinct from A and B, hopefully satisfying ! || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ). ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, B, the limits of integration. ! ! Input, real C, a parameter in the weight function, which must ! not be equal to A or B. ! ! Input, real EPSABS, EPSREL, the absolute and relative accuracy requested. ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! 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 integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more sub- ! divisions by increasing the data value of ! limit in qawc (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 appropriate integrators on the ! subranges. ! = 2 the occurrence of roundoff error is detec- ! ted, which prevents the requested ! tolerance from being achieved. ! = 3 extremely bad integrand behavior occurs ! at some points of the integration ! interval. ! = 6 the input is invalid, because ! c = a or c = b or ! epsabs < 0 and epsrel < 0, ! result, abserr, neval are set to zero. ! ! Local parameters: ! ! LIMIT is the maximum number of subintervals allowed in the ! subdivision process of qawce. take care that limit >= 1. ! implicit none integer, parameter :: limit = 500 real a real abserr real alist(limit) real b real blist(limit) real elist(limit) real c real epsabs real epsrel real, external :: f integer ier integer iord(limit) integer last integer neval real result real rlist(limit) call qawce ( f, a, b, c, epsabs, epsrel, limit, result, abserr, neval, ier, & alist, blist, rlist, elist, iord, last ) return end subroutine qawce ( f, a, b, c, epsabs, epsrel, limit, result, abserr, neval, & ier, alist, blist, rlist, elist, iord, last ) !*****************************************************************************80 ! !! QAWCE computes a Cauchy principal value. ! ! Discussion: ! ! The routine calculates an approximation RESULT to a Cauchy principal ! value ! I = integral of F*W over (A,B), ! with ! W(X) = 1 / ( X - C ), ! with C distinct from A and B, hopefully satisfying ! | I - RESULT | <= max ( EPSABS, EPSREL * |I| ). ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, B, the limits of integration. ! ! Input, real C, a parameter in the weight function, which cannot be ! equal to A or B. ! ! Input, real EPSABS, EPSREL, the absolute and relative accuracy requested. ! ! Input, integer LIMIT, the upper bound on the number of subintervals that ! will be used in the partition of [A,B]. LIMIT is typically 500. ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! 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 integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more sub- ! divisions by increasing the value of ! limit. 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 appropriate integrators ! on the subranges. ! = 2 the occurrence of roundoff error is detec- ! ted, which prevents the requested ! tolerance from being achieved. ! = 3 extremely bad integrand behavior occurs ! at some interior points of the integration ! interval. ! = 6 the input is invalid, because ! c = a or c = b or ! epsabs < 0 and epsrel < 0, ! or limit < 1. ! result, abserr, neval, rlist(1), elist(1), ! iord(1) and last are set to zero. ! alist(1) and blist(1) are set to a and b ! respectively. ! ! Workspace, real ALIST(LIMIT), BLIST(LIMIT), contains in entries 1 ! through LAST the left and right ends of the partition subintervals. ! ! Workspace, real RLIST(LIMIT), contains in entries 1 through LAST ! the integral approximations on the subintervals. ! ! Workspace, real ELIST(LIMIT), contains in entries 1 through LAST ! the absolute error estimates on the subintervals. ! ! iord - integer ! vector of dimension at least limit, the first k ! elements of which are pointers to the error ! estimates over the subintervals, so that ! elist(iord(1)), ..., elist(iord(k)) with ! k = last if last <= (limit/2+2), and ! k = limit+1-last otherwise, form a decreasing ! sequence. ! ! last - integer ! number of subintervals actually produced in ! the subdivision process ! ! Local parameters: ! ! 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)) ! elist(i) - error estimate applying to rlist(i) ! maxerr - pointer to the interval with largest error ! estimate ! errmax - elist(maxerr) ! 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 ! implicit none integer limit real a real aa real abserr real alist(limit) real area real area1 real area12 real area2 real a1 real a2 real b real bb real blist(limit) real b1 real b2 real c real elist(limit) real epsabs real epsrel real errbnd real errmax real error1 real error2 real erro12 real errsum real, external :: f integer ier integer iord(limit) integer iroff1 integer iroff2 integer krule integer last integer maxerr integer nev integer neval integer nrmax real result real rlist(limit) ! ! Test on validity of parameters. ! ier = 0 neval = 0 last = 0 alist(1) = a blist(1) = b rlist(1) = 0.0e+00 elist(1) = 0.0e+00 iord(1) = 0 result = 0.0e+00 abserr = 0.0e+00 if ( c == a ) then ier = 6 return else if ( c == b ) then ier = 6 return else if ( epsabs < 0.0e+00 .and. epsrel < 0.0e+00 ) then ier = 6 return end if ! ! First approximation to the integral. ! if ( a <= b ) then aa = a bb = b else aa = b bb = a end if krule = 1 call qc25c ( f, aa, bb, c, result, abserr, krule, neval ) last = 1 rlist(1) = result elist(1) = abserr iord(1) = 1 alist(1) = a blist(1) = b ! ! Test on accuracy. ! errbnd = max ( epsabs, epsrel * abs(result) ) if ( limit == 1 ) then ier = 1 go to 70 end if if ( abserr < min ( 1.0e-02 * abs(result), errbnd) ) then go to 70 end if ! ! Initialization ! alist(1) = aa blist(1) = bb rlist(1) = result errmax = abserr maxerr = 1 area = result errsum = abserr nrmax = 1 iroff1 = 0 iroff2 = 0 do last = 2, limit ! ! Bisect the subinterval with nrmax-th largest error estimate. ! a1 = alist(maxerr) b1 = 5.0e-01*(alist(maxerr)+blist(maxerr)) b2 = blist(maxerr) if ( c <= b1 .and. a1 < c ) then b1 = 5.0e-01*(c+b2) end if if ( b1 < c .and. c < b2 ) then b1 = 5.0e-01 * ( a1 + c ) end if a2 = b1 krule = 2 call qc25c ( f, a1, b1, c, area1, error1, krule, nev ) neval = neval+nev call qc25c ( f, a2, b2, c, area2, error2, krule, nev ) neval = neval+nev ! ! 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 ( abs ( rlist(maxerr)-area12) < 1.0e-05 * abs(area12) & .and. erro12 >= 9.9e-01 * errmax .and. krule == 0 ) & iroff1 = iroff1+1 if ( last > 10.and.erro12 > errmax .and. krule == 0 ) then iroff2 = iroff2+1 end if rlist(maxerr) = area1 rlist(last) = area2 errbnd = max ( epsabs, epsrel * abs(area) ) if ( errsum > errbnd ) then ! ! Test for roundoff error and eventually set error flag. ! if ( iroff1 >= 6 .and. iroff2 > 20 ) then ier = 2 end if ! ! Set error flag in the case that number of interval ! bisections exceeds limit. ! if ( last == limit ) then ier = 1 end if ! ! Set error flag in the case of bad integrand behavior at ! a point of the integration range. ! if ( max ( abs(a1), abs(b2) ) <= ( 1.0e+00 + 1.0e+03 * epsilon ( a1 ) ) & *( abs(a2)+1.0e+03* tiny ( a2 ) )) then ier = 3 end if end if ! ! 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 QSORT 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 qsort ( limit, last, maxerr, errmax, elist, iord, nrmax ) if ( ier /= 0 .or. errsum <= errbnd ) then exit end if end do ! ! Compute final result. ! result = sum ( rlist(1:last) ) abserr = errsum 70 continue if ( aa == b ) then result = - result end if return end subroutine qawf ( f, a, omega, integr, epsabs, result, abserr, neval, ier ) !*****************************************************************************80 ! !! QAWF computes Fourier integrals over the interval [ A, +Infinity ). ! ! Discussion: ! ! The routine calculates an approximation RESULT to a definite integral ! ! I = integral of F*COS(OMEGA*X) ! or ! I = integral of F*SIN(OMEGA*X) ! ! over the interval [A,+Infinity), hopefully satisfying ! ! || I - RESULT || <= EPSABS. ! ! If OMEGA = 0 and INTEGR = 1, the integral is calculated by means ! of QAGI, and IER has the meaning as described in the comments of QAGI. ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, the lower limit of integration. ! ! Input, real OMEGA, the parameter in the weight function. ! ! Input, integer INTEGR, indicates which weight functions is used ! = 1, w(x) = cos(omega*x) ! = 2, w(x) = sin(omega*x) ! ! Input, real EPSABS, the absolute accuracy requested. ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! 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 integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! if omega /= 0 ! ier = 6 the input is invalid because ! (integr /= 1 and integr /= 2) or ! epsabs <= 0 ! result, abserr, neval, lst are set to ! zero. ! = 7 abnormal termination of the computation ! of one or more subintegrals ! = 8 maximum number of cycles allowed ! has been achieved, i.e. of subintervals ! (a+(k-1)c,a+kc) where ! c = (2*int(abs(omega))+1)*pi/abs(omega), ! for k = 1, 2, ... ! = 9 the extrapolation table constructed for ! convergence acceleration of the series ! formed by the integral contributions ! over the cycles, does not converge to ! within the requested accuracy. ! ! Local parameters: ! ! Integer LIMLST, gives an upper bound on the number of cycles, LIMLST >= 3. ! if limlst < 3, the routine will end with ier = 6. ! ! Integer MAXP1, an upper bound on the number of Chebyshev moments which ! can be stored, i.e. for the intervals of lengths abs(b-a)*2**(-l), ! l = 0,1, ..., maxp1-2, maxp1 >= 1. if maxp1 < 1, the routine will end ! with ier = 6. ! implicit none integer, parameter :: limit = 500 integer, parameter :: limlst = 50 integer, parameter :: maxp1 = 21 real a real abserr real alist(limit) real blist(limit) real chebmo(maxp1,25) real elist(limit) real epsabs real erlst(limlst) real, external :: f integer ier integer integr integer iord(limit) integer ierlst(limlst) integer lst integer neval integer nnlog(limit) real omega real result real rlist(limit) real rslst(limlst) ier = 6 neval = 0 result = 0.0e+00 abserr = 0.0e+00 if ( limlst < 3 .or. maxp1 < 1 ) then return end if call qawfe ( f, a, omega, integr, epsabs, limlst, limit, maxp1, & result, abserr, neval, ier, rslst, erlst, ierlst, lst, alist, blist, & rlist, elist, iord, nnlog, chebmo ) return end subroutine qawfe ( f, a, omega, integr, epsabs, limlst, limit, maxp1, & result, abserr, neval, ier, rslst, erlst, ierlst, lst, alist, blist, & rlist, elist, iord, nnlog, chebmo ) !*****************************************************************************80 ! !! QAWFE computes Fourier integrals. ! ! Discussion: ! ! The routine calculates an approximation RESULT to a definite integral ! I = integral of F*COS(OMEGA*X) or F*SIN(OMEGA*X) over (A,+Infinity), ! hopefully satisfying ! || I - RESULT || <= EPSABS. ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, the lower limit of integration. ! ! Input, real OMEGA, the parameter in the weight function. ! ! Input, integer INTEGR, indicates which weight function is used ! = 1 w(x) = cos(omega*x) ! = 2 w(x) = sin(omega*x) ! ! Input, real EPSABS, the absolute accuracy requested. ! ! Input, integer LIMLST, an upper bound on the number of cycles. ! LIMLST must be at least 1. In fact, if LIMLST < 3, the routine ! will end with IER= 6. ! ! Input, integer LIMIT, an upper bound on the number of subintervals ! allowed in the partition of each cycle, limit >= 1. ! ! maxp1 - integer ! gives an upper bound on the number of ! Chebyshev moments which can be stored, i.e. ! for the intervals of lengths abs(b-a)*2**(-l), ! l=0,1, ..., maxp1-2, maxp1 >= 1 ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! ier - 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 integral and error ! are less reliable. it is assumed that ! the requested accuracy has not been ! achieved. ! if omega /= 0 ! ier = 6 the input is invalid because ! (integr /= 1 and integr /= 2) or ! epsabs <= 0 or limlst < 3. ! result, abserr, neval, lst are set ! to zero. ! = 7 bad integrand behavior occurs within ! one or more of the cycles. location ! and type of the difficulty involved ! can be determined from the vector ierlst. ! here lst is the number of cycles actually ! needed (see below). ! ierlst(k) = 1 the maximum number of ! subdivisions (= limit) ! has been achieved on the ! k th cycle. ! = 2 occurence of roundoff ! error is detected and ! prevents the tolerance ! imposed on the k th cycle ! from being acheived. ! = 3 extremely bad integrand ! behavior occurs at some ! points of the k th cycle. ! = 4 the integration procedure ! over the k th cycle does ! not converge (to within the ! required accuracy) due to ! roundoff in the ! extrapolation procedure ! invoked on this cycle. it ! is assumed that the result ! on this interval is the ! best which can be obtained. ! = 5 the integral over the k th ! cycle is probably divergent ! or slowly convergent. it ! must be noted that ! divergence can occur with ! any other value of ! ierlst(k). ! = 8 maximum number of cycles allowed ! has been achieved, i.e. of subintervals ! (a+(k-1)c,a+kc) where ! c = (2*int(abs(omega))+1)*pi/abs(omega), ! for k = 1, 2, ..., lst. ! one can allow more cycles by increasing ! the value of limlst (and taking the ! according dimension adjustments into ! account). ! examine the array iwork which contains ! the error flags over the cycles, in order ! to eventual look for local 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 appopriate integrators on the ! subranges. ! = 9 the extrapolation table constructed for ! convergence acceleration of the series ! formed by the integral contributions ! over the cycles, does not converge to ! within the required accuracy. ! as in the case of ier = 8, it is advised ! to examine the array iwork which contains ! the error flags on the cycles. ! if omega = 0 and integr = 1, ! the integral is calculated by means of qagi ! and ier = ierlst(1) (with meaning as described ! for ierlst(k), k = 1). ! ! rslst - real ! vector of dimension at least limlst ! rslst(k) contains the integral contribution ! over the interval (a+(k-1)c,a+kc) where ! c = (2*int(abs(omega))+1)*pi/abs(omega), ! k = 1, 2, ..., lst. ! note that, if omega = 0, rslst(1) contains ! the value of the integral over (a,infinity). ! ! erlst - real ! vector of dimension at least limlst ! erlst(k) contains the error estimate ! corresponding with rslst(k). ! ! ierlst - integer ! vector of dimension at least limlst ! ierlst(k) contains the error flag corresponding ! with rslst(k). for the meaning of the local error ! flags see description of output parameter ier. ! ! lst - integer ! number of subintervals needed for the integration ! if omega = 0 then lst is set to 1. ! ! alist, blist, rlist, elist - real ! vector of dimension at least limit, ! ! iord, nnlog - integer ! vector of dimension at least limit, providing ! space for the quantities needed in the ! subdivision process of each cycle ! ! chebmo - real ! array of dimension at least (maxp1,25), ! providing space for the Chebyshev moments ! needed within the cycles ! ! Local parameters: ! ! c1, c2 - end points of subinterval (of length ! cycle) ! cycle - (2*int(abs(omega))+1)*pi/abs(omega) ! psum - vector of dimension at least (limexp+2) ! (see routine qextr) ! psum contains the part of the epsilon table ! which is still needed for further computations. ! each element of psum is a partial sum of ! the series which should sum to the value of ! the integral. ! errsum - sum of error estimates over the ! subintervals, calculated cumulatively ! epsa - absolute tolerance requested over current ! subinterval ! chebmo - array containing the modified Chebyshev ! moments (see also routine qc25o) ! implicit none integer limit integer limlst integer maxp1 real a real abseps real abserr real alist(limit) real blist(limit) real chebmo(maxp1,25) real correc real cycle real c1 real c2 real dl ! real dla real drl real elist(limit) real ep real eps real epsa real epsabs real erlst(limlst) real errsum real, external :: f real fact integer ier integer ierlst(limlst) integer integr integer iord(limit) integer ktmin integer l integer ll integer lst integer momcom integer nev integer neval integer nnlog(limit) integer nres integer numrl2 real omega real, parameter :: p = 0.9E+00 real, parameter :: pi = 3.1415926535897932E+00 real p1 real psum(52) real reseps real result real res3la(3) real rlist(limit) real rslst(limlst) ! ! The dimension of psum is determined by the value of ! limexp in QEXTR (psum must be ! of dimension (limexp+2) at least). ! ! Test on validity of parameters. ! result = 0.0e+00 abserr = 0.0e+00 neval = 0 lst = 0 ier = 0 if ( (integr /= 1 .and. integr /= 2 ) .or. & epsabs <= 0.0e+00 .or. & limlst < 3 ) then ier = 6 return end if if ( omega == 0.0e+00 ) then if ( integr == 1 ) then call qagi ( f, 0.0e+00, 1, epsabs, 0.0e+00, result, abserr, neval, ier ) else result = 0.0E+00 abserr = 0.0E+00 neval = 0 ier = 0 end if rslst(1) = result erlst(1) = abserr ierlst(1) = ier lst = 1 return end if ! ! Initializations. ! l = int ( abs ( omega ) ) dl = 2 * l + 1 cycle = dl * pi / abs ( omega ) ier = 0 ktmin = 0 neval = 0 numrl2 = 0 nres = 0 c1 = a c2 = cycle+a p1 = 1.0e+00-p eps = epsabs if ( epsabs > tiny ( epsabs ) / p1 ) then eps = epsabs * p1 end if ep = eps fact = 1.0e+00 correc = 0.0e+00 abserr = 0.0e+00 errsum = 0.0e+00 do lst = 1, limlst ! ! Integrate over current subinterval. ! ! dla = lst epsa = eps * fact call qfour ( f, c1, c2, omega, integr, epsa, 0.0e+00, limit, lst, maxp1, & rslst(lst), erlst(lst), nev, ierlst(lst), alist, blist, rlist, elist, & iord, nnlog, momcom, chebmo ) neval = neval + nev fact = fact * p errsum = errsum + erlst(lst) drl = 5.0e+01 * abs(rslst(lst)) ! ! Test on accuracy with partial sum. ! if ((errsum+drl) <= epsabs.and.lst >= 6) then go to 80 end if correc = max ( correc,erlst(lst)) if ( ierlst(lst) /= 0 ) then eps = max ( ep,correc*p1) ier = 7 end if if ( ier == 7 .and. (errsum+drl) <= correc*1.0e+01.and. lst > 5) go to 80 numrl2 = numrl2+1 if ( lst <= 1 ) then psum(1) = rslst(1) go to 40 end if psum(numrl2) = psum(ll) + rslst(lst) if ( lst == 2 ) then go to 40 end if ! ! Test on maximum number of subintervals ! if ( lst == limlst ) then ier = 8 end if ! ! Perform new extrapolation ! call qextr ( numrl2, psum, reseps, abseps, res3la, nres ) ! ! Test whether extrapolated result is influenced by roundoff ! ktmin = ktmin + 1 if ( ktmin >= 15 .and. abserr <= 1.0e-03 * (errsum+drl) ) then ier = 9 end if if ( abseps <= abserr .or. lst == 3 ) then abserr = abseps result = reseps ktmin = 0 ! ! If IER is not 0, check whether direct result (partial ! sum) or extrapolated result yields the best integral ! approximation ! if ( ( abserr + 1.0e+01 * correc ) <= epsabs ) then exit end if if ( abserr <= epsabs .and. 1.0e+01 * correc >= epsabs ) then exit end if end if if ( ier /= 0 .and. ier /= 7 ) then exit end if 40 continue ll = numrl2 c1 = c2 c2 = c2+cycle end do ! ! Set final result and error estimate. ! !60 continue abserr = abserr + 1.0e+01 * correc if ( ier == 0 ) then return end if if ( result /= 0.0e+00 .and. psum(numrl2) /= 0.0e+00) go to 70 if ( abserr > errsum ) then go to 80 end if if ( psum(numrl2) == 0.0e+00 ) then return end if 70 continue if ( abserr / abs(result) <= (errsum+drl)/abs(psum(numrl2)) ) then if ( ier >= 1 .and. ier /= 7 ) then abserr = abserr + drl end if return end if 80 continue result = psum(numrl2) abserr = errsum + drl return end subroutine qawo ( f, a, b, omega, integr, epsabs, epsrel, result, abserr, & neval, ier ) !*****************************************************************************80 ! !! QAWO computes the integrals of oscillatory integrands. ! ! Discussion: ! ! The routine calculates an approximation RESULT to a given ! definite integral ! I = Integral ( A <= X <= B ) F(X) * cos ( OMEGA * X ) dx ! or ! I = Integral ( A <= X <= B ) F(X) * sin ( OMEGA * X ) dx ! hopefully satisfying following claim for accuracy ! | I - RESULT | <= max ( epsabs, epsrel * |I| ). ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, B, the limits of integration. ! ! Input, real OMEGA, the parameter in the weight function. ! ! Input, integer INTEGR, specifies the weight function: ! 1, W(X) = cos ( OMEGA * X ) ! 2, W(X) = sin ( OMEGA * X ) ! ! Input, real EPSABS, EPSREL, the absolute and relative accuracy requested. ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! 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 integral and error are ! less reliable. it is assumed that the ! requested accuracy has not been achieved. ! ier = 1 maximum number of subdivisions allowed ! (= leniw/2) has been achieved. one can ! allow more subdivisions by increasing the ! value of leniw (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 behavior occurs ! at some interior 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 due to ! roundoff in the extrapolation table, ! 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 < 0, ! result, abserr, neval are set to zero. ! ! Local parameters: ! ! limit is the maximum number of subintervals allowed in the ! subdivision process of QFOUR. take care that limit >= 1. ! ! maxp1 gives an upper bound on the number of Chebyshev moments ! which can be stored, i.e. for the intervals of lengths ! abs(b-a)*2**(-l), l = 0, 1, ... , maxp1-2. take care that ! maxp1 >= 1. implicit none integer, parameter :: limit = 500 integer, parameter :: maxp1 = 21 real a real abserr real alist(limit) real b real blist(limit) real chebmo(maxp1,25) real elist(limit) real epsabs real epsrel real, external :: f integer ier integer integr integer iord(limit) integer momcom integer neval integer nnlog(limit) real omega real result real rlist(limit) call qfour ( f, a, b, omega, integr, epsabs, epsrel, limit, 1, maxp1, & result, abserr, neval, ier, alist, blist, rlist, elist, iord, nnlog, & momcom, chebmo ) return end subroutine qaws ( f, a, b, alfa, beta, integr, epsabs, epsrel, result, & abserr, neval, ier ) !*****************************************************************************80 ! !! QAWS estimates integrals with algebraico-logarithmic endpoint singularities. ! ! Discussion: ! ! This routine calculates an approximation RESULT to a given ! definite integral ! I = integral of f*w over (a,b) ! where w shows a singular behavior at the end points, see parameter ! integr, hopefully satisfying following claim for accuracy ! abs(i-result) <= max(epsabs,epsrel*abs(i)). ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, B, the limits of integration. ! ! Input, real ALFA, BETA, parameters used in the weight function. ! ALFA and BETA should be greater than -1. ! ! Input, integer INTEGR, indicates which weight function is to be used ! = 1 (x-a)**alfa*(b-x)**beta ! = 2 (x-a)**alfa*(b-x)**beta*log(x-a) ! = 3 (x-a)**alfa*(b-x)**beta*log(b-x) ! = 4 (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x) ! ! Input, real EPSABS, EPSREL, the absolute and relative accuracy requested. ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! 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 the integral and error ! are less reliable. it is assumed that the ! requested accuracy has not been achieved. ! ier = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the data value ! of limit in qaws (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 which prevent the requested ! tolerance from being achieved. in case of ! a jump discontinuity or a local ! singularity of algebraico-logarithmic type ! at one or more interior points of the ! integration range, one should proceed by ! splitting up the interval at these points ! and calling the integrator on the ! subranges. ! = 2 the occurrence of roundoff error is ! detected, which prevents the requested ! tolerance from being achieved. ! = 3 extremely bad integrand behavior occurs ! at some points of the integration ! interval. ! = 6 the input is invalid, because ! b <= a or alfa <= (-1) or beta <= (-1) or ! integr < 1 or integr > 4 or ! epsabs < 0 and epsrel < 0, ! result, abserr, neval are set to zero. ! ! Local parameters: ! ! LIMIT is the maximum number of subintervals allowed in the ! subdivision process of qawse. take care that limit >= 2. ! implicit none integer, parameter :: limit = 500 real a real abserr real alfa real alist(limit) real b real blist(limit) real beta real elist(limit) real epsabs real epsrel real, external :: f integer ier integer integr integer iord(limit) integer last integer neval real result real rlist(limit) call qawse ( f, a, b, alfa, beta, integr, epsabs, epsrel, limit, result, & abserr, neval, ier, alist, blist, rlist, elist, iord, last ) return end subroutine qawse ( f, a, b, alfa, beta, integr, epsabs, epsrel, limit, & result, abserr, neval, ier, alist, blist, rlist, elist, iord, last ) !*****************************************************************************80 ! !! QAWSE estimates integrals with algebraico-logarithmic endpoint singularities. ! ! Discussion: ! ! This routine calculates an approximation RESULT to an integral ! I = integral of F(X) * W(X) over (a,b), ! where W(X) shows a singular behavior at the endpoints, hopefully ! satisfying: ! | I - RESULT | <= max ( epsabs, epsrel * |I| ). ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, B, the limits of integration. ! ! Input, real ALFA, BETA, parameters used in the weight function. ! ALFA and BETA should be greater than -1. ! ! Input, integer INTEGR, indicates which weight function is used: ! = 1 (x-a)**alfa*(b-x)**beta ! = 2 (x-a)**alfa*(b-x)**beta*log(x-a) ! = 3 (x-a)**alfa*(b-x)**beta*log(b-x) ! = 4 (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x) ! ! Input, real EPSABS, EPSREL, the absolute and relative accuracy requested. ! ! Input, integer LIMIT, an upper bound on the number of subintervals ! in the partition of (A,B), LIMIT >= 2. If LIMIT < 2, the routine ! will end with IER = 6. ! ! Output, real RESULT, the estimated value of the integral. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! 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 the integral and error ! are less reliable. it is assumed that the ! requested accuracy has not been achieved. ! = 1 maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the value of ! limit. however, if this yields no ! improvement it is advised to analyze the ! integrand, in order to determine the ! integration difficulties which prevent ! the requested tolerance from being ! achieved. in case of a jump discontinuity ! or a local singularity of algebraico- ! logarithmic type at one or more interior ! points of the integration range, one ! should proceed by splitting up the ! interval at these points and calling the ! integrator on the subranges. ! = 2 the occurrence of roundoff error is ! detected, which prevents the requested ! tolerance from being achieved. ! = 3 extremely bad integrand behavior occurs ! at some points of the integration ! interval. ! = 6 the input is invalid, because ! b <= a or alfa <= (-1) or beta <= (-1) or ! integr < 1 or integr > 4, or ! epsabs < 0 and epsrel < 0, ! or limit < 2. ! result, abserr, neval, rlist(1), elist(1), ! iord(1) and last are set to zero. ! alist(1) and blist(1) are set to a and b ! respectively. ! ! Workspace, real ALIST(LIMIT), BLIST(LIMIT), contains in entries 1 ! through LAST the left and right ends of the partition subintervals. ! ! Workspace, real RLIST(LIMIT), contains in entries 1 through LAST ! the integral approximations on the subintervals. ! ! Workspace, real ELIST(LIMIT), contains in entries 1 through LAST ! the absolute error estimates on the subintervals. ! ! iord - integer ! vector of dimension at least limit, the first k ! elements of which are pointers to the error ! estimates over the subintervals, so that ! elist(iord(1)), ..., elist(iord(k)) with k = last ! if last <= (limit/2+2), and k = limit+1-last ! otherwise, form a decreasing sequence. ! ! Output, integer LAST, the number of subintervals actually produced in ! the subdivision process. ! ! Local parameters: ! ! 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)) ! elist(i) - error estimate applying to rlist(i) ! maxerr - pointer to the interval with largest error ! estimate ! errmax - elist(maxerr) ! 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 ! implicit none integer limit real a real abserr real alfa real alist(limit) real area real area1 real area12 real area2 real a1 real a2 real b real beta real blist(limit) real b1 real b2 real centre real elist(limit) real epsabs real epsrel real errbnd real errmax real error1 real erro12 real error2 real errsum real, external :: f integer ier integer integr integer iord(limit) integer iroff1 integer iroff2 integer last integer maxerr integer nev integer neval integer nrmax real resas1 real resas2 real result real rg(25) real rh(25) real ri(25) real rj(25) real rlist(limit) ! ! Test on validity of parameters. ! ier = 0 neval = 0 last = 0 rlist(1) = 0.0e+00 elist(1) = 0.0e+00 iord(1) = 0 result = 0.0e+00 abserr = 0.0e+00 if ( b <= a .or. & (epsabs < 0.0e+00 .and. epsrel < 0.0e+00) .or. & alfa <= (-1.0e+00) .or. & beta <= (-1.0e+00) .or. & integr < 1 .or. & integr > 4 .or. & limit < 2 ) then ier = 6 return end if ! ! Compute the modified Chebyshev moments. ! call qmomo ( alfa, beta, ri, rj, rg, rh, integr ) ! ! Integrate over the intervals (a,(a+b)/2) and ((a+b)/2,b). ! centre = 5.0e-01 * ( b + a ) call qc25s ( f, a, b, a, centre, alfa, beta, ri, rj, rg, rh, area1, & error1, resas1, integr, nev ) neval = nev call qc25s ( f, a, b, centre, b, alfa, beta, ri, rj, rg, rh, area2, & error2, resas2, integr, nev ) last = 2 neval = neval+nev result = area1+area2 abserr = error1+error2 ! ! Test on accuracy. ! errbnd = max ( epsabs, epsrel * abs ( result ) ) ! ! Initialization. ! if ( error2 <= error1 ) then alist(1) = a alist(2) = centre blist(1) = centre blist(2) = b rlist(1) = area1 rlist(2) = area2 elist(1) = error1 elist(2) = error2 else alist(1) = centre alist(2) = a blist(1) = b blist(2) = centre rlist(1) = area2 rlist(2) = area1 elist(1) = error2 elist(2) = error1 end if iord(1) = 1 iord(2) = 2 if ( limit == 2 ) then ier = 1 return end if if ( abserr <= errbnd ) then return end if errmax = elist(1) maxerr = 1 nrmax = 1 area = result errsum = abserr iroff1 = 0 iroff2 = 0 do last = 3, limit ! ! Bisect the subinterval with largest error estimate. ! a1 = alist(maxerr) b1 = 5.0e-01 * ( alist(maxerr) + blist(maxerr) ) a2 = b1 b2 = blist(maxerr) call qc25s ( f, a, b, a1, b1, alfa, beta, ri, rj, rg, rh, area1, & error1, resas1, integr, nev ) neval = neval + nev call qc25s ( f, a, b, a2, b2, alfa, beta, ri, rj, rg, rh, area2, & error2, resas2, integr, nev ) neval = neval + nev ! ! Improve previous approximations integral and error and ! test for accuracy. ! area12 = area1+area2 erro12 = error1+error2 errsum = errsum+erro12-errmax area = area+area12-rlist(maxerr) ! ! Test for roundoff error. ! if ( a /= a1 .and. b /= b2 ) then if ( resas1 /= error1 .and. resas2 /= error2 ) then if ( abs ( rlist(maxerr) - area12 ) < 1.0e-05 * abs ( area12 ) & .and.erro12 >= 9.9e-01*errmax ) then iroff1 = iroff1 + 1 end if if ( last > 10 .and. erro12 > errmax ) then iroff2 = iroff2 + 1 end if end if end if rlist(maxerr) = area1 rlist(last) = area2 ! ! Test on accuracy. ! errbnd = max ( epsabs, epsrel * abs ( area ) ) if ( errsum > errbnd ) then ! ! Set error flag in the case that the number of interval ! bisections exceeds limit. ! if ( last == limit ) then ier = 1 end if ! ! Set error flag in the case of roundoff error. ! if ( iroff1 >= 6 .or. iroff2 >= 20 ) then ier = 2 end if ! ! Set error flag in the case of bad integrand behavior ! at interior points of integration range. ! if ( max ( abs(a1),abs(b2)) <= (1.0e+00+1.0e+03* epsilon ( a1 ) )* & ( abs(a2) + 1.0e+03* tiny ( a2) )) then ier = 3 end if end if ! ! 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 QSORT to maintain the descending ordering ! in the list of error estimates and select the subinterval ! with largest error estimate (to be bisected next). ! call qsort ( limit, last, maxerr, errmax, elist, iord, nrmax ) if ( ier /= 0 .or. errsum <= errbnd ) then exit end if end do ! ! Compute final result. ! result = sum ( rlist(1:last) ) abserr = errsum return end subroutine qc25c ( f, a, b, c, result, abserr, krul, neval ) !*****************************************************************************80 ! !! QC25C returns integration rules for Cauchy Principal Value integrals. ! ! Discussion: ! ! This routine estimates ! I = integral of F(X) * W(X) over (a,b) ! with error estimate, where ! w(x) = 1/(x-c) ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, B, the limits of integration. ! ! Input, real C, the parameter in the weight function. ! ! Output, real RESULT, the estimated value of the integral. ! RESULT is computed by using a generalized Clenshaw-Curtis method if ! C lies within ten percent of the integration interval. In the ! other case the 15-point Kronrod rule obtained by optimal addition ! of abscissae to the 7-point Gauss rule, is applied. ! ! Output, real ABSERR, an estimate of || I - RESULT ||. ! ! krul - integer ! key which is decreased by 1 if the 15-point ! Gauss-Kronrod scheme has been used ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! Local parameters: ! ! fval - value of the function f at the points ! cos(k*pi/24), k = 0, ..., 24 ! cheb12 - Chebyshev series expansion coefficients, for the ! function f, of degree 12 ! cheb24 - Chebyshev series expansion coefficients, for the ! function f, of degree 24 ! res12 - approximation to the integral corresponding to the ! use of cheb12 ! res24 - approximation to the integral corresponding to the ! use of cheb24 ! qwgtc - external function subprogram defining the weight ! function ! hlgth - half-length of the interval ! centr - mid point of the interval ! implicit none real a real abserr real ak22 real amom0 real amom1 real amom2 real b real c real cc real centr real cheb12(13) real cheb24(25) real, external :: f real fval(25) real hlgth integer i integer isym integer k integer kp integer krul integer neval real p2 real p3 real p4 real, external :: qwgtc real resabs real resasc real result real res12 real res24 real u real, parameter, dimension ( 11 ) :: x = (/ & 9.914448613738104e-01, 9.659258262890683e-01, & 9.238795325112868e-01, 8.660254037844386e-01, & 7.933533402912352e-01, 7.071067811865475e-01, & 6.087614290087206e-01, 5.000000000000000e-01, & 3.826834323650898e-01, 2.588190451025208e-01, & 1.305261922200516e-01 /) ! ! Check the position of C. ! cc = ( 2.0e+00 * c - b - a ) / ( b - a ) ! ! Apply the 15-point Gauss-Kronrod scheme. ! if ( abs ( cc ) >= 1.1e+00 ) then krul = krul - 1 call qk15w ( f, qwgtc, c, p2, p3, p4, kp, a, b, result, abserr, & resabs, resasc ) neval = 15 if ( resasc == abserr ) then krul = krul+1 end if return end if ! ! Use the generalized Clenshaw-Curtis method. ! hlgth = 5.0e-01 * ( b - a ) centr = 5.0e-01 * ( b + a ) neval = 25 fval(1) = 5.0e-01 * f(hlgth+centr) fval(13) = f(centr) fval(25) = 5.0e-01 * f(centr-hlgth) do i = 2, 12 u = hlgth * x(i-1) isym = 26 - i fval(i) = f(u+centr) fval(isym) = f(centr-u) end do ! ! Compute the Chebyshev series expansion. ! call qcheb ( x, fval, cheb12, cheb24 ) ! ! The modified Chebyshev moments are computed by forward ! recursion, using AMOM0 and AMOM1 as starting values. ! amom0 = log ( abs ( ( 1.0e+00 - cc ) / ( 1.0e+00 + cc ) ) ) amom1 = 2.0e+00 + cc * amom0 res12 = cheb12(1) * amom0 + cheb12(2) * amom1 res24 = cheb24(1) * amom0 + cheb24(2) * amom1 do k = 3, 13 amom2 = 2.0e+00 * cc * amom1 - amom0 ak22 = ( k - 2 ) * ( k - 2 ) if ( ( k / 2 ) * 2 == k ) then amom2 = amom2 - 4.0e+00 / ( ak22 - 1.0e+00 ) end if res12 = res12 + cheb12(k) * amom2 res24 = res24 + cheb24(k) * amom2 amom0 = amom1 amom1 = amom2 end do do k = 14, 25 amom2 = 2.0e+00 * cc * amom1 - amom0 ak22 = ( k - 2 ) * ( k - 2 ) if ( ( k / 2 ) * 2 == k ) then amom2 = amom2 - 4.0e+00 / ( ak22 - 1.0e+00 ) end if res24 = res24 + cheb24(k) * amom2 amom0 = amom1 amom1 = amom2 end do result = res24 abserr = abs ( res24 - res12 ) return end subroutine qc25o ( f, a, b, omega, integr, nrmom, maxp1, ksave, result, & abserr, neval, resabs, resasc, momcom, chebmo ) !*****************************************************************************80 ! !! QC25O returns integration rules for integrands with a COS or SIN factor. ! ! Discussion: ! ! This routine estimates the integral ! I = integral of f(x) * w(x) over (a,b) ! where ! w(x) = cos(omega*x) ! or ! w(x) = sin(omega*x), ! and estimates ! J = integral ( A <= X <= B ) |F(X)| dx. ! ! For small values of OMEGA or small intervals (a,b) the 15-point ! Gauss-Kronrod rule is used. In all other cases a generalized ! Clenshaw-Curtis method is used, that is, a truncated Chebyshev ! expansion of the function F is computed on (a,b), so that the ! integrand can be written as a sum of terms of the form W(X)*T(K,X), ! where T(K,X) is the Chebyshev polynomial of degree K. The Chebyshev ! moments are computed with use of a linear recurrence relation. ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real F, the name of the function routine, of the form ! function f ( x ) ! real f ! real x ! which evaluates the integrand function. ! ! Input, real A, B, the limits of integration. ! ! Input, real OMEGA, the parameter in the weight function. ! ! Input, integer INTEGR, indicates which weight function is to be used ! = 1, w(x) = cos(omega*x) ! = 2, w(x) = sin(omega*x) ! ! ?, integer NRMOM, the length of interval (a,b) is equal to the length ! of the original integration interval divided by ! 2**nrmom (we suppose that the routine is used in an ! adaptive integration process, otherwise set ! nrmom = 0). nrmom must be zero at the first call. ! ! maxp1 - integer ! gives an upper bound on the number of Chebyshev ! moments which can be stored, i.e. for the intervals ! of lengths abs(bb-aa)*2**(-l), l = 0,1,2, ..., ! maxp1-2. ! ! ksave - integer ! key which is one when the moments for the ! current interval have been computed ! ! Output, real RESULT, the estimated value of the integral. ! ! abserr - real ! estimate of the modulus of the absolute ! error, which should equal or exceed abs(i-result) ! ! Output, integer NEVAL, the number of times the integral was evaluated. ! ! Output, real RESABS, approximation to the integral J. ! ! Output, real RESASC, approximation to the integral of abs(F-I/(B-A)). ! ! on entry and return ! momcom - integer ! for each interval length we need to compute ! the Chebyshev moments. momcom counts the number ! of intervals for which these moments have already ! been computed. if nrmom < momcom or ksave = 1, ! the Chebyshev moments for the interval (a,b) !