function x = zhangSIAM00(c, A, b)

% function zhangSIAM00(c, A, b)
%
% min { <c,x>: Ax = b, x >= 0 }

A = sparse(A);
t0=cputime; 
[m,n]=size(A); 
x=sqrt(n)*ones(n,1); 
y=zeros(m,1); 
z = x; 
p = symmmd(A*A'); 
bc = 1 + max(norm(b),norm(c)); 

for iter = 1:100 
  Rd=A'*y+z-c;
  Rp=A*x-b;
  Rc=x.*z;
  residual=norm([Rd;Rp;Rc])/bc;
  fprintf('iter %2i: residual = %9.2e',iter,residual);
  fprintf('\tobj=%14.6e\n',c'*x);
  if residual<5.e-8
     break;
  end;
  gap=mean(Rc); 
  Rc=Rc-min(.1,100*gap)*gap;
  d=min(5.e+15,x./z);
  B = A*sparse(1:n,1:n,d)*A';
  R = cholinc(B(p,p),'inf');
  t1 = x.*Rd - Rc;
  t2 = -(Rp + A*(t1./z));
  dy = zeros(m,1);
  dy(p)=R\(R'\t2(p));
  dx=(x.*(A'*dy)+t1)./z;
  dz=-(z.*dx+Rc)./x;
  tau = max(.9995, 1-gap);
  ap = -1/min(min(dx./x),-1);
  ad = -1/min(min(dz./z),-1);
  ap = tau*ap;
  ad = tau*ad;
  x = x + ap*dx;
  z = z + ad*dz;
  y = y + ad*dy;
end 

fprintf('Done!\t[m n] = [%g %g]\tCPU = %g\n',m,n,cputime-t0);