function [phi, X, y] = sdps(A,TA,b,C);
% solves: max trace(CX) s.t. X psd, A(X) = b; 
%         min b'y       s.t. A^*(y) - C psd, y unconstrained,
% input:  C ... symmetric matrix
%         A ... m symmetric matrices in rows of A
%         b ... m vector
%         initx  function for initial X > 0
%         inityz  function for initial Z > 0, y
%         slv??kkt function for finding step direction by solving KKT cond.
%                 (several different choices)
% output: phi ... optimal value of primal, phi =trace(CX)
%         X   ... optimal primal matrix
%         y   ... optimal dual vector
% call:   [phi, X, y] = sdps(A,TA,b,C);

digits = 6;                     % 6 significant digits of phi
[n, n1] = size( C);             % problem size
n2=n^2;
m=size(A,1);
[X] = initx(A,b);        % initial primal matrix is pos. def.
x=reshape(X,n2,1);
[y,Z] = inityz(A,C);     % initial dual matrix is pos. def.
Fd=Z+C-reshape(A'*y,n,n);     % dual feasibility residual
Fp=b-A*x;                 % primal feasibility residual
mu = Z( :)' * X( :)/( 2*n);     % initial complementarity
XFZX=Z*X-mu*eye(n);           % complementary slackness residual
perror=norm(Fp);
derror=norm(Fd);
cerror=norm(XFZX,'fro');
phi = b'*y;                     % initial dual
psi = C(:)' * X( :);            % and primal costs
iter=0;                         % iteration count
%


disp(['      iter    alphap    alphad      gap     lower      upper perr derr']);

%  need correction for scaling feasibility errors????
% change phi-psi to cerror?? since infeasible algorithm???
while ((phi-psi)+derror+perror) > max([1,abs(phi)]) * 10^(-digits)

      iter = iter + 1;          % start a new iteration
      [dX,dy,dZ]=pslvcgkkt(mu,A,TA,b,C,X,y,Z,Fd,Fp,XFZX);
%'in sdps'
%keyboard
%pause

% line search on primal
      alphap = 1;               % initial steplength
      [dummy,posdef] = chol( X + alphap * dX );  % test if pos.def
      while posdef > 0,
              alphap = alphap * .8;
              [dummy,posdef] = chol( X + alphap * dX );
              end;
      if alphap < 1, alphap = alphap * .95; end; % stay away from boundary
% line search on dual;
      alphad = 1;
      [dummy,posdef] = chol( Z + alphad * dZ);
      while posdef > 0;
              alphad = alphad * .8;
              [dummy,posdef] = chol( Z + alphad * dZ);
              end;
      if alphad < 1, alphad = alphad * .95; end;
% update
      X = X + alphap * dX;
      y = y + alphad * dy;
      Z = Z + alphad * dZ;
      mu = X( :)' * Z( :) / (2*n);
      Fd=Z+C-reshape(A'*y,n,n);     % dual feasibility residual
      Fp=b-A*x;                 % primal feasibility residual
      perror=norm(Fp);
      derror=norm(Fd);
      if alphap + alphad > 1.8, mu = mu/2; end; % speed up for long steps
      phi = b' * y; psi = C( :)' * X( :);
% display current iteration
        disp([ iter alphap alphad (phi-psi) psi  phi perror derror]);
'in sdps'
keyboard
pause

%        test         % call test on rank of Hessians
        end;            % end of main loop
