function [R, Y, Out] = ADMM_QAP(L, Vhat, J, opts)
%function [R, Y, Out] = ADMM_QAP(L, Vhat, J, opts,kUU,kUUt,kDD)
% global Valpha rho Vv w Wbar kUU kUUt kDD Vhat W n sqrtt sqrtn onesnn

% L is a symmetric matrix
% Vhat is an orthonormal matrix
% kUU is kron(Vunscaled,Vunscaled) is 0,1 and kDD is kron for the scaling
% so kron(V,V)=KVV/kDD

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 isfield(opts,'adp_pnl'), adp_pnl = opts.adp_pnl;   else adp_pnl = 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;
    Y=sparse(Y);
    Z((abs(Z)<min([nrm_pR,nrm_dR,1e-9])))=0;
    Z=sparse(Z);

% scale matrix L for stability
n2 = size(Vhat,1)-1;
%n=round(sqrt(n2));
%nVhatcols=size(Vhat,2);
%sqrtt=sqrt(2);
%sqrtn=(sqrt(2)*n);
%onesnn=ones(n^2,1);
%Valpha=1/sqrt(2);Vv=(1/(sqrt(2)*n))*ones(n^2,1);
%Vhat0=sparse(Vhat);

%n=sqrt(n2);
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);


%obj0 = sum(L(:).*Y(:)); 
nstall = 0;

%keyboard
tic   % start time
for iter = 1:maxit
    % update R     ....    step 1 of 3
    W = Y + Z/beta;
    %%%% add variables for eigs function eval.%%%%%%%%%%%%%
    %w=W(2:end,1);
    %Wbar=W(2:end,2:end);
    %rho=W(1,1);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    WVhat = Vhat'*(W*Vhat);
  % WVhat0 = Vhat0'*W*Vhat0;
    WVhat = (WVhat'+WVhat)/2;
    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
%         tic;[v,d] =eigs(WVhat,K,'LA',opteigs);toc1=toc;  % use LA .. theory prob otherwise???
if ~issparse(WVhat),
	fprintf('WVhat is NOT sparse before eigs in iter %i \n',iter)
	keyboard
end
opteigs.v0=sum(WVhat)';
         [v,d] =eigs(WVhat,K,'LA',opteigs);  % use LA .. theory prob otherwise???
	 %%%% Use function WVhateval at bottom of file?????
        %[v,d] =eigs(@WVhateval,nVhatcols,K,'LA',opteigs);  % use LA .. theory prob otherwise???
        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,VRV-(L+Z)/beta));
    Y = min(1, max(0,Y));              % skip while testing cvx
    Y(J) = 0; Y(1,1) = 1;              % skip while testing cvx
    Y((abs(Y)<min([nrm_pR,nrm_dR,1e-9])))=0;
    Y=sparse(Y);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % update multiplier      step 3 of 3  update dual variable Z
    pR = Y-VRV;
    dR = Y-Y0; Y0 = Y;       % dR for cgnce and then update
    
    Z = Z + gamma*beta*pR;
    Z=(Z+Z')/2;
    Z((abs(Z)<min([nrm_pR,nrm_dR,1e-9])))=0;
    Z=sparse(Z);
    nrm_pR = norm(pR,'fro');
    feas(iter) = nrm_pR/norm(Y,'fro');
    %feas(iter) = norm(pR,'fro')/norm(Y,'fro');
    
    obj1 = sum(L(:).*Y(:));
    obj(iter) = obj1;
    
    %nrm_pR = norm(pR(:)); 
    %nrm_dR = beta*norm(dR(:));
    nrm_dR = beta*norm(dR,'fro');
    hist_pR(iter) = nrm_pR;
    hist_dR(iter) = nrm_dR;
    %if round(iter/100)*100==iter,
    %              fprintf('iter nrm_pr nrm_dR %i  %g  %g  \n',iter,nrm_pR,nrm_dR)
    %end
    
    if nrm_pR < tol && nrm_dR < tol 
        nstall = nstall + 1; 
    else 
        nstall = 0;
    end
    
    if nstall >=5,   
	   %fprintf('stopping break nstall=%i nrm_pR=%g  >=5 \n',...
	   %nstall,nrm_pR)
          %fprintf('iter nrm_pr nrm_dR %i  %g  %g  \n',iter,nrm_pR,nrm_dR)
	    break;
    end;
    
%    obj0 = obj1;  
    
end
Out.toc=toc;
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


	 %%K in formula is kUU/diag(kDD)
    %%%  W = Y + Z/beta;
    %%%  Vhat equiv to [Valpha zeros;Vv K]
%%	 function y=WVhateval(x)   %   for eigs command
%%% global Valpha rho Vv w Wbar kUU kUUt kDD Vhat W n sqrtt sqrtn onesnn
%%	 xbar=x(2:end);
%%	 x1=x(1);
%%	 outKxbar=kUU*(xbar./kDD);  % seems 'ok'
%%	 x1vKxbar=x1*Vv+outKxbar;
%%	 ValphwVW=(w/sqrtt)'+(Vv'*Wbar);
%%	 row1=(rho*x1)/2+ValphwVW*x1vKxbar ...
%%	            + x1*sum(w)/(2*n);
%%	 y=[row1
%%	     (kUUt*(x1*Valpha*w+Wbar*x1vKxbar))./kDD];
%%	 %   y=out;
%%	 %     y=Vhat'*(W*(Vhat*x));
%%	 %     if norm(y-out)>1e-5,
%%	 %        norm(y-out)
%%	%	      'in WV'
%%	%	      keyboard
%%	%      end
%%         end
end   % of main function ADMM_QAP
