$\begingroup$Based on some literature I tried to implement a C code based on a Matlab code for conjugate gradient solver. However, When I tried to compare the code published on Wikipedia with MATLAB pcg() function as follows:
rng default A = sprand(100,100,.5); A = A'*A; b = sum(A,2); x1 = pcg(A, b, 1e-4, 150); x2 = conj_grad(A, b, 1e-4, 150); function x_i = conj_grad(m_A, v_b, tol, iter) x_i = v_b; r = v_b - m_A * x_i; p = r; rsold = r' * r; for i = 1:iter Ap = m_A * p; alpha = rsold / (p' * Ap); x_i = x_i + alpha * p; r = r - alpha * Ap; rsnew = r' * r; if sqrt(rsnew) < tol fprintf('cg converged at iteration %d to a solution with relative residual %f', i, sqrt(rsnew)); break; end p = r + (rsnew / rsold) * p; rsold = rsnew; end end I discovered that the custom made function does not converge as pcg() and the results are somehow different, for example:
conj_grad() pcg() ---------------------- 1.057764 1.050005 0.959692 0.957766 1.055046 1.017192 0.990600 1.012397 0.945796 0.966080 1.017147 0.989121 0.957176 0.975864 1.084515 1.087173 1.057552 1.047540 1.003535 1.005227 1.010863 0.995646 Could someone give some hints, please? Thank you
$\endgroup$ Reset to default
