This article is concerned with a matrix splitting preconditioning technique with two selective relaxations and algebraic multigrid subsolves for $(G+2)×(G+2)$ block-structured sparse linear systems derived from the three-dimensional flux-limited multi-group radiation diffusion equations, where $G$ is the number of photon energy groups. We introduce an easy-to-implement algebraic selection strategy for the sole contributing parameter, report a spectral analysis and investigate the degree of the minimal polynomial of its left and right preconditioned matrices, and discuss its sequential practical implementation together with the two-level parallelization. Experiments are run with the representative real-world unstructured capsule implosion test cases and it is found that the numerical robustness, computational efficiency and parallel scalability of the proposed preconditioner evaluated on the Tianhe-2A supercomputer with up to 2,816 processor cores are superior to some existing popular monolithic and block preconditioning approaches.