+ n = length(x0);
+
+ if (isequal(M, eye(n)))
+ [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations);
+ return;
+ end
+
+ zero_vector = zeros(n, 1);
+
+ k = 0;
+ x = x0; % Eschew the 'k' suffix on 'x' for simplicity.
+ rk = A*x - b; % The first residual must be computed the hard way.
+ zk = M \ rk;
+ dk = -zk;