function output = run_admm(output,C_tilde,Bhat_blk,U,V,opts)

% this function runs the alternating direction method of multipliers(ADMM)

%

% timer

fprintf('admm starts here: \n') 

admm_t = tic; 

m = length(C_tilde)-1; 

maxit = opts.maxit; 

beta = opts.beta; 

gamma = opts.gamma; 

scalar = opts.scalar; 

tol = opts.tol; 

blk_sizes = output.blk_sizes; 

blk_nz = output.blk_nz; 

B_mat_mean = output.B_mat_mean; 

B_matrix = output.B_mat; 

orbit_size = size(B_matrix,2); 

nr_blocks = length(blk_sizes); 

bt = cumsum(blk_sizes);          % the end index of each blocks 

bs = [1; bt(1:end-1)+1];           % the begin index of each blocks 

% Initialization of the alternating direction method of multipliers(ADMM)

% preallocate variables

fvals = zeros(maxit,1); 

p_res = zeros(maxit,1); 

d_res = zeros(maxit,1); 

Z = cell(nr_blocks,1); 

R = cell(nr_blocks,1); 

VRV = cell(nr_blocks,1); 

U_blk = cell(nr_blocks,1);        % columns of U associated to the k-th block of Y 

C_blk = cell(nr_blocks,1);        % k-th block of U'A0U 

Bhat_x = cell(nr_blocks,1);       % Bhat_x contains k-th block of sum \hat{B}^{*}_{i}x_{i} 

block_z = cell(nr_blocks,1);      % the zero pattern of k-th block 

% initial value of x

% x = rand(orbit_size,1);

x = zeros(orbit_size,1);          % initialize with zeros has superior convergence 

% x(1) = 1;

for i = 1:nr_blocks 

    % initial value of Z

    Z{i} = zeros(blk_sizes(i)); 

    Z{i}(:) = Bhat_blk{i}*x; 

    % initial value of R

    R{i} = zeros(size(V{i},2)); 

    % initialize variables

    idx = bs(i):bt(i);             % idx for the i-th block of Y 

    U_blk{i} = U(:,idx); 

    C_blk{i} = C_tilde(idx,idx); 

    Bhat_x{i} = zeros(blk_sizes(i)); 

    block_z{i} = ~blk_nz(idx,idx); 

end 

%%

VRV_CZ = cell(nr_blocks,1); 

%% execute ADMM

for k = 1:maxit 

    % ==========solve the R-subproblem=========

    Rold=R; 

    Rdiff = 0; 

    for i = 1:nr_blocks 

        VBZV = V{i}'*(Bhat_x{i}+(1/beta)*Z{i})*V{i}; 

        VBZV = (VBZV+VBZV')/2; 

        %       R{i} = proj_psd(VBZV);

        if size(V{i},1) > 1 

            R{i} = proj_psd(VBZV); 

        else 

            if VBZV < 0 

                R{i} = 0;

            else 

                R{i} = VBZV; 

            end 

        end 

        %         if norm(R2{i}-R{i}) > 1e-3

        %             keyboard

        %         end

        Rdiff = Rdiff + norm(Rold{i}-R{i},'fro'); 

    end 

    % ==========solve the x-subproblem=========

    xold=x; 

    %%%%% (*) is slower than (**)???????

    % (*)

    %     M = zeros(m+1);

    %     for i = 1:nr_blocks

    %         VRV{i} = V{i}*R{i}*V{i}';

    %         M = M + U_blk{i}*(VRV{i} -(C_blk{i}+Z{i})/beta)*U_blk{i}';

    %     end

    %     M = (M+M')/2;

    % (*)

    for i = 1:nr_blocks 

        VRV{i} = V{i}*R{i}*V{i}'; 

        VRV_CZ{i} = VRV{i}-(C_blk{i}+Z{i})/beta; 

    end 

    M1 = U*blkdiag(VRV_CZ{:})*U'; 

    M = (M1+M1')/2; 

    x = B_mat_mean*M(:); 

    x = max(min(x,1),0); 

    x(1) = 1; 

    xdiff=norm(xold-x); 

    %==========dual variable Z-update=========

    Z_old = Z; 

    for i = 1:nr_blocks 

        Bhat_x{i}(:) = Bhat_blk{i}*x;                   % precompute 

        Z{i} = Z{i} + gamma*beta*(Bhat_x{i}-VRV{i});    % update the dual variable 

        Z{i} = (Z{i}+Z{i}')/2; 

        Z{i}(block_z{i}) = 0;                           % set entries to zero based on zero-pattern 

    end 

    %==========update obj&res values==========

    %%%%% no need for f_vals(k)????????????????????

    for i = 1:nr_blocks 

        fvals(k) = fvals(k) + trace(C_blk{i}*Bhat_x{i})/scalar; 

        p_res(k) = p_res(k) + norm(Bhat_x{i}-VRV{i},'fro'); 

        d_res(k) = d_res(k) + norm(Z_old{i}-Z{i},'fro'); 

    end 

    %==========print the information==========

    print_admm(k, fvals(k), p_res(k),d_res(k),Rdiff,xdiff, 0); 

    % ==========stopping criterions==========

    % STOP when the step-size is less than tol

    if (xdiff < tol)&&(Rdiff < tol) 

        output.exitflag = 'converged'; 

        print_admm(k, fvals(k), p_res(k),d_res(k),Rdiff,xdiff, 1); 

        break 

    end

end 

% timer and the optimal solution/value

admm_t = toc(admm_t); 

Y = reshape(B_matrix*x,m+1,m+1); 

fval = fvals(k); 

% save the output

output = save_output(output,x,R,Z,Y,admm_t,k,fvals,p_res,d_res); 

% ==========print the outputs==========

print_output(output) 

end