  This paper is concerned with developing an efficient numerical algorithm for
fast implementation of the sparse grid method for computing the $d$-dimensional
integral of a given function. The new algorithm, called the MDI-SG ({\em
multilevel dimension iteration sparse grid}) method, implements the sparse grid
method based on a dimension iteration/reduction procedure, it does not need to
store the integration points, neither does it compute the function values
independently at each integration point, instead, it re-uses the computation
for function evaluations as much as possible by performing the function
evaluations at all integration points in a cluster and iteratively along
coordinate directions. It is shown numerically that the computational
complexity (in terms of CPU time) of the proposed MDI-SG method is of
polynomial order $O(Nd^3 )$ or better, compared to the exponential order
$O(N(\log N)^{d-1})$ for the standard sparse grid method, where $N$ denotes the
maximum number of integration points in each coordinate direction. As a result,
the proposed MDI-SG method effectively circumvents the curse of dimensionality
suffered by the standard sparse grid method for high-dimensional numerical
integration.
