To reduce computational complexity and memory requirement for 3-D elastodynamics using the boundary element method (BEM), a multi-level fast multipole BEM (FM-BEM) is proposed. The diagonal form for the expansion of the elastodynamic fundamental solution is used, with a truncation parameter adjusted to the subdivision level, a feature necessary for achieving optimal computational efficiency. Both the single-level and multi-level forms of the elastodynamic FM-BEM are considered, with emphasis on the latter. Crucial implementation issues, including the truncation of the multipole expansion, the optimal number of levels, the direct and inverse extrapolation steps are examined in detail with the backing of numerical experiments. A complexity analysis for both the single-level and multi-level versions is conducted. The correctness and computational performances of the proposed elastodynamic FMM are demonstrated on numerical examples featuring up to O(106) DOFs run on a single-processor PC.