gap     = z - b'*y
S       = C-A'*diag(y)*A;
potenf  = rho*log(gap)-log(det(S))
R       = chol(S);
R       = inv(R);
B       = A*R;
M       = B*B';
(M.*M)\[b diag(M)];
dy1     = ans(:,1);dy2=ans(:,2);
dy      = (rho/gap)*dy1-dy2;
Pk      = B'*diag(dy)*B;
normpk  = norm(Pk,'fro')
while normpk < .75, 
  Xk     = (gap/rho)*(R*(eye(n)+Pk)*R');
  z      = trace(C*Xk);
  gap    = z-b'*y;
  dy     = (rho/gap)*dy1-dy2;
  Pk     = B'*diag(dy)*B;
  normpk =norm(Pk,'fro')
end;
%
%  Update the dual solution
%
%z      = z
%normpk = normpk 
y       = y+(0.5/normpk)*dy;



