Robust absorbing boundary conditions are central to the utility and advancement of 3-D numerical wave propagation methods. It is in general preferred that an absorbing boundary method be capable of broadband absorption, be efficient in terms of memory and computation time, and be widely stable in connection with sophisticated numerical schemes. Here we discuss these issues for a promising absorbing boundary method, perfectly matched layers (PML), as implemented in the widely used fourth-order accurate three-dimensional (3-D) staggered-grid velocity-stress finite difference (FD) scheme. Numerical results for point (explosive and double couple) and extended sources, velocity structures (homogeneous, 1-D and 3-D), and different thickness PML zones are excellent, in general, leaving no observable reflections in PML seismograms compared to the amplitudes of the primary phases. For both homogeneous half-space and 1-D models, typical amplitude reduction factors (with respect to the maximum trace amplitude) range between 1/100 and 1/625 for PML thicknesses of 5--20 nodes. A PML region of thickness 5 outperforms a simple exponential damping region of thickness 20 in a homogeneous half-space model by a factor of 3. We find that PML is effective across the simulation bandwidth. For example, permanent offset artifacts due to particularly poor absorption of long-period energy by the simple exponential damping are effectively absent when PML is used. The computational efficiency and storage requirements of PML, compared to the simple exponential damping, are reduced due to the need for only narrow absorbing regions. We also discuss stability and present the complete PML model for the 3-D velocity-stress system.