The self-consistent evaluation of Hubbard parameters using linear-response theory is crucial for quantitatively predictive calculations based on Hubbard-corrected density-functional theory. Here, we extend a recently introduced approach based on density-functional perturbation theory (DFPT) for the calculation of the onsite Hubbard U to also compute the intersite Hubbard V. DFPT allows us to reduce significantly computational costs, improve numerical accuracy, and fully automate the calculation of the Hubbard parameters by recasting the linear response of a localized perturbation into an array of monochromatic perturbations that can be calculated in the primitive cell. In addition, here we generalize the entire formalism from norm-conserving to ultrasoft and projector-augmented wave formulations, and to metallic ground states. After benchmarking DFPT against the conventional real-space Hubbard linear response in a supercell, we demonstrate the effectiveness of the present extended Hubbard formulation in determining the equilibrium crystal structure of LixMnPO4 (x = 0, 1) and the subtle energetics of Li intercalation.