This paper presents a problem-oriented approach, designed for the numerical simulation of seismic wave propagation in models containing geological formations with complex properties such as anisotropy, attenuation, and small-scale heterogeneities. Each of the named property requires a special treatment that increases the computational complexity of an algorithm in comparison with ideally elastic isotropic media. At the same time, such formations are typically relatively small, filling about 25% of the model, thus the local use of computationally expensive approaches can speed-up the simulation essentially. In this paper we discuss both mathematical and numerical aspects of the hybrid algorithm paying most attention to its parallel implementation. At the same time essential efforts are spent to couple different equations and, hence, different finite-difference stencils to describe properly the different nature of seismic wave propagation in different areas. The main issue in the coupling is to suppress numerical artifacts down to the acceptable level, usually a few tenth of the percent.