Large-scale topology optimization problems demand the solution of a large number of linear systems arising in the finite element analysis. These systems can be solved efficiently by special iterative solvers. Because the linear systems in the sequence of optimization steps change slowly from one step to the next, we can significantly reduce the number of iterations and the runtime of the linear solver by recycling selected search spaces from previous linear systems, and by using preconditioning and scaling techniques. We also provide a new implementation of the 8-node brick (B8) element for the continuous approximation of material distribution (CAMD) approach to improve designs of functionally graded materials. Specifically, we develop a B8/B8 implementation in which the element shape functions are used for the approximation of both displacements and material density at nodal locations. Finally, we evaluate the effectiveness of several solver and preconditioning strategies, and we investigate large-scale examples, including functionally graded materials, which are solved with a special version of the SIMP (solid isotropic material with penalization) model. The effectiveness of the solver is demonstrated by means of a topology optimization problem in a functionally graded material with 1.6 million unknowns on a fast PC.