mitsuba-af602c6fd98a/src/libcore/quad.cpp
/* adaptive integration, analogous to adaptintegrator.cpp in HIntLib */
adaptive integration, analogous to adaptintegrator.cpp in HIntLib
adaptive integration analogous to adaptintegrator . cpp in HIntLib
static NDIntegrator::EResult ruleadapt_integrate(rule *r, unsigned int fdim,
		const VectorizedIntegrand & f, const hypercube *h, size_t maxEval,
		Float reqAbsError, Float reqRelError, Float *val, Float *err, size_t &numEval, int parallel) {
	heap regions;
	unsigned int i, j;
	region *R = NULL; 
	unsigned int 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, 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; 
		if (parallel) {
			unsigned int 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; 
			} while (regions.n > 0 && (numEval < maxEval || !maxEval));
			if (eval_regions(nR, R, f, r) || heap_push_many(&regions, nR, R))
				goto bad;
		} else { 
			R[0] = heap_pop(&regions); 
			if (cut_region(R, R+1) || eval_regions(2, R, f, r)
				|| heap_push_many(&regions, 2, R))
				goto bad;
			numEval += r->num_points * 2;
		}
	}
	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]);
	}
	free(ee);
	heap_free(&regions);
	free(R);
	return NDIntegrator::ESuccess;
bad:
	free(ee);
	heap_free(&regions);
	free(R);
	return NDIntegrator::EFailure;
}
static NDIntegrator :: EResult ruleadapt_integrate ( rule * r , unsigned int fdim , const VectorizedIntegrand & f , const hypercube * h , size_t maxEval , Float reqAbsError , Float reqRelError , Float * val , Float * err , size_t & numEval , int parallel ) { heap regions ; unsigned int i , j ; region * R = NULL ; unsigned int 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 , 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 ; if ( parallel ) { unsigned int 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 ; } while ( regions . n > 0 && ( numEval < maxEval || ! maxEval ) ) ; if ( eval_regions ( nR , R , f , r ) || heap_push_many ( & regions , nR , R ) ) goto bad ; } else { R [ 0 ] = heap_pop ( & regions ) ; if ( cut_region ( R , R + 1 ) || eval_regions ( 2 , R , f , r ) || heap_push_many ( & regions , 2 , R ) ) goto bad ; numEval += r -> num_points * 2 ; } } 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 ] ) ; } free ( ee ) ; heap_free ( & regions ) ; free ( R ) ; return NDIntegrator :: ESuccess ; bad : free ( ee ) ; heap_free ( & regions ) ; free ( R ) ; return NDIntegrator :: EFailure ; }
adaptive integration, analogous to adaptintegrator.cpp in HIntLib
adaptive integration analogous to adaptintegrator . cpp in HIntLib
/* array of regions to evaluate */ /* convergence */ /* 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

			   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.] */ /* other regions have small errs */ /* minimize number of function evaluations */ /* get worst region */ /* re-sum integral and errors */ /* printf("regions.nalloc = %d\n", regions.nalloc); */