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.qc25f; 9 10 11 import std.math; 12 13 import scid.core.fortran; 14 import scid.ports.linpack.gtsl; 15 import scid.ports.quadpack.qk15w; 16 import scid.ports.quadpack.qwgtf; 17 import scid.ports.quadpack.qcheb; 18 19 version(unittest) import scid.core.testing; 20 21 22 23 24 /// 25 void qc25f(Real, Func)(Func f, Real a, Real b, Real omega, int integr, 26 int nrmom, int maxp1, int ksave, out Real result, 27 out Real abserr, out int neval, out Real resabs, out Real resasc, 28 out int momcom, Real* chebmo_) 29 { 30 //***begin prologue dqc25f 31 //***date written 810101 (yymmdd) 32 //***revision date 830518 (yymmdd) 33 //***category no. h2a2a2 34 //***keywords integration rules for functions with cos or sin 35 // factor, clenshaw-curtis, gauss-kronrod 36 //***author piessens,robert,appl. math. & progr. div. - k.u.leuven 37 // de doncker,elise,appl. math. & progr. div. - k.u.leuven 38 //***purpose to compute the integral i=integral of f(x) over (a,b) 39 // where w(x) = cos(omega*x) or w(x)=sin(omega*x) and to 40 // compute j = integral of abs(f) over (a,b). for small value 41 // of omega or small intervals (a,b) the 15-point gauss-kronro 42 // rule is used. otherwise a generalized clenshaw-curtis 43 // method is used. 44 //***description 45 // 46 // integration rules for functions with cos or sin factor 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 55 // be declared e x t e r n a l in the calling program. 56 // 57 // a - double precision 58 // lower limit of integration 59 // 60 // b - double precision 61 // upper limit of integration 62 // 63 // omega - double precision 64 // parameter in the weight function 65 // 66 // integr - integer 67 // indicates which weight function is to be used 68 // integr = 1 w(x) = cos(omega*x) 69 // integr = 2 w(x) = sin(omega*x) 70 // 71 // nrmom - integer 72 // the length of interval (a,b) is equal to the length 73 // of the original integration interval divided by 74 // 2**nrmom (we suppose that the routine is used in an 75 // adaptive integration process, otherwise set 76 // nrmom = 0). nrmom must be zero at the first call. 77 // 78 // maxp1 - integer 79 // gives an upper bound on the number of chebyshev 80 // moments which can be stored, i.e. for the 81 // intervals of lengths abs(bb-aa)*2**(-l), 82 // l = 0,1,2, ..., maxp1-2. 83 // 84 // ksave - integer 85 // key which is one when the moments for the 86 // current interval have been computed 87 // 88 // on return 89 // result - double precision 90 // approximation to the integral i 91 // 92 // abserr - double precision 93 // estimate of the modulus of the absolute 94 // error, which should equal or exceed abs(i-result) 95 // 96 // neval - integer 97 // number of integrand evaluations 98 // 99 // resabs - double precision 100 // approximation to the integral j 101 // 102 // resasc - double precision 103 // approximation to the integral of abs(f-i/(b-a)) 104 // 105 // on entry and return 106 // momcom - integer 107 // for each interval length we need to compute the 108 // chebyshev moments. momcom counts the number of 109 // intervals for which these moments have already been 110 // computed. if nrmom.lt.momcom or ksave = 1, the 111 // chebyshev moments for the interval (a,b) have 112 // already been computed and stored, otherwise we 113 // compute them and we increase momcom. 114 // 115 // chebmo - double precision 116 // array of dimension at least (maxp1,25) containing 117 // the modified chebyshev moments for the first momcom 118 // momcom interval lengths 119 // 120 // ...................................................................... 121 //***references (none) 122 //***routines called d1mach,dgtsl,dqcheb,dqk15w,dqwgtf 123 //***end prologue dqc25f 124 // 125 Real ac,an,an2,as,asap,ass,centr, 126 conc,cons,cospar, 127 estc,ests,hlgth,oflow,parint,par2,par22, 128 p2,p3,p4,resc12,resc24,ress12,ress24, 129 sinpar; 130 int i=1,iers=1,isym=1,j=1,k=1,m=1, 131 noequ=1,noeq1=1; 132 // 133 Real[13] cheb12_; 134 Real[25] cheb24_, d_, d1_, d2_, fval_; 135 Real[28] v_; 136 auto chebmo = dimension(chebmo_, maxp1, 25); 137 auto cheb12 = dimension(cheb12_.ptr, 13); 138 auto cheb24 = dimension(cheb24_.ptr, 25); 139 auto d = dimension(d_.ptr, 25); 140 auto d1 = dimension(d1_.ptr, 25); 141 auto d2 = dimension(d2_.ptr, 25); 142 auto fval = dimension(fval_.ptr, 25); 143 auto v = dimension(v_.ptr, 28); 144 // 145 // the vector x contains the values cos(k*pi/24) 146 // k = 1, ...,11, to be used for the chebyshev expansion of f 147 // 148 static immutable Real[11] x_ = [ 149 0.9914448613_7381041114_4557526928_563, 150 0.9659258262_8906828674_9743199728_897, 151 0.9238795325_1128675612_8183189396_788, 152 0.8660254037_8443864676_3723170752_936, 153 0.7933533402_9123516457_9776961501_299, 154 0.7071067811_8654752440_0844362104_849, 155 0.6087614290_0872063941_6097542898_164, 156 0.5000000000_0000000000_0000000000_000, 157 0.3826834323_6508977172_8459984030_399, 158 0.2588190451_0252076234_8898837624_048, 159 0.1305261922_2005159154_8406227895_489 160 ]; 161 auto x = dimension(x_.ptr, 11); 162 // 163 // list of major variables 164 // ----------------------- 165 // 166 // centr - mid point of the integration interval 167 // hlgth - half-length of the integration interval 168 // fval - value of the function f at the points 169 // (b-a)*0.5*cos(k*pi/12) + (b+a)*0.5, k = 0, ..., 24 170 // cheb12 - coefficients of the chebyshev series expansion 171 // of degree 12, for the function f, in the 172 // interval (a,b) 173 // cheb24 - coefficients of the chebyshev series expansion 174 // of degree 24, for the function f, in the 175 // interval (a,b) 176 // resc12 - approximation to the integral of 177 // cos(0.5*(b-a)*omega*x)*f(0.5*(b-a)*x+0.5*(b+a)) 178 // over (-1,+1), using the chebyshev series 179 // expansion of degree 12 180 // resc24 - approximation to the same integral, using the 181 // chebyshev series expansion of degree 24 182 // ress12 - the analogue of resc12 for the sine 183 // ress24 - the analogue of resc24 for the sine 184 // 185 // 186 // machine dependent constant 187 // -------------------------- 188 // 189 // oflow is the largest positive magnitude. 190 // 191 //***first executable statement dqc25f 192 oflow = Real.max; 193 // 194 centr = 0.5*(b+a); 195 hlgth = 0.5*(b-a); 196 parint = omega*hlgth; 197 // 198 // compute the integral using the 15-point gauss-kronrod 199 // formula if the value of the parameter in the integrand 200 // is small. 201 // 202 if (fabs(parint) > 0.2e1) goto l10; 203 qk15w!(Real, Func)(f,&qwgtf!Real,omega,p2,p3,p4,integr,a,b,result, 204 abserr,resabs,resasc); 205 neval = 15; 206 goto l170; 207 // 208 // compute the integral using the generalized clenshaw- 209 // curtis method. 210 // 211 l10: conc = hlgth*cos(centr*omega); 212 cons = hlgth*sin(centr*omega); 213 resasc = oflow; 214 neval = 25; 215 // 216 // check whether the chebyshev moments for this interval 217 // have already been computed. 218 // 219 if (nrmom < momcom || ksave == 1) goto l120; 220 // 221 // compute a new set of chebyshev moments. 222 // 223 m = momcom+1; 224 par2 = parint*parint; 225 par22 = par2+0.2e1; 226 sinpar = sin(parint); 227 cospar = cos(parint); 228 // 229 // compute the chebyshev moments with respect to cosine. 230 // 231 v[1] = 0.2e1*sinpar/parint; 232 v[2] = (0.8e1*cospar+(par2+par2-0.8e1)*sinpar/parint)/par2; 233 v[3] = (0.32e2*(par2-0.12e2)*cospar+(0.2e1* 234 ((par2-0.80e2)*par2+0.192e3)*sinpar)/parint)/(par2*par2); 235 ac = 0.8e1*cospar; 236 as = 0.24e2*parint*sinpar; 237 if (fabs(parint) > 0.24e2) goto l30; 238 // 239 // compute the chebyshev moments as the solutions of a 240 // boundary value problem with 1 initial value (v(3)) and 1 241 // end value (computed using an asymptotic formula). 242 // 243 noequ = 25; 244 noeq1 = noequ-1; 245 an = 0.6e1; 246 for (k=1; k<=noeq1; k++) { // end: l20 247 an2 = an*an; 248 d[k] = -0.2e1*(an2-0.4e1)*(par22-an2-an2); 249 d2[k] = (an-0.1e1)*(an-0.2e1)*par2; 250 d1[k+1] = (an+0.3e1)*(an+0.4e1)*par2; 251 v[k+3] = as-(an2-0.4e1)*ac; 252 an = an+0.2e1; 253 l20:;} 254 an2 = an*an; 255 d[noequ] = -0.2e1*(an2-0.4e1)*(par22-an2-an2); 256 v[noequ+3] = as-(an2-0.4e1)*ac; 257 v[4] = v[4]-0.56e2*par2*v[3]; 258 ass = parint*sinpar; 259 asap = (((((0.210e3*par2-0.1e1)*cospar-(0.105e3*par2 260 -0.63e2)*ass)/an2-(0.1e1-0.15e2*par2)*cospar 261 +0.15e2*ass)/an2-cospar+0.3e1*ass)/an2-cospar)/an2; 262 v[noequ+3] = v[noequ+3]-0.2e1*asap*par2*(an-0.1e1)* 263 (an-0.2e1); 264 // 265 // solve the tridiagonal system by means of gaussian 266 // elimination with partial pivoting. 267 // 268 //*** call to dgtsl must be replaced by call to 269 //*** double precision version of linpack routine sgtsl 270 // 271 gtsl!Real(noequ,d1.ptr,d.ptr,d2.ptr,v.ptr+3,iers); 272 goto l50; 273 // 274 // compute the chebyshev moments by means of forward 275 // recursion. 276 // 277 l30: an = 0.4e1; 278 for (i=4; i<=13; i++) { // end: l40 279 an2 = an*an; 280 v[i] = ((an2-0.4e1)*(0.2e1*(par22-an2-an2)*v[i-1]-ac) 281 +as-par2*(an+0.1e1)*(an+0.2e1)*v[i-2])/ 282 (par2*(an-0.1e1)*(an-0.2e1)); 283 an = an+0.2e1; 284 l40:;} 285 l50: for (j=1; j<=13; j++) { // end: l60 286 chebmo[m,2*j-1] = v[j]; 287 l60:;} 288 // 289 // compute the chebyshev moments with respect to sine. 290 // 291 v[1] = 0.2e1*(sinpar-parint*cospar)/par2; 292 v[2] = (0.18e2-0.48e2/par2)*sinpar/par2 293 +(-0.2e1+0.48e2/par2)*cospar/parint; 294 ac = -0.24e2*parint*cospar; 295 as = -0.8e1*sinpar; 296 if (fabs(parint) > 0.24e2) goto l80; 297 // 298 // compute the chebyshev moments as the solutions of a boundary 299 // value problem with 1 initial value (v(2)) and 1 end value 300 // (computed using an asymptotic formula). 301 // 302 an = 0.5e1; 303 for (k=1 ; k<=noeq1; k++) { // end: l70 304 an2 = an*an; 305 d[k] = -0.2e1*(an2-0.4e1)*(par22-an2-an2); 306 d2[k] = (an-0.1e1)*(an-0.2e1)*par2; 307 d1[k+1] = (an+0.3e1)*(an+0.4e1)*par2; 308 v[k+2] = ac+(an2-0.4e1)*as; 309 an = an+0.2e1; 310 l70:;} 311 an2 = an*an; 312 d[noequ] = -0.2e1*(an2-0.4e1)*(par22-an2-an2); 313 v[noequ+2] = ac+(an2-0.4e1)*as; 314 v[3] = v[3]-0.42e2*par2*v[2]; 315 ass = parint*cospar; 316 asap = (((((0.105e3*par2-0.63e2)*ass+(0.210e3*par2 317 -0.1e1)*sinpar)/an2+(0.15e2*par2-0.1e1)*sinpar- 318 0.15e2*ass)/an2-0.3e1*ass-sinpar)/an2-sinpar)/an2; 319 v[noequ+2] = v[noequ+2]-0.2e1*asap*par2*(an-0.1e1) 320 *(an-0.2e1); 321 // 322 // solve the tridiagonal system by means of gaussian 323 // elimination with partial pivoting. 324 // 325 //*** call to dgtsl must be replaced by call to 326 //*** double precision version of linpack routine sgtsl 327 // 328 gtsl!Real(noequ,d1.ptr,d.ptr,d2.ptr,v.ptr+2,iers); 329 goto l100; 330 // 331 // compute the chebyshev moments by means of forward recursion. 332 // 333 l80: an = 0.3e1; 334 for (i=3; i<=12; i++) { // end: l90 335 an2 = an*an; 336 v[i] = ((an2-0.4e1)*(0.2e1*(par22-an2-an2)*v[i-1]+as) 337 +ac-par2*(an+0.1e1)*(an+0.2e1)*v[i-2]) 338 /(par2*(an-0.1e1)*(an-0.2e1)); 339 an = an+0.2e1; 340 l90:;} 341 l100: for (j=1; j<=12; j++) { // end: l110 342 chebmo[m,2*j] = v[j]; 343 l110:;} 344 l120: if (nrmom < momcom) m = nrmom+1; 345 if (momcom < (maxp1-1) && nrmom >= momcom) momcom = momcom+1; 346 // 347 // compute the coefficients of the chebyshev expansions 348 // of degrees 12 and 24 of the function f. 349 // 350 fval[1] = 0.5*f(centr+hlgth); 351 fval[13] = f(centr); 352 fval[25] = 0.5*f(centr-hlgth); 353 for(i=2; i<=12; i++) { // end: l130 354 isym = 26-i; 355 fval[i] = f(hlgth*x[i-1]+centr); 356 fval[isym] = f(centr-hlgth*x[i-1]); 357 l130:;} 358 qcheb!Real(x.ptr,fval.ptr,cheb12.ptr,cheb24.ptr); 359 // 360 // compute the integral and error estimates. 361 // 362 resc12 = cheb12[13]*chebmo[m,13]; 363 ress12 = 0.0; 364 k = 11; 365 for (j=1; j<=6; j++) { // end: l140 366 resc12 = resc12+cheb12[k]*chebmo[m,k]; 367 ress12 = ress12+cheb12[k+1]*chebmo[m,k+1]; 368 k = k-2; 369 l140:;} 370 resc24 = cheb24[25]*chebmo[m,25]; 371 ress24 = 0.0; 372 resabs = fabs(cheb24[25]); 373 k = 23; 374 for (j=1; j<=12; j++) { // end: l150 375 resc24 = resc24+cheb24[k]*chebmo[m,k]; 376 ress24 = ress24+cheb24[k+1]*chebmo[m,k+1]; 377 resabs = fabs(cheb24[k])+fabs(cheb24[k+1]); 378 k = k-2; 379 l150:;} 380 estc = fabs(resc24-resc12); 381 ests = fabs(ress24-ress12); 382 resabs = resabs*fabs(hlgth); 383 if(integr == 2) goto l160; 384 result = conc*resc24-cons*ress24; 385 abserr = fabs(conc*estc)+fabs(cons*ests); 386 goto l170; 387 l160: result = conc*ress24+cons*resc24; 388 abserr = fabs(conc*ests)+fabs(cons*estc); 389 l170: return; 390 } 391 392 393 unittest 394 { 395 alias qc25f!(float, float delegate(float)) fqc25f; 396 alias qc25f!(double, double delegate(double)) dqc25f; 397 alias qc25f!(double, double function(double)) dfqc25f; 398 alias qc25f!(real, real delegate(real)) rqc25f; 399 } 400 401 402 unittest 403 { 404 double f(double x) { return x*x; } 405 406 double a = 0.0, b = 2.0; 407 int integr = 1, nrmom = 0, maxp1 = 21, ksave = 0, momcom = 0; 408 double result, abserr, resabs, resasc; 409 int neval; 410 double[] chebmo = new double[maxp1*25]; 411 412 // Compute using 15-point Gauss-Kronrod formula. 413 double omega = 1.0; 414 double ans = 2*(2*cos(2.0) + sin(2.0)); 415 qc25f(&f, a, b, omega, integr, nrmom, maxp1, ksave, result, abserr, 416 neval, resabs, resasc, momcom, chebmo.ptr); 417 assert (isAccurate(result, abserr, ans, 1e-8)); 418 419 420 // Compute using Clenshaw-Curtis method. 421 omega = 3.0; 422 ans = 2*(6*cos(6.0) + 17*sin(6.0))/27; 423 qc25f(&f, a, b, omega, integr, nrmom, maxp1, ksave, result, abserr, 424 neval, resabs, resasc, momcom, chebmo.ptr); 425 assert (isAccurate(result, 1e-10, ans, 1e-8)); 426 427 // Compute again, using the stored moments. 428 ksave = 1; 429 qc25f(&f, a, b, omega, integr, nrmom, maxp1, ksave, result, abserr, 430 neval, resabs, resasc, momcom, chebmo.ptr); 431 assert (isAccurate(result, 1e-10, ans, 1e-8)); 432 } 433 434 435 unittest 436 { 437 double alpha = 3.0; 438 double f(double x) { return x <= 0.0 ? 0.0 : exp(-(2.0^^(-alpha))*x)/sqrt(x); } 439 440 double a = 1.0, omega = 1.0; 441 int integr = 1, nrmom = 0, maxp1 = 21, ksave = 0, momcom = 0; 442 double result, abserr, resabs, resasc; 443 int neval; 444 double[] chebmo = new double[maxp1*25]; 445 446 // Compute using 15-point Gauss-Kronrod formula. 447 double b = 2.0; 448 double ans = 0.0729155596656107239492740504503422453; 449 qc25f(&f, a, b, omega, integr, nrmom, maxp1, ksave, result, abserr, 450 neval, resabs, resasc, momcom, chebmo.ptr); 451 assert (approxEqual(result, ans, 1e-8)); 452 453 454 // Compute using Clenshaw-Curtis method. 455 b = 6.0; 456 ans = -0.50789465569969289260403404655063469; 457 qc25f(&f, a, b, omega, integr, nrmom, maxp1, ksave, result, abserr, 458 neval, resabs, resasc, momcom, chebmo.ptr); 459 assert (approxEqual(result, ans, 1e-8)); 460 461 // Compute again, using the stored moments. 462 ksave = 1; 463 qc25f(&f, a, b, omega, integr, nrmom, maxp1, ksave, result, abserr, 464 neval, resabs, resasc, momcom, chebmo.ptr); 465 assert (approxEqual(result, ans, 1e-8)); 466 } 467