Why is my custom Conjugate gradient solver and MATLAB pcg() do not have the same result?

10 views (last 30 days)
Hello
Based on some literature I tried to implement a C code based on this 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(400,400,.5);
A = A'*A;
b = sum(A,2);
x1 = pcg(A, b, 1e-4, 150);
x2 = conj_grad(A, b, 1e-4, 150);
cmp = abs(x1-x2) < 1e6*eps(min(abs(x1),abs(x2))); % to compare
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 different for example:
pcg() conj_grad()
---------------------------------------------
0.9906 0.9959
1.0118 1.0056
0.9989 0.9999
1.0201 0.9970
1.0049 1.0139
......
Could someone give some hints please?
Thank you

Answers (0)

Categories

Find more on Particle & Nuclear Physics in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!