% preconditioned conjugate gradient
% with initialisation vector
function x = CGxPre(A,b,x,P)

  r = b - A*x;
  p = P * r;
  a = r' * p;
  for k=1:size(A,1)
    v = A * p;
    lambda = a / (v' * p);
    x = x + lambda * p;
    r = r - lambda * v;

    z = P * r;
    anew = r' * z;
    p = z + anew / a * p;
    a = anew;
  end
