function [Z,u,psi,iter]= specialized_ipa(G,rho,eta,Z,u,tol,max_iter) 
% == Nov. 2008 ==
% By Wing-Kin Ma, CUHK
% A dedicated primal-dual interior-point algorithm for
% a special diagonally constrained SDP
%     max tr(GZ) s.t. Z>=0
%                     \sum_{i=1}^{n-1} Z_{ii} + \rho Z_{nn}= 1
%                     Z_{nn} <= Z_{ii} <= \eta Z_{nn}, i=1,...,n-1
% or its dual
%     min v s.t. Y>=0, t>=0
%                Y=D([ 1- t(1:n-1) +t(n:2n-2); rho+ \sum_{i=1}^{n-1} 
%                      t_i - \eta t_{i+n-1} ])- G
%
% This is a M-file prototype. The syntax:
% 
% [Z,u,psi,iter]= specialized_ipa(G,rho,eta,Z,u,tol,max_iter)
%
% For the inputs:
%   - Z is the initial primal point (must be strictly feasible)
%   - u= [ v', t' ]' is the initial dual point (also strictly feasible)
%   - tol is the solution accuracy (say, 1e-3 or 1e-4 is reasonable)
%   - max_iter is the maximum number of iteration, beyond which the
%     algorithm will terminate even if the sol. accuracy is not achieved
%
% For the outputs:
%   - Z is the primal solution
%   - u= [ v', t' ]' is the dual solution
%   - psi is the optimal value
%   - iter is the no. of iterations run
%
% Note: this specialized IPA MATLAB code is for demonstration purposes.
%   The running speed is further improved if the code is written in C.
%   However C codes are generally not friendly to read as MATLAB codes. 
%   In fact, the running speed reported in the paper is based on a C 
%   code that was written by an operation-optimized approach (instead of
%   using friendly-to-read operations).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

SIGMA= 0.5; %centrality parameter
BETA= 0.8; %backoff parameter in line search
GAMMA= 0.95; %backoff parameter to stay away from boundary
n= size(G,1);
k= 1; % one equality
m= 2*(n-1); % 2(n-1) inequalities
kpm= m+k;
iter= 0;    %iteration count

%phi= q'*u; %initial dual objective
phi= u(1);
psi= G(:)'*Z(:); %initial primal objective
mu= (phi-psi)*SIGMA/(n+m); %initial central trajectory
%Y= diag(P'*u)- G; %initial Y= Diag(A'*v)+ Diag(B'*t)- G
Y= diag( [ ones(n-1,1)*u(1)- u(2:n)+ u(n+1:(2*n-1)); rho*u(1) + sum( u(2:n)- eta*u(n+1:(2*n-1)) ) ] )- G;

while phi-psi > tol & iter <= max_iter,
    iter= iter+1;
    %computations of search directions
    Yi= inv(Y); %inverse of Y explicitly needed
    Yi= (Yi+Yi')/2; %symmetrize (optional?)
    W= (Yi.*Z); %W= Y^{-1} o Z
    %F= P*W*P';
    F= zeros(2*n-1,2*n-1);
    h= sum(W(n,1:(n-1))); %h= ones(n-1,1)'*w_{12}
    iota= W(1:n-1,1:n-1)*ones(n-1,1); %iota= W_{11}*ones(n-1,1)
    F(1,1)= sum(iota)+ 2*rho*sum(h)+ rho^2*W(n,n);
    F(2:n,1)= -iota- rho*W(n,1:(n-1))'+ (h+ rho*W(n,n))*ones(n-1,1);
    F((n+1):(2*n-1),1)= iota+ rho*W(n,1:(n-1))'- eta*(h+ rho*W(n,n))*ones(n-1,1);
    F(1,2:(2*n-1))= F(2:(2*n-1),1)';
    F(2:n,2:n)= W(1:n-1,1:n-1)- W(n,1:(n-1))'*ones(n-1,1)'- ones(n-1,1)*W(n,1:(n-1)) + W(n,n)*ones(n-1,n-1);
    F(n+1:(2*n-1),2:n)= -W(1:n-1,1:n-1)+ W(n,1:(n-1))'*ones(n-1,1)'+ eta*ones(n-1,1)*W(n,1:(n-1))- eta*W(n,n)*ones(n-1,n-1);
    F(2:n,n+1:(2*n-1))= F(n+1:(2*n-1),2:n)';
    F(n+1:(2*n-1),n+1:(2*n-1))= W(1:n-1,1:n-1)- eta*W(n,1:(n-1))'*ones(n-1,1)'- eta*ones(n-1,1)*W(n,1:(n-1))+ eta^2*W(n,n)*ones(n-1,n-1);
    %r= (q(k+1:kpm)- P(k+1:kpm,:)*diag(Z))./u(k+1:kpm);
    diagZ= diag(Z); r= [ diagZ(1:n-1)- diagZ(n); eta*diagZ(n)- diagZ(1:n-1) ]./u(k+1:kpm);
    F(k+1:kpm,k+1:kpm)= F(k+1:kpm,k+1:kpm)+ diag(r);
    %g= mu*( P*diag(Yi)+ [ zeros(k,1); 1./u(k+1:kpm) ] )- q;
    diagYi= diag(Yi);
    g= mu*[ sum(diagYi(1:n-1))+ rho*diagYi(n); -diagYi(1:n-1)+ diagYi(n); diagYi(1:n-1)- eta*diagYi(n) ]+ ...
        [ -1; mu./u(k+1:kpm) ];
    du= F\g;
    %dY= diag(P'*du);
    dY= diag( [ ones(n-1,1)*du(1)- du(2:n)+ du(n+1:(2*n-1)); rho*du(1) + sum( du(2:n)- eta*du(n+1:(2*n-1)) ) ] );
    dZ= mu*Yi- Z- Yi*dY*Z;
    dZ= (dZ+dZ')/2; %symmetrize  
    %line search on primal
    alphap= 1;
    Zt= Z+alphap*dZ;
    [dummy,posdef]= chol(Zt); %see if pos def
    diagZt= diag(Zt);
    while posdef > 0 | any(diagZt(1:n-1)- ones(n-1,1)*diagZt(n) <=0) | any(diagZt(1:n-1)- eta*ones(n-1,1)*diagZt(n) >=0),
        alphap= alphap*BETA;
        Zt= Z+alphap*dZ;
        [dummy,posdef]= chol(Zt);
        diagZt= diag(Zt);
    end;
    if alphap < 1, alphap= alphap*GAMMA; end; %stay away from boundary
    %line search on dual
    alphad= 1;
    [dummy,posdef]= chol(Y+alphad*dY);
    while posdef> 0 | any(u(k+1:kpm)+alphad*du(k+1:kpm)<= 0),
        alphad= alphad*BETA;
        [dummy,posdef]= chol(Y+alphad*dY);
    end;
    if alphad < 1, alphad= alphad*GAMMA; end; %stay away from boundary
    
    %update
    Z= Z+alphap*dZ; Y= Y+alphad*dY;
    u= u+alphad*du;
    %phi= q'*u; %current dual objective
    phi= u(1);
    psi= G(:)'*Z(:); %current primal objective
    mu= (phi-psi)*SIGMA/(n+m); %update central trajectory
    if alphap + alphad > 1.8, mu= mu*SIGMA; end; % speed up for long steps
    %fprintf('%d  Duality gap= %4.4f, tp= %4.4f, td= %4.4f\n',iter,phi-psi,alphap,alphad);
end;
    

