  Anomalous diffusion is often modelled in terms of the subdiffusion equation,
which can involve a weakly singular source term. For this case, many
predominant time stepping methods, including the correction of high-order BDF
schemes [{\sc Jin, Li, and Zhou}, SIAM J. Sci. Comput., 39 (2017),
A3129--A3152], may suffer from a severe order reduction.
  To fill in this gap, we propose a smoothing method for time stepping schemes,
where the singular term is regularized by using a $m$-fold
integral-differential calculus and the equation is discretized by the $k$-step
BDF convolution quadrature, called ID$m$-BDF$k$ method. We prove that the
desired $k$th-order convergence can be recovered even if the source term is a
weakly singular and the initial data is not compatible. Numerical experiments
illustrate the theoretical results.
