function [x,kk,err]=conjgra(A,b,x0,kmax,tol)
%
% Conjugate Gradient method for linear systems
%
% Input
%  A     Matrix (square matrix)
%  b     r.h.s (column vector)
%  x0    initial guess (column vector)
%  kmax  maximum number of iterations
%  tol  Tolerance for the stopping test
%
% Output
%  x    solution (column vector)
%  kk  Iterations
%  err   array with errors at each iteration
%
%    [x,kk,err]=conjgra(A,b,x0,kmax,tol)
%
x=x0;
r=b-A*x;
nrk=r'*r;
p=r;
bb=norm(b,2);
%
k=0;err1=sqrt(nrk)/bb;
err=[err1];
while k<kmax & err> tol
q=A*p;
alpha=nrk/(p'*q);
%
x=x+alpha*p;
%
r=r-alpha*q;
nrk1=r'*r;
beta=nrk1/nrk;
p=r+beta*p;
res=norm(r,2);
err1=res/bb;
err=[err;err1];
k=k+1;
nrk=nrk1;
end
kk=k;
fprintf('\n The CG has satisfied the stopping test in')
fprintf('%2g iterations,\n with tolerance = %6.3g\n', kk, tol)

return
