  We investigate the matter power spectra in the power law and exponential
types of viable $f(R)$ theories along with massive neutrinos. The enhancement
of the matter power spectrum is found to be a generic feature in these models.
In particular, we show that in the former type, such as the Starobinsky model,
the spectrum is magnified much larger than the latter one, such as the
exponential model. A greater scale of the total neutrino mass, $\Sigma
m_{\nu}$, is allowed in the viable $f(R)$ models than that in the $\Lambda$CDM
one. We obtain the constraints on the neutrino masses by using the CosmoMC
package with the modified MGCAMB. Explicitly, we get $\Sigma m_{\nu} < 0.451 \
(0.214)\ \mathrm{eV}$ at 95% C.L. in the Starobinsky (exponential) model, while
the corresponding one for the $\Lambda$CDM model is $\Sigma m_{\nu} < 0.200\
\mathrm{eV}$. Furthermore, by treating the effective number of neutrino species
$N_{\mathrm{eff}}$ as a free parameter along with $\Sigma m_{\nu}$, we find
that $N_{\mathrm{eff}} = 3.78^{+0.64}_{-0.84} (3.47^{+0.74}_{-0.60})$ and
$\Sigma m_{\nu} = 0.533^{+0.254}_{-0.411}$ ($< 0.386) \ \mathrm{eV}$ at 95%
C.L. in the Starobinsky (exponential) model.
