| /* |
| * R : A Computer Language for Statistical Data Analysis |
| * Copyright (C) 2001-2014 The R Core Team |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| * |
| * Most of this file is C translations of Fortran routines in |
| * QUADPACK: the latter is part of SLATEC 'and therefore in the public |
| * domain' (https://en.wikipedia.org/wiki/QUADPACK). |
| * |
| * |
| */ |
| |
| #ifdef HAVE_CONFIG_H |
| #include <config.h> |
| #endif |
| |
| #include <math.h> |
| #include <float.h> |
| #include <Rmath.h> /* for fmax2, fmin2, imin2 */ |
| #include <R_ext/Applic.h> /* exporting the API , particularly */ |
| /*--- typedef void integr_fn(double *x, int n, void *ex) --- |
| * vectorizing function f(x[1:n], ...) -> x[] {overwriting x[]}. |
| * Vectorization can be used to speed up the integrand |
| * instead of calling it n times. |
| */ |
| |
| /* f2c-ed translations + modifications of QUADPACK functions */ |
| |
| static void rdqagie(integr_fn f, void *ex, |
| double *, int *, double * , double *, int *, |
| double *, double *, int *, |
| int *, double *, double *, double *, double *, |
| int *, int *); |
| |
| static void rdqk15i(integr_fn f, void *ex, |
| double *, int *, double * , double *, |
| double *, double *, double *, double *); |
| |
| static void rdqagse(integr_fn f, void *ex, double *, double *, |
| double *, double *, int *, double *, double *, |
| int *, int *, double *, double *, double *, |
| double *, int *, int *); |
| |
| static void rdqk21(integr_fn f, void *ex, |
| double *, double *, double *, double *, double *, double *); |
| |
| static void rdqpsrt(int *, int *, int *, double *, double *, int *, int *); |
| |
| static void rdqelg(int *, double *, double *, double *, double *, int *); |
| |
| /* Table of constant values */ |
| |
| static double c_b6 = 0.; |
| static double c_b7 = 1.; |
| |
| void Rdqagi(integr_fn f, void *ex, double *bound, int *inf, |
| double *epsabs, double *epsrel, |
| double *result, double *abserr, int *neval, int *ier, |
| int *limit, int *lenw, int *last, |
| int *iwork, double *work) |
| { |
| int l1, l2, l3; |
| |
| /* |
| ***begin prologue dqagi |
| ***date written 800101 (yymmdd) |
| ***revision date 830518 (yymmdd) |
| ***category no. h2a3a1,h2a4a1 |
| ***keywords automatic integrator, infinite intervals, |
| general-purpose, transformation, extrapolation, |
| globally adaptive |
| ***author piessens,robert,appl. math. & progr. div. - k.u.leuven |
| de doncker,elise,appl. math. & progr. div. -k.u.leuven |
| ***purpose the routine calculates an approximation result to a given |
| integral i = integral of f over (bound,+infinity) |
| or i = integral of f over (-infinity,bound) |
| or i = integral of f over (-infinity,+infinity) |
| hopefully satisfying following claim for accuracy |
| abs(i-result) <= max(epsabs,epsrel*abs(i)). |
| ***description |
| |
| integration over infinite intervals |
| standard fortran subroutine |
| |
| parameters |
| on entry |
| f - double precision |
| function subprogram defining the integrand |
| function f(x). the actual name for f needs to be |
| declared e x t e r n a l in the driver program. |
| |
| bound - double precision |
| finite bound of integration range |
| (has no meaning if interval is doubly-infinite) |
| |
| inf - int |
| indicating the kind of integration range involved |
| inf = 1 corresponds to (bound,+infinity), |
| inf = -1 to (-infinity,bound), |
| inf = 2 to (-infinity,+infinity). |
| |
| epsabs - double precision |
| absolute accuracy requested |
| epsrel - double precision |
| relative accuracy requested |
| if epsabs <= 0 |
| and epsrel < max(50*rel.mach.acc.,0.5d-28), |
| the routine will end with ier = 6. |
| |
| |
| on return |
| result - double precision |
| approximation to the integral |
| |
| abserr - double precision |
| estimate of the modulus of the absolute error, |
| which should equal or exceed abs(i-result) |
| |
| neval - int |
| number of integrand evaluations |
| |
| ier - int |
| ier = 0 normal and reliable termination of the |
| routine. it is assumed that the requested |
| accuracy has been achieved. |
| - ier > 0 abnormal termination of the routine. the |
| estimates for result and error are less |
| reliable. it is assumed that the requested |
| accuracy has not been achieved. |
| error messages |
| ier = 1 maximum number of subdivisions allowed |
| has been achieved. one can allow more |
| subdivisions by increasing the value of |
| limit (and taking the according dimension |
| adjustments into account). however, if |
| this yields no improvement it is advised |
| to analyze the integrand in order to |
| determine the integration difficulties. if |
| the position of a local difficulty can be |
| determined (e.g. singularity, |
| discontinuity within the interval) one |
| will probably gain from splitting up the |
| interval at this point and calling the |
| integrator on the subranges. if possible, |
| an appropriate special-purpose integrator |
| should be used, which is designed for |
| handling the type of difficulty involved. |
| = 2 the occurrence of roundoff error is |
| detected, which prevents the requested |
| tolerance from being achieved. |
| the error may be under-estimated. |
| = 3 extremely bad integrand behaviour occurs |
| at some points of the integration |
| interval. |
| = 4 the algorithm does not converge. |
| roundoff error is detected in the |
| extrapolation table. |
| it is assumed that the requested tolerance |
| cannot be achieved, and that the returned |
| result is the best which can be obtained. |
| = 5 the integral is probably divergent, or |
| slowly convergent. it must be noted that |
| divergence can occur with any other value |
| of ier. |
| = 6 the input is invalid, because |
| (epsabs <= 0 and |
| epsrel < max(50*rel.mach.acc.,0.5d-28)) |
| or limit < 1 or leniw < limit*4. |
| result, abserr, neval, last are set to |
| zero. exept when limit or leniw is |
| invalid, iwork(1), work(limit*2+1) and |
| work(limit*3+1) are set to zero, work(1) |
| is set to a and work(limit+1) to b. |
| |
| dimensioning parameters |
| limit - int |
| dimensioning parameter for iwork |
| limit determines the maximum number of subintervals |
| in the partition of the given integration interval |
| (a,b), limit >= 1. |
| if limit < 1, the routine will end with ier = 6. |
| |
| lenw - int |
| dimensioning parameter for work |
| lenw must be at least limit*4. |
| if lenw < limit*4, the routine will end |
| with ier = 6. |
| |
| last - int |
| on return, last equals the number of subintervals |
| produced in the subdivision process, which |
| determines the number of significant elements |
| actually in the work arrays. |
| |
| work arrays |
| iwork - int |
| vector of dimension at least limit, the first |
| k elements of which contain pointers |
| to the error estimates over the subintervals, |
| such that work(limit*3+iwork(1)),... , |
| work(limit*3+iwork(k)) form a decreasing |
| sequence, with k = last if last <= (limit/2+2), and |
| k = limit+1-last otherwise |
| |
| work - double precision |
| vector of dimension at least lenw |
| on return |
| work(1), ..., work(last) contain the left |
| end points of the subintervals in the |
| partition of (a,b), |
| work(limit+1), ..., work(limit+last) contain |
| the right end points, |
| work(limit*2+1), ...,work(limit*2+last) contain the |
| integral approximations over the subintervals, |
| work(limit*3+1), ..., work(limit*3) |
| contain the error estimates. |
| |
| ***routines called dqagie |
| ***end prologue dqagi */ |
| |
| *ier = 6; |
| *neval = 0; |
| *last = 0; |
| *result = 0.; |
| *abserr = 0.; |
| if (*limit < 1 || *lenw < *limit << 2) return; |
| |
| l1 = *limit; |
| l2 = *limit + l1; |
| l3 = *limit + l2; |
| |
| rdqagie(f, ex, bound, inf, epsabs, epsrel, limit, result, abserr, neval, ier, |
| work, &work[l1], &work[l2], &work[l3], iwork, last); |
| |
| return; |
| } /* Rdqagi */ |
| |
| static |
| void rdqagie(integr_fn f, void *ex, double *bound, int *inf, double * |
| epsabs, double *epsrel, int *limit, double *result, |
| double *abserr, int *neval, int *ier, double *alist, |
| double *blist, double *rlist, double *elist, int * |
| iord, int *last) |
| { |
| /* Local variables */ |
| double area, dres; |
| int ksgn; |
| double boun; |
| int nres; |
| double area1, area2, area12; |
| int k; |
| double small = 0.0, erro12; |
| int ierro; |
| double a1, a2, b1, b2, defab1, defab2, oflow; |
| int ktmin, nrmax; |
| double uflow; |
| Rboolean noext; |
| int iroff1, iroff2, iroff3; |
| double res3la[3], error1, error2; |
| int id; |
| double rlist2[52]; |
| int numrl2; |
| double defabs, epmach, erlarg = 0.0, abseps, correc = 0.0, errbnd, resabs; |
| int jupbnd; |
| double erlast, errmax; |
| int maxerr; |
| double reseps; |
| Rboolean extrap; |
| double ertest = 0.0, errsum; |
| |
| /**begin prologue dqagie |
| ***date written 800101 (yymmdd) |
| ***revision date 830518 (yymmdd) |
| ***category no. h2a3a1,h2a4a1 |
| ***keywords automatic integrator, infinite intervals, |
| general-purpose, transformation, extrapolation, |
| globally adaptive |
| ***author piessens,robert,appl. math & progr. div - k.u.leuven |
| de doncker,elise,appl. math & progr. div - k.u.leuven |
| ***purpose the routine calculates an approximation result to a given |
| integral i = integral of f over (bound,+infinity) |
| or i = integral of f over (-infinity,bound) |
| or i = integral of f over (-infinity,+infinity), |
| hopefully satisfying following claim for accuracy |
| abs(i-result) <= max(epsabs,epsrel*abs(i)) |
| ***description |
| |
| integration over infinite intervals |
| standard fortran subroutine |
| |
| f - double precision |
| function subprogram defining the integrand |
| function f(x). the actual name for f needs to be |
| declared e x t e r n a l in the driver program. |
| |
| bound - double precision |
| finite bound of integration range |
| (has no meaning if interval is doubly-infinite) |
| |
| inf - double precision |
| indicating the kind of integration range involved |
| inf = 1 corresponds to (bound,+infinity), |
| inf = -1 to (-infinity,bound), |
| inf = 2 to (-infinity,+infinity). |
| |
| epsabs - double precision |
| absolute accuracy requested |
| epsrel - double precision |
| relative accuracy requested |
| if epsabs <= 0 |
| and epsrel < max(50*rel.mach.acc.,0.5d-28), |
| the routine will end with ier = 6. |
| |
| limit - int |
| gives an upper bound on the number of subintervals |
| in the partition of (a,b), limit >= 1 |
| |
| on return |
| result - double precision |
| approximation to the integral |
| |
| abserr - double precision |
| estimate of the modulus of the absolute error, |
| which should equal or exceed abs(i-result) |
| |
| neval - int |
| number of integrand evaluations |
| |
| ier - int |
| ier = 0 normal and reliable termination of the |
| routine. it is assumed that the requested |
| accuracy has been achieved. |
| - ier > 0 abnormal termination of the routine. the |
| estimates for result and error are less |
| reliable. it is assumed that the requested |
| accuracy has not been achieved. |
| error messages |
| ier = 1 maximum number of subdivisions allowed |
| has been achieved. one can allow more |
| subdivisions by increasing the value of |
| limit (and taking the according dimension |
| adjustments into account). however,if |
| this yields no improvement it is advised |
| to analyze the integrand in order to |
| determine the integration difficulties. |
| if the position of a local difficulty can |
| be determined (e.g. singularity, |
| discontinuity within the interval) one |
| will probably gain from splitting up the |
| interval at this point and calling the |
| integrator on the subranges. if possible, |
| an appropriate special-purpose integrator |
| should be used, which is designed for |
| handling the type of difficulty involved. |
| = 2 the occurrence of roundoff error is |
| detected, which prevents the requested |
| tolerance from being achieved. |
| the error may be under-estimated. |
| = 3 extremely bad integrand behaviour occurs |
| at some points of the integration |
| interval. |
| = 4 the algorithm does not converge. |
| roundoff error is detected in the |
| extrapolation table. |
| it is assumed that the requested tolerance |
| cannot be achieved, and that the returned |
| result is the best which can be obtained. |
| = 5 the integral is probably divergent, or |
| slowly convergent. it must be noted that |
| divergence can occur with any other value |
| of ier. |
| = 6 the input is invalid, because |
| (epsabs <= 0 and |
| epsrel < max(50*rel.mach.acc.,0.5d-28), |
| result, abserr, neval, last, rlist(1), |
| elist(1) and iord(1) are set to zero. |
| alist(1) and blist(1) are set to 0 |
| and 1 respectively. |
| |
| alist - double precision |
| vector of dimension at least limit, the first |
| last elements of which are the left |
| end points of the subintervals in the partition |
| of the transformed integration range (0,1). |
| |
| blist - double precision |
| vector of dimension at least limit, the first |
| last elements of which are the right |
| end points of the subintervals in the partition |
| of the transformed integration range (0,1). |
| |
| rlist - double precision |
| vector of dimension at least limit, the first |
| last elements of which are the integral |
| approximations on the subintervals |
| |
| elist - double precision |
| vector of dimension at least limit, the first |
| last elements of which are the moduli of the |
| absolute error estimates on the subintervals |
| |
| iord - int |
| vector of dimension 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 |
| |
| last - int |
| number of subintervals actually produced |
| in the subdivision process |
| |
| ***routines called dqelg,dqk15i,dqpsrt |
| ***end prologue dqagie |
| |
| |
| the dimension of rlist2 is determined by the value of |
| limexp in subroutine dqelg. |
| |
| list of major variables |
| ----------------------- |
| |
| alist - list of left end points of all subintervals |
| considered up to now |
| blist - list of right end points of all subintervals |
| considered up to now |
| rlist(i) - approximation to the integral over |
| (alist(i),blist(i)) |
| rlist2 - array of dimension at least (limexp+2), |
| containing the part of the epsilon table |
| wich is still needed for further computations |
| elist(i) - error estimate applying to rlist(i) |
| maxerr - pointer to the interval with largest error |
| estimate |
| errmax - elist(maxerr) |
| erlast - error on the interval currently subdivided |
| (before that subdivision has taken place) |
| area - sum of the integrals over the subintervals |
| errsum - sum of the errors over the subintervals |
| errbnd - requested accuracy max(epsabs,epsrel* |
| abs(result)) |
| *****1 - variable for the left subinterval |
| *****2 - variable for the right subinterval |
| last - index for subdivision |
| nres - number of calls to the extrapolation routine |
| numrl2 - number of elements currently in rlist2. if an |
| appropriate approximation to the compounded |
| integral has been obtained, it is put in |
| rlist2(numrl2) after numrl2 has been increased |
| by one. |
| small - length of the smallest interval considered up |
| to now, multiplied by 1.5 |
| erlarg - sum of the errors over the intervals larger |
| than the smallest interval considered up to now |
| extrap - logical variable denoting that the routine |
| is attempting to perform extrapolation. i.e. |
| before subdividing the smallest interval we |
| try to decrease the value of erlarg. |
| noext - logical variable denoting that extrapolation |
| is no longer allowed (true-value) |
| |
| machine dependent constants |
| --------------------------- |
| |
| epmach is the largest relative spacing. |
| uflow is the smallest positive magnitude. |
| oflow is the largest positive magnitude. */ |
| |
| /* ***first executable statement dqagie */ |
| /* Parameter adjustments */ |
| --iord; |
| --elist; |
| --rlist; |
| --blist; |
| --alist; |
| |
| /* Function Body */ |
| epmach = DBL_EPSILON; |
| |
| /* test on validity of parameters */ |
| /* ----------------------------- */ |
| |
| *ier = 0; |
| *neval = 0; |
| *last = 0; |
| *result = 0.; |
| *abserr = 0.; |
| alist[1] = 0.; |
| blist[1] = 1.; |
| rlist[1] = 0.; |
| elist[1] = 0.; |
| iord[1] = 0; |
| if (*epsabs <= 0. && (*epsrel < fmax2(epmach * 50., 5e-29))) *ier = 6; |
| if (*ier == 6) return; |
| |
| /* 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). */ |
| |
| boun = *bound; |
| if (*inf == 2) { |
| boun = 0.; |
| } |
| rdqk15i(f, ex, &boun, inf, &c_b6, &c_b7, result, abserr, &defabs, &resabs); |
| |
| /* test on accuracy */ |
| |
| *last = 1; |
| rlist[1] = *result; |
| elist[1] = *abserr; |
| iord[1] = 1; |
| dres = fabs(*result); |
| errbnd = fmax2(*epsabs, *epsrel * dres); |
| if (*abserr <= epmach * 100. * defabs && *abserr > errbnd) *ier = 2; |
| if (*limit == 1) *ier = 1; |
| if (*ier != 0 || (*abserr <= errbnd && *abserr != resabs) |
| || *abserr == 0.) goto L130; |
| |
| /* initialization */ |
| /* -------------- */ |
| |
| uflow = DBL_MIN; |
| oflow = DBL_MAX; |
| rlist2[0] = *result; |
| errmax = *abserr; |
| maxerr = 1; |
| area = *result; |
| errsum = *abserr; |
| *abserr = oflow; |
| nrmax = 1; |
| nres = 0; |
| ktmin = 0; |
| numrl2 = 2; |
| extrap = FALSE; |
| noext = FALSE; |
| ierro = 0; |
| iroff1 = 0; |
| iroff2 = 0; |
| iroff3 = 0; |
| ksgn = -1; |
| if (dres >= (1. - epmach * 50.) * defabs) { |
| ksgn = 1; |
| } |
| |
| /* main do-loop */ |
| /* ------------ */ |
| |
| for (*last = 2; *last <= *limit; ++(*last)) { |
| |
| /* bisect the subinterval with nrmax-th largest error estimate. */ |
| |
| a1 = alist[maxerr]; |
| b1 = (alist[maxerr] + blist[maxerr]) * .5; |
| a2 = b1; |
| b2 = blist[maxerr]; |
| erlast = errmax; |
| rdqk15i(f, ex, &boun, inf, &a1, &b1, &area1, &error1, &resabs, &defab1); |
| rdqk15i(f, ex, &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 || defab2 == error2)) { |
| if (fabs(rlist[maxerr] - area12) <= fabs(area12) * 1e-5 && |
| erro12 >= errmax * .99) { |
| if (extrap) |
| ++iroff2; |
| else /* if (! extrap) */ |
| ++iroff1; |
| } |
| if (*last > 10 && erro12 > errmax) |
| ++iroff3; |
| } |
| |
| rlist[maxerr] = area1; |
| rlist[*last] = area2; |
| errbnd = fmax2(*epsabs, *epsrel * fabs(area)); |
| |
| /* test for roundoff error and eventually set error flag. */ |
| |
| if (iroff1 + iroff2 >= 10 || iroff3 >= 20) |
| *ier = 2; |
| if (iroff2 >= 5) |
| ierro = 3; |
| |
| /* set error flag in the case that the number of |
| subintervals equals limit. */ |
| |
| if (*last == *limit) |
| *ier = 1; |
| |
| /* set error flag in the case of bad integrand behaviour |
| at some points of the integration range. */ |
| |
| if (fmax2(fabs(a1), fabs(b2)) <= |
| (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3)) |
| { |
| *ier = 4; |
| } |
| |
| /* append the newly-created intervals to the list. */ |
| |
| if (error2 <= error1) { |
| 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; |
| } |
| |
| /* call subroutine dqpsrt to maintain the descending ordering |
| in the list of error estimates and select the subinterval |
| with nrmax-th largest error estimate (to be bisected next). */ |
| |
| rdqpsrt(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax); |
| if (errsum <= errbnd) { |
| goto L115; |
| } |
| if (*ier != 0) break; |
| if (*last == 2) { /* L80: */ |
| small = .375; |
| erlarg = errsum; |
| ertest = errbnd; |
| rlist2[1] = area; continue; |
| } |
| if (noext) continue; |
| |
| erlarg -= erlast; |
| if (fabs(b1 - a1) > small) { |
| erlarg += erro12; |
| } |
| if (!extrap) { |
| |
| /* test whether the interval to be bisected next is the |
| smallest interval. */ |
| |
| if (fabs(blist[maxerr] - alist[maxerr]) > small) { |
| continue; |
| } |
| extrap = TRUE; |
| nrmax = 2; |
| } |
| |
| if (ierro != 3 && erlarg > ertest) { |
| |
| /* the smallest interval has the largest error. |
| before bisecting decrease the sum of the errors over the |
| larger intervals (erlarg) and perform extrapolation. */ |
| |
| id = nrmax; |
| jupbnd = *last; |
| if (*last > *limit / 2 + 2) { |
| jupbnd = *limit + 3 - *last; |
| } |
| for (k = id; k <= jupbnd; ++k) { |
| maxerr = iord[nrmax]; |
| errmax = elist[maxerr]; |
| if (fabs(blist[maxerr] - alist[maxerr]) > small) { |
| goto L90; |
| } |
| ++nrmax; |
| /* L50: */ |
| } |
| } |
| /* perform extrapolation. L60: */ |
| ++numrl2; |
| rlist2[numrl2 - 1] = area; |
| rdqelg(&numrl2, rlist2, &reseps, &abseps, res3la, &nres); |
| ++ktmin; |
| if (ktmin > 5 && *abserr < errsum * .001) { |
| *ier = 5; |
| } |
| if (abseps >= *abserr) { |
| goto L70; |
| } |
| ktmin = 0; |
| *abserr = abseps; |
| *result = reseps; |
| correc = erlarg; |
| ertest = fmax2(*epsabs, *epsrel * fabs(reseps)); |
| if (*abserr <= ertest) { |
| break; |
| } |
| |
| /* prepare bisection of the smallest interval. */ |
| |
| L70: |
| if (numrl2 == 1) { |
| noext = TRUE; |
| } |
| if (*ier == 5) { |
| break; |
| } |
| maxerr = iord[1]; |
| errmax = elist[maxerr]; |
| nrmax = 1; |
| extrap = FALSE; |
| small *= .5; |
| erlarg = errsum; |
| L90: |
| ; |
| } |
| |
| /* L100: set final result and error estimate. */ |
| /* ------------------------------------ */ |
| |
| if (*abserr == oflow) { |
| goto L115; |
| } |
| if (*ier + ierro == 0) { |
| goto L110; |
| } |
| if (ierro == 3) { |
| *abserr += correc; |
| } |
| if (*ier == 0) { |
| *ier = 3; |
| } |
| if (*result == 0. || area == 0.) { |
| if (*abserr > errsum) |
| goto L115; |
| |
| if (area == 0.) |
| goto L130; |
| } |
| else { /* L105: */ |
| if (*abserr / fabs(*result) > errsum / fabs(area)) { |
| goto L115; |
| } |
| } |
| |
| /* test on divergence */ |
| L110: |
| if (ksgn == -1 && fmax2(fabs(*result), fabs(area)) <= defabs * .01) { |
| goto L130; |
| } |
| if (.01 > *result / area || *result / area > 100. || errsum > fabs(area)) { |
| *ier = 6; |
| } |
| goto L130; |
| |
| /* compute global integral sum. */ |
| |
| L115: |
| *result = 0.; |
| for (k = 1; k <= *last; ++k) |
| *result += rlist[k]; |
| |
| *abserr = errsum; |
| L130: |
| *neval = *last * 30 - 15; |
| if (*inf == 2) { |
| *neval <<= 1; |
| } |
| if (*ier > 2) { |
| --(*ier); |
| } |
| return; |
| } /* rdqagie_ */ |
| |
| void Rdqags(integr_fn f, void *ex, double *a, double *b, |
| double *epsabs, double *epsrel, |
| double *result, double *abserr, int *neval, int *ier, |
| int *limit, int *lenw, int *last, int *iwork, double *work) |
| { |
| int l1, l2, l3; |
| |
| /* |
| ***begin prologue dqags |
| ***date written 800101 (yymmdd) |
| ***revision date 830518 (yymmdd) |
| ***category no. h2a1a1 |
| ***keywords automatic integrator, general-purpose, |
| (end-point) singularities, extrapolation, |
| globally adaptive |
| ***author piessens,robert,appl. math. & progr. div. - k.u.leuven |
| de doncker,elise,appl. math. & prog. div. - k.u.leuven |
| ***purpose the routine calculates an approximation result to a given |
| definite integral i = integral of f over (a,b), |
| hopefully satisfying following claim for accuracy |
| abs(i-result) <= max(epsabs,epsrel*abs(i)). |
| ***description |
| |
| computation of a definite integral |
| standard fortran subroutine |
| double precision version |
| |
| |
| parameters |
| on entry |
| f - double precision |
| function subprogram defining the integrand |
| function f(x). the actual name for f needs to be |
| declared e x t e r n a l in the driver program. |
| |
| a - double precision |
| lower limit of integration |
| |
| b - double precision |
| upper limit of integration |
| |
| epsabs - double precision |
| absolute accuracy requested |
| epsrel - double precision |
| relative accuracy requested |
| if epsabs <= 0 |
| and epsrel < max(50*rel.mach.acc.,0.5d-28), |
| the routine will end with ier = 6. |
| |
| on return |
| result - double precision |
| approximation to the integral |
| |
| abserr - double precision |
| estimate of the modulus of the absolute error, |
| which should equal or exceed abs(i-result) |
| |
| neval - int |
| number of integrand evaluations |
| |
| ier - int |
| 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. |
| error messages |
| ier = 1 maximum number of subdivisions allowed |
| has been achieved. one can allow more sub- |
| divisions by increasing the value of limit |
| (and taking the according dimension |
| adjustments into account. however, if |
| this yields no improvement it is advised |
| to analyze the integrand in order to |
| determine the integration difficulties. if |
| the position of a local difficulty can be |
| determined (e.g. singularity, |
| discontinuity within the interval) one |
| will probably gain from splitting up the |
| interval at this point and calling the |
| integrator on the subranges. if possible, |
| an appropriate special-purpose integrator |
| should be used, which is designed for |
| handling the type of difficulty involved. |
| = 2 the occurrence of roundoff error is detec- |
| ted, which prevents the requested |
| tolerance from being achieved. |
| the error may be under-estimated. |
| = 3 extremely bad integrand behaviour |
| occurs at some points of the integration |
| interval. |
| = 4 the algorithm does not converge. |
| roundoff error is detected in the |
| extrapolation table. it is presumed that |
| the requested tolerance cannot be |
| achieved, and that the returned result is |
| the best which can be obtained. |
| = 5 the integral is probably divergent, or |
| slowly convergent. it must be noted that |
| divergence can occur with any other value |
| of ier. |
| = 6 the input is invalid, because |
| (epsabs <= 0 and |
| epsrel < max(50*rel.mach.acc.,0.5d-28) |
| or limit < 1 or lenw < limit*4. |
| result, abserr, neval, last are set to |
| zero.except when limit or lenw is invalid, |
| iwork(1), work(limit*2+1) and |
| work(limit*3+1) are set to zero, work(1) |
| is set to a and work(limit+1) to b. |
| |
| dimensioning parameters |
| limit - int |
| dimensioning parameter for iwork |
| limit determines the maximum number of subintervals |
| in the partition of the given integration interval |
| (a,b), limit >= 1. |
| if limit < 1, the routine will end with ier = 6. |
| |
| lenw - int |
| dimensioning parameter for work |
| lenw must be at least limit*4. |
| if lenw < limit*4, the routine will end |
| with ier = 6. |
| |
| last - int |
| on return, last equals the number of subintervals |
| produced in the subdivision process, detemines the |
| number of significant elements actually in the work |
| arrays. |
| |
| work arrays |
| iwork - int |
| vector of dimension at least limit, the first k |
| elements of which contain pointers |
| to the error estimates over the subintervals |
| such that work(limit*3+iwork(1)),... , |
| work(limit*3+iwork(k)) form a decreasing |
| sequence, with k = last if last <= (limit/2+2), |
| and k = limit+1-last otherwise |
| |
| work - double precision |
| vector of dimension at least lenw |
| on return |
| work(1), ..., work(last) contain the left |
| end-points of the subintervals in the |
| partition of (a,b), |
| work(limit+1), ..., work(limit+last) contain |
| the right end-points, |
| work(limit*2+1), ..., work(limit*2+last) contain |
| the integral approximations over the subintervals, |
| work(limit*3+1), ..., work(limit*3+last) |
| contain the error estimates. |
| |
| ***routines called dqagse |
| ***end prologue dqags */ |
| |
| /* check validity of limit and lenw. */ |
| |
| *ier = 6; |
| *neval = 0; |
| *last = 0; |
| *result = 0.; |
| *abserr = 0.; |
| if (*limit < 1 || *lenw < *limit *4) return; |
| |
| /* prepare call for dqagse. */ |
| |
| l1 = *limit; |
| l2 = *limit + l1; |
| l3 = *limit + l2; |
| |
| rdqagse(f, ex, a, b, epsabs, epsrel, limit, result, abserr, neval, ier, |
| work, &work[l1], &work[l2], &work[l3], iwork, last); |
| |
| return; |
| } /* rdqags_ */ |
| |
| static |
| void rdqagse(integr_fn f, void *ex, double *a, double *b, double * |
| epsabs, double *epsrel, int *limit, double *result, |
| double *abserr, int *neval, int *ier, double *alist, |
| double *blist, double *rlist, double *elist, int * |
| iord, int *last) |
| { |
| /* Local variables */ |
| Rboolean noext, extrap; |
| int k,ksgn, nres; |
| int ierro; |
| int ktmin, nrmax; |
| int iroff1, iroff2, iroff3; |
| int id; |
| int numrl2; |
| int jupbnd; |
| int maxerr; |
| double res3la[3]; |
| double rlist2[52]; |
| double abseps, area, area1, area2, area12, dres, epmach; |
| double a1, a2, b1, b2, defabs, defab1, defab2, oflow, uflow, resabs, reseps; |
| double error1, error2, erro12, errbnd, erlast, errmax, errsum; |
| |
| double correc = 0.0, erlarg = 0.0, ertest = 0.0, small = 0.0; |
| /* |
| ***begin prologue dqagse |
| ***date written 800101 (yymmdd) |
| ***revision date 830518 (yymmdd) |
| ***category no. h2a1a1 |
| ***keywords automatic integrator, general-purpose, |
| (end point) singularities, extrapolation, |
| globally adaptive |
| ***author piessens,robert,appl. math. & progr. div. - k.u.leuven |
| de doncker,elise,appl. math. & progr. div. - k.u.leuven |
| ***purpose the routine calculates an approximation result to a given |
| definite integral i = integral of f over (a,b), |
| hopefully satisfying following claim for accuracy |
| abs(i-result) <= max(epsabs,epsrel*abs(i)). |
| ***description |
| |
| computation of a definite integral |
| standard fortran subroutine |
| double precision version |
| |
| parameters |
| on entry |
| f - double precision |
| function subprogram defining the integrand |
| function f(x). the actual name for f needs to be |
| declared e x t e r n a l in the driver program. |
| |
| a - double precision |
| lower limit of integration |
| |
| b - double precision |
| upper limit of integration |
| |
| epsabs - double precision |
| absolute accuracy requested |
| epsrel - double precision |
| relative accuracy requested |
| if epsabs <= 0 |
| and epsrel < max(50*rel.mach.acc.,0.5d-28), |
| the routine will end with ier = 6. |
| |
| limit - int |
| gives an upperbound on the number of subintervals |
| in the partition of (a,b) |
| |
| on return |
| result - double precision |
| approximation to the integral |
| |
| abserr - double precision |
| estimate of the modulus of the absolute error, |
| which should equal or exceed abs(i-result) |
| |
| neval - int |
| number of integrand evaluations |
| |
| ier - int |
| 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. |
| error messages |
| = 1 maximum number of subdivisions allowed |
| has been achieved. one can allow more sub- |
| divisions by increasing the value of limit |
| (and taking the according dimension |
| adjustments into account). however, if |
| this yields no improvement it is advised |
| to analyze the integrand in order to |
| determine the integration difficulties. if |
| the position of a local difficulty can be |
| determined (e.g. singularity, |
| discontinuity within the interval) one |
| will probably gain from splitting up the |
| interval at this point and calling the |
| integrator on the subranges. if possible, |
| an appropriate special-purpose integrator |
| should be used, which is designed for |
| handling the type of difficulty involved. |
| = 2 the occurrence of roundoff error is detec- |
| ted, which prevents the requested |
| tolerance from being achieved. |
| the error may be under-estimated. |
| = 3 extremely bad integrand behaviour |
| occurs at some points of the integration |
| interval. |
| = 4 the algorithm does not converge. |
| roundoff error is detected in the |
| extrapolation table. |
| it is presumed that the requested |
| tolerance cannot be achieved, and that the |
| returned result is the best which can be |
| obtained. |
| = 5 the integral is probably divergent, or |
| slowly convergent. it must be noted that |
| divergence can occur with any other value |
| of ier. |
| = 6 the input is invalid, because |
| epsabs <= 0 and |
| epsrel < max(50*rel.mach.acc.,0.5d-28). |
| result, abserr, neval, last, rlist(1), |
| iord(1) and elist(1) are set to zero. |
| alist(1) and blist(1) are set to a and b |
| respectively. |
| |
| alist - double precision |
| vector of dimension at least limit, the first |
| last elements of which are the left end points |
| of the subintervals in the partition of the |
| given integration range (a,b) |
| |
| blist - double precision |
| vector of dimension at least limit, the first |
| last elements of which are the right end points |
| of the subintervals in the partition of the given |
| integration range (a,b) |
| |
| rlist - double precision |
| vector of dimension at least limit, the first |
| last elements of which are the integral |
| approximations on the subintervals |
| |
| elist - double precision |
| vector of dimension at least limit, the first |
| last elements of which are the moduli of the |
| absolute error estimates on the subintervals |
| |
| iord - int |
| vector of dimension at least limit, the first k |
| elements of which are pointers to the |
| error estimates over the subintervals, |
| such that elist(iord(1)), ..., elist(iord(k)) |
| form a decreasing sequence, with k = last |
| if last <= (limit/2+2), and k = limit+1-last |
| otherwise |
| |
| last - int |
| number of subintervals actually produced in the |
| subdivision process |
| |
| ***references (none) |
| ***routines called dqelg,dqk21,dqpsrt |
| ***end prologue dqagse |
| |
| |
| |
| the dimension of rlist2 is determined by the value of |
| limexp in subroutine dqelg (rlist2 should be of dimension |
| (limexp+2) at least). |
| |
| list of major variables |
| ----------------------- |
| |
| alist - list of left end points of all subintervals |
| considered up to now |
| blist - list of right end points of all subintervals |
| considered up to now |
| rlist(i) - approximation to the integral over |
| (alist(i),blist(i)) |
| rlist2 - array of dimension at least limexp+2 containing |
| the part of the epsilon table which is still |
| needed for further computations |
| elist(i) - error estimate applying to rlist(i) |
| maxerr - pointer to the interval with largest error |
| estimate |
| errmax - elist(maxerr) |
| erlast - error on the interval currently subdivided |
| (before that subdivision has taken place) |
| area - sum of the integrals over the subintervals |
| errsum - sum of the errors over the subintervals |
| errbnd - requested accuracy max(epsabs,epsrel* |
| abs(result)) |
| *****1 - variable for the left interval |
| *****2 - variable for the right interval |
| last - index for subdivision |
| nres - number of calls to the extrapolation routine |
| numrl2 - number of elements currently in rlist2. if an |
| appropriate approximation to the compounded |
| integral has been obtained it is put in |
| rlist2(numrl2) after numrl2 has been increased |
| by one. |
| small - length of the smallest interval considered up |
| to now, multiplied by 1.5 |
| erlarg - sum of the errors over the intervals larger |
| than the smallest interval considered up to now |
| extrap - logical variable denoting that the routine is |
| attempting to perform extrapolation i.e. before |
| subdividing the smallest interval we try to |
| decrease the value of erlarg. |
| noext - logical variable denoting that extrapolation |
| is no longer allowed (true value) |
| |
| machine dependent constants |
| --------------------------- |
| |
| epmach is the largest relative spacing. |
| uflow is the smallest positive magnitude. |
| oflow is the largest positive magnitude. */ |
| |
| /* ***first executable statement dqagse */ |
| /* Parameter adjustments */ |
| --iord; |
| --elist; |
| --rlist; |
| --blist; |
| --alist; |
| |
| /* Function Body */ |
| epmach = DBL_EPSILON; |
| |
| /* test on validity of parameters */ |
| /* ------------------------------ */ |
| *ier = 0; |
| *neval = 0; |
| *last = 0; |
| *result = 0.; |
| *abserr = 0.; |
| alist[1] = *a; |
| blist[1] = *b; |
| rlist[1] = 0.; |
| elist[1] = 0.; |
| if (*epsabs <= 0. && *epsrel < fmax2(epmach * 50., 5e-29)) { |
| *ier = 6; |
| return; |
| } |
| |
| /* first approximation to the integral */ |
| /* ----------------------------------- */ |
| |
| uflow = DBL_MIN; |
| oflow = DBL_MAX; |
| ierro = 0; |
| rdqk21(f, ex, a, b, result, abserr, &defabs, &resabs); |
| |
| /* test on accuracy. */ |
| |
| dres = fabs(*result); |
| errbnd = fmax2(*epsabs, *epsrel * dres); |
| *last = 1; |
| rlist[1] = *result; |
| elist[1] = *abserr; |
| iord[1] = 1; |
| if (*abserr <= epmach * 100. * defabs && *abserr > errbnd) |
| *ier = 2; |
| if (*limit == 1) |
| *ier = 1; |
| if (*ier != 0 || (*abserr <= errbnd && *abserr != resabs) |
| || *abserr == 0.) goto L140; |
| |
| /* initialization */ |
| /* -------------- */ |
| |
| rlist2[0] = *result; |
| errmax = *abserr; |
| maxerr = 1; |
| area = *result; |
| errsum = *abserr; |
| *abserr = oflow; |
| nrmax = 1; |
| nres = 0; |
| numrl2 = 2; |
| ktmin = 0; |
| extrap = FALSE; |
| noext = FALSE; |
| iroff1 = 0; |
| iroff2 = 0; |
| iroff3 = 0; |
| ksgn = -1; |
| if (dres >= (1. - epmach * 50.) * defabs) { |
| ksgn = 1; |
| } |
| |
| /* main do-loop */ |
| /* ------------ */ |
| |
| for (*last = 2; *last <= *limit; ++(*last)) { |
| |
| /* bisect the subinterval with the nrmax-th largest error estimate. */ |
| |
| a1 = alist[maxerr]; |
| b1 = (alist[maxerr] + blist[maxerr]) * .5; |
| a2 = b1; |
| b2 = blist[maxerr]; |
| erlast = errmax; |
| rdqk21(f, ex, &a1, &b1, &area1, &error1, &resabs, &defab1); |
| rdqk21(f, ex, &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 || defab2 == error2)) { |
| |
| if (fabs(rlist[maxerr] - area12) <= fabs(area12) * 1e-5 && |
| erro12 >= errmax * .99) { |
| if (extrap) |
| ++iroff2; |
| else /* if(! extrap) */ |
| ++iroff1; |
| } |
| if (*last > 10 && erro12 > errmax) |
| ++iroff3; |
| } |
| rlist[maxerr] = area1; |
| rlist[*last] = area2; |
| errbnd = fmax2(*epsabs, *epsrel * fabs(area)); |
| |
| /* test for roundoff error and eventually set error flag. */ |
| |
| if (iroff1 + iroff2 >= 10 || iroff3 >= 20) |
| *ier = 2; |
| if (iroff2 >= 5) |
| ierro = 3; |
| |
| /* set error flag in the case that the number of subintervals equals limit. */ |
| if (*last == *limit) |
| *ier = 1; |
| |
| /* set error flag in the case of bad integrand behaviour |
| at a point of the integration range. */ |
| |
| if (fmax2(fabs(a1), fabs(b2)) <= |
| (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3)) { |
| *ier = 4; |
| } |
| |
| /* append the newly-created intervals to the list. */ |
| |
| if (error2 > error1) { |
| alist[maxerr] = a2; |
| alist[*last] = a1; |
| blist[*last] = b1; |
| rlist[maxerr] = area2; |
| rlist[*last] = area1; |
| elist[maxerr] = error2; |
| elist[*last] = error1; |
| } else { |
| alist[*last] = a2; |
| blist[maxerr] = b1; |
| blist[*last] = b2; |
| elist[maxerr] = error1; |
| elist[*last] = error2; |
| } |
| |
| /* call subroutine dqpsrt to maintain the descending ordering |
| in the list of error estimates and select the subinterval |
| with nrmax-th largest error estimate (to be bisected next). */ |
| |
| /*L30:*/ |
| rdqpsrt(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax); |
| |
| if (errsum <= errbnd) goto L115;/* ***jump out of do-loop */ |
| if (*ier != 0) break; |
| if (*last == 2) { /* L80: */ |
| small = fabs(*b - *a) * .375; |
| erlarg = errsum; |
| ertest = errbnd; |
| rlist2[1] = area; continue; |
| } |
| if (noext) continue; |
| |
| erlarg -= erlast; |
| if (fabs(b1 - a1) > small) { |
| erlarg += erro12; |
| } |
| if (!extrap) { |
| |
| /* test whether the interval to be bisected next is the |
| smallest interval. */ |
| |
| if (fabs(blist[maxerr] - alist[maxerr]) > small) { |
| continue; |
| } |
| extrap = TRUE; |
| nrmax = 2; |
| } |
| |
| if (ierro != 3 && erlarg > ertest) { |
| |
| /* the smallest interval has the largest error. |
| before bisecting decrease the sum of the errors over the |
| larger intervals (erlarg) and perform extrapolation. */ |
| |
| id = nrmax; |
| jupbnd = *last; |
| if (*last > *limit / 2 + 2) { |
| jupbnd = *limit + 3 - *last; |
| } |
| for (k = id; k <= jupbnd; ++k) { |
| maxerr = iord[nrmax]; |
| errmax = elist[maxerr]; |
| if (fabs(blist[maxerr] - alist[maxerr]) > small) { |
| goto L90; |
| } |
| ++nrmax; |
| /* L50: */ |
| } |
| } |
| /* perform extrapolation. L60: */ |
| |
| ++numrl2; |
| rlist2[numrl2 - 1] = area; |
| rdqelg(&numrl2, rlist2, &reseps, &abseps, res3la, &nres); |
| ++ktmin; |
| if (ktmin > 5 && *abserr < errsum * .001) { |
| *ier = 5; |
| } |
| if (abseps < *abserr) { |
| ktmin = 0; |
| *abserr = abseps; |
| *result = reseps; |
| correc = erlarg; |
| ertest = fmax2(*epsabs, *epsrel * fabs(reseps)); |
| if (*abserr <= ertest) { |
| break; |
| } |
| } |
| |
| /* prepare bisection of the smallest interval. L70: */ |
| |
| if (numrl2 == 1) { |
| noext = TRUE; |
| } |
| if (*ier == 5) { |
| break; |
| } |
| maxerr = iord[1]; |
| errmax = elist[maxerr]; |
| nrmax = 1; |
| extrap = FALSE; |
| small *= .5; |
| erlarg = errsum; |
| L90: |
| ; |
| } |
| |
| |
| /* L100: set final result and error estimate. */ |
| /* ------------------------------------ */ |
| |
| if (*abserr == oflow) goto L115; |
| if (*ier + ierro == 0) goto L110; |
| if (ierro == 3) |
| *abserr += correc; |
| if (*ier == 0) |
| *ier = 3; |
| if (*result == 0. || area == 0.) { |
| if (*abserr > errsum) goto L115; |
| if (area == 0.) goto L130; |
| } |
| else { /* L105:*/ |
| if (*abserr / fabs(*result) > errsum / fabs(area)) |
| goto L115; |
| } |
| |
| L110:/* test on divergence. */ |
| if (ksgn == -1 && fmax2(fabs(*result), fabs(area)) <= defabs * .01) { |
| goto L130; |
| } |
| if (.01 > *result / area || *result / area > 100. || errsum > fabs(area)) { |
| *ier = 5; |
| } |
| goto L130; |
| |
| L115:/* compute global integral sum. */ |
| *result = 0.; |
| for (k = 1; k <= *last; ++k) |
| *result += rlist[k]; |
| *abserr = errsum; |
| L130: |
| if (*ier > 2) |
| L140: |
| *neval = *last * 42 - 21; |
| return; |
| } /* rdqagse_ */ |
| |
| |
| static void rdqk15i(integr_fn f, void *ex, |
| double *boun, int *inf, double *a, double *b, |
| double *result, |
| double *abserr, double *resabs, double *resasc) |
| { |
| /* Initialized data */ |
| |
| static double wg[8] = { |
| 0., .129484966168869693270611432679082, |
| 0., .27970539148927666790146777142378, |
| 0., .381830050505118944950369775488975, |
| 0., .417959183673469387755102040816327 }; |
| static double xgk[8] = { |
| .991455371120812639206854697526329, |
| .949107912342758524526189684047851, |
| .864864423359769072789712788640926, |
| .741531185599394439863864773280788, |
| .58608723546769113029414483825873, |
| .405845151377397166906606412076961, |
| .207784955007898467600689403773245, 0. }; |
| static double wgk[8] = { |
| .02293532201052922496373200805897, |
| .063092092629978553290700663189204, |
| .104790010322250183839876322541518, |
| .140653259715525918745189590510238, |
| .16900472663926790282658342659855, |
| .190350578064785409913256402421014, |
| .204432940075298892414161999234649, |
| .209482141084727828012999174891714 }; |
| |
| /* Local variables */ |
| double absc, dinf, resg, resk, fsum, absc1, absc2, fval1, fval2; |
| int j; |
| double hlgth, centr, reskh, uflow; |
| double tabsc1, tabsc2, fc, epmach; |
| double fv1[7], fv2[7], vec[15], vec2[15]; |
| |
| /* |
| ***begin prologue dqk15i |
| ***date written 800101 (yymmdd) |
| ***revision date 830518 (yymmdd) |
| ***category no. h2a3a2,h2a4a2 |
| ***keywords 15-point transformed gauss-kronrod rules |
| ***author piessens,robert,appl. math. & progr. div. - k.u.leuven |
| de doncker,elise,appl. math. & progr. div. - k.u.leuven |
| ***purpose the original (infinite integration range is mapped |
| onto the interval (0,1) and (a,b) is a part of (0,1). |
| it is the purpose to compute |
| i = integral of transformed integrand over (a,b), |
| j = integral of abs(transformed integrand) over (a,b). |
| ***description |
| |
| integration rule |
| standard fortran subroutine |
| double precision version |
| |
| parameters |
| on entry |
| f - double precision |
| fuction subprogram defining the integrand |
| function f(x). the actual name for f needs to be |
| declared e x t e r n a l in the calling program. |
| |
| boun - double precision |
| finite bound of original integration |
| range (set to zero if inf = +2) |
| |
| inf - int |
| if inf = -1, the original interval is |
| (-infinity,bound), |
| if inf = +1, the original interval is |
| (bound,+infinity), |
| if inf = +2, the original interval is |
| (-infinity,+infinity) and |
| the integral is computed as the sum of two |
| integrals, one over (-infinity,0) and one over |
| (0,+infinity). |
| |
| a - double precision |
| lower limit for integration over subrange |
| of (0,1) |
| |
| b - double precision |
| upper limit for integration over subrange |
| of (0,1) |
| |
| on return |
| result - double precision |
| approximation to the integral i |
| result is computed by applying the 15-point |
| kronrod rule(resk) obtained by optimal addition |
| of abscissae to the 7-point gauss rule(resg). |
| |
| abserr - double precision |
| estimate of the modulus of the absolute error, |
| which should equal or exceed abs(i-result) |
| |
| resabs - double precision |
| approximation to the integral j |
| |
| resasc - double precision |
| approximation to the integral of |
| abs((transformed integrand)-i/(b-a)) over (a,b) |
| |
| ***references (none) |
| ***end prologue dqk15i |
| |
| |
| the abscissae and weights are supplied for the interval |
| (-1,1). because of symmetry only the positive abscissae and |
| their corresponding weights are given. |
| |
| xgk - abscissae of the 15-point kronrod rule |
| xgk(2), xgk(4), ... abscissae of the 7-point |
| gauss rule |
| xgk(1), xgk(3), ... abscissae which are optimally |
| added to the 7-point gauss rule |
| |
| wgk - weights of the 15-point kronrod rule |
| |
| wg - weights of the 7-point gauss rule, corresponding |
| to the abscissae xgk(2), xgk(4), ... |
| wg(1), wg(3), ... are set to zero. |
| |
| |
| |
| |
| |
| list of major variables |
| ----------------------- |
| |
| centr - mid point of the interval |
| hlgth - half-length of the interval |
| absc* - abscissa |
| tabsc* - transformed abscissa |
| fval* - function value |
| resg - result of the 7-point gauss formula |
| resk - result of the 15-point kronrod formula |
| reskh - approximation to the mean value of the transformed |
| integrand over (a,b), i.e. to i/(b-a) |
| |
| machine dependent constants |
| --------------------------- |
| |
| epmach is the largest relative spacing. |
| uflow is the smallest positive magnitude. |
| */ |
| |
| /* ***first executable statement dqk15i */ |
| epmach = DBL_EPSILON; |
| uflow = DBL_MIN; |
| dinf = (double) imin2(1, *inf); |
| |
| centr = (*a + *b) * .5; |
| hlgth = (*b - *a) * .5; |
| tabsc1 = *boun + dinf * (1. - centr) / centr; |
| vec[0] = tabsc1; |
| if (*inf == 2) { |
| vec2[0] = -tabsc1; |
| } |
| for (j = 1; j <= 7; ++j) { |
| absc = hlgth * xgk[j - 1]; |
| absc1 = centr - absc; |
| absc2 = centr + absc; |
| tabsc1 = *boun + dinf * (1. - absc1) / absc1; |
| tabsc2 = *boun + dinf * (1. - absc2) / absc2; |
| vec[(j << 1) - 1] = tabsc1; |
| vec[j * 2] = tabsc2; |
| if (*inf == 2) { |
| vec2[(j << 1) - 1] = -tabsc1; |
| vec2[j * 2] = -tabsc2; |
| } |
| /* L5: */ |
| } |
| f(vec, 15, ex); /* -> new vec[] overwriting old vec[] */ |
| if (*inf == 2) f(vec2, 15, ex); |
| fval1 = vec[0]; |
| if (*inf == 2) fval1 += vec2[0]; |
| fc = fval1 / centr / centr; |
| |
| /* compute the 15-point kronrod approximation to |
| the integral, and estimate the error. */ |
| |
| resg = wg[7] * fc; |
| resk = wgk[7] * fc; |
| *resabs = fabs(resk); |
| for (j = 1; j <= 7; ++j) { |
| absc = hlgth * xgk[j - 1]; |
| absc1 = centr - absc; |
| absc2 = centr + absc; |
| tabsc1 = *boun + dinf * (1. - absc1) / absc1; |
| tabsc2 = *boun + dinf * (1. - absc2) / absc2; |
| fval1 = vec[(j << 1) - 1]; |
| fval2 = vec[j * 2]; |
| if (*inf == 2) { |
| fval1 += vec2[(j << 1) - 1]; |
| } |
| if (*inf == 2) { |
| fval2 += vec2[j * 2]; |
| } |
| fval1 = fval1 / absc1 / absc1; |
| fval2 = fval2 / absc2 / absc2; |
| fv1[j - 1] = fval1; |
| fv2[j - 1] = fval2; |
| fsum = fval1 + fval2; |
| resg += wg[j - 1] * fsum; |
| resk += wgk[j - 1] * fsum; |
| *resabs += wgk[j - 1] * (fabs(fval1) + fabs(fval2)); |
| /* L10: */ |
| } |
| reskh = resk * .5; |
| *resasc = wgk[7] * fabs(fc - reskh); |
| for (j = 1; j <= 7; ++j) { |
| *resasc += wgk[j - 1] * (fabs(fv1[j - 1] - reskh) + |
| fabs(fv2[j - 1] - reskh)); |
| /* L20: */ |
| } |
| *result = resk * hlgth; |
| *resasc *= hlgth; |
| *resabs *= hlgth; |
| *abserr = fabs((resk - resg) * hlgth); |
| if (*resasc != 0. && *abserr != 0.) { |
| *abserr = *resasc * fmin2(1., pow(*abserr * 200. / *resasc, 1.5)); |
| } |
| if (*resabs > uflow / (epmach * 50.)) { |
| *abserr = fmax2(epmach * 50. * *resabs, *abserr); |
| } |
| return; |
| } /* rdqk15i_ */ |
| |
| static void rdqelg(int *n, double *epstab, double * |
| result, double *abserr, double *res3la, int *nres) |
| { |
| /* Local variables */ |
| int i__, indx, ib, ib2, ie, k1, k2, k3, num, newelm, limexp; |
| double delta1, delta2, delta3, e0, e1, e1abs, e2, e3, epmach, epsinf; |
| double oflow, ss, res; |
| double errA, err1, err2, err3, tol1, tol2, tol3; |
| |
| /* ***begin prologue dqelg |
| ***refer to dqagie,dqagoe,dqagpe,dqagse |
| ***revision date 830518 (yymmdd) |
| ***keywords epsilon algorithm, convergence acceleration, |
| extrapolation |
| ***author piessens,robert,appl. math. & progr. div. - k.u.leuven |
| de doncker,elise,appl. math & progr. div. - k.u.leuven |
| ***purpose the routine determines the limit of a given sequence of |
| approximations, by means of the epsilon algorithm of |
| p.wynn. an estimate of the absolute error is also given. |
| the condensed epsilon table is computed. only those |
| elements needed for the computation of the next diagonal |
| are preserved. |
| ***description |
| |
| epsilon algorithm |
| standard fortran subroutine |
| double precision version |
| |
| parameters |
| n - int |
| epstab(n) contains the new element in the |
| first column of the epsilon table. |
| |
| epstab - double precision |
| vector of dimension 52 containing the elements |
| of the two lower diagonals of the triangular |
| epsilon table. the elements are numbered |
| starting at the right-hand corner of the |
| triangle. |
| |
| result - double precision |
| resulting approximation to the integral |
| |
| abserr - double precision |
| estimate of the absolute error computed from |
| result and the 3 previous results |
| |
| res3la - double precision |
| vector of dimension 3 containing the last 3 |
| results |
| |
| nres - int |
| number of calls to the routine |
| (should be zero at first call) |
| |
| ***end prologue dqelg |
| |
| |
| list of major variables |
| ----------------------- |
| |
| e0 - the 4 elements on which the computation of a new |
| e1 element in the epsilon table is based |
| e2 |
| e3 e0 |
| e3 e1 new |
| e2 |
| |
| newelm - number of elements to be computed in the new diagonal |
| errA - errA = abs(e1-e0)+abs(e2-e1)+abs(new-e2) |
| result - the element in the new diagonal with least value of errA |
| |
| machine dependent constants |
| --------------------------- |
| |
| epmach is the largest relative spacing. |
| oflow is the largest positive magnitude. |
| limexp is the maximum number of elements the epsilon |
| table can contain. if this number is reached, the upper |
| diagonal of the epsilon table is deleted. */ |
| |
| /* ***first executable statement dqelg */ |
| /* Parameter adjustments */ |
| --res3la; |
| --epstab; |
| |
| /* Function Body */ |
| epmach = DBL_EPSILON; |
| oflow = DBL_MAX; |
| ++(*nres); |
| *abserr = oflow; |
| *result = epstab[*n]; |
| if (*n < 3) { |
| goto L100; |
| } |
| limexp = 50; |
| epstab[*n + 2] = epstab[*n]; |
| newelm = (*n - 1) / 2; |
| epstab[*n] = oflow; |
| num = *n; |
| k1 = *n; |
| for (i__ = 1; i__ <= newelm; ++i__) { |
| k2 = k1 - 1; |
| k3 = k1 - 2; |
| res = epstab[k1 + 2]; |
| e0 = epstab[k3]; |
| e1 = epstab[k2]; |
| e2 = res; |
| e1abs = fabs(e1); |
| delta2 = e2 - e1; |
| err2 = fabs(delta2); |
| tol2 = fmax2(fabs(e2), e1abs) * epmach; |
| delta3 = e1 - e0; |
| err3 = fabs(delta3); |
| tol3 = fmax2(e1abs, fabs(e0)) * epmach; |
| if (err2 <= tol2 && err3 <= tol3) { |
| /* if e0, e1 and e2 are equal to within machine |
| accuracy, convergence is assumed. */ |
| *result = res;/* result = e2 */ |
| *abserr = err2 + err3;/* abserr = fabs(e1-e0)+fabs(e2-e1) */ |
| |
| goto L100; /* ***jump out of do-loop */ |
| } |
| |
| e3 = epstab[k1]; |
| epstab[k1] = e1; |
| delta1 = e1 - e3; |
| err1 = fabs(delta1); |
| tol1 = fmax2(e1abs, fabs(e3)) * epmach; |
| |
| /* if two elements are very close to each other, omit |
| a part of the table by adjusting the value of n */ |
| |
| if (err1 > tol1 && err2 > tol2 && err3 > tol3) { |
| ss = 1. / delta1 + 1. / delta2 - 1. / delta3; |
| epsinf = fabs(ss * e1); |
| |
| /* test to detect irregular behaviour in the table, and |
| eventually omit a part of the table adjusting the value of n. */ |
| |
| if (epsinf > 1e-4) { |
| goto L30; |
| } |
| } |
| |
| *n = i__ + i__ - 1; |
| goto L50;/* ***jump out of do-loop */ |
| |
| |
| L30:/* compute a new element and eventually adjust the value of result. */ |
| |
| res = e1 + 1. / ss; |
| epstab[k1] = res; |
| k1 += -2; |
| errA = err2 + fabs(res - e2) + err3; |
| if (errA <= *abserr) { |
| *abserr = errA; |
| *result = res; |
| } |
| } |
| |
| /* shift the table. */ |
| |
| L50: |
| if (*n == limexp) { |
| *n = (limexp / 2 << 1) - 1; |
| } |
| |
| if (num / 2 << 1 == num) ib = 2; else ib = 1; |
| ie = newelm + 1; |
| for (i__ = 1; i__ <= ie; ++i__) { |
| ib2 = ib + 2; |
| epstab[ib] = epstab[ib2]; |
| ib = ib2; |
| } |
| if (num != *n) { |
| indx = num - *n + 1; |
| for (i__ = 1; i__ <= *n; ++i__) { |
| epstab[i__] = epstab[indx]; |
| ++indx; |
| } |
| } |
| /*L80:*/ |
| if (*nres >= 4) { |
| /* L90: */ |
| *abserr = fabs(*result - res3la[3]) + |
| fabs(*result - res3la[2]) + |
| fabs(*result - res3la[1]); |
| res3la[1] = res3la[2]; |
| res3la[2] = res3la[3]; |
| res3la[3] = *result; |
| } else { |
| res3la[*nres] = *result; |
| *abserr = oflow; |
| } |
| |
| L100:/* compute error estimate */ |
| *abserr = fmax2(*abserr, epmach * 5. * fabs(*result)); |
| return; |
| } /* rdqelg_ */ |
| |
| static void rdqk21(integr_fn f, void *ex, double *a, double *b, double *result, |
| double *abserr, double *resabs, double *resasc) |
| { |
| /* Initialized data */ |
| |
| static double wg[5] = { .066671344308688137593568809893332, |
| .149451349150580593145776339657697, |
| .219086362515982043995534934228163, |
| .269266719309996355091226921569469, |
| .295524224714752870173892994651338 }; |
| static double xgk[11] = { .995657163025808080735527280689003, |
| .973906528517171720077964012084452, |
| .930157491355708226001207180059508, |
| .865063366688984510732096688423493, |
| .780817726586416897063717578345042, |
| .679409568299024406234327365114874, |
| .562757134668604683339000099272694, |
| .433395394129247190799265943165784, |
| .294392862701460198131126603103866, |
| .14887433898163121088482600112972,0. }; |
| static double wgk[11] = { .011694638867371874278064396062192, |
| .03255816230796472747881897245939, |
| .05475589657435199603138130024458, |
| .07503967481091995276704314091619, |
| .093125454583697605535065465083366, |
| .109387158802297641899210590325805, |
| .123491976262065851077958109831074, |
| .134709217311473325928054001771707, |
| .142775938577060080797094273138717, |
| .147739104901338491374841515972068, |
| .149445554002916905664936468389821 }; |
| |
| |
| /* Local variables */ |
| double fv1[10], fv2[10], vec[21]; |
| double absc, resg, resk, fsum, fval1, fval2; |
| double hlgth, centr, reskh, uflow; |
| double fc, epmach, dhlgth; |
| int j, jtw, jtwm1; |
| |
| /* ***begin prologue dqk21 |
| ***date written 800101 (yymmdd) |
| ***revision date 830518 (yymmdd) |
| ***category no. h2a1a2 |
| ***keywords 21-point gauss-kronrod rules |
| ***author piessens,robert,appl. math. & progr. div. - k.u.leuven |
| de doncker,elise,appl. math. & progr. div. - k.u.leuven |
| ***purpose to compute i = integral of f over (a,b), with error |
| estimate |
| j = integral of abs(f) over (a,b) |
| ***description |
| |
| integration rules |
| standard fortran subroutine |
| double precision version |
| |
| parameters |
| on entry |
| f - double precision |
| function subprogram defining the integrand |
| function f(x). the actual name for f needs to be |
| declared e x t e r n a l in the driver program. |
| |
| a - double precision |
| lower limit of integration |
| |
| b - double precision |
| upper limit of integration |
| |
| on return |
| result - double precision |
| approximation to the integral i |
| result is computed by applying the 21-point |
| kronrod rule (resk) obtained by optimal addition |
| of abscissae to the 10-point gauss rule (resg). |
| |
| abserr - double precision |
| estimate of the modulus of the absolute error, |
| which should not exceed abs(i-result) |
| |
| resabs - double precision |
| approximation to the integral j |
| |
| resasc - double precision |
| approximation to the integral of abs(f-i/(b-a)) |
| over (a,b) |
| |
| ***references (none) |
| ***end prologue dqk21 |
| |
| |
| |
| the abscissae and weights are given for the interval (-1,1). |
| because of symmetry only the positive abscissae and their |
| corresponding weights are given. |
| |
| xgk - abscissae of the 21-point kronrod rule |
| xgk(2), xgk(4), ... abscissae of the 10-point |
| gauss rule |
| xgk(1), xgk(3), ... abscissae which are optimally |
| added to the 10-point gauss rule |
| |
| wgk - weights of the 21-point kronrod rule |
| |
| wg - weights of the 10-point gauss rule |
| |
| |
| gauss quadrature weights and kronron quadrature abscissae and weights |
| as evaluated with 80 decimal digit arithmetic by l. w. fullerton, |
| bell labs, nov. 1981. |
| |
| |
| |
| |
| |
| list of major variables |
| ----------------------- |
| |
| centr - mid point of the interval |
| hlgth - half-length of the interval |
| absc - abscissa |
| fval* - function value |
| resg - result of the 10-point gauss formula |
| resk - result of the 21-point kronrod formula |
| reskh - approximation to the mean value of f over (a,b), |
| i.e. to i/(b-a) |
| |
| |
| machine dependent constants |
| --------------------------- |
| |
| epmach is the largest relative spacing. |
| uflow is the smallest positive magnitude. */ |
| |
| /* ***first executable statement dqk21 */ |
| epmach = DBL_EPSILON; |
| uflow = DBL_MIN; |
| |
| centr = (*a + *b) * .5; |
| hlgth = (*b - *a) * .5; |
| dhlgth = fabs(hlgth); |
| |
| /* compute the 21-point kronrod approximation to |
| the integral, and estimate the absolute error. */ |
| |
| resg = 0.; |
| vec[0] = centr; |
| for (j = 1; j <= 5; ++j) { |
| jtw = j << 1; |
| absc = hlgth * xgk[jtw - 1]; |
| vec[(j << 1) - 1] = centr - absc; |
| /* L5: */ |
| vec[j * 2] = centr + absc; |
| } |
| for (j = 1; j <= 5; ++j) { |
| jtwm1 = (j << 1) - 1; |
| absc = hlgth * xgk[jtwm1 - 1]; |
| vec[(j << 1) + 9] = centr - absc; |
| vec[(j << 1) + 10] = centr + absc; |
| } |
| f(vec, 21, ex); |
| fc = vec[0]; |
| resk = wgk[10] * fc; |
| *resabs = fabs(resk); |
| for (j = 1; j <= 5; ++j) { |
| jtw = j << 1; |
| absc = hlgth * xgk[jtw - 1]; |
| fval1 = vec[(j << 1) - 1]; |
| fval2 = vec[j * 2]; |
| fv1[jtw - 1] = fval1; |
| fv2[jtw - 1] = fval2; |
| fsum = fval1 + fval2; |
| resg += wg[j - 1] * fsum; |
| resk += wgk[jtw - 1] * fsum; |
| *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2)); |
| /* L10: */ |
| } |
| for (j = 1; j <= 5; ++j) { |
| jtwm1 = (j << 1) - 1; |
| absc = hlgth * xgk[jtwm1 - 1]; |
| fval1 = vec[(j << 1) + 9]; |
| fval2 = vec[(j << 1) + 10]; |
| fv1[jtwm1 - 1] = fval1; |
| fv2[jtwm1 - 1] = fval2; |
| fsum = fval1 + fval2; |
| resk += wgk[jtwm1 - 1] * fsum; |
| *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2)); |
| /* L15: */ |
| } |
| reskh = resk * .5; |
| *resasc = wgk[10] * fabs(fc - reskh); |
| for (j = 1; j <= 10; ++j) { |
| *resasc += wgk[j - 1] * (fabs(fv1[j - 1] - reskh) + |
| fabs(fv2[j - 1] - reskh)); |
| /* L20: */ |
| } |
| *result = resk * hlgth; |
| *resabs *= dhlgth; |
| *resasc *= dhlgth; |
| *abserr = fabs((resk - resg) * hlgth); |
| if (*resasc != 0. && *abserr != 0.) { |
| *abserr = *resasc * fmin2(1., pow(*abserr * 200. / *resasc, 1.5)); |
| } |
| if (*resabs > uflow / (epmach * 50.)) { |
| *abserr = fmax2(epmach * 50. * *resabs, *abserr); |
| } |
| return; |
| } /* rdqk21_ */ |
| |
| static void rdqpsrt(int *limit, int *last, int *maxerr, |
| double *ermax, double *elist, int *iord, int *nrmax) |
| { |
| /* Local variables */ |
| int i, j, k, ido, jbnd, isucc, jupbn; |
| double errmin, errmax; |
| |
| /* ***begin prologue dqpsrt |
| ***refer to dqage,dqagie,dqagpe,dqawse |
| ***routines called (none) |
| ***revision date 810101 (yymmdd) |
| ***keywords sequential sorting |
| ***author piessens,robert,appl. math. & progr. div. - k.u.leuven |
| de doncker,elise,appl. math. & progr. div. - k.u.leuven |
| ***purpose this routine maintains the descending ordering in the |
| list of the local error estimated resulting from the |
| interval subdivision process. at each call two error |
| estimates are inserted using the sequential search |
| method, top-down for the largest error estimate and |
| bottom-up for the smallest error estimate. |
| ***description |
| |
| ordering routine |
| standard fortran subroutine |
| double precision version |
| |
| parameters (meaning at output) |
| limit - int |
| maximum number of error estimates the list |
| can contain |
| |
| last - int |
| number of error estimates currently in the list |
| |
| maxerr - int |
| maxerr points to the nrmax-th largest error |
| estimate currently in the list |
| |
| ermax - double precision |
| nrmax-th largest error estimate |
| ermax = elist(maxerr) |
| |
| elist - double precision |
| vector of dimension last containing |
| the error estimates |
| |
| iord - int |
| vector of dimension last, the first k elements |
| of which contain pointers to the error |
| estimates, such that |
| elist(iord(1)),..., elist(iord(k)) |
| form a decreasing sequence, with |
| k = last if last <= (limit/2+2), and |
| k = limit+1-last otherwise |
| |
| nrmax - int |
| maxerr = iord(nrmax) |
| |
| ***end prologue dqpsrt |
| */ |
| |
| |
| /* Parameter adjustments */ |
| --iord; |
| --elist; |
| |
| /* Function Body */ |
| |
| /* check whether the list contains more than |
| two error estimates. */ |
| if (*last <= 2) { |
| iord[1] = 1; |
| iord[2] = 2; |
| goto Last; |
| } |
| /* this part of the routine is only executed if, due to a |
| difficult integrand, subdivision increased the error |
| estimate. in the normal case the insert procedure should |
| start after the nrmax-th largest error estimate. */ |
| |
| errmax = elist[*maxerr]; |
| if (*nrmax > 1) { |
| ido = *nrmax - 1; |
| for (i = 1; i <= ido; ++i) { |
| isucc = iord[*nrmax - 1]; |
| if (errmax <= elist[isucc]) |
| break; /* out of for-loop */ |
| iord[*nrmax] = isucc; |
| --(*nrmax); |
| /* L20: */ |
| } |
| } |
| |
| /*L30: compute the number of elements in the list to be maintained |
| in descending order. this number depends on the number of |
| subdivisions still allowed. */ |
| if (*last > *limit / 2 + 2) |
| jupbn = *limit + 3 - *last; |
| else |
| jupbn = *last; |
| |
| errmin = elist[*last]; |
| |
| /* insert errmax by traversing the list top-down, |
| starting comparison from the element elist(iord(nrmax+1)). */ |
| |
| jbnd = jupbn - 1; |
| for (i = *nrmax + 1; i <= jbnd; ++i) { |
| isucc = iord[i]; |
| if (errmax >= elist[isucc]) {/* ***jump out of do-loop */ |
| /* L60: insert errmin by traversing the list bottom-up. */ |
| iord[i - 1] = *maxerr; |
| for (j = i, k = jbnd; j <= jbnd; j++, k--) { |
| isucc = iord[k]; |
| if (errmin < elist[isucc]) { |
| /* goto L80; ***jump out of do-loop */ |
| iord[k + 1] = *last; |
| goto Last; |
| } |
| iord[k + 1] = isucc; |
| } |
| iord[i] = *last; |
| goto Last; |
| } |
| iord[i - 1] = isucc; |
| } |
| |
| iord[jbnd] = *maxerr; |
| iord[jupbn] = *last; |
| |
| Last:/* set maxerr and ermax. */ |
| |
| *maxerr = iord[*nrmax]; |
| *ermax = elist[*maxerr]; |
| return; |
| } /* rdqpsrt_ */ |
