ruleadapt_integrate()

static int ruleadapt_integrate ( rule *  r, unsigned  fdim, integrand_v  f, void *  fdata, const hypercube *  h, unsigned  maxEval, double  reqAbsError, double  reqRelError, double *  val, double *  err, int  parallel  )

static

{
     unsigned numEval = 0;
     heap regions;
     unsigned i, j;
     region *R = NULL; /* array of regions to evaluate */
     unsigned nR_alloc = 0;
     esterr *ee = NULL;

     regions = heap_alloc(1, fdim);
     if (!regions.ee || !regions.items) goto bad;

     ee = (esterr *) malloc(sizeof(esterr) * fdim);
     if (!ee) goto bad;
     
     nR_alloc = 2;
     R = (region *) malloc(sizeof(region) * nR_alloc);
     if (!R) goto bad;
     R[0] = make_region(h, fdim);
     if (!R[0].ee
         || eval_regions(1, R, f, fdata, r)
         || heap_push(&regions, R[0]))
               goto bad;
     numEval += r->num_points;
     
     while (numEval < maxEval || !maxEval) {
          for (j = 0; j < fdim && (regions.ee[j].err <= reqAbsError
                                   || relError(regions.ee[j]) <= reqRelError);
               ++j) ;
          if (j == fdim)
               break; /* convergence */

          if (parallel) { /* maximize potential parallelism */
               /* adapted from I. Gladwell, "Vectorization of one
                  dimensional quadrature codes," pp. 230--238 in
                  _Numerical Integration. Recent Developments,
                  Software and Applications_, G. Fairweather and
                  P. M. Keast, eds., NATO ASI Series C203, Dordrecht
                  (1987), as described in J. M. Bull and
                  T. L. Freeman, "Parallel Globally Adaptive
                  Algorithms for Multi-dimensional Integration,"
                  http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.6638
                  (1994). 

                  Basically, this evaluates in one shot all regions
                  that *must* be evaluated in order to reduce the
                  error to the requested bound: the minimum set of
                  largest-error regions whose errors push the total
                  error over the bound.

                  [Note: Bull and Freeman claim that the Gladwell
                  approach is intrinsically inefficent because it
                  "requires sorting", and propose an alternative
                  algorithm that "only" requires three passes over the
                  entire set of regions.  Apparently, they didn't
                  realize that one could use a heap data structure, in
                  which case the time to pop K biggest-error regions
                  out of N is only O(K log N), much better than the
                  O(N) cost of the Bull and Freeman algorithm if K <<
                  N, and it is also much simpler.] */
               unsigned nR = 0;
               for (j = 0; j < fdim; ++j) ee[j] = regions.ee[j];
               do {
                    if (nR + 2 > nR_alloc) {
                         nR_alloc = (nR + 2) * 2;
                         R = (region *) realloc(R, nR_alloc * sizeof(region));
                         if (!R) goto bad;
                    }
                    R[nR] = heap_pop(&regions);
                    for (j = 0; j < fdim; ++j) ee[j].err -= R[nR].ee[j].err;
                    if (cut_region(R+nR, R+nR+1)) goto bad;
                    numEval += r->num_points * 2;
                    nR += 2;
                    for (j = 0; j < fdim 
                              && (ee[j].err <= reqAbsError
                                  || relError(ee[j]) <= reqRelError); ++j) ;
                    if (j == fdim) break; /* other regions have small errs */
               } while (regions.n > 0 && (numEval < maxEval || !maxEval));
               if (eval_regions(nR, R, f, fdata, r)
                   || heap_push_many(&regions, nR, R))
                    goto bad;
          }
          else { /* minimize number of function evaluations */
               R[0] = heap_pop(&regions); /* get worst region */
               if (cut_region(R, R+1)
                   || eval_regions(2, R, f, fdata, r)
                   || heap_push_many(&regions, 2, R))
                    goto bad;
               numEval += r->num_points * 2;
          }
     }

     /* re-sum integral and errors */
     for (j = 0; j < fdim; ++j) val[j] = err[j] = 0;  
     for (i = 0; i < regions.n; ++i) {
          for (j = 0; j < fdim; ++j) { 
               val[j] += regions.items[i].ee[j].val;
               err[j] += regions.items[i].ee[j].err;
          }
          destroy_region(&regions.items[i]);
     }

     /* printf("regions.nalloc = %d\n", regions.nalloc); */
     free(ee);
     heap_free(&regions);
     free(R);
     return SUCCESS;

bad:
     free(ee);
     heap_free(&regions);
     free(R);
     return FAILURE;
}