% conjugate gradient with initialisation vector
function x = CGx(A,b,x)
  r = b - A*x;
  p = r;
  a = r' * r;
  for k=1:size(A,1)
    v = A * p;
    lambda = a / (v' * p);
    x = x + lambda * p;
    r = r - lambda * v;
    anew = r' * r;
    p = r + anew / a * p;
    a = anew;
  end
