In this article, we concentrate on the fast numerical computation of the radiation energy densities together with electron and ion temperatures of three-dimensional multi-group radiation diffusion equations, which is temporally discretized with the adaptive backward Eulerian scheme, linearized iteratively via the method of frozen coefficients and spatially approximated through a cell-centered finite volume discretization on the adaptive unstructured computational meshes. We present, analyze and implement an alternating positive semidefinite splitting preconditioning technique with two selective relaxations and algebraic multigrid subsolves, and provide an algebraic quasi-optimal selection approach to determine the involved parameters. Our parallel implementation is based on the software package jxpamg and the preconditioned flexible restarted generalized minimal residual solver has been examined by running realistic simulations of hydrodynamic instability on the Tianhe-2A supercomputer to demonstrate its numerical robustness, computational efficiency, parallel strong and weak scalabilities, and the competitiveness with some existing popular monolithic and block preconditioning strategies.