1 // This code has been mechanically translated from the original FORTRAN 2 // code at http://netlib.org/quadpack. 3 4 /** Authors: Lars Tandle Kyllingstad 5 Copyright: Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved. 6 License: Boost License 1.0 7 */ 8 module scid.ports.quadpack.qage; 9 10 11 import std.algorithm: max; 12 import std.math: fabs; 13 14 import scid.ports.quadpack.qk15; 15 import scid.ports.quadpack.qk21; 16 import scid.ports.quadpack.qk31; 17 import scid.ports.quadpack.qk41; 18 import scid.ports.quadpack.qk51; 19 import scid.ports.quadpack.qk61; 20 import scid.ports.quadpack.qpsrt; 21 22 import scid.core.fortran; 23 24 25 26 27 /// 28 void qage(Real, Func)(Func f, Real a, Real b, Real epsabs, Real epsrel, 29 int key, int limit, out Real result, out Real abserr, out int neval, 30 out int ier, Real* alist_, Real* blist_, Real* rlist_, Real* elist_, 31 int* iord_, out int last) 32 { 33 //***begin prologue dqage 34 //***date written 800101 (yymmdd) 35 //***revision date 830518 (yymmdd) 36 //***category no. h2a1a1 37 //***keywords automatic integrator, general-purpose, 38 // integrand examinator, globally adaptive, 39 // gauss-kronrod 40 //***author piessens,robert,appl. math. & progr. div. - k.u.leuven 41 // de doncker,elise,appl. math. & progr. div. - k.u.leuven 42 //***purpose the routine calculates an approximation result to a given 43 // definite integral i = integral of f over (a,b), 44 // hopefully satisfying following claim for accuracy 45 // abs(i-reslt).le.max(epsabs,epsrel*abs(i)). 46 //***description 47 // 48 // computation of a definite integral 49 // standard fortran subroutine 50 // double precision version 51 // 52 // parameters 53 // on entry 54 // f - double precision 55 // function subprogram defining the integrand 56 // function f(x). the actual name for f needs to be 57 // declared e x t e r n a l in the driver program. 58 // 59 // a - double precision 60 // lower limit of integration 61 // 62 // b - double precision 63 // upper limit of integration 64 // 65 // epsabs - double precision 66 // absolute accuracy requested 67 // epsrel - double precision 68 // relative accuracy requested 69 // if epsabs.le.0 70 // and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 71 // the routine will end with ier = 6. 72 // 73 // key - integer 74 // key for choice of local integration rule 75 // a gauss-kronrod pair is used with 76 // 7 - 15 points if key.lt.2, 77 // 10 - 21 points if key = 2, 78 // 15 - 31 points if key = 3, 79 // 20 - 41 points if key = 4, 80 // 25 - 51 points if key = 5, 81 // 30 - 61 points if key.gt.5. 82 // 83 // limit - integer 84 // gives an upperbound on the number of subintervals 85 // in the partition of (a,b), limit.ge.1. 86 // 87 // on return 88 // result - double precision 89 // approximation to the integral 90 // 91 // abserr - double precision 92 // estimate of the modulus of the absolute error, 93 // which should equal or exceed abs(i-result) 94 // 95 // neval - integer 96 // number of integrand evaluations 97 // 98 // ier - integer 99 // ier = 0 normal and reliable termination of the 100 // routine. it is assumed that the requested 101 // accuracy has been achieved. 102 // ier.gt.0 abnormal termination of the routine 103 // the estimates for result and error are 104 // less reliable. it is assumed that the 105 // requested accuracy has not been achieved. 106 // error messages 107 // ier = 1 maximum number of subdivisions allowed 108 // has been achieved. one can allow more 109 // subdivisions by increasing the value 110 // of limit. 111 // however, if this yields no improvement it 112 // is rather advised to analyze the integrand 113 // in order to determine the integration 114 // difficulties. if the position of a local 115 // difficulty can be determined(e.g. 116 // singularity, discontinuity within the 117 // interval) one will probably gain from 118 // splitting up the interval at this point 119 // and calling the integrator on the 120 // subranges. if possible, an appropriate 121 // special-purpose integrator should be used 122 // which is designed for handling the type of 123 // difficulty involved. 124 // = 2 the occurrence of roundoff error is 125 // detected, which prevents the requested 126 // tolerance from being achieved. 127 // = 3 extremely bad integrand behaviour occurs 128 // at some points of the integration 129 // interval. 130 // = 6 the input is invalid, because 131 // (epsabs.le.0 and 132 // epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 133 // result, abserr, neval, last, rlist(1) , 134 // elist(1) and iord(1) are set to zero. 135 // alist(1) and blist(1) are set to a and b 136 // respectively. 137 // 138 // alist - double precision 139 // vector of dimension at least limit, the first 140 // last elements of which are the left 141 // end points of the subintervals in the partition 142 // of the given integration range (a,b) 143 // 144 // blist - double precision 145 // vector of dimension at least limit, the first 146 // last elements of which are the right 147 // end points of the subintervals in the partition 148 // of the given integration range (a,b) 149 // 150 // rlist - double precision 151 // vector of dimension at least limit, the first 152 // last elements of which are the 153 // integral approximations on the subintervals 154 // 155 // elist - double precision 156 // vector of dimension at least limit, the first 157 // last elements of which are the moduli of the 158 // absolute error estimates on the subintervals 159 // 160 // iord - integer 161 // vector of dimension at least limit, the first k 162 // elements of which are pointers to the 163 // error estimates over the subintervals, 164 // such that elist(iord(1)), ..., 165 // elist(iord(k)) form a decreasing sequence, 166 // with k = last if last.le.(limit/2+2), and 167 // k = limit+1-last otherwise 168 // 169 // last - integer 170 // number of subintervals actually produced in the 171 // subdivision process 172 // 173 //***references (none) 174 //***routines called d1mach,dqk15,dqk21,dqk31, 175 // dqk41,dqk51,dqk61,dqpsrt 176 //***end prologue dqage 177 // 178 Real area,area1,area12,area2,a1,a2, 179 b1,b2,defabs,defab1,defab2,epmach, 180 errbnd,errmax,error1,error2,erro12,errsum, 181 resabs,uflow; 182 int iroff1=1,iroff2=1,k=1,keyf=1,maxerr=1, nrmax; 183 // 184 auto alist = dimension(alist_, limit); 185 auto blist = dimension(blist_, limit); 186 auto elist = dimension(elist_, limit); 187 auto iord = dimension(iord_, limit); 188 auto rlist = dimension(rlist_, limit); 189 // 190 // list of major variables 191 // ----------------------- 192 // 193 // alist - list of left end points of all subintervals 194 // considered up to now 195 // blist - list of right end points of all subintervals 196 // considered up to now 197 // rlist(i) - approximation to the integral over 198 // (alist(i),blist(i)) 199 // elist(i) - error estimate applying to rlist(i) 200 // maxerr - pointer to the interval with largest 201 // error estimate 202 // errmax - elist(maxerr) 203 // area - sum of the integrals over the subintervals 204 // errsum - sum of the errors over the subintervals 205 // errbnd - requested accuracy max(epsabs,epsrel* 206 // abs(result)) 207 // *****1 - variable for the left subinterval 208 // *****2 - variable for the right subinterval 209 // last - index for subdivision 210 // 211 // 212 // machine dependent constants 213 // --------------------------- 214 // 215 // epmach is the largest relative spacing. 216 // uflow is the smallest positive magnitude. 217 // 218 //***first executable statement dqage 219 epmach = Real.epsilon; 220 uflow = Real.epsilon; 221 // 222 // test on validity of parameters 223 // ------------------------------ 224 // 225 ier = 0; 226 neval = 0; 227 last = 0; 228 result = 0.0; 229 abserr = 0.0; 230 alist[1] = a; 231 blist[1] = b; 232 rlist[1] = 0.0; 233 elist[1] = 0.0; 234 iord[1] = 0; 235 if(epsabs <= 0.0 && 236 epsrel < max(0.5e2*epmach,0.5e-28)) ier = 6; 237 if(ier == 6) goto l999; 238 // 239 // first approximation to the integral 240 // ----------------------------------- 241 // 242 keyf = key; 243 if(key <= 0) keyf = 1; 244 if(key >= 7) keyf = 6; 245 neval = 0; 246 if(keyf == 1) qk15!(Real, Func)(f,a,b,result,abserr,defabs,resabs); 247 else if(keyf == 2) qk21!(Real, Func)(f,a,b,result,abserr,defabs,resabs); 248 else if(keyf == 3) qk31!(Real, Func)(f,a,b,result,abserr,defabs,resabs); 249 else if(keyf == 4) qk41!(Real, Func)(f,a,b,result,abserr,defabs,resabs); 250 else if(keyf == 5) qk51!(Real, Func)(f,a,b,result,abserr,defabs,resabs); 251 else if(keyf == 6) qk61!(Real, Func)(f,a,b,result,abserr,defabs,resabs); 252 last = 1; 253 rlist[1] = result; 254 elist[1] = abserr; 255 iord[1] = 1; 256 // 257 // test on accuracy. 258 // 259 errbnd = max(epsabs,epsrel*fabs(result)); 260 if(abserr <= 0.5e2*epmach*defabs && abserr > errbnd) ier = 2; 261 if(limit == 1) ier = 1; 262 if(ier != 0 || (abserr <= errbnd && abserr != resabs) 263 || abserr == 0.0) goto l60; 264 // 265 // initialization 266 // -------------- 267 // 268 // 269 errmax = abserr; 270 maxerr = 1; 271 area = result; 272 errsum = abserr; 273 nrmax = 1; 274 iroff1 = 0; 275 iroff2 = 0; 276 // 277 // main do-loop 278 // ------------ 279 // 280 for (last=2; last<=limit; last++) { //do 30 last = 2,limit 281 // 282 // bisect the subinterval with the largest error estimate. 283 // 284 a1 = alist[maxerr]; 285 b1 = 0.5*(alist[maxerr]+blist[maxerr]); 286 a2 = b1; 287 b2 = blist[maxerr]; 288 if(keyf == 1) qk15!(Real, Func)(f,a1,b1,area1,error1,resabs,defab1); 289 else if(keyf == 2) qk21!(Real, Func)(f,a1,b1,area1,error1,resabs,defab1); 290 else if(keyf == 3) qk31!(Real, Func)(f,a1,b1,area1,error1,resabs,defab1); 291 else if(keyf == 4) qk41!(Real, Func)(f,a1,b1,area1,error1,resabs,defab1); 292 else if(keyf == 5) qk51!(Real, Func)(f,a1,b1,area1,error1,resabs,defab1); 293 else if(keyf == 6) qk61!(Real, Func)(f,a1,b1,area1,error1,resabs,defab1); 294 if(keyf == 1) qk15!(Real, Func)(f,a2,b2,area2,error2,resabs,defab2); 295 else if(keyf == 2) qk21!(Real, Func)(f,a2,b2,area2,error2,resabs,defab2); 296 else if(keyf == 3) qk31!(Real, Func)(f,a2,b2,area2,error2,resabs,defab2); 297 else if(keyf == 4) qk41!(Real, Func)(f,a2,b2,area2,error2,resabs,defab2); 298 else if(keyf == 5) qk51!(Real, Func)(f,a2,b2,area2,error2,resabs,defab2); 299 else if(keyf == 6) qk61!(Real, Func)(f,a2,b2,area2,error2,resabs,defab2); 300 // 301 // improve previous approximations to integral 302 // and error and test for accuracy. 303 // 304 neval = neval+1; 305 area12 = area1+area2; 306 erro12 = error1+error2; 307 errsum = errsum+erro12-errmax; 308 area = area+area12-rlist[maxerr]; 309 if(defab1 == error1 || defab2 == error2) goto l5; 310 if(fabs(rlist[maxerr]-area12) <= 0.1e-4*fabs(area12) 311 && erro12 >= 0.99*errmax) iroff1 = iroff1+1; 312 if(last > 10 && erro12 > errmax) iroff2 = iroff2+1; 313 l5: rlist[maxerr] = area1; 314 rlist[last] = area2; 315 errbnd = max(epsabs,epsrel*fabs(area)); 316 if(errsum <= errbnd) goto l8; 317 // 318 // test for roundoff error and eventually set error flag. 319 // 320 if(iroff1 >= 6 || iroff2 >= 20) ier = 2; 321 // 322 // set error flag in the case that the number of subintervals 323 // equals limit. 324 // 325 if(last == limit) ier = 1; 326 // 327 // set error flag in the case of bad integrand behaviour 328 // at a point of the integration range. 329 // 330 if(max(fabs(a1),fabs(b2)) <= (0.1e1+0.1e3* 331 epmach)*(fabs(a2)+0.1e4*uflow)) ier = 3; 332 // 333 // append the newly-created intervals to the list. 334 // 335 l8: if(error2 > error1) goto l10; 336 alist[last] = a2; 337 blist[maxerr] = b1; 338 blist[last] = b2; 339 elist[maxerr] = error1; 340 elist[last] = error2; 341 goto l20; 342 l10: alist[maxerr] = a2; 343 alist[last] = a1; 344 blist[last] = b1; 345 rlist[maxerr] = area2; 346 rlist[last] = area1; 347 elist[maxerr] = error2; 348 elist[last] = error1; 349 // 350 // call subroutine dqpsrt to maintain the descending ordering 351 // in the list of error estimates and select the subinterval 352 // with the largest error estimate (to be bisected next). 353 // 354 l20: qpsrt!Real(limit,last,maxerr,errmax,elist.ptr,iord.ptr,nrmax); 355 // ***jump out of do-loop 356 if(ier != 0 || errsum <= errbnd) goto l40; 357 l30: ;} 358 // 359 // compute final result. 360 // --------------------- 361 // 362 l40: result = 0.0; 363 for (k=1; k<=last; k++) { //do 50 k=1,last 364 result = result+rlist[k]; 365 l50: ;} 366 abserr = errsum; 367 l60: if(keyf != 1) neval = (10*keyf+1)*(2*neval+1); 368 else neval = 30*neval+15; 369 l999: return; 370 } 371 372 unittest 373 { 374 alias qage!(float, float delegate(float)) fqage; 375 alias qage!(double, double delegate(double)) dqage; 376 alias qage!(double, double function(double)) dfqage; 377 alias qage!(real, real delegate(real)) rqage; 378 }