Driven by the challenging task of pursuing the robust and accurate iterative numerical solution of the three-dimensional flux-limited multi-group radiation diffusion equations in an efficient and scalable manner, we propose and analyze a generalized matrix splitting preconditioning scheme with two selective relaxations and algebraic multigrid subsolves, introduce an algebraic quasi-optimal choice strategy to determine the involved parameters and consider its sequential implementation and two-level parallelization. A great deal of numerical results for typical unstructured twenty-group problems arising from realistic simulations of the hydrodynamic instability are presented and discussed to demonstrate the robustness, efficiency, strong and weak parallel scaling properties with up to 2,816 parallel processor cores together with the competitiveness of the proposed preconditioner when compared with several state-of-the-art monolithic and block preconditioning approaches.