function [u,uk,res]=CG(A,f,u0,tol,m)
% CG conjugate gradient method
%   [u,uk,res]=CG(A,f,u0,tol,m) solves Au=f using the conjugate gradient 
%   method starting at the initial guess u0 up to a tolerance tol using 
%   at most m iterations. The matrix A has to be symmetric positive 
%   definite. CG returns in the matrix uk the iterates, in u the solution
%   computed, and in res the history of the norm of the residuals.

if nargin<5, m=100; end             % default values
if nargin<4, tol=1e-6; end
uk(:,1)=u0; 
r=f-A*uk(:,1); 
p=r;
oldrho=r'*r; 
res(1)=sqrt(oldrho);
k=0;
while sqrt(oldrho)/norm(f)>tol && k<=m
  k=k+1;
  Ap=A*p;
  alpha=oldrho/(p'*Ap);
  uk(:,k+1)=uk(:,k)+alpha*p; 
  r=r-alpha*Ap;
  rho=r'*r;
  beta=rho/oldrho; oldrho=rho;
  res(k+1)=sqrt(oldrho);
  p=r+beta*p;
end
u=uk(:,k+1);
