function [normpf, normdf, X, y, theta, iter] = fthetaimp( adj,n,indsE,indsEc,indsEAL);
%        [normpf, normdf, X, y, theta, iter] = fthetaimp( adj,n,indsE,indsEc,indsEAL);
% changed rendl's program ..
% input:  w ... node weight vector (of size n) (w > 0 is assumed)
%         adj ... ( list of nonedges of graph (no loops allowed)
% output: theta ... weighted theta function of graph
%         X   ... primal variables
%         y   ... dual variables
% solves: max trace(WX) 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) - W psd
%         where W=sqrt(w)sqrt(w'), and A(ij)= e_i*e_j' + e_je_i' ij nonedge.
% call:   [X, y, theta] = fthetaimp( adj, n,indsE,indsEc);
%         if only adj is provided, then  w = 1, and n = max(max(adj))
%         predictor - corrector version


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

% setting some parameters:
mc = length(indsEc);
m = length(indsE);
w = ones(n,1); 
t0 = clock;                     % take the time
[n, n1] = size( w);             % problem size
[edges, n1] = size( adj);
I = zeros( edges, 1);           % index array (points to matrix entry (i,j)
index = 1:edges;                % auxiliary array
for i = 1:edges;
I( i) = (adj( i,1) - 1) * n + adj( i,2);
end;
e1 = edges + 1;                 % dimension of dual
w1 = sqrt( w);                  % square root of weights
a = zeros( e1,1); a(e1,1) = 1;  % right hand side vector
digits = 8;                     % stop if required accuracy reached

% starting triple (X,y,Z) is feasible and in interior
X = eye(n) / n;
y = zeros( e1, 1); y( e1) = 1.1 * w1'*w1;
z = zeros( n*n,1);
z( I) = y( index);
Z = reshape( z, n, n);
Z = Z + Z' + y(e1) * eye( n) - w1 * w1';

phi = a' * y;           % initial dual cost
psi = w1'* X * w1;      % intial prinak cost
delta = phi-psi;        % delta controls stopping condition
mu = Z( :)' * X( :) /( 2*n);
iter = 0;
text = ['      iter    alphap    alphad    log(gap)  lower      upper'];
disp( text); disp([ iter 1 1 log10(delta) psi  phi ]);

% main iteration
while delta > phi * 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

	i1= adj(:,1); i2=adj(:,2); xx=X(:); zz=Zi(:);
	A = th_system1(n,edges,i1,i2,xx,zz);
	A=reshape(A,edges+1,edges+1);

	for i1=1 : edges;
	  i = adj(i1,1); j=adj(i1,2);
	  tmp=Zi(j,:)*X(:,i) + Zi(i,:)*X(:,j);
	  A(i1,e1)=tmp; A(e1,i1)=tmp;
	end;
	A(e1,e1)= Zi(:)'*X(:);

%norm(A-A')
        A=(A+A')/2;
        [lowtri,p] = chol( A); %clear A;
if p> 0, 
disp(['A not pos. def. in fthetaimp ']);
%keyboard 
end;
	lowtri = sparse(lowtri);
        uptri = lowtri'; 
        rhs = -a;
        ytmp = uptri \rhs; 
	dy = lowtri \  ytmp;
        z = zeros( n*n,1);
        z( I) = dy( index);
        dZ = reshape( z, n, n);

        dZt = dZ'; 
        ddy = find(dZt);   % OURS
        dZ = dZ + dZ' + dy( e1) * eye( n);  %rendk's dZ
        ddz0 = dy( e1);  
  
        tmp1 = Zi * dZ;
        dX = - tmp1 * X - X;
        dX = ( dX + dX')/2;     % symmetrise

         
% OURS:  calculating dx
%        usMatEdy = usMatE(ddy,n,m,indsE);
%        iZusMatX = Zi*(usMatEdy*X);
%        dz0iZX = ddz0*Zi*X;
%        iZusMatX = (iZusMatX+ iZusMatX')/2;
%        muiZmX = mu*Zi-X;
%        rhsx = muiZmX - dz0iZX - iZusMatX;
%        %rhsx = ( rhsx + rhsx')/2;       % symmetrize
%        dx = sqrt(2)*rhsx(indsEc);  
%        usMatEcdx = usMatEc(dx,n,mc,indsEc);
%keyboard   

% corrector step
        tmp = mu * Zi - tmp1 * dX; tmp=(tmp+tmp')/2;
        rhs = zeros( e1,1);
        z = reshape( tmp, n*n, 1);
        rhs( index) = 2 * z( I);
        rhs( e1) = trace( tmp);
%        ytmp = rhs'/ lowtri; ytmp= ytmp'; 
        ytmp = uptri \ rhs; 
	dy1 = lowtri \ ytmp;
        clear lowtri uptri
        z = zeros( n*n,1);
        z( I) = dy1( index);
        dZ1 = reshape( z, n, n);
        dZ1 = dZ1 + dZ1' + dy1( e1) * eye( n);
        dX1 = tmp - Zi * dZ1 * X;
        dX1 = ( dX1 + dX1')/2;     % symmetrise

% the actual steps
        dy = dy + dy1;
        dX = dX + dX1;
        dZ = dZ + dZ1;

%disp([sum(abs(dX(indsE)))]);
%disp([sum(diag(dX))]);
dX(indsEAL) = zeros(2*m,1);
dX(1:n+1:end) = diag(dX) - sum(diag(dX))/n;
%trace(dX)    %   for echo only ????


% 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;
        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 = a' * y;
        psi = w1' * X * w1;
        delta = phi-psi;

% display current iteration
        disp([ iter alphap alphad log10(delta) psi  phi]);
        end;            % end of main loop

theta = phi;
disp('elapsed time:'); disp(etime(clock, t0));
% evaluate feasibilities
normpf=norm(X(indsEc))+abs(trace(X)-1);
normdf=norm(Z);
disp('at end of thetaimp - using improved HKM direction');
disp(['norm(X(indsEc))+abs(trace(X)-1)) is: ',num2str(normpf)]);
