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