function [Q, blk_nz, blk_sizes] = get_Q(basis)

% INPUT:

% Given A basis b = {b_1,...,b_k} of the commutant of a permutation

% group G, i.e., b_i is a matrix such that P*b_{i}*P' = b_{i}

% OUTPUT:

% 1. Q ==> the orthogonal matrix such that it simulatenously block diagonalize

% all the elements in the basis, i.e., Q'*bi*Q is a  block-diagonal matrix.

% 2. blk_nz ==> the non-zero pattern of the block-diagonal matrices

% 3. blk_sizes ==> the size of each block

tic;

fprintf('get_Q ')

% sample a symmetric matrix from the basis

B_sample = sample_basis(basis);

% perform a spectral decomposition of B_sample

[Q,~] = eig(B_sample); % Q is orthogonal s.t. Q'BiQ is block-diagonal

% the default tolerance 

if nargin < 2 

    tol = 1e-10;

end

% block-diagonalize each element B_{i} in the basis 

m = length(basis{1}); 

nr_B = length(basis); 

blk_nz = ~logical(m); 

for i = 1:nr_B 

    basis_hat = Q'*basis{i}*Q;                  % block diagonalize the basis matrices B 

    basis_hat(abs(basis_hat)<tol) = 0;          % remove zero-entires 

    blk_nz = or(blk_nz,basis_hat);              % remember non-zero entires.

end

% find the ordering s.t. Q'BQ is block diagonal for every B in G

r = symrcm(blk_nz);

% order Q and blk_nz correspondingly 

Q = Q(:,r);

blk_nz = blk_nz(r,r);

% for convienience, we sort the blocks based on its size

% get the size of each block 

blk_sizes = get_blk_sizes(blk_nz); 

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

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

% sort the blocks in descending order

[blk_sizes,idx] = sort(blk_sizes,'descend'); 

k = 0; 

r = zeros(m,1); 

for i = 1:length(blk_sizes) 

        % consider the b-th block

        b = idx(i);

        % map the entries b(s):s(t) to k:k+blk_sizes(i)

        r(k+1:k+blk_sizes(i)) = s(b):t(b);

        % update the starting index 

        k = k+blk_sizes(i); 

end 

Q = Q(:,r);

blk_nz = blk_nz(r,r);

% save the output 

G.Q = Q; 

G.blk_nz = blk_nz;

G.blk_sizes = blk_sizes;

fprintf('  %g \n', toc)

end