Finite-difference (FD) method for 3D simulation of sonic waves propagation in a borehole and surrounding 3D heterogeneous viscoelastic medium is presented. It is based on explicit second-order staggered grid FD scheme that solves the first-order viscoelastic wave equations (velocity-stress formulation) in cylindrical coordinates. Attenuation is implemented via Generalize Standard Linear Solid model. In order to provide desired Quality Factors τ - method is applied (Blanch et al., 1995). In order to bound area of computations special modification of Perfectly Matched Layer for cylindrical coordinate system is developed and implemented (Kostin et al., 2007). Essential 3D nature of the waves processes for realistic models claims necessity to use parallel computations. Parallelization is performed on the base of domain decomposition approach and implemented under Message Passing Interface. Results of numerical experiments are presented.