function [R, Y, Out] = ADMM_QAP(L, Vhat, J, opts)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Alternating direction method of multiplier (ADMM) for the SDP relaxation
% of quadratic assignment problem (QAP)
% min_{R,Y} trace(L*Y)
% subject to: Y(id) = 0 for all id in J\{(1,1)}, and Y(1,1) = 1;
%             0 <= Y <= 1,
%             Vhat*R*Vhat' = Y,
%             R is a psd matrix
%%%%%%%%%%%==================================================%%%%%%%%%%%%%%
%%% Input:
%           L: a symmetric matrix
%           Vhat: an orthonormal matrix
%           J: index set in the gangster constraint
%           opts. (default value)
%               maxit: maximum number of iterations (500)
%               tol: tolerance for stopping (0.1)
%               beta: penalty parameter (50)
%               low_rank: if there is low-rank constraint (FALSE)
%               K: rank that Y is constrained to, valid only if low_rank = TRUE
%               gamma: stepsize for multiplier update (1)
%
%%% Output:
%           R, Y: solutions
%           Out.
%               obj: history of objective values
%               iter: number of iterations
%               pr: primal residuals
%               dr: dual residuals
%               time: total running time
%               Z: dual variable for computing lower bound
%               pos_eig: track the number of positive eigenvalues
%%%%%%%%%%%%================================================%%%%%%%%%%%%%%
%%%
% Copyright: Danilo E. Oliveira, Henry Wolkowicz, Yangyang Xu, 2016
% Please cite to use the code:
% "Danilo E. Oliveira, Henry Wolkowicz, Yangyang Xu. ADMM for the SDP
% relaxation of the QAP. arXiv:1512.05448, 2015."
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% parameter setting
if isfield(opts,'maxit'),   maxit = opts.maxit;       else maxit = 500; end;
if isfield(opts,'tol'),     tol = opts.tol;           else tol = 1e-1;  end;
if isfield(opts,'beta'),    beta = opts.beta;         else beta = 50;   end;
if isfield(opts,'low_rank'),low_rank = opts.low_rank; else low_rank = 0;end;
if isfield(opts,'K'),       K = opts.K;               else K = 1;       end;
if isfield(opts,'gamma'),   gamma = opts.gamma;       else gamma = 1;   end;
if norm(L-L','fro') > 10*eps
    error('The matrix L is not symmetric! so symmetrizing here');
    L=(L+L')/2;
end
%%%%%%%%%%%%%%%%%%opteigs
opteigs.issym=1;
opteigs.real=1;
%%%%%%%%%%%%%%%%%%
R0 = opts.R0;
Y0 = opts.Y0;
Z0 = opts.Z0;
R = R0; Y = Y0; Z = Z0;
nrm_pR=1;
nrm_dR=1;
Y((abs(Y)<min([nrm_pR,nrm_dR,1e-9]))) = 0;
Z((abs(Z)<min([nrm_pR,nrm_dR,1e-9]))) = 0;


% scale matrix L for stability
n2 = size(Vhat,1)-1;
nrmL = norm(L,'fro');
L = L/nrmL*n2;

feas = zeros(maxit,1);
obj = zeros(maxit,1);
pos_eig = zeros(maxit,1);
max_eig = zeros(maxit,K);
hist_pR = zeros(maxit,1);
hist_dR = zeros(maxit,1);

nstall = 0;

time0 = tic;
for iter = 1:maxit
    % update R     ....    step 1 of 3
    W = Y + Z/beta;
    
    % usually Z is not sparse, and thus neither W
    % sparse matrix multiplies dense matrix is expensive
    WVhat = Vhat'*(W*Vhat);
    WVhat = (WVhat'+WVhat)/2; % make sure WVhat symmetric
    if ~low_rank
        [U,S] = eig(full(WVhat));
        s = diag(S);
        id = s>0;
        R = U(:,id)*diag(s(id))*U(:,id)';
        pos_eig(iter) = sum(id);
    else        

        opteigs.v0 = sum(WVhat)';
        [v,d] = eigs(sparse(WVhat),K,'LA',opteigs);        
        d = diag(d);
        id = d>0;
        if sum(id) > 0
            R = v(:,id)*diag(d(id))*v(:,id)';
        else
            R = zeros(size(R));
        end
        max_eig(iter,:) = d';
    end
    R=(R+R')/2;
    
    % update Y    .....      step 2 of 3 with projections >=0 and <=1
    VRV = Vhat*R*Vhat';
    Y = VRV-(L+Z)/beta;
    Y=(Y+Y')/2;
    Y(J) = 0; Y(1,1) = 1;
    
    Y = min(1, max(0,Y));
    Y(J) = 0; Y(1,1) = 1;
    Y((abs(Y)<min([nrm_pR,nrm_dR,1e-9]))) = 0;

    
    % update multiplier      step 3 of 3  update dual variable Z
    pR = Y-VRV;
    dR = Y-Y0; Y0 = Y;       
    
    Z = Z + gamma*beta*pR;
    Z=(Z+Z')/2;
    Z((abs(Z)<min([nrm_pR,nrm_dR,1e-9])))=0;

    nrm_pR = norm(pR,'fro');
    feas(iter) = nrm_pR/norm(Y,'fro');
    
    % evaluate current objective
    obj1 = sum(L(:).*Y(:));
    obj(iter) = obj1;
    
    nrm_dR = beta*norm(dR,'fro');
    hist_pR(iter) = nrm_pR;
    hist_dR(iter) = nrm_dR;
    
    % if both primal and dual residual are small, increase nstall
    if nrm_pR < tol && nrm_dR < tol
        nstall = nstall + 1;
    else
        nstall = 0;
    end
    % only if primal and dual residual keep smaller than tol in a few steps
    if nstall >=5,
        break;
    end
    
end
Out.time=toc(time0);
Out.obj = obj(1:iter)*nrmL/n2;
Out.iter = iter;
Out.feas = feas(1:iter);
Out.pr = hist_pR(1:iter);
Out.dr = hist_dR(1:iter);
Out.Z = Z*nrmL/n2;
if ~low_rank
    Out.pos_eig = pos_eig(1:iter);
else
    Out.max_eig = max_eig(1:iter);
end
end   % of main function ADMM_QAP
