The basic theory of surface-related multiple elimination (SRME) can be formulated easily for 3D seismic data. However, because standard 3D seismic acquisition geometries violate the requirements of the method, the practical implementation for 3D seismic data is far from trivial. A major problem is to perform the crossline-summation step of 3D SRME, which becomes aliased because of the large separation between receiver cables and between source lines. A solution to this problem, based on hyperbolic sparse inversion, has been presented previously. This method is an alternative to extensive interpolation and extrapolation of data. The hyperbolic sparse inversion is formulated in the time domain and leads to few, but large, systems of equations. In this paper, we propose an alternative formulation using parabolic sparse inversion based on an efficient weighted minimum-norm solution that can be computed in the angular frequency domain. The main advantage of the new method is numerical efficiency because solving many small systems of equations often is faster than solving a few big ones. The method is demonstrated on 3D synthetic and real data with reflected and diffracted multiples. Numerical results show that the proposed method gives improved results compared to 2D SRME. For typical seismic acquisition geometries, the numerical cost running on 50 processors is 0.05 s per output trace. This makes production-scale processing of 3D seismic data feasible on current Linux clusters.