x      = (n/(c'*x))*x;
gk     = c - ones(n,1)./x;
Xgk    = diag(x)*gk;
Ak     = A*diag(x);
Mk     = Ak*Ak';
yk     = Mk\(Ak*Xgk);
pk     = Xgk-Ak'*yk;
normpk = norm(pk)
y      = yk;
x      = x.*(ones(n,1)-pk);


