  We calculate the integrated galactic initial stellar mass function (IGIMF) in
the presence of IMF variations in clusters. IMF Variations for a population of
clusters are taken into account in the form of Gaussian distributions of the
IMF parameters. For the tapered power law function used here, these are the
slopes at the high and low mass ends, $\Gamma$ and $\gamma$, and the
characteristic mass $M_{ch}$. Variations are modeled by varying the width of
the Gaussian distributions. The reference values are the standard deviations of
the parameters observed for young clusters in the present-day Milky Way
$\sigma_{\Gamma}=0.6$, $\sigma_{\gamma}=0.25$, and $\sigma_{M_{ch}}=0.27$
M$_{\odot}$. Increasing the dispersions of $\gamma$ and $\Gamma$ moderately
flattens the IGIMF at the low and high mass ends. Increasing $\sigma_{M_{ch}}$
shifts the peak of the IGIMF to lower masses, rendering the IGIMF more bottom
heavy. This can explain the bottom heavy stellar mass function of Early-type
galaxies as they are the result of the merger of disk galaxies where the
physical conditions of the star forming gas vary significantly both in time and
space. The effect of IMF variations is compared to that due to other effects
such as variations in the shape of the initial cluster mass function,
metallicity, and galactic SFR. We find that the effect of IMF variations is a
dominant factor that always affects the characteristic mass of the IGIMF. We
compare our results to a sample of ultra-faint dwarf satellite galaxies (UFDs).
Their present-day stellar mass function is an analog to their IGIMF at the time
their stellar populations have formed. We show that the slope of the IGIMF of
the UFDs can only be reproduced when IMF variations of the same order as those
measured in the present-day Milky Way are included. (Abridged)
