function [R,Z,VT,t] = mainHKM(A,B,C);

% min tr AXBX^t + CX^t
% 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') = E_oo; R >= 0       
% (D) max -t_0 
%     s.t. L+ VT'*G(t)*VT - Z = 0; Z >= 0,
% 14.6.2002, wolkowicz & sotirov, waterloo 
% gangster indices: [Jb,i,j] = Jbar(n);
% call: [R,Z,VT,t] = mainHKM(A,B,C);


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

[Jb,P,Q] = Jbar(n);
P = [1;P];
Q = [1;Q];
eq = length(P); 

% 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 = [ 1, zeros(1,(n-1)*(n-1));  zeros((n-1)*(n-1),1),kron(n*eye(n-1)-ones(n-1),n*eye(n-1) - ones(n-1))/(n*n*(n-1))];   %clanka
Di=inv(MD);
R = Di*R*Di';

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)


[Y,t] = initial(L,VT,P,Q,eq,n,Jb);
Z =  L + VT'*Y*VT; 
t(2:end)=t(2:end)*2; % just fixing for the first step


done = 0;
k1 = 0;

disp(' ');
disp('------------------------------------------------------------------------- ');
disp('      min tr AXBX^t + CX^t,  G_g( VT*R*VT^t)= E_oo, R \succ 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); 

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;
dZ = (dZ + dZ')/2;
dR = invZ*(mu*eye((n-1)*(n-1)+1) - dZ*R ) - R; 
VdRV = VT*dR*VT';

% FOR TESTING:
% 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


dR = ( dR + dR' )/2;                       % delta R
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;
t = t + t1*dt;  
[GTt] = GT(t,P,Q,eq,n);

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 %9.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;  
%*********************************

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

end;   % while .... done


 disp(' ');
disp(['norm dual residual is: ',...
   num2str(norm(  (L + VT'*GTt*VT)-Z))]);
[Jbr,ir,jr] = gangster(n-1); 
indsG = find(Jbr);
disp(['norm primal residuals, Gangster,Arrow are: ',...
   num2str([norm(R(indsG)) norm(diag(R)-R(:,1))])]);



