<s>[INST] How could we define the factorial of a matrix? Suppose I have a square matrix $\mathsf{A}$ with $\det \mathsf{A}\neq 0$.
How could we define the following operation? $$\mathsf{A}!$$
Maybe we could make some simple example, admitted it makes any sense, with 
$$\mathsf{A} =
\left(\begin{matrix}
1 & 3 \\ 
2 & 1 
\end{matrix}
\right)
$$
 [/INST] I use the well known (and simple) definitions 
$$n!=\Gamma (n+1)$$
and
$$\Gamma (A+1)=\int_0^{\infty } \exp (-t) \exp (A \log (t)) \, dt$$
Now, if A is a (square) matrix, all we need is to define the exponential function for a matrix.
This can always be done (in principle) via the power series which only requires to calculate powers of the matrix and adding matrices.
We skip here a possible diagonalization procedure (which was shown by others before) and use the function MatrixExp[] of Mathematica. For the Matrix given in the OP
$$A=\left(
\begin{array}{cc}
 1 & 3 \\
 2 & 1 \\
\end{array}
\right);$$
we have
Ax = MatrixExp[A Log[t]]

which gives
$$\left(
\begin{array}{cc}
 \frac{t^{1-\sqrt{6}}}{2}+\frac{t^{1+\sqrt{6}}}{2} & \frac{1}{2} \sqrt{\frac{3}{2}} t^{1+\sqrt{6}}-\frac{1}{2} \sqrt{\frac{3}{2}} t^{1-\sqrt{6}} \\
 \frac{t^{1+\sqrt{6}}}{\sqrt{6}}-\frac{t^{1-\sqrt{6}}}{\sqrt{6}} & \frac{t^{1-\sqrt{6}}}{2}+\frac{t^{1+\sqrt{6}}}{2} \\
\end{array}
\right)$$
We observe that in some cases the exponent of t is less than -1 ($1-\sqrt{6}=-1.44949$). This leads to a divergent integral which will then be understood as the analytic continuation.
This is accomplished simply by replacing each $t^q$ by $\Gamma (q+1)$.
fA = Ax /. t^q_ -> Gamma[q + 1]

giving
$$\left(
\begin{array}{cc}
 \frac{\Gamma \left(2-\sqrt{6}\right)}{2}+\frac{\Gamma \left(2+\sqrt{6}\right)}{2} & \frac{1}{2} \sqrt{\frac{3}{2}} \Gamma \left(2+\sqrt{6}\right)-\frac{1}{2} \sqrt{\frac{3}{2}} \Gamma \left(2-\sqrt{6}\right) \\
 \frac{\Gamma \left(2+\sqrt{6}\right)}{\sqrt{6}}-\frac{\Gamma \left(2-\sqrt{6}\right)}{\sqrt{6}} & \frac{\Gamma \left(2-\sqrt{6}\right)}{2}+\frac{\Gamma \left(2+\sqrt{6}\right)}{2} \\
\end{array}
\right)$$
Numerically this is
$$\left(
\begin{array}{cc}
 3.62744 & 8.84231 \\
 5.89488 & 3.62744 \\
\end{array}
\right)$$
This is in agreement with the result of hans and travis.
Discussion
(1) Let me point out that the definitions presented here do not need any specific property of the matrix. For example, is does not matter whether the eigenvalues and eigenvectors are singular or not, the matrix might well be deficient.
(2) I have used Mathematica here just to facilitate things. In the end we all use some tools at a certain stage to come to terms. The main ideas are idependent.
(3) The procedure desribed here obviously generalizes to other more or less complicated analytic functions.
As a more exotic example let us take the harmonic number H(A) of a matrix A. This function can be defined using the integral representation (see e.g. Relation between binomial coefficients and harmonic numbers)
$$H(\text{A})\text{=}\int_0^1 \frac{1-(1-x)^A}{x} \, dx$$
This definition also needs only the exponential function of the matrix.
The result for our matrix A is (after some analytic continuation)
$$\left(
\begin{array}{cc}
 \frac{1}{2} \left(H_{-\sqrt{6}}+H_{\sqrt{6}}+\frac{1}{1-\sqrt{6}}+\frac{1}{1+\sqrt{6}}\right) & \frac{1}{2} \sqrt{\frac{3}{2}} \left(-H_{-\sqrt{6}}+H_{\sqrt{6}}-\frac{1}{1-\sqrt{6}}+\frac{1}{1+\sqrt{6}}\right) \\
 \frac{-H_{-\sqrt{6}}+H_{\sqrt{6}}-\frac{1}{1-\sqrt{6}}+\frac{1}{1+\sqrt{6}}}{\sqrt{6}} & \frac{1}{2} \left(H_{-\sqrt{6}}+H_{\sqrt{6}}+\frac{1}{1-\sqrt{6}}+\frac{1}{1+\sqrt{6}}\right) \\
\end{array}
\right)$$
Numerically,
$$\left(
\begin{array}{cc}
 1.51079 & 0.542134 \\
 0.361423 & 1.51079 \\
\end{array}
\right)$$
</s>