The main advantages of the Boundary Element Method (BEM) to model frequency-domain elastic waves in large scale domains with complex topographies are that only the domain boundaries are discretized and no artificial truncation of semi-infinite domains is needed. However, the standard BEM leads to high computational costs. Using a hierarchical representation of the matrix and low-rank approximations, we derive a fast direct solver for elastodynamics. We assess the numerical efficiency and accuracy of the method.