In the example, there are two CG systems. In one, solve() is called with b as parameter. The default initialization is zero. In the other system, solveWithGuess is called with b and x0 (which is a vector of zeros). Hence, both systems start with zero vector as initial point. So, I expected both to take the same step after one iteration.

The default initialization is not zero.

Thanks for the clarification.

Please, change the documentation here:

"By default the iterations start with x=0 as an initial guess of the solution."

It looks like zero vector is going to be the default in upcoming release:

"ConjugateGradient and BiCGSTAB now properly use a zero vector as the default guess."