%
% This is a Matlab program for decomposing an n-by-n unitary matrix U as the
% product of N = n(n-1)/2 two-level matrices according to the permutation
% P = [i_1, ..., i_n] of [1, ..., n]. See arXiv:1210.7366.
%
% The Matlab command: A = pud(U,P) will generate the two-level matrices
% A{1}, ..., A{N} such that A{1} ... A{N}U = I_N
% One also needs the auxiliary Matlab program: genDeletePos.m to do the
% decomposition.
%
function finalMatrix = pud(matrixGiven, permutation)
count = 0;
for k = 1:numel(matrixGiven)
count = count + 1;
end
count = sqrt(count);
N = (count * (count-1)) / 2 ; %N is the total number of matrices
% Create a cell array to hold all of the outputs
finalMatrix = cell(N, 1);
j = count;
count = count - 1;
%while loop goes through each column in the matrix
while count > 0
col = zeros(j, 1);
for i = 1:j
col(i,1) = matrixGiven(i, permutation(j-count));
end
%for loop zeros all positions in a column but 1
for idx = 1:count
[col, unitary] = genDeletePos(col, permutation(j-idx+1), permutation(j-idx));
matrixGiven = unitary * matrixGiven;
finalMatrix{N} = unitary;
N = N - 1;
end
count = count - 1;
end