This has a very straightforward answer using the Burnside lemma. With
$n$ rows, $m$ columns and $q$ possible values we simply compute the
cycle index of the cartesian product group ($S_n \times S_m$, consult
Harary and Palmer, Graphical Enumeration, section 4.3) and evaluate
it at $a[p]=q$ as we have $q$ possibilities for an assignment that is
constant on the cycle. The cycle index is easy too -- for two cycles
of length $p_1$ and $p_2$ that originate in a permutation $\alpha$
from $S_n$ and $\beta$ from $S_2$ the contribution is
$a[\mathrm{lcm}(p_1, p_2)]^{\gcd(p_1, p_2)}.$
We get for a $3\times3$ the following colorings of at most $q$ colors:
$$1, 36, 738, 8240, 57675, 289716, 1144836, 3780288,\ldots$$
which points us to OEIS A058001 where
these values are confirmed.
We get for a $4\times 4$ the following colorings of at most $q$ colors:
$$1, 317, 90492, 7880456, 270656150, 4947097821,
\\ 58002778967, 490172624992,\ldots$$
which points us to OEIS A058002 where
again these values are confirmed.
We get for a $5\times 5$ the following colorings of at most $q$ colors:
$$1, 5624, 64796982, 79846389608, 20834113243925, 1979525296377132,
\\ 93242242505023122, 2625154125717590496,\ldots$$
which points us to OEIS A058003 where
here too these values are confirmed.
This was the Maple code.
with(combinat);
pet_cycleind_symm :=
proc(n)
option remember;
if n=0 then return 1; fi;
expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;
pet_cycleind_symmNM :=
proc(n, m)
local indA, indB, res, termA, termB, varA, varB,
lenA, lenB, instA, instB, p, lcmv;
option remember;
if n=1 then
return pet_cycleind_symm(m);
else
indA := pet_cycleind_symm(n);
fi;
if m=1 then
return pet_cycleind_symm(n);
else
indB := pet_cycleind_symm(m);
fi;
res := 0;
for termA in indA do
for termB in indB do
p := 1;
for varA in indets(termA) do
lenA := op(1, varA);
instA := degree(termA, varA);
for varB in indets(termB) do
lenB := op(1, varB);
instB := degree(termB, varB);
lcmv := lcm(lenA, lenB);
p :=
p*a[lcmv]^(instA*instB*lenA*lenB/lcmv);
od;
od;
res := res + lcoeff(termA)*lcoeff(termB)*p;
od;
od;
res;
end;
mat_count :=
proc(n, m, q)
subs([seq(a[p]=q, p=1..n*m)],
pet_cycleind_symmNM(n, m));
end;
Addendum Nov 17 2018. Note that a product of powers of variables
implements the multiset of cycles concept through indets (distinct
elements) and degree (number of occurrences). Here we iterate
over pairs of monomials representing a conjugacy class from $Z(S_n)$
and $Z(S_m)$ and compute $a[\mathrm{lcm}(p_1, p_2)]^{\gcd(p_1, p_2)}$
for pairs of cycles $a_{p_1}$ and $a_{p_2}.$ This makes for the highly
compact algorithm shown above, which will produce e.g. for a three by
four,
$${\frac {{a_{{1}}}^{12}}{144}}+1/24\,{a_{{1}}}^{6}{a_{{2}}}^{3}
+1/18\,{a_{{1}}}^{3}{a_{{3}}}^{3}+1/12\,{a_{{2}}}^{6}
\\+1/6\,{a_{{4}}}^{3}+1/48\,{a_{{1}}}^{4}{a_{{2}}}^{4}
+1/8\,{a_{{2}}}^{5}{a_{{1}}}^{2}+1/6\,a_{{1}}a_{{2}}a_{{3}}a_{{6}}
\\+1/8\,{a_{{3}}}^{4}+1/12\,{a_{{3}}}^{2}a_{{6}}
+1/24\,{a_{{6}}}^{2}+1/12\,a_{{12}}.$$
Let me just observe that while rigorous notation is a must for any
serious mathematician it can sometimes block the path of the beginning
reader. What is happening here is very simple. We need the cycle index
of the permutations from $S_2\times S_3$ acting on the slots of the
matrix. While the present case can be computed by inspection I will
explain the general method.
We require the cycle index $Z(M_{2,3})$ of the pairs $(\sigma,\tau)$
of row and column permutations acting simultaneously on the slots of
the matrix ($6$ of them, and with $12$ permutations total). Start
from the two cycle indices for the row and column permutations
$$Z(S_2) = \frac{1}{2} a_1^2 + \frac{1}{2} a_2$$
and
$$Z(S_3) =\frac{1}{6} b_1^3 + \frac{1}{2} b_1 b_2 + \frac{1}{3} b_3.$$
The method here is as follows. We draw a diagram of the cycle types of
$\sigma$ and $\tau,$ perhaps one beneath the other. Now for the cycle
index of the cartesian product of $S_2$ and $S_3$ we must factor the
combined action of the two permutations on the slots into cycles. We
represent row-column pairs i.e. slots by marking say the row and the
column and connecting them with a edge, not to be confused with the
directed edges of the two cycles. This edge travels in parallel along
the two cycles it is on and returns to its initial position after
$\mathrm{lcm}(l_1, l_2)$ steps, where $l_{1,2}$ are the lengths of the
cycles. As the pair of cycles contributes $l_1\times l_2$ slot
identifiers we get for the contribition to the combined cycle index
the term
$$a_{\mathrm{lcm}(l_1, l_2)}^{l_1 l_2/\mathrm{lcm}(l_1, l_2)}.$$
We now do the computation. There are thus six possible combinations of
cycle types that combine to form $Z(M_{2,3})$. We process these in
turn, leaving the most difficult part for last.
First, combining $a_1^2$ and $b_1^3.$ This fixes all six slots for a
contribution of
$$\frac{1}{12} a_1^6.$$
Second, combining $a_1^2$ and $b_3.$ This partitions the pairs into
three-cycles for a contribution of
$$\frac{1}{6} a_3^2.$$
Third, combining $a_2$ and $b_1^3.$ Here we have everything on
two-cycles to get
$$\frac{1}{12} a_2^3.$$
Fourth, combining $a_2$ and $b_3$. This will produce a six-cycle to get
$$\frac{1}{6} a_6.$$
Now for the tricky part. Fifth, combining $a_1^2$ and $b_1 b_2.$ There
are two pairs that are fixed by these, and two pairs on two-cycles
and we obtain
$$\frac{1}{4} a_1^2 a_2^2.$$
Finally, sixth, combining $a_2$ and $b_1 b_2.$ We have everything on
two-cycles and obtain
$$\frac{1}{4} a_2^3.$$
Adding everything we now have the cycle index
$$Z(M_{2,3})
= \frac{1}{12} a_1^6
+ \frac{1}{3} a_2^3
+ \frac{1}{6} a_3^2
+ \frac{1}{4} a_1^2 a_2^2
+ \frac{1}{6} a_6.$$
In order to apply Burnside and ask about colorings with at most $N$
colors we have that the assignment of the colors must be constant on
each cycle and we obtain
$$\frac{1}{12} N^6
+ \frac{1}{3} N^3
+ \frac{1}{6} N^2
+ \frac{1}{4} N^4
+ \frac{1}{6} N.$$
This yields the sequence
$$M_n = 1, 13, 92, 430, 1505, 4291, 10528, 23052, \ldots$$
which is OEIS A027670. In particular we
find there are $92$ colorings using at most three distinct colors. We
could apply PET at this point since we have the cycle index. E.g. for
three colors we obtain
$$1/12\, \left( R+G+B \right) ^{6}
\\ +1/4\, \left( R+G+B \right) ^{2} \left( {B}^{2}+{G}^{2}+{R
}^{2} \right) ^{2}+1/6\, \left( {B}^{3}+{G}^{3}+{R}^{3} \right) ^{2}
\\ +1/3\, \left( {B}^{2}+
{G}^{2}+{R}^{2} \right) ^{3}+1/6\,{B}^{6}
+1/6\,{G}^{6}+1/6\,{R}^{6}$$
which expands to
$${B}^{6}+{B}^{5}G+{B}^{5}R+3\,{B}^{4}{G}^{2}+3\,{B}^{4}GR
\\ +3\,{B}^{4}{R}^{2}+3\,{B}^{3}{G}^{3}+6\,{B}^{3}{G}^{2}R
\\ +6\,{B}^{3}G{R}^{2}+3\,{B}^{3}{R}^{3}+3\,{B}^{2}{G}^{4}
\\ +6\,{B}^{2}{G}^{3}R+11\,{B}^{2}{G}^{2}{R}^{2}+6\,{B}^{2}G{R}^{3}
\\ +3\,{B}^{2}{R}^{4}+B{G}^{5}+3\,B{G}^{4}R+6\,B{G}^{3}{R}^{2}
\\ +6\,B{G}^{2}{R}^{3}+3\,BG{R}^{4}+B{R}^{5}+{G}^{6}+{G}^{5}R
\\ +3\,{G}^{4}{R}^{2}+3\,{G}^{3}{R}^{3}+3\,{G}^{2}{R}^{4}
\\+G{R}^{5}+{R}^{6}.$$
The reader might want to verify some of these with pen and paper.
We may also apply inclusion-exclusion to obtain the count of colorings
of the matrix with exactly $N$ colors. The nodes of the poset
correspond to sets $P\subseteq [N]$ of colors which include colorings
that use some subset of these colors, with the top node representing
at most $N$ colors, which is the only $P$ that includes colorings with
exactly $N$ colors. Colorings with exactly $p \lt N$ colors where
$p\ge 1$ are included at all nodes that are supersets of the $p$
colors. We thus obtain for the total weight
$$\sum_{q=p}^N {N-p\choose q-p} (-1)^{N-q}
= (-1)^{N-p} \sum_{q=0}^{N-p} {N-p\choose q} (-1)^q = 0$$
since $N-p\ge 1.$ Colorings with less than $N$ colors have total
weight of zero in the poset. We thus obtain
$$\sum_{q=1}^N {N\choose q} (-1)^{N-q} M_q.$$
We get a finite sequence since with six slots it is not possible to
have a coloring with more than six distinct colors:
$$1, 11, 56, 136, 150, 60$$
The last term represents colorings with exactly six colors. This means
all slots in the matrix are distinct. Therefore all orbits have the
same size, the number of permutations in the group, which is twelve,
and indeed $6!/12 = 60.$ The term for two colors indicates that the
two monochrome colorings have been excluded.
This MSE link
has the Maple code for the general case.
Addendum. Here is what we mean when we say in the introduction
that the cycle index can be computed by inspection. This refers to the
isomorphism between $M_{2,3} = S_2\times S_3$ and $D_6$, the dihedral
group (reflections and rotations of regular polygons) acting on six
slots in this case. The cycle indices $Z(D_p)$ are tabulated and have
simple closed forms, consult
e.g. Wikipedia. Label the
vertices of a hexagon in clockwise order with the labels $(0,0),
(1,1), (0,2), (1,0), (0,1)$ and $(1,2).$ Then it is not difficult to
see that the rotations of the hexagon are in a bijection with the
pairs of cycles from $C_2 \times C_3$ embedded in $M_{2,3}.$ E.g. the
rotation that takes the top vertex to its clockwise neighbor
corresponds to the two cycles $(0,1)$ and $(0,1,2)$. The reflections
in an axis passing through opposite vertices preserve parity
(permutation $(0)(1)$ from $S_2$) and fix one element from $0,1,2$
while permuting the other two in two-cycles. The reflections in an
axis passing through opposite edges flip parity (permutation $(0,1)$
from $S_2$) and fix one of three elements from $0,1,2$ while
exchanging the other two. In this way we have bjectively accounted for
all permutations and the proof of the isomorphism is complete.
Best Answer
Basically this is about bipartite graphs: graphs where the vertex set is split into two parts that I will call type $X$ vertices and type $Y$ vertices, and where all edges have one end of type $X$ and another of type $Y$. These graphs represent a relation between elements of type $X$ and $Y$, the relation holding between $x$ and $y$ if there is an edge between $x$ and $y$. In your example type $X$ elements would be row indices and type $Y$ elements would be column indices (so a number that is both a row index and a column index gives rise both to a type $X$ vertex and an unrelated type $Y$ vertex), and the edges of the graph are between pairs $(x,y)$ where the matrix has a nonzero entry in position $(x,y)$. Permuting rows respectively columns now just means renumbering the type $X$ respectively type $Y$ vertices, so you can think of these sets as being unlabelled, which helps you abstract away from the order in which rows and columns are initially given.
Now what you want amounts to partitioning the whole graph into chunks such that there are no edges between different chunks (and moreover you want each chunk to have specified numbers of type $X$ and type $Y$ vertices). The finest partitioning with this property is obtained by computing the connected components of you bipartite graph. There is a straightforward and efficient algorithm to compute connected components. Even for sparse graphs, it is quite possible that there is only one connected component, in which case you are out of luck. Supposing however that there are many connected components, then you next task is to group together connected components into chunks so that the number of type $X$ and type $Y$ vertices matches your specification for each chunk. This is a kind of bin packing problem, but with an extra twist that sizes of items to be packed are pairs of numbers rather than single numbers (namely the numbers of type $X$ and type $Y$ vertices), which sizes are added together as $2$-component vectors.
This does not look as an easy problem, neither with respect to it being likely that there is a solution at all, nor with respect to allowing efficient algorithms for solving the general solvable case. I would personally be somewhat demotivated to even spend much effort on this problem.