Why won't my conjugate gradient algorithm work?

Question

I made this simple Conjugate Algorithm on Matlab

n = length(b);

r0 = b - A*x0;

p0=r0;

k=1;

n0=(r0')*r0; while n0 >= eps && k <= n

interm = A*p0;

alpha = n0/((p0')*interm);

x0=x0+alpha*p0;

r0=r0-alpha*interm;

n1=(r0')*r0;

if n1 <=eps

    y=x0;

    break

end
betta = n1/n0;

p0=r0+betta*p0;

k=k+1;

n0=n1;

end

I checked it dozen of times and saw through a couple of sources but couldn't find any mistake, and either way it doesn't work when I test it on the command window. Anyone knows what is wrong?

Thanks in advance

Answer 1

I think one problem is that eps is machine epsilon in Matlab, and you can't expect the residual to be that small. You could try introducing a variable $\text{epsilon} = 10^{-8}$ or something, and using that instead of eps.