function [X,Z,theta] = pcthetaimp( n,indsE,indsEc);

%%   THETA solves the SDP version of the Lovacz-theta function problem
%%      max trace EX    s.t.  trace(X)=1, X_ij=0 for all ij in EE, X psd
%%  where E is the matrix of all ones, EE is the set of edges in the graph
%%   graph is  G=(V,EE) with n=|V| and m=|EE| and mc=(n(n-1)/2)-m;
%  input: n=|V|
%         indsE ... edge set
%         indsEc ... edge set complement
% output: theta ... weighted theta function of graph
%         X   ... primal variables
%         y   ... dual variables
% solves: max trace(EX) s.t. X psd, X(i,j) =0, (ij) nonedge, trace(X)=1
%         min y0 s.t. y0*eye(n) + sum_ij y(ij)A(ij) - E psd
%         E(ij)= (e_i*e_j' + e_je_i')/sqrt(2);  (ij nonedge)
% call:   [X,Z,theta] = pcthetaimp( n,indsE,indsEc);


[iE,jE] = ind2sub([n,n],indsE);
[iEc,jEc] = ind2sub([n,n],indsEc);

mc = length(indsEc);
m = length(indsE);
index = 1:m;
t0 = clock;                     % take the time
w = ones(n,1);
digits = 8;                     % stop if required accuracy reached

% starting triple (d,x,y,Z) is feasible and in interior
X = eye(n)/n;
x = zeros(mc,1);
d = zeros(n-1,1);
z0 = 2*n;
%z0=1.1*w'*w;
y = zeros(m,1);
Z = z0*speye(n)-ones(n); 

phi = z0;               % initial dual cost
psi = sum(sum( X));        % intial primal cost
delta = phi-psi;        % duality gap delta controls stopping condition
mu = Z(:)' * X(:) /(2*n);

iter = 0;
text = ['iter        alphap    alphad    log(gap)     lower                upper'];
disp( text); 
disp([sprintf(' %3g      %6.2e    %3.1e    %5.3e     %13.9e      %13.9e', ...
        iter,1,1,-log10(delta),psi,phi)]);
%disp([ iter 1 1 -log10(delta) psi  phi ]);

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

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

% form the matrix of the system! 
        iZX = Zi*X;
        %alf2 = 1/trace(iZX);
        alf2 = trace(iZX);
	A = []; 
	A1 = []; 
	A2 = []; 
        iZXs =(iZX+iZX')/2;
        alf4 = 1/trace(iZX);

    for cnt=1:m, 
       p=iE(cnt);q=jE(cnt);
       iZMatEX = (   Zi(:,p)*X(q,:)+Zi(:,q)*X(p,:)  )/sqrt(2); 
       %[ iZMatEX ] = iZu2sMatEX( n,Zi,X,iE(cnt),jE(cnt) );   
       [ iZMatEX ] = (iZMatEX+iZMatEX')/2;
       alf1 = trace(iZMatEX)/alf2;
       ls1 = - alf1*iZXs; % left hand side without projection
       ls2 =  iZMatEX; % left hand side without projection
       lhs1 = sqrt(2)*ls1(indsE);   
       lhs2 = sqrt(2)*ls2(indsE);   
       A1 = [A1 lhs1];
       A2 = [A2 lhs2];

       alf3 = trace(iZMatEX);
       ls = - alf3*alf4*(iZX+iZX')/2 + iZMatEX; % left hand side without projection
       lhs = sqrt(2)*ls(indsE);   
       A = [A lhs];
    end;
    Asum=A1+A2;
    % temp. test
    if norm(A-A')>1e-12,
       ['error in A not symmetric  ' num2str(norm(A-A')) ]
       A=(A+A')/2;
       if min(eig(A))<0,
          'error A not p.d.'
       end
     else
       A=(A+A')/2;
       if min(eig(A))<0,
          'error A not p.d.'
       end
     end

        muiZmX = mu*Zi-X;
        trmuiZmX = trace(muiZmX);
        rs = muiZmX - trmuiZmX*alf2*iZX;
        rhs = sqrt(2)*rs(indsE);

        lowtri = chol(A); 
        %clear A;
	lowtri = sparse(lowtri);
        uptri = lowtri'; 
        ytmp = uptri\rhs; 
	dy = lowtri\ytmp;   % dy !!

        usMatEdy = usMatE(dy,n,m,indsE);
        iZusMatX = Zi*usMatEdy*X;

        dz0 = alf2*(- trace(iZusMatX) + trmuiZmX );


        % calculating dx:
        dz0iZX = dz0*iZX;
        iZusMatX= (iZusMatX+ iZusMatX')/2;
        rhsx = muiZmX - dz0iZX - iZusMatX;
        %rhsx = ( rhsx + rhsx')/2;       % symmetrize
        dx = sqrt(2)*rhsx(indsEc);  
        usMatEcdx = usMatEc(dx,n,mc,indsEc);

        % calculating dd:
        rhsdd = muiZmX - dz0iZX -iZusMatX; %change this because it is dag(Vdd) 
        ddiag = diag(rhsdd); 
        dd = ddiag(2:n); 

        dX = usMatEcdx + diag(ddiag); 
        dZ = dz0*eye(n) + usMatEdy;


% find steplengths alphap and alphad
        alphap =  1;

        [dummy,posdef] = chol( X + alphap * dX);
        while posdef > 0,
                alphap = alphap * .8;
		[dummy,posdef] = chol( X + alphap * dX);
               end;
        if alphap < 1, alphap = alphap * .95; end;

        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;
        x = x +  alphap *dx;
        z0 = z0 + alphad *dz0;
        y = y + alphad * dy;
        Z = Z + alphad * dZ;
        mu = .6 * X( :)' * Z( :) / n;
        if alphap + alphad > 1.5, mu = mu * .6; end;
        if alphap + alphad > 1.7, mu = mu * .4; end;

        phi = z0;
        psi = sum(sum(X));
        delta = phi-psi;

% display current iteration
disp([sprintf(' %3g      %6.2e    %3.1e    %5.3e     %13.9e      %13.9e', ...
        iter,alphap,alphad,-log10(delta),psi,phi)]);
        %disp([ iter alphap alphad -log10(delta) psi  phi]);
        end;            % end of main loop

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