function [ z, rnorm ] = cgit ( A, x )

  z = zeros ( size ( x ) );
  r = x;
  rho = r' * r;
  p = r;

  for i = 1 : 15
    q = A * p;
    alpha = rho / ( p' * q );
    z = z + alpha * p;
    rho0 = rho;
    r = r - alpha * q;
    rho = r' * r;
    beta = rho / rho0;
    p = r + beta * p;
  end

  rnorm =  norm ( x - A * z );

  return
end
