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