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
