- zero_vector = zeros(length(x0), 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.
- pk = -rk;
-
- for k = [ 0 : max_iterations ]
- if (norm(rk) < tolerance)
- % Success.
- return;
- end
-
- alpha_k = step_length_cgm(rk, A, pk);
- x_next = x + alpha_k*pk;
- r_next = rk + alpha_k*A*pk;
- beta_next = (r_next' * r_next)/(rk' * rk);
- p_next = -r_next + beta_next*pk;
+ % The rather verbose name of this function was chosen to avoid
+ % conflicts with other implementations.
+ %
+ n = length(x0);
+ M = eye(n);