A highly parallel algorithm has been developed and exploited to compute the planetary normal modes of the elastic-gravitational system, which is approximated via the mixed finite element method on unstructured tetrahedral meshes. The eigenmodes of the relevant generalized eigenvalue problem were extracted by a Lanczos approach combined with polynomial filtering. In contrast with the standard shift-and-invert and the full-mode coupling algorithms, the polynomial filtering technique is ideally suited for solving large-scale 3-D interior eigenvalue problems since it significantly enhances the memory and computational efficiency without loss of accuracy. The parallel efficiency and scalability of this approach are demonstrated on Stampede2 at the Texas Advanced Computing Center. To our knowledge, this is the first time that the direct calculation of the normal modes of 3-D strongly heterogeneous planets, in particular, Earth and Mars, is made feasible via a combination of multiple matrix-free methods and a separation of the essential spectra.