  State-of-the-art computations of the gravitational self-force (GSF) on
massive particles in black hole spacetimes involve numerical evolution of the
metric perturbation equations in the time-domain, which is computationally very
costly. We present here a new strategy, based on a frequency-domain treatment
of the perturbation equations, which offers considerable computational saving.
The essential ingredients of our method are (i) a Fourier-harmonic
decomposition of the Lorenz-gauge metric perturbation equations and a numerical
solution of the resulting coupled set of ordinary equations with suitable
boundary conditions; (ii) a generalized version of the method of extended
homogeneous solutions [Phys. Rev. D {\bf 78}, 084021 (2008)] used to circumvent
the Gibbs phenomenon that would otherwise hamper the convergence of the Fourier
mode-sum at the particle's location; and (iii) standard mode-sum
regularization, which finally yields the physical GSF as a sum over regularized
modal contributions. We present a working code that implements this strategy to
calculate the Lorenz-gauge GSF along eccentric geodesic orbits around a
Schwarzschild black hole. The code is far more efficient than existing
time-domain methods; the gain in computation speed (at a given precision) is
about an order of magnitude at an eccentricity of 0.2, and up to three orders
of magnitude for circular or nearly circular orbits. This increased efficiency
was crucial in enabling the recently reported calculation of the long-term
orbital evolution of an extreme mass ratio inspiral [Phys. Rev. D {\bf 85},
061501(R) (2012)]. Here we provide full technical details of our method to
complement the above report.
