function [x_k,rnormvec] = cg(A,b,P,max_iter);

n = length(b);
x_k = zeros(n,1);

r_k = A*x_k - b;
% Preconditioning step
y_k = P\r_k;
p_k = -y_k;

% Save residual norm
rnormvec = sqrt(r_k'*r_k);

%PCG iteration

for k = 1:max_iter
  Ap_k = A*p_k;
  rty = r_k'*y_k;
  alpha_k = rty / (p_k'*Ap_k);
  x_k = x_k + alpha_k * p_k;
  r_k = r_k + alpha_k*Ap_k;
  % Preconditioning step
  y_k = P\r_k;
  rty_new = r_k'*y_k;
  beta_k = rty_new / rty;
  rty = rty_new;
  p_k = -y_k + beta_k*p_k;
  
  % Save residual norms
  rnormvec = [rnormvec; sqrt(r_k'*r_k)];

end

figure(1)
   semilogy(rnormvec,'o')
   title('residual conjugate gradient')
   xlabel('iteration')