00001 subroutine qagp(f,a,b,npts2,points,epsabs,epsrel,result,abserr, 00002 * neval,ier,leniw,lenw,last,iwork,work) 00003 c***begin prologue qagp 00004 c***date written 800101 (yymmdd) 00005 c***revision date 830518 (yymmdd) 00006 c***category no. h2a2a1 00007 c***keywords automatic integrator, general-purpose, 00008 c singularities at user specified points, 00009 c extrapolation, globally adaptive 00010 c***author piessens,robert,appl. math. & progr. div - k.u.leuven 00011 c de doncker,elise,appl. math. & progr. div. - k.u.leuven 00012 c***purpose the routine calculates an approximation result to a given 00013 c definite integral i = integral of f over (a,b), 00014 c hopefully satisfying following claim for accuracy 00015 c break points of the integration interval, where local 00016 c difficulties of the integrand may occur(e.g. singularities, 00017 c discontinuities), are provided by the user. 00018 c***description 00019 c 00020 c computation of a definite integral 00021 c standard fortran subroutine 00022 c real version 00023 c 00024 c parameters 00025 c on entry 00026 c f - subroutine f(x,ierr,result) defining the integrand 00027 c function f(x). the actual name for f needs to be 00028 c declared e x t e r n a l in the driver program. 00029 c 00030 c a - real 00031 c lower limit of integration 00032 c 00033 c b - real 00034 c upper limit of integration 00035 c 00036 c npts2 - integer 00037 c number equal to two more than the number of 00038 c user-supplied break points within the integration 00039 c range, npts.ge.2. 00040 c if npts2.lt.2, the routine will end with ier = 6. 00041 c 00042 c points - real 00043 c vector of dimension npts2, the first (npts2-2) 00044 c elements of which are the user provided break 00045 c points. if these points do not constitute an 00046 c ascending sequence there will be an automatic 00047 c sorting. 00048 c 00049 c epsabs - real 00050 c absolute accuracy requested 00051 c epsrel - real 00052 c relative accuracy requested 00053 c if epsabs.le.0 00054 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 00055 c the routine will end with ier = 6. 00056 c 00057 c on return 00058 c result - real 00059 c approximation to the integral 00060 c 00061 c abserr - real 00062 c estimate of the modulus of the absolute error, 00063 c which should equal or exceed abs(i-result) 00064 c 00065 c neval - integer 00066 c number of integrand evaluations 00067 c 00068 c ier - integer 00069 c ier = 0 normal and reliable termination of the 00070 c routine. it is assumed that the requested 00071 c accuracy has been achieved. 00072 c ier.gt.0 abnormal termination of the routine. 00073 c the estimates for integral and error are 00074 c less reliable. it is assumed that the 00075 c requested accuracy has not been achieved. 00076 c error messages 00077 c ier = 1 maximum number of subdivisions allowed 00078 c has been achieved. one can allow more 00079 c subdivisions by increasing the value of 00080 c limit (and taking the according dimension 00081 c adjustments into account). however, if 00082 c this yields no improvement it is advised 00083 c to analyze the integrand in order to 00084 c determine the integration difficulties. if 00085 c the position of a local difficulty can be 00086 c determined (i.e. singularity, 00087 c discontinuity within the interval), it 00088 c should be supplied to the routine as an 00089 c element of the vector points. if necessary 00090 c an appropriate special-purpose integrator 00091 c must be used, which is designed for 00092 c handling the type of difficulty involved. 00093 c = 2 the occurrence of roundoff error is 00094 c detected, which prevents the requested 00095 c tolerance from being achieved. 00096 c the error may be under-estimated. 00097 c = 3 extremely bad integrand behaviour occurs 00098 c at some points of the integration 00099 c interval. 00100 c = 4 the algorithm does not converge. 00101 c roundoff error is detected in the 00102 c extrapolation table. 00103 c it is presumed that the requested 00104 c tolerance cannot be achieved, and that 00105 c the returned result is the best which 00106 c can be obtained. 00107 c = 5 the integral is probably divergent, or 00108 c slowly convergent. it must be noted that 00109 c divergence can occur with any other value 00110 c of ier.gt.0. 00111 c = 6 the input is invalid because 00112 c npts2.lt.2 or 00113 c break points are specified outside 00114 c the integration range or 00115 c (epsabs.le.0 and 00116 c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 00117 c result, abserr, neval, last are set to 00118 c zero. exept when leniw or lenw or npts2 is 00119 c invalid, iwork(1), iwork(limit+1), 00120 c work(limit*2+1) and work(limit*3+1) 00121 c are set to zero. 00122 c work(1) is set to a and work(limit+1) 00123 c to b (where limit = (leniw-npts2)/2). 00124 c 00125 c dimensioning parameters 00126 c leniw - integer 00127 c dimensioning parameter for iwork 00128 c leniw determines limit = (leniw-npts2)/2, 00129 c which is the maximum number of subintervals in the 00130 c partition of the given integration interval (a,b), 00131 c leniw.ge.(3*npts2-2). 00132 c if leniw.lt.(3*npts2-2), the routine will end with 00133 c ier = 6. 00134 c 00135 c lenw - integer 00136 c dimensioning parameter for work 00137 c lenw must be at least leniw*2-npts2. 00138 c if lenw.lt.leniw*2-npts2, the routine will end 00139 c with ier = 6. 00140 c 00141 c last - integer 00142 c on return, last equals the number of subintervals 00143 c produced in the subdivision process, which 00144 c determines the number of significant elements 00145 c actually in the work arrays. 00146 c 00147 c work arrays 00148 c iwork - integer 00149 c vector of dimension at least leniw. on return, 00150 c the first k elements of which contain 00151 c pointers to the error estimates over the 00152 c subintervals, such that work(limit*3+iwork(1)),..., 00153 c work(limit*3+iwork(k)) form a decreasing 00154 c sequence, with k = last if last.le.(limit/2+2), and 00155 c k = limit+1-last otherwise 00156 c iwork(limit+1), ...,iwork(limit+last) contain the 00157 c subdivision levels of the subintervals, i.e. 00158 c if (aa,bb) is a subinterval of (p1,p2) 00159 c where p1 as well as p2 is a user-provided 00160 c break point or integration limit, then (aa,bb) has 00161 c level l if abs(bb-aa) = abs(p2-p1)*2**(-l), 00162 c iwork(limit*2+1), ..., iwork(limit*2+npts2) have 00163 c no significance for the user, 00164 c note that limit = (leniw-npts2)/2. 00165 c 00166 c work - real 00167 c vector of dimension at least lenw 00168 c on return 00169 c work(1), ..., work(last) contain the left 00170 c end points of the subintervals in the 00171 c partition of (a,b), 00172 c work(limit+1), ..., work(limit+last) contain 00173 c the right end points, 00174 c work(limit*2+1), ..., work(limit*2+last) contain 00175 c the integral approximations over the subintervals, 00176 c work(limit*3+1), ..., work(limit*3+last) 00177 c contain the corresponding error estimates, 00178 c work(limit*4+1), ..., work(limit*4+npts2) 00179 c contain the integration limits and the 00180 c break points sorted in an ascending sequence. 00181 c note that limit = (leniw-npts2)/2. 00182 c 00183 c***references (none) 00184 c***routines called qagpe,xerror 00185 c***end prologue qagp 00186 c 00187 real a,abserr,b,epsabs,epsrel,points,result,work 00188 integer ier,iwork,leniw,lenw,limit,lvl,l1,l2,l3,neval,npts2 00189 c 00190 dimension iwork(leniw),points(npts2),work(lenw) 00191 c 00192 external f 00193 c 00194 c check validity of limit and lenw. 00195 c 00196 c***first executable statement qagp 00197 ier = 6 00198 neval = 0 00199 last = 0 00200 result = 0.0e+00 00201 abserr = 0.0e+00 00202 if(leniw.lt.(3*npts2-2).or.lenw.lt.(leniw*2-npts2).or.npts2.lt.2) 00203 * go to 10 00204 c 00205 c prepare call for qagpe. 00206 c 00207 limit = (leniw-npts2)/2 00208 l1 = limit+1 00209 l2 = limit+l1 00210 l3 = limit+l2 00211 l4 = limit+l3 00212 c 00213 call qagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,abserr, 00214 * neval,ier,work(1),work(l1),work(l2),work(l3),work(l4), 00215 * iwork(1),iwork(l1),iwork(l2),last) 00216 c 00217 c call error handler if necessary. 00218 c 00219 lvl = 0 00220 10 if(ier.eq.6) lvl = 1 00221 if(ier.ne.0) call xerror('abnormal return from qagp',26,ier,lvl) 00222 return 00223 end
1.7.1
