A simple but robust method for joint (multiple) earthquake location in arbitrary 3D velocity media is presented in this article. It exploits the advantages of the irregular shortest-path method (irregular spm) for the raytracing, and uses a damped minimum norm, constrained least-squares solution for the source parameter (hypocenter and origin time) update. Special features of the scheme, made possible by the irregular spm, are that the raytracing needs to be carried out only once and the elements of the Jacobian matrix are calculated directly. Six supplementary sources are added around a trial (to be updated) source and trilinear interpolation used to form the travel times and the derivatives. Comparison tests against the regular spm and the finite-difference method show that the newly developed procedure is able to achieve comparable accuracy in both the travel-time calculation and in the hypocentral location solution, but at considerable computational economy. Further numerical tests for a complex 3D model indicate that the method is accurate, with a location error in depth comparable to that in the horizontal plane, at least for the station distribution used. In a practical sense, there is no requirement to have a close initial guess of the source coordinates and origin time before the location process is undertaken if the sources are not too far from the seismic network. The results show that this approach is fairly insensitive to modest levels of random noise, but like all hypocenter determination algorithms, it strongly depends on having the correct velocity model.