  We introduce a new implementation of time-dependent density-functional theory
which allows the \emph{entire} spectrum of a molecule or extended system to be
computed with a numerical effort comparable to that of a \emph{single} standard
ground-state calculation. This method is particularly well suited for large
systems and/or large basis sets, such as plane waves or real-space grids. By
using a super-operator formulation of linearized time-dependent
density-functional theory, we first represent the dynamical polarizability of
an interacting-electron system as an off-diagonal matrix element of the
resolvent of the Liouvillian super-operator. One-electron operators and density
matrices are treated using a representation borrowed from time-independent
density-functional perturbation theory, which permits to avoid the calculation
of unoccupied Kohn-Sham orbitals. The resolvent of the Liouvillian is evaluated
through a newly developed algorithm based on the non-symmetric Lanczos method.
Each step of the Lanczos recursion essentially requires twice as many
operations as a single step of the iterative diagonalization of the unperturbed
Kohn-Sham Hamiltonian. Suitable extrapolation of the Lanczos coefficients
allows for a dramatic reduction of the number of Lanczos steps necessary to
obtain well converged spectra, bringing such number down to hundreds (or a few
thousands, at worst) in typical plane-wave pseudopotential applications. The
resulting numerical workload is only a few times larger than that needed by a
ground-state Kohn-Sham calculation for a same system. Our method is
demonstrated with the calculation of the spectra of benzene, C$_{60}$
fullerene, and of chlorofyll a.
