  Using a super-operator formulation of linearized time-dependent
density-functional theory, the dynamical polarizability of a system of
interacting electrons is given a matrix continued-fraction representation whose
coefficients can be obtained from the non-symmetric block-Lanczos method. The
resulting algorithm allows for the calculation of the {\em full spectrum} of a
system with a computational workload which is only a few times larger than that
needed for {\em static} polarizabilities within time-independent
density-functional perturbation theory. The method is demonstrated with the
calculation of the spectrum of benzene, and prospects for its application to
the large-scale calculation of optical spectra are discussed.
