function [X, y, psi] = psd( C,a,store_a,nonz);
% general purpose ip-program for semidefinite optimization
% input:  C ... primal matrix ( sym)
%         a ... right hand side vector
% output: phi ... optimal dual cost
%         X   ... primal matrix
%         y   ... dual variables
% solves: max trace(CX) s.t. X psd, A(X) = a;
%         min y'a       s.t. A'(y) - C psd;
% uses :  pd_check  (for pos.def test)
%         pd_solve  (for solution of pos.def. system)
%         op_a     ... evaluates A(X)
%         op_at    ... evalueates A'(y)
%         op_biga  ... forms system matrix
%         psd_strt ... generates a feasible interior point [X,y,Z]
%         psd_feas ... checks feasibility of starting point
% call:   [X, y, psi] = psd( C,a,store_a, nonz);
%         predictor - corrector version
% written by: Franz Rendl, TU Graz, Mathematics, A -8010 Graz, Austria
t0 = clock;
[n, n1] = size( C);             % problem size
[m, n1] = size( a);             % size of dual
digits = 7;                     % accuracy

% starting has to be feasible and in interior
[X, y,Z] = psd_strt( C, a, store_a, nonz);
okay = psd_feas( X,y, Z, C, a, store_a, nonz);
if okay > 0, break, end;        % check interiority and feasibility
phi = a' * y;
psi = trace(C * X);
delta = phi-psi;                % stopping criterion
mu = trace( Z * X) /( 2*n);     % mu drives iterations
iter = 0;
tx=['      iter    alphap    alphad   log(gap)  primal     dual     unsym'];
disp( tx); disp([ iter 1 1 log10(delta) psi  phi ]);

while delta > (abs(phi)+1) * 10^(-digits)  % while duality gap is too large

        iter = iter + 1;        % start new iteration
        Zi = inv( Z);
        Zi = (Zi + Zi')/2;      % Zi is symmetric

% form system matrix
        A = op_biga( X, Zi, store_a, nonz);
        unsym = norm(A-A',1);
        unsym = max([ unsym, 1e-99]);
        if unsym > 1e-7,
                disp('numerical instability.');
                break;
                end;
%        A = ( A + A')/2;        % A should be pos.def. and sym.

% predictor step (mu = 0 )
        uptri = chol ( A);
        uptri = uptri';
        rhs = -a;
        dy = pd_solve( uptri, rhs);
        dZ = op_at( dy, store_a, nonz, n);
        tmp1 = Zi * dZ;
        dX = - tmp1 * X - X;
        dX = ( dX + dX')/2;     % symmetrise
% corrector step
        tmp = mu * Zi - tmp1 * dX;
        rhs = op_a( tmp, store_a, nonz);
        dy1 = pd_solve( uptri, rhs);
        dZ1 = op_at( dy1, store_a, nonz, n);
        dX1 = tmp - Zi * dZ1 * X;
        dX1 = ( dX1 + dX1')/2;     % symmetrise
% the actual steps
        dy = dy + dy1;
        dX = dX + dX1;
        dZ = dZ + dZ1;
        okay = 0;
% find steplengths alphap and alphad
        alphap =  1;
        posdef = pd_check( X + alphap * dX);
        while posdef == 0,
                alphap = alphap * .8;
                if alphap < 1e-4, okay = 1; break; end;
                posdef = pd_check( X + alphap * dX);
                end;
% stay away from boundary
        if alphap < 1, alphap = alphap * .95; end;

        alphad = 1;
        posdef = pd_check( Z + alphad * dZ);
        while posdef == 0;
                alphad = alphad * .8;
                if alphad < 1e-4, okay = 1; break; end;
                posdef = pd_check( Z + alphad * dZ);
                end;
% stay away from boundary
        if alphad < 1, alphad = alphad * .95; end;

% update
        if okay == 1; disp('convergence problems.'); break; end;
        X = X + alphap * dX;
        y = y + alphad * dy;
        Z = Z + alphad * dZ;
        mu = trace( X * Z) / (2*n);
        if alphap + alphad > 1.5, mu = mu/2; end;
        if alphap + alphad > 1.8, mu = mu/2; end;
        phi = a' * y;
        psi = trace( C * X);
        delta = phi-psi;

% display current iteration
        disp([ iter alphap alphad log10(delta) psi  phi log10(unsym)]);

        end;            % end of main loop

disp('elapsed time:'); disp(etime(clock, t0));
