Iterative image reconstruction algorithms can model complicated imaging physics, compensate for imperfect data acquisition systems, and exploit prior information regarding the object. Hence, they produce higher quality images than do analytical image reconstruction algorithms. However, three-dimensional (3D) iterative image reconstruction is computationally burdensome, which greatly hinders its use with applications requiring a large field-of-view (FOV), such as breast imaging. In this study, an improved GPU-based implementation of a numerical imaging model and its adjoint have been developed for use with general gradient-based iterative image reconstruction algorithms. Both computer simulations and experimental studies are conducted to investigate the efficiency and accuracy of the proposed implementation for optoacoustic tomography (OAT). The results suggest that the proposed implementation is more than five times faster than the previous implementation.