In some of my numerical programming using conjugate gradient solvers, I noticed an alarming problem: The residuals were not monotonically decreasing to zero, but were sometimes increasing. In this document, the author writes
About fifteen years ago enough experience had been gained with matrix preconditoning that variants on the Conjugate Gradient method had come into wide use. Evolution of this methodology has continued with the introduction of several variations on the basic algorithm. The most popular of these is currently Sonneveld's conjugate gradient squared (CGS) algorithm. This class of methods has a rate of convergence that is generally very good, but is not monotonic. Plots of residual versus iteration count can show oscillations.
So I guess it might not be an error in my code. However, given the structure of the CG algorithm and the principles on which it was based, I cannot understand why the residuals don't decrease monotonically with iteration number. Is it finite arithmetic?