1 subroutine dqagp(f,a,b,npts2,points,epsabs,epsrel,result,abserr, 2 * neval,ier,leniw,lenw,last,iwork,work) 3c***begin prologue dqagp 4c***date written 800101 (yymmdd) 5c***revision date 830518 (yymmdd) 6c***category no. h2a2a1 7c***keywords automatic integrator, general-purpose, 8c singularities at user specified points, 9c extrapolation, globally adaptive 10c***author piessens,robert,appl. math. & progr. div - k.u.leuven 11c de doncker,elise,appl. math. & progr. div. - k.u.leuven 12c***purpose the routine calculates an approximation result to a given 13c definite integral i = integral of f over (a,b), 14c hopefully satisfying following claim for accuracy 15c break points of the integration interval, where local 16c difficulties of the integrand may occur (e.g. 17c singularities, discontinuities), are provided by the user. 18c***description 19c 20c computation of a definite integral 21c standard fortran subroutine 22c double precision version 23c 24c parameters 25c on entry 26c f - double precision 27c function subprogram defining the integrand 28c function f(x). the actual name for f needs to be 29c declared e x t e r n a l in the driver program. 30c 31c a - double precision 32c lower limit of integration 33c 34c b - double precision 35c upper limit of integration 36c 37c npts2 - integer 38c number equal to two more than the number of 39c user-supplied break points within the integration 40c range, npts.ge.2. 41c if npts2.lt.2, the routine will end with ier = 6. 42c 43c points - double precision 44c vector of dimension npts2, the first (npts2-2) 45c elements of which are the user provided break 46c points. if these points do not constitute an 47c ascending sequence there will be an automatic 48c sorting. 49c 50c epsabs - double precision 51c absolute accuracy requested 52c epsrel - double precision 53c relative accuracy requested 54c if epsabs.le.0 55c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 56c the routine will end with ier = 6. 57c 58c on return 59c result - double precision 60c approximation to the integral 61c 62c abserr - double precision 63c estimate of the modulus of the absolute error, 64c which should equal or exceed abs(i-result) 65c 66c neval - integer 67c number of integrand evaluations 68c 69c ier - integer 70c ier = 0 normal and reliable termination of the 71c routine. it is assumed that the requested 72c accuracy has been achieved. 73c ier.gt.0 abnormal termination of the routine. 74c the estimates for integral and error are 75c less reliable. it is assumed that the 76c requested accuracy has not been achieved. 77c error messages 78c ier = 1 maximum number of subdivisions allowed 79c has been achieved. one can allow more 80c subdivisions by increasing the value of 81c limit (and taking the according dimension 82c adjustments into account). however, if 83c this yields no improvement it is advised 84c to analyze the integrand in order to 85c determine the integration difficulties. if 86c the position of a local difficulty can be 87c determined (i.e. singularity, 88c discontinuity within the interval), it 89c should be supplied to the routine as an 90c element of the vector points. if necessary 91c an appropriate special-purpose integrator 92c must be used, which is designed for 93c handling the type of difficulty involved. 94c = 2 the occurrence of roundoff error is 95c detected, which prevents the requested 96c tolerance from being achieved. 97c the error may be under-estimated. 98c = 3 extremely bad integrand behaviour occurs 99c at some points of the integration 100c interval. 101c = 4 the algorithm does not converge. 102c roundoff error is detected in the 103c extrapolation table. 104c it is presumed that the requested 105c tolerance cannot be achieved, and that 106c the returned result is the best which 107c can be obtained. 108c = 5 the integral is probably divergent, or 109c slowly convergent. it must be noted that 110c divergence can occur with any other value 111c of ier.gt.0. 112c = 6 the input is invalid because 113c npts2.lt.2 or 114c break points are specified outside 115c the integration range or 116c (epsabs.le.0 and 117c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 118c result, abserr, neval, last are set to 119c zero. except when leniw or lenw or npts2 is 120c invalid, iwork(1), iwork(limit+1), 121c work(limit*2+1) and work(limit*3+1) 122c are set to zero. 123c work(1) is set to a and work(limit+1) 124c to b (where limit = (leniw-npts2)/2). 125c 126c dimensioning parameters 127c leniw - integer 128c dimensioning parameter for iwork 129c leniw determines limit = (leniw-npts2)/2, 130c which is the maximum number of subintervals in the 131c partition of the given integration interval (a,b), 132c leniw.ge.(3*npts2-2). 133c if leniw.lt.(3*npts2-2), the routine will end with 134c ier = 6. 135c 136c lenw - integer 137c dimensioning parameter for work 138c lenw must be at least leniw*2-npts2. 139c if lenw.lt.leniw*2-npts2, the routine will end 140c with ier = 6. 141c 142c last - integer 143c on return, last equals the number of subintervals 144c produced in the subdivision process, which 145c determines the number of significant elements 146c actually in the work arrays. 147c 148c work arrays 149c iwork - integer 150c vector of dimension at least leniw. on return, 151c the first k elements of which contain 152c pointers to the error estimates over the 153c subintervals, such that work(limit*3+iwork(1)),..., 154c work(limit*3+iwork(k)) form a decreasing 155c sequence, with k = last if last.le.(limit/2+2), and 156c k = limit+1-last otherwise 157c iwork(limit+1), ...,iwork(limit+last) contain the 158c subdivision levels of the subintervals, i.e. 159c if (aa,bb) is a subinterval of (p1,p2) 160c where p1 as well as p2 is a user-provided 161c break point or integration limit, then (aa,bb) has 162c level l if abs(bb-aa) = abs(p2-p1)*2**(-l), 163c iwork(limit*2+1), ..., iwork(limit*2+npts2) have 164c no significance for the user, 165c note that limit = (leniw-npts2)/2. 166c 167c work - double precision 168c vector of dimension at least lenw 169c on return 170c work(1), ..., work(last) contain the left 171c end points of the subintervals in the 172c partition of (a,b), 173c work(limit+1), ..., work(limit+last) contain 174c the right end points, 175c work(limit*2+1), ..., work(limit*2+last) contain 176c the integral approximations over the subintervals, 177c work(limit*3+1), ..., work(limit*3+last) 178c contain the corresponding error estimates, 179c work(limit*4+1), ..., work(limit*4+npts2) 180c contain the integration limits and the 181c break points sorted in an ascending sequence. 182c note that limit = (leniw-npts2)/2. 183c 184c***references (none) 185c***routines called dqagpe,xerror 186c***end prologue dqagp 187c 188 double precision a,abserr,b,epsabs,epsrel,f,points,result,work 189 integer ier,iwork,last,leniw,lenw,limit,lvl,l1,l2,l3,l4,neval, 190 * npts2 191c 192 dimension iwork(leniw),points(npts2),work(lenw) 193c 194 external f 195c 196c check validity of limit and lenw. 197c 198c***first executable statement dqagp 199 ier = 6 200 neval = 0 201 last = 0 202 result = 0.0d+00 203 abserr = 0.0d+00 204 if(leniw.lt.(3*npts2-2).or.lenw.lt.(leniw*2-npts2).or.npts2.lt.2) 205 * go to 10 206c 207c prepare call for dqagpe. 208c 209 limit = (leniw-npts2)/2 210 l1 = limit+1 211 l2 = limit+l1 212 l3 = limit+l2 213 l4 = limit+l3 214c 215 call dqagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,abserr, 216 * neval,ier,work(1),work(l1),work(l2),work(l3),work(l4), 217 * iwork(1),iwork(l1),iwork(l2),last) 218c 219c call error handler if necessary. 220c 221 lvl = 0 22210 if(ier.eq.6) lvl = 1 223 if(ier.ne.0) call xerror('abnormal return from dqagp',26,ier,lvl) 224 return 225 end 226