function [x, j] = PCG_newton(Amap,b,D,toler,maxiter)
%Solve for Ax = b, A is symmetric and PD
%Stopping follows LIBLINEAR to use the obj value of Q
%Meaning that we consider
%min_x Q(x) := -b'*x + x'*A*x/2

if (size(D,1) == 1)
	D = D';
end
if (isscalar(D))
	Dmap = @(x) D*x;
elseif (size(D,2) == 1)
	Dmap = @(x) D.*x;
else
	Dmap = @(x) D*x;
end

if (~isa(Amap,'function_handle'))
	Amap = @(x) Amap*x;
end
x = zeros(size(b));
r = b;
%r = -(Ax-b)
z = Dmap(r);
p = z;
newQ = 0;
Q = 0;%Q = (-(-(Ax-b))'*x - b'*x)/2 = -(-r'*x - b'*x)/2
rz = r' * z;
if (nargin >= 4)
	cgtol = min(toler, sqrt(rz));
else
	cgtol = sqrt(rz);
end
if (nargin >= 5)
	cg_maxiter = min(maxiter, length(x));
else
	cg_maxiter = length(x);
end
	
for j = 1:cg_maxiter
	Ap = Amap(p);
	kap = p' * Ap;
	alpha = rz / kap;
	x = x + alpha * p;
	r = r - alpha * Ap;
	newQ = (-r'*x - b'*x) / 2;
	Qdiff = newQ - Q;
	if (newQ <= 0 && Qdiff <= 0)
		if (j * Qdiff >= cgtol * newQ)
			break;
		end
	else
		fprintf('WARNING: CG increased subproblem obj, exit\n');
		break;
	end
	Q = newQ;
	z = Dmap(r);
	newrz = r' * z;
	beta = newrz / rz;
	rz = newrz;
	p = z + beta * p;
end
end

