We complete the study of two-loop infrared singularities of scattering amplitudes with an arbitrary number of massive and massless partons in non-abelian gauge theories. To this end, we calculate the universal functions F_1 and f_2, which completely specify the structure of three-parton correlations in the soft anomalous-dimension matrix, at two-loop order in closed analytic form. Both functions are found to be suppressed like O(m^4/s^2) in the limit of small parton masses, in accordance with mass factorization theorems proposed in the literature. On the other hand, they are unsuppressed and diverge logarithmically near the threshold for pair production of two heavy particles. As an application, we calculate the two-loop anomalous-dimension matrix for q q_bar --> t t_bar near threshold and show that it is not diagonal in the s-channel singlet-octet basis.