+function PMs = permutation_matrices(integerN)
+ ## Generate all permutation matrices of size ``integerN``.
+ ##
+ ## INPUT:
+ ##
+ ## - ``integerN`` -- The dimension of the resulting matrices.
+ ##
+ ## OUTPUT:
+ ##
+ ## - ``PMs`` -- A cell array of permutation matrices.
+ ##
+
+ if (integerN < 1)
+ PMs = NA;
+ return;
+ end
+
+ ## Append to this as we generate them.
+ PMs = {};
+
+ ## Generate all permutations of [1,2,...,integerN].
+ permutations = perms([1:integerN]);
+
+ for idx = [ 1 : factorial(integerN) ]
+ sigma = permutations(idx,:);
+ ## Create a permutation matrix from the permutation, sigma.
+ P = eye(integerN) (sigma,:);
+ PMs{end+1} = P;
+ end
+
+end