I tested my conjugate gradients implementation with float and double precision and contrary to my guess the double code was twice faster than the single precision code. The reason is that I need many more iterations to achieve the same residual. Mind you, that was already on a pretty small system of about just 2000 unknowns (though the matrix is ill-conditioned). Considering this, are there any remedies to this? For example I looked here and they simply implement the dot product in double precision. Are there any approaches that are accurate and efficient in practice (not just theoretically) that work with single precision (e.g. if I want to use this on the GPU where double comes at a premium)? Or maybe there are some reformulations of the conjugate gradients algorithm that make it more stable (reorthogonalization?).