We investigate the solution of linear systems of equations that arise when singularly perturbed reaction–diffusion partial differential equations are solved using a standard finite difference method on layer adapted grids. It is known that there are difficulties in solving such systems by direct methods when the perturbation parameter, $$\varepsilon $$ , is small (MacLachlan and Madden in SIAM J Sci Comput 35(5):A2225–A2254, 2013). Therefore, iterative methods are natural choices. However, we show that, in two dimensions, the condition number of the coefficient matrix grows unboundedly when $$\varepsilon $$ tends to zero, and so unpreconditioned iterative schemes, such as the conjugate gradient algorithm, perform poorly with respect to $$\varepsilon $$ . We provide a careful analysis of diagonal and incomplete Cholesky preconditioning methods, and show that the condition number of the preconditioned linear system is independent of the perturbation parameter. We demonstrate numerically the surprising fact that these schemes are more efficient when $$\varepsilon $$ is small, than when $$\varepsilon $$ is $$\mathcal {O}(1)$$ . Furthermore, our analysis shows that when the singularly perturbed problem features no corner layers, an incomplete Cholesky preconditioner performs extremely well when $$\varepsilon \ll 1$$ . We provide numerical evidence that our findings extend to three-dimensional problems.