A model of the evolution of an ensemble of magnetic massive stars on the main sequence is developed. We use our own population synthesis code, which allows us to obtain distributions of stars by radii, ages, masses, temperatures, effective magnetic fields, and magnetic fluxes from the pre-main sequence up to the TAMS stages. We assume that magnetic fields in massive stars decrease with time. The rate of magnetic field dissipation may depend on the mass of a star on ZAMS. The distribution of magnetic fluxes of the ZAMS stars is assumed to be log-normal. We show that such kind of distribution may be a result of the dynamo action occurring at the pre-MS evolutionary stage of magnetic stars. Our model also includes capabilities for statistical simulations and parameter estimation necessary for the analysis of real data. Comparison of model magnetic field distributions with those obtained from recent measurements of stellar magnetic fields allows us to conclude that the evolution of magnetic fields of massive stars is very slow if not absent. The shape of the real magnetic field distribution has no indications of the "magnetic desert," previously suggested by Lingieres et al. (2014). Based on those findings we argue that the observed fraction of magnetic stars is determined only by physical conditions at early stages of stellar evolution.