Experimental data suggest that the viscoelastic behavior of rocks is more easily and accurately described in the frequency domain than in the time domain, supporting the idea of simulating seismic wave propagation in the frequency domain. We evaluated weighted-averaged 27-point finite-difference operators for 3D viscoelastic wave modeling in the frequency domain. Within the proposed framework, we developed general equations for normalized phase velocities that can be used with arbitrary finite-difference operators. Three sets of weighting coefficients for second-order central finite-difference operators that minimize the numerical dispersion for up to five grid points per wavelength were found using a damped least-squares (LS) criterion as well as a global optimization scheme based on - and -norm criteria. The three sets produced very similar dispersion curves, and improvement provided by global optimization appeared marginal in this respect. We also evaluated a discrete form for the heterogeneous formulation of the 3D viscoelastic equations with a perfectly match layer (PML). Heuristic performance assessment of frequency-dependent PML absorption coefficients provided a simple rule giving good results for eight PMLs at all frequencies. The proposed formalism was implemented with a massively parallel direct solver. Modeling results were compared with an analytic solution and a time-domain finite-difference code, and they gave good agreement when using LS and -norm optimal coefficients. On the other hand, -norm coefficients produced noisy results, indicating that minimizing the difference between analytic and numerical phase velocities, although necessary, is not a sufficient condition to guarantee low-numerical noise. Finally, analysis of the computational resources required to factorize the impedance matrix revealed that the memory complexity of the factorization is for an grid, compared to for the viscoacoustic case.