We considered the discretization and approximate solutions of equations describing time-harmonic qP-polarized waves in 3D inhomogeneous anisotropic media. The anisotropy comprises general (tilted) transversely isotropic symmetries. We are concerned with solving these equations for a large number of different sources. We considered higher-order partial differential equations and variable-order finite-difference schemes to accommodate anisotropy on the one hand and allow higher-order accuracy — to control sampling rates for relatively high frequencies — on the other hand. We made use of a nested dissection based domain decomposition in a massively parallel multifrontal solver combined with hierarchically semiseparable matrix compression techniques. The higher-order partial differential operators and the variable-order finite-difference schemes require the introduction of separators with variable thickness in the nested dissection; the development of these and their integration with the multifrontal solver is the main focus of our study. The algorithm that we developed is a powerful tool for anisotropic full-waveform inversion.