conjgrad.m

function[x]=conjgrad(A,b,x0,ni)
%%
%ni: number of iteration
%%
x=[];
r=[];
p=[];
alfa=[];
beta=[];
n=size(A);
n=n(1,1);
k=1;
for nn=1:n
    x(nn,k)=x0(nn,1);
end
r(:,k)=b-A*x(:,k);
p(:,k)=r(:,k);

for k=1:ni
    alfa(:,k)=r(:,k)'*r(:,k)/(p(:,k)'*A*p(:,k));
    x(:,(k+1))=x(:,k)+alfa(:,k)*p(:,k);
    r(:,(k+1))=r(:,k)-alfa(:,k)*A*p(:,k);
    if (norm(r(:,(k+1)))<1e-10)
        break;
    end
    beta(:,k)=r(:,(k+1))'*r(:,(k+1))/(r(:,k)'*r(:,k));
    p(:,(k+1))=r(:,(k+1))+beta(:,k)*p(:,k);
end
x=x(:,(k+1));