- x0 = zeros(5, 1);
- b = unifrnd(-1000, 1000, 5, 1);
- [o_x, o_flag, o_relres, o_iter] = pcg(A, b, tolerance, max_iterations, C, C');
+ x0 = zeros(n, 1);
+ b = unifrnd(-100, 100, n, 1);
+ g = @(x) A*x - b;
+
+ ## pcg() stops when the /relative/ norm falls below tolerance. To
+ ## eliminate the relativity, we divide the tolerance by the
+ ## quantity that pcg() will divide by.
+ [o_x, o_flag, o_relres, o_iter] = pcg(A, b, tolerance/norm(g(x0)), ...
+ max_iterations, C, C');