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); s1=A'*dy1; x1=d.*s1;
dy2=ans(:,2); s2=A'*dy2; x2=ones(n,1)-d.*s2;

%
%  Update the primal solution
%
z1=rho*sum(x1)/(rho-sum(x2));z2=0;z3=z1+1.e10;
for i=1:n,
  if x2(i)>0, 
     z1=max(z1,-x1(i)*rho/x2(i));
  elseif x2(i)==0, 
     z2=min(z2,rho*x1(i));
  else   
     z3=min(z3,-x1(i)*rho/x2(i)); 
  end;
end;
if min([z3-z1,z2,gap-z1])>=0, 
   x=d.*(x1+(z1/rho)*x2);
   z=c'*x
   gap=s'*x;
end;
%
%  Update the dual solution
%  
nora=max( d.*((rho/gap)*s1-s2) );
dy  =(0.5/nora)*((rho/gap)*dy1-dy2);
y   =y+dy;



