gap      =z-b'*y
s     =c-A'*y;
potenf   =rho*log(gap)-sum(log(s))
d        =ones(n,1)./s;
(A*diag(d.*d)*A')\[b A*d];
dy1=ans(:,1);dy2=ans(:,2);
dy       =(rho/gap)*dy1-dy2;
pk       =d.*(A'*dy);
normpk   =norm(pk)
while normpk < .75,
  xk     =(gap/rho)*d.*(ones(n,1)+pk);
  z      =c'*xk;
  gap    =z-b'*y;
   dy    =(rho/gap)*dy1-dy2;
   pk    =d.*(A'*dy);
   normpk=norm(pk);
end;
%
%  Update the dual solution
%
z = z
normpk = normpk 
y     =y+(0.5/normpk)*dy;



