One key feature of massive multiple-input multiple-output systems is the large number of antennas and users. As a result, reducing the computational complexity of beamforming design becomes imperative. To this end, the goal of this paper is to achieve a lower complexity order than that of existing beamforming methods, via the parallel accelerated random coordinate descent (ARCD). However, it is known that ARCD is only applicable when the problem is convex, smooth, and separable. In contrast, the beamforming design problem is nonconvex, nonsmooth, and nonseparable. Despite these challenges, this paper shows that it is possible to incorporate ARCD for multicast beamforming by leveraging majorization minimization and strong duality. Numerical results show that the proposed method reduces the execution time by one order of magnitude compared to state-of-the-art methods.