  (Abridged) We selected high-z massive galaxies at 5.8 microns, in the SWIRE
ELAIS-S1 field (1 sq. deg.). Galaxies with the 1.6 microns stellar peak
redshifted into the IRAC bands (z~1-3, called ``IR-peakers'') were identified.
Stellar masses were derived by means of spectro-photometric fitting and used to
compute the stellar mass function (MF) at z=1-2 and 2-3. A parametric fit to
the MF was performed, based on a Bayesian formalism, and the stellar mass
density of massive galaxies above z=2 determined. We present the first
systematic study of the very-massive tail of the galaxy stellar mass function
at high redshift. A total of 326 sources were selected. The majority of these
galaxies have stellar masses in excess of 1e11 Msun and lie at z>1.5. The
availability of mid-IR data turned out to be a valuable tool to constrain the
contribution of young stars to galaxy SEDs, and thus their M(stars)/L ratio.
The influence of near-IR data and of the chosen stellar library on the SED
fitting are also discussed. A significant evolution is found not only for
galaxies with M~1e11 Msun, but also in the highest mass bins considered. The
comoving number density of these galaxies was lower by more than a factor of 10
at z=2-3, with respect to the local estimate. SWIRE 5.8 micron peakers more
massive than 1.6x1e11 Msun provide 30-50% of the total stellar mass density in
galaxies at z=2-3.
