There are a few parts of the code where it's not clear if we're tracking residual of the original Ax=b system, or the residual of the preconditioned system M^{-1}Ax=M^{-1}b. If it isn't too costly (e.g. no extra residual computations), I would prefer to track the residual of the original system. @BallasticMonkey Would you like to have a look into this? Maybe some tests with larger matrices could help verify this.