We have developed a fast direct solver for numerical simulation of acoustic waves in 3D heterogeneous media. The Helmholtz equation is approximated by a 27-point finite-difference stencil of second-order accuracy that is optimized to reduce the numerical dispersion. Due to the optimization, dispersion errors less than 1% are achieved with model discretization of only five points per wavelength. Wave-propagation problem requires solving a large system of linear equations with complex sparse symmetric coefficient matrix with seismic shots representing the right-hand side. The first step is the triangular factorization of the coefficient matrix followed by solving systems of linear equations with triangular coefficient matrices as a second step. Having defined the triangular factors, the second step is very cheap, and its linear scaling with respect to the number of shots is the main advantage of direct methods. To reduce memory consumption and computational time at the factorization step, the lower triangular factors are compressed using a low-rank approximation of their nonzero blocks. The compression enables rapid solving of systems of more than equations corresponding to realistic geophysical models. Accuracy and performance comparison of our solver with a highly optimized time-domain solver proves that these approaches complement each other — depending on the problem size and computing configuration, either solver may be preferable.