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

%  Use Richardson iteration for the solution of the lss problem
%  
%  Call:  [dX,dy,dZ]=pslvkkt(mu,A,TA,b,C,X,y,Z,Fd,Fp,ZFZX)

digitsiter= 3;                     % 6 significant digits of iteration
alpha=.01;   % Richardson parameter - guess for now!!!
m=length(b);n=size(X,2);n2=size(A,2);tn=size(TA,2);
[ii,jj]=find(triu(ones(n)));
trind=(jj-1)*n+ii;       % column vector with indices for triu
tFd=Fd(trind);          % vector of upper triang part
rhsnorm=norm(Fd,'fro')+norm(Fp)+norm(ZFZX,'fro');   % norm of rhs
%  initial estimates use range of operator K'
dx=-(A'*ones(m,1),n,n) + Z;
dy=A*reshape(eye(n),n^2,1);
dZ=eye(n) + X;
% 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-AdXn;
rdy=-Fp-Adyn;
rdZ=-ZFZX-AdZn;
% now updates
dXn=dX+alpha*(-reshape(A'*rdy,n,n) + (Z*rdZ+rdZ'*Z)/2   );
dyn=dy+alpha*(-A*reshape(rdX,n^2,1));
dZn=dZ+alpha*(rdX+(rdZ*X + X*rdZ')/2  );
change=norm(dX-dXn,'fro')+norm(dy-dyn)+norm(dZ-dZn,'fro'); 
dX=dXn;dy=dyn;dZ=dZn;
iter=0;
while change > max(1,rhsnorm)*10^(-digitsiter),
iter=iter+1;
% 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-AdXn;
rdy=-Fp-Adyn;
rdZ=-ZFZX-AdZn;
% now updates
dXn=dX+alpha*(-reshape(A'*rdy,n,n) + (Z*rdZ+rdZ'*Z)/2   );
dyn=dy+alpha*(-A*reshape(rdX,n^2,1));
dZn=dZ+alpha*(rdX+(rdZ*X + X*rdZ')/2  );
change=norm(dX-dXn,'fro')+norm(dy-dyn)+norm(dZ-dZn,'fro');
dX=dXn;dy=dyn;dZ=dZn;


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

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