We investigate the solution of linear systems of equations that arise when a two-dimensional singularly perturbed reaction-diffusion problem is solved using a standard finite difference method on a layer adapted, piecewise uniform (``Shishkin'') mesh. It is known that %, due to the discretization on layer adapted mesh, there are difficulties in solving such systems by direct methods when the perturbation parameter, $\varepsilon$, is small [MacLachlan and Madden, 2013]. Therefore, iterative methods are natural choices. However, we show that 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 preconditionings, and show that the condition number of the preconditioned linear system is independent of 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, the incomplete Cholesky preconditioner performs extremely well when $\varepsilon\ll 1$.