Double parton distribution functions (DPDFs) are used in the QCD description of double parton scattering. The DPDFs evolve with hard scales through relatively new QCD evolution equations which obey nontrivial momentum and valence quark number sum rules. Based on the constructed numerical program, we present results on the QCD evolution of the DPDFs. In particular, we discuss the problem how to specify initial conditions for the evolution equations which exactly fulfill the sum rules.