  The Bethe-Salpeter formalism represents the most accurate method available
nowadays for computing neutral excitation energies and optical spectra of
crystalline systems from first principles. Bethe-Salpeter calculations yield
very good agreement with experiment but are notoriously difficult to converge
with respect to the sampling of the electronic wavevectors. Well-converged
spectra therefore require significant computational and memory resources, even
by today's standards. These bottlenecks hinder the investigation of systems of
great technological interest. They are also barriers to the study of derived
quantities like piezoreflectance, thermoreflectance or resonant Raman
intensities.
  We present a new methodology that decreases the workload needed to reach a
given accuracy. It is based on a double-grid on-the-fly interpolation within
the Brillouin zone, combined with the Lanczos algorithm. It achieves
significant speed-up and reduction of memory requirements. The technique is
benchmarked in terms of accuracy on silicon, gallium arsenide and lithium
fluoride. The scaling of the performance of the method as a function of the
Brillouin Zone point density is much better than a conventional implementation.
We also compare our method with other similar techniques proposed in the
literature.
