function [dX,dy,dZ]=slvcgkkt(mu,A,TA,b,C,X,y,Z,Fd,Fp,ZFZX)
% Solve the KKT conditions for the step direction

%  Use conjugate gradients approach
%  
%  Call:  [dX,dy,dZ]=slvcgkkt(mu,A,TA,b,C,X,y,Z,Fd,Fp,ZFZX)

digitsiter= 3;                     %  significant digits of iteration
m=length(b);n=size(X,2);n2=size(A,2);tn=size(TA,2);
rhsnorm=norm(Fd,'fro')+norm(Fp)+norm(ZFZX,'fro');   % norm of rhs
%%  initial estimates use range of operator K'
%dX=-reshape(A'*ones(m,1),n,n) + Z;
%dy=A*reshape(eye(n),n^2,1);
%dZ=eye(n) + X;
% try initial estimate is 0
dX=zeros(n);
dy=zeros(m,1);
dZ=zeros(n);
% calculate update
%AdXn=-reshape(A'*dy,n,n) + dZ;
%Adyn=-A*reshape(dX,n^2,1);
%AdZn=Z*dX+dZ*X;
% residuals first
rdX=-Fd;
rdy=-Fp;
rdZ=-ZFZX;
% now p and s
pX=(-reshape(A'*rdy,n,n) + (Z*rdZ+rdZ'*Z)/2   );
py=(-A*reshape(rdX,n^2,1));
pZ=(rdX+(rdZ*X + X*rdZ')/2  );
sX=pX;
sy=py;
sZ=pZ;
gamma=norm(sX,'fro')+norm(sy)+norm(sZ,'fro');   % norm of rhs
change=2*max(1,rhsnorm)*10^(-digitsiter);
iter=0;
while change > max(1,rhsnorm)*10^(-digitsiter),
iter=iter+1;
qX=-reshape(A'*py,n,n) + dZ;
qy=-A*reshape(pX,n^2,1);
qZ=Z*pX+pZ*X;
normq=norm(qX,'fro')+norm(qy)+norm(qZ,'fro');   
alpha=gamma/(normq^2);
dXn=dX+alpha*pX;
dyn=dy+alpha*py;
dZn=dZ+alpha*pZ;
change=norm(dX-dXn,'fro')+norm(dy-dyn)+norm(dZ-dZn,'fro'); 
rdXn=rdX-alpha*qX;
rdyn=rdy-alpha*qy;
rdZn=rdZ-alpha*qZ;
sXn=(-reshape(A'*rdyn,n,n) + (Z*rdZ+rdZ'*Z)/2   );
syn=(-A*reshape(rdXn,n^2,1));
sZn=(rdX+(rdZn*X + X*rdZn')/2  );
gamman=norm(sXn,'fro')^2+norm(syn)^2+norm(sZn,'fro')^2; 
beta=gamman/gamma;
pXn=sXn+beta*pX;
pyn=syn+beta*py;
pZn=sZn+beta*pZ;

% update
dX=dXn;
dy=dyn;
dZ=dZn;
pX=pXn;
py=pyn;
pZ=pZn;
sX=sXn;
sy=syn;
sZ=sZn;
rdX=rdXn;
rdy=rdyn;
rdZ=rdZn;
gamma=gamman;

disp(['iteration number and change norm in the CG iteration scheme'])
disp([iter change])

%keyboard
%pause
end   % end of while loop
disp(['number of iterations in the CG iteration scheme'])
disp([iter])
keyboard
pause
