function [R,Z,VT,t] = impmain1(A,B,C);
% for improved  HKM!
% min tr AXBX^t + CX^tal
% input: A = A^t ... n*n
%        B = B^t ... n*n
%        C       ... n*n      
%        matrix R: ((n-1)^2+1)*((n-1)^2+1)
%        matrix Z: ((n-1)^2+1)*((n-1)^2+1)
% (P) min LR  s.t.  G_{P,Q}(VT*R*VT') = 0; R >= 0       
% (D) max w_0 
%     s.t. L+ G(VT*R*VT) - Z = 0; Z >= 0,  w... n*n+1
% 14.6.2002, waterloo  only gangster part
% gangster indices: [Jb,i,j] = Jbar(n);
% version that works!
% call: [R,Z,VT,t] = impmain1(A,B,C);


n = size(A,1);   %size of the problem
stoptoler=1e-7;

% indices for both (n^2+1) \times (n^2+1) and ((n-1)^2+1) \times ((n-1)^2+1)
[Jb,P,Q] = Jbar(n);            % indices for RT ...  n^2+1
P = [1;P];
Q = [1;Q];
eq = length(P);

[Jbr,ir,jr] = gangster(n-1);   Jbr=Jbr+Jbr';
%[Jbr,ir,jr] = Jbar(n-1);        % indices for R ... (n-1)^2+1
indsG = find(Jbr);
indsGc = find(Jbr==0);
nnzdR=size(indsGc,1);
brr = size(indsG,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initialization
t = zeros(eq,1);           % belongs to gangster part
a = zeros(n*n+1,1); 
a(1,1) = -1;
b = zeros(eq,1); b(1,1) = -1;

 
% old definition of V
V = [ eye( n-1 ); - ones( n-1,1)'];         
e = ones(n,1);
VT = [1, zeros(1,(n-1)*(n-1));kron(e,e)/n, kron(V,V) ];   
VT = sparse(VT);

% new definition of V to guarantee zeros in R:
nm1=(n-1)^2;
MD=[1 spalloc(1,nm1,0);-ones(nm1,1)/n speye(nm1)];
VT=VT*MD;

R=spalloc((n-1)^2+1,(n-1)^2+1,n*(n-1)^2+1);
R = [ 1, spalloc(1,(n-1)*(n-1),0);...
    spalloc((n-1)*(n-1),1,0),kron(n*speye(n-1)-ones(n-1),...
                       n*speye(n-1) - ones(n-1))/(n*n*(n-1))];   %clanka
Di=inv(MD);
R = Di*R*Di';
R(indsG)=zeros(brr,1);

c = C(:);
LT = [ 0, c'/2; c/2, kron(B,A)];   %dimension (n^2+1)*(n^2+1)
L = VT'*LT*VT;                     %dimension ((n - 1)^2 + 1)*((n - 1)^2 + 1)


[Z,t] = initial(L,VT,P,Q,eq,n,Jb);


done = 0;
k1 = 0;

disp(' ');
disp('-------------------------------------------------------------- ');
disp('min tr AXBX^t + CX^t,  arrow(R) = a, G_g( VT*R*VT^t)=0 ');
disp('-------------------------------------------------------------- ');
disp(' num      L(:)^t*R(:)      dual    mu        t1       t2     relres')
disp('-------------------------------------------------------------------- ');

RVT = R*VT';
RT = VT*RVT;       % RT = VT*R*VT'

% calculating initial mu:
mu = 1.5*Z(:)'*R(:)/((n-1)*(n -1) + 1); 

dR=spalloc((n-1)^2+1,(n-1)^2+1,nnzdR);
while done < 1;             % while done < 1
  k1 = k1 + 1;
 
% left side
invZ = inv(Z);
VTiZ = VT*invZ;
ZT = VTiZ * VT';

[M] = left(ZT,RT,P,Q,eq,n);    % forming the matrix of op.al

% right side
[D,GTt,U1] = right(mu,n,VT,Z,invZ,R,L,VTiZ,t,RT,a,P,Q,eq,b);

[Mpom,p] = chol(M);

if p> 0, 
disp('problem with cho');
keyboard
end;

wttemp = (Mpom')\ D;
rez = Mpom \ wttemp;                        % delta omega 
dt = rez;

if norm(M*rez-D)> stoptoler; 
disp(['relaccuracy: ',num2str(norm(M*rez-D)/norm(D))]); end;

[GTt] = GT(dt,P,Q,eq,n);
U2 =  VT'*GTt*VT;
dZ = U1 + U2;     %????huh? not perfectly feasible????
dZ = (dZ + dZ')/2;
tempdR = invZ*(mu*eye((n-1)*(n-1)+1) - dZ*R ) - R; 

% FOR TESTING:
%VdRV = VT*dR*VT';
% pom = 2:n*n+1;
% L1 = diag(VdRV) - [0; (VdRV(pom,1) + VdRV(1,pom)')/2];
% R1 = - (diag(RT) - [0; (RT(pom,1) + RT(1,pom)')/2]) - a; nn1=norm(L1-R1);
%%--------------------------------------------------------------------
% L2 = VT'*(ATw+GTt)*VT-dZ; R2 = - U1; nn2=norm(L2 - R2);
%%--------------------------------------------------------------------
% L3=  G(VdRV,P,Q,eq,n); R3 =  G(RT,P,Q,eq,n); nn3=norm(L3-R3);
%%--------------------------------------------------------------------
% L4 = dZ*R + Z*dR; R4 = mu*eye((n-1)*(n-1) + 1) - Z*R; nn4=norm(R4-L4);

%nnn=[nn1 nn2 nn3 nn4];
%if max(nnn) >1e-3; disp('error'); keyboard; end


tempdR = ( tempdR + tempdR' )/2;                       % delta R

% fixing for the improved version:
dR(indsGc) = tempdR(indsGc);  %setting  nonzeros
%dR(indsG) = zeros(brr,1);  %setting  zeros on GANG positions
dim = (n-1)^2+1;
d1 = dR(1:dim+1:end)'; 
d = d1(2:dim);
c = dR(2:dim,1); 
pom = (d+c)/2;
dR(1:dim+1:end) = [0; pom];
dR(2:dim,1) = pom;
dR(1,2:dim) = pom';


VdRV = VT*dR*VT';
t1 = 1;
stop1 = 0;
while stop1 < 1
   Ztmp = Z + t1*dZ;
   [Z1,p] = chol(Ztmp);
  
   if p == 0, stop1 = 1;
      else
     t1 = t1 * 0.8; 
      end; 
   if t1 < 1e-7 stop1 = 1; disp('problem with stepsize...');  end;
end;
if t1 <1; t1 = t1 * 0.95; end;
                   % to make sure the next point is not to close to boundary
%Z = Z + t1*dZ;    % old way
t = t + t1*dt;  
[GTt] = GT(t,P,Q,eq,n);
Z = L + VT'*GTt*VT;      % added to get exact dual feasibility by Henry

t2 = 1;
stop2 = 0;
Rpom = VT*R*VT';
while stop2 < 1
  Rtmp = R + t2*dR;
  [R1,p] = chol(Rtmp);  % first condition
  
  RTtmp = Rpom + t2*VdRV;
  
  if p == 0, stop2 = 1;
     else
   t2 = t2 * 0.8; 
  end; 
end; 
if t2 <1; t2 = t2 * 0.95; end;
                     %to make sure the next point is not to close to baundary
                   
R = R + t2*dR;
RT = Rpom + t2*VdRV;

if  k1 >= 100;   done = 1; disp(' max iter.');  end;

pom = 2:n*n+1;
arrow1 = diag(RT) - [0; (RT(pom,1) + RT(1,pom)')/2];
proba = norm(L + VT'*GTt*VT - Z,1) + norm( arrow1+a,1);

if proba >0.1; mu_fact = .95; 
elseif min(t1,t2) < .2;
  mu_fact=1.05;
elseif min(t1,t2) > .6;
  mu_fact = .25;
  else
  mu_fact = .8;
end;

[GRT] = G(RT,P,Q,eq,n);
mu = mu * mu_fact;
dual = -t(1);

%****************************************
primal = L(:)'*R(:);

fprintf(' %3.0d %12.5f %12.5f %12.7f %8.4f %8.4f %8.4e \n',[k1, L(:)'*R(:) dual mu t1  t2 (norm(M*rez-D)/norm(D))]);


if (proba < stoptoler) & (abs(primal-dual)/max(1,abs(primal))<stoptoler);  
  done = 1;  
end;  

if  (mu<10^-9);  
 done = 1;
end;  
%*********************************

%fprintf(' %3.0d %12.5f %12.5f %12.7f %8.4f %8.4f %8.4f\n',[k1, L(:)'*R(:) dual mu t1  t2 proba]);

%if  (proba < stoptoler) & (mu<stoptoler) ;  done = 1;  end;  

end;   % while .... done


%disp(' ');
%disp('<t,b> =');
%fprintf('              %12.6f\n',t'*b);


