This is an implementation of the conjugate gradient method for
RealLinearOperator. It follows closely the template by Barrett et al. (1994) (figure 2.5). The linear system at
hand is A · x = b, and the residual is r = b - A · x.

A default stopping criterion is implemented. The iterations stop when || r ||
≤ δ || b ||, where b is the right-hand side vector, r the current
estimate of the residual, and δ a user-specified tolerance. It should
be noted that r is the so-called updated residual, which might
differ from the true residual due to rounding-off errors (see e.g. Strakos and Tichy, 2002).

Iteration count

In the present context, an iteration should be understood as one evaluation
of the matrix-vector product A · x. The initialization phase therefore
counts as one iteration.