In this paper, we propose a three-dimensional (3D) beamforming scheme for the massive multiple-input multiple-output (MIMO) system where the base station (BS) employs a uniform rectangular array (URA). In order to avoid the high computational complexity involving large-dimensional channel matrices, a two-stage beamforming method is applied where the second-stage beamforming is a Kronecker product of azimuth and elevation discrete Fourier transform (DFT) beamforming. These DFT prebeamformers are used for cell splitting and form effective channels with lower dimension for first-stage precoding. We develop a low-complexity user grouping algorithm based on the statistical channel state information at the transmitter (CSIT) to partition users. Each group of users is served by the signal-to-leakage-and-noise ratio (SLNR) precoding aiming at suppressing the intra-group and adjacent-group interferences, which is a good balance between performance and complexity. We derive the approximate signal-to-interference-plus-noise ratio (SINR) of our proposed scheme. Numerical results validate that the SINR approximations are tight and indicate the significance of the proposed 3D beamforming scheme.