function [Y, fval, output] = qap_admm_blkdiag(qap_A,qap_B,G,opts)

% load the group information

Q = G.Q; 

B = G.B; 

Bhat = G.Bhat; 

blk_nz = G.blk_nz; 

blk_sizes = G.blk_sizes; 

% the objective kron(B,A)

BXA = kron(qap_B,qap_A); 

%% load the setting

[opts,output] = load_setting(opts,BXA); 

b_tol = 1e-10; % the tolerance for the block-diagonalization 

%% load the instance parameter

n = length(qap_A); 

% m = n^2;

% print information

    fprintf('--------------------------------------------\n') 

    fprintf('Instance information: number of facilitys/locations %d \n', n) 

%% obtain the data and datamatrices for the sdp below

% min{ <A0,Y> | <Ai,Y>==0 for i = 1,2 }

scalar = opts.scalar; 

% cost matrix Q in the objective

A0 = blkdiag(0,BXA)*scalar;    

% A1 is the gangster constraint <A1,Y> = 0

A1 = kron(eye(n),ones(n)-eye(n)) + kron(ones(n)-eye(n),eye(n)); 

A1 = blkdiag(0,A1); 

% T_hat is in the null-space of Y,i.e., T_hatY = 0

% Thus T_hat'T_hat is an exposing vector

T = [kron(eye(n),ones(1,n)); kron(ones(1,n),eye(n))]; 

T_hat = [-ones(2*n,1) T]; 

% ev = [2*n -2*ones(1,n^2); -2*ones(n^2,1) kron(eye(n),ones(n))+kron(ones(n),eye(n))];

% if norm(ev - T_hat'*T_hat) > 0

%     keyboard

% end

%% group symmetry - preprocessing

% get the exposing vector of the symmetry reduced problem

Wbar = get_blkdiag(T_hat'*T_hat,Q,blk_nz,b_tol); 

C_tilde = get_blkdiag(A0,Q,blk_nz,b_tol); 

output.Wbar = Wbar; 

% remove the base B_{i} if the associated entry is 0 

idx = false(length(B),1); 

for i = 1:length(B) 

   if norm(double(B{i}(logical(A1)))) == 0 

       idx(i) = true; 

   end     

end 

B = B(idx);         % only keep base B_{i} for non-zeros 

Bhat = Bhat(idx); % only keep base Bhat_{i} for non-zeros 

output.B_hat = Bhat; 

%% ADMM

% preprocess the data-matrices for admm

[B_mat, B_mat_mean, Bhat_blk, V] = get_admm_matrix(B, Bhat, Wbar, blk_sizes); 

% check if the dimension of the problem is correct

r = sum(cellfun(@(x) size(x,2), V));        % the dim. of facially reduced program 

if r ~= (n-1)^2+1                               % the strictly feasible problem has rank m-n+2 

    error('the dimension is NOT correct.')

end

% save the data-matrices for admm

output.blk_sizes = blk_sizes; 

output.blk_nz = blk_nz; 

output.B_mat_mean = B_mat_mean; 

output.B_mat = B_mat; 

% run admm to solve the problem

output = run_admm(output,C_tilde,Bhat_blk,Q,V,opts); 

% obtain the optimal solution&value

Y = output.Y; 

fval = output.fval; 

% compute the upper&lower bound

% output.lb = getlb(C_tilde,blkdiag(output.Z{:}),blkdiag(V{:}),U)/scalar;

output.lb = -1; 

% [output.ub, output.ub_x] = getub(I,Q,s,t,Y(m+3:m+2:end));

output.ub =-1; 

output.ub_x  =-1; 

fprintf('lower bound      : %.3f (*)\n', output.lb) 

fprintf('upper bound      : %.3f (*)\n', output.ub) 

end