subroutine qagp(f,a,b,npts2,points,epsabs,epsrel,result,abserr, * neval,ier,leniw,lenw,last,iwork,work) c***begin prologue qagp c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a1 c***keywords automatic integrator, general-purpose, c singularities at user specified points, c extrapolation, globally adaptive c***author piessens,robert,appl. math. & progr. div - k.u.leuven c de doncker,elise,appl. math. & progr. div. - k.u.leuven c***purpose the routine calculates an approximation result to a given c definite integral i = integral of f over (a,b), c hopefully satisfying following claim for accuracy c break points of the integration interval, where local c difficulties of the integrand may occur(e.g. singularities, c discontinuities), are provided by the user. c***description c c computation of a definite integral c standard fortran subroutine c real version c c parameters c on entry c f - real c function subprogram defining the integrand c function f(x). the actual name for f needs to be c declared e x t e r n a l in the driver program. c c a - real c lower limit of integration c c b - real c upper limit of integration c c npts2 - integer c number equal to two more than the number of c user-supplied break points within the integration c range, npts.ge.2. c if npts2.lt.2, the routine will end with ier = 6. c c points - real c vector of dimension npts2, the first (npts2-2) c elements of which are the user provided break c points. if these points do not constitute an c ascending sequence there will be an automatic c sorting. c c epsabs - real c absolute accuracy requested c epsrel - real c relative accuracy requested c if epsabs.le.0 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), c the routine will end with ier = 6. c c on return c result - real c approximation to the integral c c abserr - real c estimate of the modulus of the absolute error, c which should equal or exceed abs(i-result) c c neval - integer c number of integrand evaluations c c ier - integer c ier = 0 normal and reliable termination of the c routine. it is assumed that the requested c accuracy has been achieved. c ier.gt.0 abnormal termination of the routine. c the estimates for integral and error are c less reliable. it is assumed that the c requested accuracy has not been achieved. c error messages c ier = 1 maximum number of subdivisions allowed c has been achieved. one can allow more c subdivisions by increasing the value of c limit (and taking the according dimension c adjustments into account). however, if c this yields no improvement it is advised c to analyze the integrand in order to c determine the integration difficulties. if c the position of a local difficulty can be c determined (i.e. singularity, c discontinuity within the interval), it c should be supplied to the routine as an c element of the vector points. if necessary c an appropriate special-purpose integrator c must be used, which is designed for c handling the type of difficulty involved. c = 2 the occurrence of roundoff error is c detected, which prevents the requested c tolerance from being achieved. c the error may be under-estimated. c = 3 extremely bad integrand behaviour occurs c at some points of the integration c interval. c = 4 the algorithm does not converge. c roundoff error is detected in the c extrapolation table. c it is presumed that the requested c tolerance cannot be achieved, and that c the returned result is the best which c can be obtained. c = 5 the integral is probably divergent, or c slowly convergent. it must be noted that c divergence can occur with any other value c of ier.gt.0. c = 6 the input is invalid because c npts2.lt.2 or c break points are specified outside c the integration range or c (epsabs.le.0 and c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) c result, abserr, neval, last are set to c zero. exept when leniw or lenw or npts2 is c invalid, iwork(1), iwork(limit+1), c work(limit*2+1) and work(limit*3+1) c are set to zero. c work(1) is set to a and work(limit+1) c to b (where limit = (leniw-npts2)/2). c c dimensioning parameters c leniw - integer c dimensioning parameter for iwork c leniw determines limit = (leniw-npts2)/2, c which is the maximum number of subintervals in the c partition of the given integration interval (a,b), c leniw.ge.(3*npts2-2). c if leniw.lt.(3*npts2-2), the routine will end with c ier = 6. c c lenw - integer c dimensioning parameter for work c lenw must be at least leniw*2-npts2. c if lenw.lt.leniw*2-npts2, the routine will end c with ier = 6. c c last - integer c on return, last equals the number of subintervals c produced in the subdivision process, which c determines the number of significant elements c actually in the work arrays. c c work arrays c iwork - integer c vector of dimension at least leniw. on return, c the first k elements of which contain c pointers to the error estimates over the c subintervals, such that work(limit*3+iwork(1)),..., c work(limit*3+iwork(k)) form a decreasing c sequence, with k = last if last.le.(limit/2+2), and c k = limit+1-last otherwise c iwork(limit+1), ...,iwork(limit+last) contain the c subdivision levels of the subintervals, i.e. c if (aa,bb) is a subinterval of (p1,p2) c where p1 as well as p2 is a user-provided c break point or integration limit, then (aa,bb) has c level l if abs(bb-aa) = abs(p2-p1)*2**(-l), c iwork(limit*2+1), ..., iwork(limit*2+npts2) have c no significance for the user, c note that limit = (leniw-npts2)/2. c c work - real c vector of dimension at least lenw c on return c work(1), ..., work(last) contain the left c end points of the subintervals in the c partition of (a,b), c work(limit+1), ..., work(limit+last) contain c the right end points, c work(limit*2+1), ..., work(limit*2+last) contain c the integral approximations over the subintervals, c work(limit*3+1), ..., work(limit*3+last) c contain the corresponding error estimates, c work(limit*4+1), ..., work(limit*4+npts2) c contain the integration limits and the c break points sorted in an ascending sequence. c note that limit = (leniw-npts2)/2. c c***references (none) c***routines called qagpe,xerror c***end prologue qagp c real a,abserr,b,epsabs,epsrel,f,points,result,work integer ier,iwork,leniw,lenw,limit,lvl,l1,l2,l3,neval,npts2 c dimension iwork(leniw),points(npts2),work(lenw) c external f c c check validity of limit and lenw. c c***first executable statement qagp ier = 6 neval = 0 last = 0 result = 0.0e+00 abserr = 0.0e+00 if(leniw.lt.(3*npts2-2).or.lenw.lt.(leniw*2-npts2).or.npts2.lt.2) * go to 10 c c prepare call for qagpe. c limit = (leniw-npts2)/2 l1 = limit+1 l2 = limit+l1 l3 = limit+l2 l4 = limit+l3 c call qagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,abserr, * neval,ier,work(1),work(l1),work(l2),work(l3),work(l4), * iwork(1),iwork(l1),iwork(l2),last) c c call error handler if necessary. c lvl = 0 10 if(ier.eq.6) lvl = 1 if(ier.ne.0) call xerror(26habnormal return from qagp, * 26,ier,lvl) return end