The Globally Adaptive Algorithm

In this section we describe the adaptive integration algorithm in detail. All of the adaptive algorithms for numerical integration have a general structure that consists of the following four main steps:
i) Choose some subregion(s) from a set of subregions.
ii) Subdivide the chosen subregion(s).
iii) Apply an integration rule to the new subregions; update the subre-
gion set.
iv) Update global integral and error estimates; check for convergence.

The corresponding components that distinguish these adaptive algorithms are i) a data structure for the subregion set, ii) a subdivision strategy and iii) an integration rule and an error estimator.

The original van Dooren and de Ridder algorithm was initialized by applying the integration rule to the integrand on the whole integration region. The integration rule provided an estimate for the integral, as well as an error estimate and a coordinate axis for the next subdivision. In the original algorithm, the subregion set was an ordered list with the subregions ordered according to the corresponding error estimates, and the algorithm proceeded from one stage to the next by always choosing the subregion with largest estimated error for subdivision. This subregion was divided in half along one coordinate axis to produce two new subregions, and then the algorithm continued with steps iii) and iv) above. Following the recommendations made by Malcolm and Simpson [21] for one dimensional globally adaptive algorithms, Genz and Malik changed this subregion set data structure to a heap, in order to reduce the overheads associated with maintenance of the set. The algorithm in this paper uses a heap data structure for the subregion set.

The new algorithm uses the same integration rule for all of the integrands; in a particular subregion the use of this rule gives an s-vector of integration rule results for that subregion, an s-vector of error estimates and a coordinate axis k for a subsequent subdivision of that subregion. The original and improved algorithms divided one selected subregion into two pieces along the coordinate axis where the integrand had largest fourth divided difference. The new algorithm uses this same strategy, with two modifications. The first modification is introduced because the new algorithm applies to a vector of integrands. The same subdivision is used for all of the integrands, so the subdivision axis for the next subdivision of a particular subregion should be determined using information from all of the integrands. Let be the center of a chosen subregion with dimensions given by the components of the width vector , and let be the vector obtained from by adding to the ith component of . Given some positive parameters and , define a fourth difference operator along coordinate axis i by

The subdivision axis chosen for that subregion is a coordinate axis k where = . The integrand values used for these difference operators are also used in the integration rule, so the calculation of the differences does not significantly contribute to the computation time. The computation of these differences is done at the same time as the integration rule and error estimate computations are done. In order to ensure that several different axes are chosen for subsequent subdivision when the differences are near zero, any term in the sum for is set to zero if that term is less than four times the machine epsilon, relative to the absolute value of the integrand at the center of the subregion. If there are several axes where the maximum occurs, k is chosen from this group to be the axis where the subregion width () is largest.

A second modification to the previous algorithms is introduced to allow the new algorithm to run more efficiently on a parallel computer with p processors. At each stage the p/2 subregions with the largest errors are divided into two pieces along a (possibly different) selected coordinate axis. This produces p subregions for the parallel computation of the integration rule results.

The error assigned to each subregion for the purpose of determining the position of the subregion in the subregion heap is the maximum of the estimated errors from the integration rule results from that subregion for the different components of the vector integrand. In addition to this maximum error, the information that is maintained for a subregion in the heap is the subregion center (n components), the subregion width vector (n components), the subregion integration rule results (s components), the subregion error estimates (s components) and the coordinate axis k chosen for the next subdivision. A heap of M subregions therefore requires M(2n+2s+2) storage locations. In our implementation of the algorithm there must be a user specified limit on the working storage space; the algorithm described here does not continue if the next step requires more than this amount of heap storage. The global variables that need to be maintained by the algorithm are the number of subregions M, the number of integrand evaluations , the integral estimates and the global error estimates .

We finish this section with a more complete description of the algorithm steps. The integration rule and error estimation are assumed to require L integrand evaluations.

 if t#tex2html_wrap_inline683# then

 		 use the integration rule in the whole integration region,

 		 		 returning ,  and a subdivision axis k

 		initialize the heap, M, ,  and 

 		 set P = 1

 		 do while ( and 
      and )

  		 		 take P subregions from the heap

  		 		 divide the P subregions in half along the chosen axes

  		 		 use the integration rule in 2P new subregions (in parallel),

  		 		 		 returning ,  and k for each new subregion

  		 		 update the heap, M, ,  and 

  		 		 set P = max(1, min(p/2, M, -M))

  		 endo

 endif.

The new algorithm retains the features that made the original and improved algorithms efficient, while allowing for parallelism at both the subregion and integrand levels. The heap maintenance could also be parallelized, but test results from implementations of earlier versions of this algorithm have shown that the time for heap maintenance is usually insignificant compared to the integrand evaluation time, even when the integrand evaluations are done in parallel. The other significant changes to the previous algorithms occur in the integration rule and error estimation components, and these changes are discussed in the next two sections.