Many scientific applications require accurate modeling of seismic wave propagation in complex media. These objectives can include fundamental understanding of seismic wave propagation of the Earth on a global scale, including fluid envelopes, mitigation of seismic risk with better quantitative estimates of seismic hazard, and improved exploitation of the natural resources in the crust of the Earth. Accurate quantification is a continual quest, and benchmark protocols have been designed for model definitions and comparison of solutions. Through this quest, an impressive number of numerical tools have been developed, ranging from efficient finite-difference methods to more sophisticated finite-element methods, and including the so-called pseudospectral methods (see Wu and Maupin for a review). The main motivation behind these permanent developments has been to improve the efficiency and accuracy of forward modeling. To achieve this, one systematic choice for both finite-difference and finite-element methods has involved explicit time-stepping integration to avoid matrix inversion.