The iterative solvers have a tendency to 'stall', that is to stop removing errors after some number of iterations. The problem is more prominent when the meshes are refined or if the meshes are very very fine. For many iterative methods the number of iteration required to converge the solution are linearly proportional to number of nodes in one direction. This behavior could be rooted out to the fact that during the iteration process, the information travels only one grid size per iteration. And for proper convergence the information has to travel back and forward several times.
Since direct matrix inversion is out of the question for realistic problems and the solvers that rely on conjugate-gradient methods have robustness problems associated with them. The methods of choice are point implicit solvers like Gauss-Seidel and Jacobi. Although the Gauss-Seidel scheme rapidly removes local or high-frequency errors in the solution, global or low-frequency errors are reduced at prohitively slow rates, which is inversely related to the grid size.
Multistage solution schemes reduce low frequency errors nicely and this make them ideal solvers when used with combination of good smoothers as Gauss-Seidel solvers.