We used the Green's function method to describe the multiple scattering of surface waves generated by interactions with the complex 3-D earth structures. The basic process involves an integral equation of convolutional type. An efficient multilevel fast multipole method is utilized to accelerate the matrix-vector products that correspond to the integrals on the horizontal plane. This new algorithm is shown to be quite successful when compared to its more traditional complement (e.g. direct integration method, DIM) with particularly the large number of data points. The fast execution is achieved through well-organized truncated multipole expansions, functions grouping and translation operators. In general, the new algorithm has a fruitful logarithmic time complexity as opposed to an uncomfortable exponential time complexity attainable with the traditional algorithms. This algorithm requires the user provide some important parameters to operate. One of them is the truncation number of the infinite series associated with the multipole expansions for which two linear relationships are derived for an automatic determination. The other central parameter is the clustering number that is used to group data points in the data structure. We showed that the clustering number around 5 is mostly an optimum value providing a minimum run time. However, greater clustering numbers (i.e. ∼10) become necessary for an optimum operation when the number of data points gets real large.: In order to test the convergence of the current algorithm we compared our numerical results with analytical solutions provided for cylindrical obstacles by other researchers. A half-space model in which the total wavefield is represented by a finite number of propagating modes describes the reference structure. Not all waves are trapped in this representation since leaking modes and body waves downward radiating into the half-space can exist. The comparisons with the exact solutions revealed that the quality of vertical component amplitudes is reasonably well while some amplitude discrepancies on the horizontal components particularly near and inside the heterogeneities exist. Extra fast half-space velocities complementing the deficiencies on the modal structure greatly help reduce the discrepancies on the horizontal components. © 2006 The Author Journal compilation © 2006 RAS.