We briefly review the basis of the indirect boundary element method (IBEM) and apply it to compute the seismic response of three-dimensional topographic features for incident P and S waves. The method is based on the integral representation for scattered and diffracted waves using single-layer boundary sources. This approach is called indirect because source strengths should be obtained as an intermediate step. Boundary conditions lead to a system of integral equations for surface's sources. The discretization is based on the approximate rectification of involved surfaces using circles. This allows simple numerical and analytical integration of the exact Green's function for an unbounded, homogeneous, isotropic elastic space. In previous work, the accuracy of this approach has been verified. However, a serious limitation comes from the maximum frequency of interest because the size of the coefficient matrix grows as the square of frequency, and it can be easily too large for problems of practical interest. To overcome this difficulty, various threshold criteria were designed taking advantage from the significant spatial decay exhibited by Green's functions. Fast approximate solutions were obtained by using sparse matrix computations. Encouraging results are presented.