Over the last few decades, developing efficient iterative methods for solving discretized partial differential equations (PDEs) has been a topic of intensive research. Though these efforts have yielded many mathematically optimal solvers, such as the multigrid method, the unfortunate reality is that multigrid methods have not been used much in practical applications. This marked gap between theory and practice is mainly due to the fragility of traditional multigrid methodology and the complexity of its implementation. This paper aims to develop theories and techniques that will narrow this gap. Specifically, its aim is to develop mathematically optimal solvers that are robust and easy to use for a variety of problems in practice. One central mathematical technique for reaching this goal is a general framework called the Fast Auxiliary Space Preconditioning (FASP) method. FASP methodology represents a class of methods that (1) transform a complicated system into a sequence of simpler systems by using auxiliary spaces and (2) produces an efficient and robust preconditioner (to be used with Krylov space methods such as CG and GMRes) in terms of efficient solvers for these simpler systems. By carefully making use of the special features of each problem, the FASP method can be efficiently applied to a large class of commonly used partial differential equations including equations of Poisson, diffusion-convection-reaction, linear elasticity, Stokes, Brinkman, Navier-Stokes, complex fluids models, and magnetohydrodynamics. This paper will give a summary of results that have been obtained mostly by the author and his collaborators on this topic in recent years.