This chapter provides an overview of state-of-the-art adaptive finite element methods (AFEMs) for the numerical solution of second-order elliptic partial differential equations (PDEs), where the primary focus is on the optimal interplay of local mesh refinement and iterative solution of the arising discrete systems. Particular emphasis is placed on the thorough description of the essential ingredients necessary to design adaptive algorithms of optimal complexity, i.e., algorithms that mathematically guarantee the optimal rate of convergence with respect to the overall computational cost and, hence, time. Crucially, adaptivity induces reliability of the computed numerical approximations by means of a-posteriori error control. This ensures that the error committed by the numerical scheme is bounded from above by computable quantities. The analysis of the adaptive algorithms is based on the study of appropriate quasi-error quantities that include and balance different components of the overall error. Importantly, the quasi-errors stemming from an adaptive algorithm with contractive iterative solver satisfy a centerpiece concept, namely, full R-linear convergence. This guarantees that the adaptive algorithm is essentially contracting this quasi-error at each step and it turns out to be the cornerstone for the optimal complexity of AFEM. The unified analysis of the adaptive algorithms is presented in the context of symmetric linear PDEs. Extensions to goal-oriented, non-symmetric, as well as non-linear PDEs are presented with suitable nested iterative solvers fitting into the general analytical framework of the linear symmetric case. Numerical experiments highlight the theoretical results and emphasize the practical relevance and gain of adaptivity with iterative solvers for numerical simulations with optimal complexity.