function [R,Z,VT,w,t, primal, dual] = main2(A,B,C,P,Q);

% 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.  arrow(VT*R*VT') + a=0, G_{P,Q}(VT*R*VT') = 0; R >= 0       
% (D) max w_0 
%     s.t. L+ Arrow(w) + G(VT*R*VT) - Z = 0; Z >= 0,  w... n*n+1
% 12.1.2002   ARROW PART and EQUALITIES (gangster part)
% [Jb,i,j] = Jbar_ind(n); .. call for lin-indep gagn. indices
% isto kao i main1 ali bez chol - iskljuce se komentari pa radi bez frke
% call: [R,Z,VT,w,t, primal, dual] = main2(A,B,C,P,Q);
  
  
n = size(A,1);   %size of the problem

eq = length(P); 

% initialization

w = ones(n*n+1,1); w(1)=2*n; % belongs to arrow part
t = 10*ones(eq,1);           % belongs to gangster part
a = zeros(n*n+1,1); 
a(1,1) = -1;
b = zeros(eq,1);

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

% 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) ];   %clanak
VT = sparse(VT);

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] = initial(L,VT,w,t,P,Q,eq,n);
     
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   ')
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,ATw,BTt,U1] = right(w,mu,n,VT,Z,invZ,R,L,VTiZ,t,RT,a,P,Q,eq);

test = 1;
if test>0;
  rez = M\D;
  else
[Mpom,p] = chol(M);
wttemp = (Mpom')\ D;
rez = Mpom \ wttemp;   % delta omega
end
dw = rez(1:(n*n+1));
dt = rez((n*n+1)+1 : (n*n+1) + eq);  % delta t

if norm(M*rez-D)> 1e-5; disp('accuracy ...'); end;

[ATw] = AT(dw,n);
[GTt] = GT(dt,P,Q,eq,n);
U2 =  VT'*(ATw + 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-6; 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;
w = w + t1*dw;
t = t + t1*dt;  
[GTt] = GT(t,P,Q,eq,n);
[ATw] = AT(w,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'*( ATw + 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 = w'*a-t'*b;
primal = L(:)'*R(:);

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


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

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


end;   % while .... done


disp(' ');

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

%'at end of main2'
keyboard
%primal: primfiz = diag(RT) - [0; RT(pom,1)];
%dual: [ATw] = AT(w,n); [GTt] = GT(t,P,Q,eq,n);
%      dualfiz = L + VT'*( ATw + GTt )*VT; norm(dualfiz-Z)
