Computational acoustics offers a powerful tool for investigating the interaction between ultrasound waves and the human body. However, in many common ultrasound settings, performing realistic simulations is computationally difficult due to the large size of the tissue volume compared to the size of the acoustic wavelength. This is particularly true in the case of high-intensity focused ultrasound where large diameter transducers are used to tightly focus ultrasound waves, often deep within the tissue. Here, an efficient model for large-scale ultrasound simulation is presented. The model is based on coupled first-order acoustic equations valid for nonlinear wave propagation in heterogeneous media with power law absorption. The equations are discretized using the k-space pseudospectral method and encoded using advanced programming techniques for parallel computer architectures. This allows the efficient simulation of nonlinear ultrasound propagation in three-dimensions over hundreds of wavelengths. Applications to both diagnostic and therapeutic ultrasound are discussed, and the results from several simulation examples are presented.