Since the permuted matrices are all symmetric, you really only need to permute the lower (or upper) triangle of the matrix, excluding the diagonal, and then reflect the values. Since your matrix is 4x4, there are 6 values in the lower triangle excluding the diagonal. That results in 6! permutations (720). But there are also 16 variations of the diagonal (4 choices from 2 values [0,2] with replacement; if that's what you're looking for since it was unlcear).
so in total, there are 720 * 16 = 11520 permutations of matrix A.
Some of the lines below could be combined but I chose to keep them separated to increase readability. There may be additional shortcuts that I didn't think of but the code completes in less than 150 milliseconds.
The code requires Jos' permn() function from the file exchange in order to compute the permutations of diagonal elements with replacement.
See the inline comments for detials.
A = [0 2 0 1; 2 0 1 0; 0 1 0 1; 1 0 1 0];
trilMat = tril(reshape(1:numel(A),size(A)),-1);
trilIdx = trilMat(trilMat>0);
trilPerms = perms(trilIdx);
diagPermitted = [0,2];
diagPerms = permn(diagPermitted,numel(diag(trilMat)));
diagMask = logical(eye(size(A)));
nPerms = size(trilPerms,1) * size(diagPerms,1);
A_PERM = nan([size(A),nPerms]);
for i = 1:size(trilPerms,1)
for j = 1:size(diagPerms,1)
c = (i-1)*size(diagPerms,1)+j;
base = zeros(size(A));
base(trilIdx) = A(trilPerms(i,:));
base = base + tril(base,-1).';
base(diagMask) = diagPerms(j,:);
error('Matrix not symmetic. i = %d j = %d',i,j)
A_PERM(:,:,c) = base;
A_PERM is a [4 x 4 x 11520] array containing all 11520 permutations of the 4x4 matrix A.
"I would like to restrict the resultant, permuted matrices to also have rows and columns that sum to 3,3,2,2, more specifically two of the rows (or columns) in the "accepted" permuted matrices must sum to 2 and the other two must sum to 3, the order of the row (column) sums does not matter."
It would be complicated to work that restriction into the permutation rules. Instead, keep all of the code above, produce all permutations, and then eliminate any that do no follow that rule. That just requires adding the the following two lines below at the end of the code above.
goodRows = ismember(squeeze(sort(sum(A_PERM,2))).',[2 2 3 3],'rows');
A_PERM(:,:,~goodRows) = ;
To see the sum of each row,
Extract the unique matrices (see the discussion below for details).
megaMat = squeeze(reshape(A_PERM,1,prod(size(A_PERM,[1,2])),size(A_PERM,3))).';
[~, unqRowIdx] = unique(megaMat,'rows');
A_PERM_UNQ = A_PERM(:,:,unqRowIdx);