The dynamic response of a viscoelastic suspension of spheres under small and large amplitude oscillatory shear is investigated by three-dimensional direct numerical simulations. A sliding triperiodic domain is implemented whereby the computational domain is regarded as the bulk of an infinite suspension. A fictitious domain method is used to manage the particle motion. After the stress field is computed, the bulk properties are recovered by an averaging procedure. The numerical method is validated by comparing the computed linear viscoelastic response of Newtonian and non-Newtonian suspensions with previous theories and simulations. The numerical predictions are in very good quantitative agreement with experimental data for the Newtonian case, whereas deviations are found with respect to some sets of experiments for semidilute and concentrated viscoelastic suspensions. To investigate on such discrepancies, the effect of aggregates in the bulk of the suspension is examined. The simulations show that the presence of structures significantly alters the loss modulus. Such an effect is more pronounced as the volume fraction increases. In this light, the above mentioned disagreement between simulations and data (and among experimental data themselves) can be rationalized, as its origin can be attributed to inhomogeneous particle configurations. For increasing strain amplitudes, both loss and storage moduli depart from the linear viscoelastic values. Although the deviations are qualitatively similar to the large amplitude response of the unfilled suspending matrix, our results for dilute and semidilute suspensions show that the decrease of the moduli is more and more pronounced as the volume fraction is higher. Furthermore, a higher concentration of solid particles reduces the value of strain amplitude such that the nonlinear behavior is observed. Simulations at higher frequencies also correctly capture the overshoot in the loss modulus for intermediate strain amplitudes. Finally, the effect of fluid elasticity on the particle motion is analyzed. The particles are found to move away from their starting positions and the average distance, computed at the beginning of each cycle with respect to the initial configuration, linearly increases with the number of cycles. The change in the microstructure is attributed to the long-range hydrodynamic interactions mediated by fluid viscoelasticity. (C) 2013 The Society of Rheology. [http://dx.doi.org/10.1122/1.4798626]