Many studies for nonequilibrium systems, e.g., the pre-equilibration puzzle in heavy-ion collisions, require solving the relativistic Boltzmann equation (BE) with the full collisional kernel to high precision. It is challenging to solve relativistic BE due to its high dimensional phase-space integrals and limited computing resources. We have developed a numerical framework for a full solution of the relativistic Boltzmann equations for the quark-gluon matter using the multiple graphics processing units (GPUs) on distributed clusters. Including all the 2 → 2 scattering processes of 3-flavor quarks and gluons, we compute the time evolution of distribution functions in both coordinate and momentum spaces for the cases of pure gluons, quarks, and the mixture of quarks and gluons. By introducing a symmetrical sampling method on GPUs which ensures the particle number conservation, our framework is able to perform the space-time evolution of quark-gluon system toward thermal equilibrium with high performance. We also observe that the gluons naturally accumulate in the soft region at the early time, which may indicate the gluon condensation.