tcubature.c - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tcubature.c (22028B)
---
1 /*
2 * Copyright (c) 2005 Steven G. Johnson
3 *
4 * Portions (see comments) based on HIntLib (also distributed under
5 * the GNU GPL), copyright (c) 2002-2005 Rudolf Schuerer.
6 *
7 * Portions (see comments) based on GNU GSL (also distributed under
8 * the GNU GPL), copyright (c) 1996-2000 Brian Gough.
9 *
10 * This program is free software; you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation; either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with this program; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 *
24 */
25
26 /* CHANGELOG:
27
28 15jan09 E. Bueler: initialization values for esterr struct; stops warnings
29 from -Wall
30
31 16Sep10 C. Khroulev: very minor changes to get rid of *very* pedantic warnings
32
33 23Jan15 C. Khroulev: avoid calling realloc(ptr, 0) (see heap_free())
34 */
35
36
37 #include <stdio.h>
38 #include <stdlib.h>
39 #include <math.h>
40 #include <limits.h>
41 #include <float.h>
42
43 #include "cubature.h"
44
45 /* Basic datatypes */
46
47 typedef struct {
48 double val, err;
49 } esterr;
50
51 static double relError(esterr ee)
52 {
53 return (ee.val == 0.0 ? HUGE_VAL : fabs(ee.err / ee.val));
54 }
55
56 typedef struct {
57 unsigned dim;
58 double *data; /* length 2*dim = center followed by half-width */
59 double vol; /* cache volume = product of widths */
60 } hypercube;
61
62 static double compute_vol(const hypercube *h)
63 {
64 unsigned i;
65 double vol = 1;
66 for (i = 0; i < h->dim; ++i)
67 vol *= 2 * h->data[i + h->dim];
68 return vol;
69 }
70
71 static hypercube make_hypercube(unsigned dim, const double *center, const double *width)
72 {
73 unsigned i;
74 hypercube h;
75 h.dim = dim;
76 h.data = (double *) malloc(sizeof(double) * dim * 2);
77 for (i = 0; i < dim; ++i) {
78 h.data[i] = center[i];
79 h.data[i + dim] = width[i];
80 }
81 h.vol = compute_vol(&h);
82 return h;
83 }
84
85 static hypercube make_hypercube_range(unsigned dim, const double *xmin, const double *xmax)
86 {
87 hypercube h = make_hypercube(dim, xmin, xmax);
88 unsigned i;
89 for (i = 0; i < dim; ++i) {
90 h.data[i] = 0.5 * (xmin[i] + xmax[i]);
91 h.data[i + dim] = 0.5 * (xmax[i] - xmin[i]);
92 }
93 h.vol = compute_vol(&h);
94 return h;
95 }
96
97 static void destroy_hypercube(hypercube *h)
98 {
99 free(h->data);
100 h->dim = 0;
101 }
102
103 typedef struct {
104 hypercube h;
105 esterr ee;
106 unsigned splitDim;
107 } region;
108
109 static region make_region(const hypercube *h)
110 {
111 region R;
112 /* initialization values by E. Bueler 1/15/09 */
113 R.ee.err = 1.0e30; R.ee.val = 1.0e30;
114 R.h = make_hypercube(h->dim, h->data, h->data + h->dim);
115 R.splitDim = 0;
116 return R;
117 }
118
119 static void destroy_region(region *R)
120 {
121 destroy_hypercube(&R->h);
122 }
123
124 static void cut_region(region *R, region *R2)
125 {
126 unsigned d = R->splitDim, dim = R->h.dim;
127 *R2 = *R;
128 R->h.data[d + dim] *= 0.5;
129 R->h.vol *= 0.5;
130 R2->h = make_hypercube(dim, R->h.data, R->h.data + dim);
131 R->h.data[d] -= R->h.data[d + dim];
132 R2->h.data[d] += R->h.data[d + dim];
133 }
134
135 typedef struct rule_s {
136 unsigned dim; /* the dimensionality */
137 unsigned num_points; /* number of evaluation points */
138 unsigned (*evalError)(struct rule_s *r, integrand f, void *fdata,
139 const hypercube *h, esterr *ee);
140 void (*destroy)(struct rule_s *r);
141 } rule;
142
143 static void destroy_rule(rule *r)
144 {
145 if (r->destroy) r->destroy(r);
146 free(r);
147 }
148
149 static region eval_region(region R, integrand f, void *fdata, rule *r)
150 {
151 R.splitDim = r->evalError(r, f, fdata, &R.h, &R.ee);
152 return R;
153 }
154
155 /***************************************************************************/
156 /* Functions to loop over points in a hypercube. */
157
158 /* Based on orbitrule.cpp in HIntLib-0.0.10 */
159
160 /* ls0 returns the least-significant 0 bit of n (e.g. it returns
161 0 if the LSB is 0, it returns 1 if the 2 LSBs are 01, etcetera). */
162
163 #if (defined(__GNUC__) || defined(__ICC)) && (defined(__i386__) || defined (__x86_64__))
164 /* use x86 bit-scan instruction, based on count_trailing_zeros()
165 macro in GNU GMP's longlong.h. */
166 static unsigned ls0(unsigned n)
167 {
168 unsigned count;
169 n = ~n;
170 __asm__("bsfl %1,%0": "=r"(count):"rm"(n));
171 return count;
172 }
173 #else
174 static unsigned ls0(unsigned n)
175 {
176 const unsigned bits[] = {
177 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
178 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
179 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
180 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
181 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
182 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
183 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
184 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 7,
185 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
186 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
187 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
188 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
189 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
190 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
191 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
192 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 8,
193 };
194 unsigned bit;
195 bit = 0;
196 while ((n & 0xff) == 0xff) {
197 n >>= 8;
198 bit += 8;
199 }
200 return bit + bits[n & 0xff];
201 }
202 #endif
203
204 /**
205 * Evaluate the integral on all 2^n points (+/-r,...+/-r)
206 *
207 * A Gray-code ordering is used to minimize the number of coordinate updates
208 * in p.
209 */
210 static double evalR_Rfs(integrand f, void *fdata, unsigned dim, double *p, const double *c, const
211 double *r)
212 {
213 double sum = 0;
214 unsigned i;
215 unsigned signs = 0; /* 0/1 bit = +/- for corresponding element of r[] */
216
217 /* We start with the point where r is ADDed in every coordinate (This implies signs=0) */
218 for (i = 0; i < dim; ++i)
219 p[i] = c[i] + r[i];
220
221 /* Loop through the points in gray-code ordering */
222 for (i = 0;; ++i) {
223 unsigned mask, d;
224
225 sum += f(dim, p, fdata);
226
227 d = ls0(i); /* which coordinate to flip */
228 if (d >= dim)
229 break;
230
231 /* flip the d-th bit and add/subtract r[d] */
232 mask = 1U << d;
233 signs ^= mask;
234 p[d] = (signs & mask) ? c[d] - r[d] : c[d] + r[d];
235 }
236 return sum;
237 }
238
239 static double evalRR0_0fs(integrand f, void *fdata, unsigned dim, double *p, const double *c,
240 const double *r)
241 {
242 unsigned i, j;
243 double sum = 0;
244
245 for (i = 0; i < dim - 1; ++i) {
246 p[i] = c[i] - r[i];
247 for (j = i + 1; j < dim; ++j) {
248 p[j] = c[j] - r[j];
249 sum += f(dim, p, fdata);
250 p[i] = c[i] + r[i];
251 sum += f(dim, p, fdata);
252 p[j] = c[j] + r[j];
253 sum += f(dim, p, fdata);
254 p[i] = c[i] - r[i];
255 sum += f(dim, p, fdata);
256
257 p[j] = c[j]; /* Done with j -> Restore p[j] */
258 }
259 p[i] = c[i]; /* Done with i -> Restore p[i] */
260 }
261 return sum;
262 }
263
264 /* static double evalR0_0fs4d(integrand f, void *fdata, unsigned dim, double *p, const double *c,
265 double *sum0_, const double *r1, double *sum1_, const double *r2, double *sum2_) */
266 static unsigned evalR0_0fs4d(integrand f, void *fdata, unsigned dim, double *p, const double *c,
267 double *sum0_, const double *r1, double *sum1_, const double *r2, double *sum2_)
268 {
269 double maxdiff = 0;
270 unsigned i, dimDiffMax = 0;
271 double sum0, sum1 = 0, sum2 = 0; /* copies for aliasing, performance */
272
273 double ratio = r1[0] / r2[0];
274
275 ratio *= ratio;
276 sum0 = f(dim, p, fdata);
277
278 for (i = 0; i < dim; i++) {
279 double f1a, f1b, f2a, f2b, diff;
280
281 p[i] = c[i] - r1[i];
282 sum1 += (f1a = f(dim, p, fdata));
283 p[i] = c[i] + r1[i];
284 sum1 += (f1b = f(dim, p, fdata));
285 p[i] = c[i] - r2[i];
286 sum2 += (f2a = f(dim, p, fdata));
287 p[i] = c[i] + r2[i];
288 sum2 += (f2b = f(dim, p, fdata));
289 p[i] = c[i];
290
291 diff = fabs(f1a + f1b - 2 * sum0 - ratio * (f2a + f2b - 2 * sum0));
292
293 if (diff > maxdiff) {
294 maxdiff = diff;
295 dimDiffMax = i;
296 }
297 }
298
299 *sum0_ += sum0;
300 *sum1_ += sum1;
301 *sum2_ += sum2;
302
303 return dimDiffMax;
304 }
305
306 #define num0_0(dim) (1U)
307 #define numR0_0fs(dim) (2 * (dim))
308 #define numRR0_0fs(dim) (2 * (dim) * (dim-1))
309 #define numR_Rfs(dim) (1U << (dim))
310
311 /***************************************************************************/
312 /* Based on rule75genzmalik.cpp in HIntLib-0.0.10: An embedded
313 cubature rule of degree 7 (embedded rule degree 5) due to A. C. Genz
314 and A. A. Malik. See:
315
316 A. C. Genz and A. A. Malik, "An imbedded [sic] family of fully
317 symmetric numerical integration rules," SIAM
318 J. Numer. Anal. 20 (3), 580-588 (1983).
319 */
320
321 typedef struct {
322 rule parent;
323
324 /* temporary arrays of length dim */
325 double *widthLambda, *widthLambda2, *p;
326
327 /* dimension-dependent constants */
328 double weight1, weight3, weight5;
329 double weightE1, weightE3;
330 } rule75genzmalik;
331
332 #define real(x) ((double)(x))
333 #define to_int(n) ((int)(n))
334
335 static int isqr(int x)
336 {
337 return x * x;
338 }
339
340 static void destroy_rule75genzmalik(rule *r_)
341 {
342 rule75genzmalik *r = (rule75genzmalik *) r_;
343 free(r->p);
344 }
345
346 static unsigned rule75genzmalik_evalError(rule *r_, integrand f, void *fdata, const hypercube *h,
347 esterr *ee)
348 {
349 /* lambda2 = sqrt(9/70), lambda4 = sqrt(9/10), lambda5 = sqrt(9/19) */
350 const double lambda2 = 0.3585685828003180919906451539079374954541;
351 const double lambda4 = 0.9486832980505137995996680633298155601160;
352 const double lambda5 = 0.6882472016116852977216287342936235251269;
353 const double weight2 = 980. / 6561.;
354 const double weight4 = 200. / 19683.;
355 const double weightE2 = 245. / 486.;
356 const double weightE4 = 25. / 729.;
357
358 rule75genzmalik *r = (rule75genzmalik *) r_;
359 unsigned i, dimDiffMax, dim = r_->dim;
360 double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0, sum4, sum5, result, res5th;
361 const double *center = h->data;
362 const double *width = h->data + dim;
363
364 for (i = 0; i < dim; ++i)
365 r->p[i] = center[i];
366
367 for (i = 0; i < dim; ++i)
368 r->widthLambda2[i] = width[i] * lambda2;
369 for (i = 0; i < dim; ++i)
370 r->widthLambda[i] = width[i] * lambda4;
371
372 /* Evaluate function in the center, in f(lambda2,0,...,0) and
373 f(lambda3=lambda4, 0,...,0). Estimate dimension with largest error */
374 dimDiffMax = evalR0_0fs4d(f, fdata, dim, r->p, center, &sum1, r->widthLambda2, &sum2,
375 r->widthLambda, &sum3);
376
377 /* Calculate sum4 for f(lambda4, lambda4, 0, ...,0) */
378 sum4 = evalRR0_0fs(f, fdata, dim, r->p, center, r->widthLambda);
379
380 /* Calculate sum5 for f(lambda5, lambda5, ..., lambda5) */
381 for (i = 0; i < dim; ++i)
382 r->widthLambda[i] = width[i] * lambda5;
383 sum5 = evalR_Rfs(f, fdata, dim, r->p, center, r->widthLambda);
384
385 /* Calculate fifth and seventh order results */
386
387 result = h->vol * (r->weight1 * sum1 + weight2 * sum2 + r->weight3 * sum3 + weight4 * sum4 +
388 r->weight5 * sum5);
389 res5th = h->vol * (r->weightE1 * sum1 + weightE2 * sum2 + r->weightE3 * sum3 + weightE4 *
390 sum4);
391
392 ee->val = result;
393 ee->err = fabs(res5th - result);
394
395 return dimDiffMax;
396 }
397
398 static rule *make_rule75genzmalik(unsigned dim)
399 {
400 rule75genzmalik *r;
401
402 if (dim < 2) return 0; /* this rule does not support 1d integrals */
403
404 r = (rule75genzmalik *) malloc(sizeof(rule75genzmalik));
405 r->parent.dim = dim;
406
407 r->weight1 = (real(12824 - 9120 * to_int(dim) + 400 * isqr(to_int(dim)))
408 / real(19683));
409 r->weight3 = real(1820 - 400 * to_int(dim)) / real(19683);
410 r->weight5 = real(6859) / real(19683) / real(1U << dim);
411 r->weightE1 = (real(729 - 950 * to_int(dim) + 50 * isqr(to_int(dim)))
412 / real(729));
413 r->weightE3 = real(265 - 100 * to_int(dim)) / real(1458);
414
415 r->p = (double *) malloc(sizeof(double) * dim * 3);
416 r->widthLambda = r->p + dim;
417 r->widthLambda2 = r->p + 2 * dim;
418
419 r->parent.num_points = num0_0(dim) + 2 * numR0_0fs(dim)
420 + numRR0_0fs(dim) + numR_Rfs(dim);
421
422 r->parent.evalError = rule75genzmalik_evalError;
423 r->parent.destroy = destroy_rule75genzmalik;
424
425 return (rule *) r;
426 }
427
428 /***************************************************************************/
429 /* 1d 15-point Gaussian quadrature rule, based on qk15.c and qk.c in
430 GNU GSL (which in turn is based on QUADPACK). */
431
432 static unsigned rule15gauss_evalError(rule *r, integrand f, void *fdata,
433 const hypercube *h, esterr *ee)
434 {
435 /* Gauss quadrature weights and kronrod quadrature abscissae and
436 weights as evaluated with 80 decimal digit arithmetic by
437 L. W. Fullerton, Bell Labs, Nov. 1981. */
438 const unsigned n = 8;
439 const double xgk[8] = { /* abscissae of the 15-point kronrod rule */
440 0.991455371120812639206854697526329,
441 0.949107912342758524526189684047851,
442 0.864864423359769072789712788640926,
443 0.741531185599394439863864773280788,
444 0.586087235467691130294144838258730,
445 0.405845151377397166906606412076961,
446 0.207784955007898467600689403773245,
447 0.000000000000000000000000000000000
448 /* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule.
449 xgk[0], xgk[2], ... to optimally extend the 7-point gauss rule */
450 };
451 static const double wg[4] = { /* weights of the 7-point gauss rule */
452 0.129484966168869693270611432679082,
453 0.279705391489276667901467771423780,
454 0.381830050505118944950369775488975,
455 0.417959183673469387755102040816327
456 };
457 static const double wgk[8] = { /* weights of the 15-point kronrod rule */
458 0.022935322010529224963732008058970,
459 0.063092092629978553290700663189204,
460 0.104790010322250183839876322541518,
461 0.140653259715525918745189590510238,
462 0.169004726639267902826583426598550,
463 0.190350578064785409913256402421014,
464 0.204432940075298892414161999234649,
465 0.209482141084727828012999174891714
466 };
467
468 const double center = h->data[0];
469 const double width = h->data[1];
470 double fv1[7], fv2[7]; /* C. Khroulev, 16Sep10: hard-wired 7 = n - 1 */
471 const double f_center = f(1, ¢er, fdata);
472 double result_gauss = f_center * wg[n/2 - 1];
473 double result_kronrod = f_center * wgk[n - 1];
474 double result_abs = fabs(result_kronrod);
475 double result_asc, mean, err;
476 unsigned j;
477
478 if (r == NULL) {} /* should quash unused parameter pedantic msg */
479
480 for (j = 0; j < (n - 1) / 2; ++j) {
481 unsigned int j2 = 2*j + 1;
482 double x, f1, f2, fsum, w = width * xgk[j2];
483 x = center - w; fv1[j2] = f1 = f(1, &x, fdata);
484 x = center + w; fv2[j2] = f2 = f(1, &x, fdata);
485 fsum = f1 + f2;
486 result_gauss += wg[j] * fsum;
487 result_kronrod += wgk[j2] * fsum;
488 result_abs += wgk[j2] * (fabs(f1) + fabs(f2));
489 }
490
491 for (j = 0; j < n/2; ++j) {
492 unsigned int j2 = 2*j;
493 double x, f1, f2, w = width * xgk[j2];
494 x = center - w; fv1[j2] = f1 = f(1, &x, fdata);
495 x = center + w; fv2[j2] = f2 = f(1, &x, fdata);
496 result_kronrod += wgk[j2] * (f1 + f2);
497 result_abs += wgk[j2] * (fabs(f1) + fabs(f2));
498
499 }
500
501 ee->val = result_kronrod * width;
502
503 /* compute error estimate: */
504 mean = result_kronrod * 0.5;
505 result_asc = wgk[n - 1] * fabs(f_center - mean);
506 for (j = 0; j < n - 1; ++j)
507 result_asc += wgk[j] * (fabs(fv1[j]-mean) + fabs(fv2[j]-mean));
508 err = (result_kronrod - result_gauss) * width;
509 result_abs *= width;
510 result_asc *= width;
511 if (result_asc != 0 && err != 0) {
512 double scale = pow((200 * err / result_asc), 1.5);
513 if (scale < 1)
514 err = result_asc * scale;
515 else
516 err = result_asc;
517 }
518 if (result_abs > DBL_MIN / (50 * DBL_EPSILON)) {
519 double min_err = 50 * DBL_EPSILON * result_abs;
520 if (min_err > err)
521 err = min_err;
522 }
523 ee->err = err;
524
525 return 0; /* no choice but to divide 0th dimension */
526 }
527
528 static rule *make_rule15gauss(unsigned dim)
529 {
530 rule *r;
531 if (dim != 1) return 0; /* this rule is only for 1d integrals */
532 r = (rule *) malloc(sizeof(rule));
533 r->dim = dim;
534 r->num_points = 15;
535 r->evalError = rule15gauss_evalError;
536 r->destroy = 0;
537 return r;
538 }
539
540 /***************************************************************************/
541 /* binary heap implementation (ala _Introduction to Algorithms_ by
542 Cormen, Leiserson, and Rivest), for use as a priority queue of
543 regions to integrate. */
544
545 typedef region heap_item;
546 #define KEY(hi) ((hi).ee.err)
547
548 typedef struct {
549 unsigned n, nalloc;
550 heap_item *items;
551 esterr ee;
552 } heap;
553
554 static void heap_resize(heap *h, unsigned nalloc)
555 {
556 if (nalloc != 0) {
557 h->nalloc = nalloc;
558 h->items = (heap_item *) realloc(h->items, sizeof(heap_item) * nalloc);
559 } else {
560 h->nalloc = 0;
561 free(h->items);
562 h->items = NULL;
563 }
564 }
565
566 static heap heap_alloc(unsigned nalloc)
567 {
568 heap h;
569 h.n = 0;
570 h.nalloc = 0;
571 h.items = 0;
572 h.ee.val = h.ee.err = 0;
573 heap_resize(&h, nalloc);
574 return h;
575 }
576
577 /* note that heap_free does not deallocate anything referenced by the items */
578 static void heap_free(heap *h)
579 {
580 h->n = 0;
581 heap_resize(h, 0);
582 }
583
584 static void heap_push(heap *h, heap_item hi)
585 {
586 int insert;
587
588 h->ee.val += hi.ee.val;
589 h->ee.err += hi.ee.err;
590 insert = h->n;
591 if (++(h->n) > h->nalloc)
592 heap_resize(h, h->n * 2);
593
594 while (insert) {
595 int parent = (insert - 1) / 2;
596 if (KEY(hi) <= KEY(h->items[parent]))
597 break;
598 h->items[insert] = h->items[parent];
599 insert = parent;
600 }
601 h->items[insert] = hi;
602 }
603
604 static heap_item heap_pop(heap *h)
605 {
606 heap_item ret;
607 int i, n, child;
608
609 if (!(h->n)) {
610 fprintf(stderr, "attempted to pop an empty heap\n");
611 exit(EXIT_FAILURE);
612 }
613
614 ret = h->items[0];
615 h->items[i = 0] = h->items[n = --(h->n)];
616 while ((child = i * 2 + 1) < n) {
617 int largest;
618 heap_item swap;
619
620 if (KEY(h->items[child]) <= KEY(h->items[i]))
621 largest = i;
622 else
623 largest = child;
624 if (++child < n && KEY(h->items[largest]) < KEY(h->items[child]))
625 largest = child;
626 if (largest == i)
627 break;
628 swap = h->items[i];
629 h->items[i] = h->items[largest];
630 h->items[i = largest] = swap;
631 }
632
633 h->ee.val -= ret.ee.val;
634 h->ee.err -= ret.ee.err;
635 return ret;
636 }
637
638 /***************************************************************************/
639
640 /* adaptive integration, analogous to adaptintegrator.cpp in HIntLib */
641
642 static int ruleadapt_integrate(rule *r, integrand f, void *fdata, const hypercube *h, unsigned
643 maxEval, double reqAbsError, double reqRelError, esterr *ee)
644 {
645 unsigned initialRegions; /* number of initial regions (non-adaptive) */
646 unsigned minIter; /* minimum number of adaptive subdivisions */
647 unsigned maxIter; /* maximum number of adaptive subdivisions */
648 unsigned initialPoints;
649 heap regions;
650 unsigned i;
651 int status = -1; /* = ERROR */
652
653 initialRegions = 1; /* or: use some percentage of maxEval/r->num_points */
654 initialPoints = initialRegions * r->num_points;
655 minIter = 0;
656 if (maxEval) {
657 if (initialPoints > maxEval) {
658 initialRegions = maxEval / r->num_points;
659 initialPoints = initialRegions * r->num_points;
660 }
661 maxEval -= initialPoints;
662 maxIter = maxEval / (2 * r->num_points);
663 } else
664 maxIter = UINT_MAX;
665
666 if (initialRegions == 0)
667 return status; /* ERROR */
668
669 regions = heap_alloc(initialRegions);
670
671 heap_push(®ions, eval_region(make_region(h), f, fdata, r));
672 /* or:
673 if (initialRegions != 1)
674 partition(h, initialRegions, EQUIDISTANT, ®ions, f,fdata, r); */
675
676 for (i = 0; i < maxIter; ++i) {
677 region R, R2;
678 if (i >= minIter && (regions.ee.err <= reqAbsError
679 || relError(regions.ee) <= reqRelError)) {
680 status = 0; /* converged! */
681 break;
682 }
683 R = heap_pop(®ions); /* get worst region */
684 cut_region(&R, &R2);
685 heap_push(®ions, eval_region(R, f, fdata, r));
686 heap_push(®ions, eval_region(R2, f, fdata, r));
687 }
688
689 ee->val = ee->err = 0; /* re-sum integral and errors */
690 for (i = 0; i < regions.n; ++i) {
691 ee->val += regions.items[i].ee.val;
692 ee->err += regions.items[i].ee.err;
693 destroy_region(®ions.items[i]);
694 }
695 /* printf("regions.nalloc = %d\n", regions.nalloc); */
696 heap_free(®ions);
697
698 return status;
699 }
700
701 int adapt_integrate(integrand f, void *fdata,
702 unsigned dim, const double *xmin, const double *xmax,
703 unsigned maxEval, double reqAbsError, double reqRelError,
704 double *val, double *estimated_error)
705 {
706 rule *r;
707 hypercube h;
708 esterr ee;
709 int status;
710 /* initialization values by E. Bueler 1/15/09 */
711 ee.err = 1.0e30; ee.val = 1.0e30;
712
713 if (dim == 0) { /* trivial integration */
714 ee.val = f(0, xmin, fdata);
715 ee.err = 0;
716 return 0;
717 }
718 r = dim == 1 ? make_rule15gauss(dim) : make_rule75genzmalik(dim);
719 h = make_hypercube_range(dim, xmin, xmax);
720 status = ruleadapt_integrate(r, f, fdata, &h,
721 maxEval, reqAbsError, reqRelError,
722 &ee);
723 *val = ee.val;
724 *estimated_error = ee.err;
725 destroy_hypercube(&h);
726 destroy_rule(r);
727 return status;
728 }
729