This is an attempt at an analogy with prime numbers. Let's consider only square matrices with positive integer entries. Which of them are 'prime' and how to decompose such a matrix in general?
To illustrate, there is a product of two general $2 \times 2$ matrices:
$$AB=\left[ \begin{matrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix} \right] \left[ \begin{matrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{matrix} \right]=\left[ \begin{matrix} a_{11} b_{11}+a_{12} b_{21} & a_{11} b_{12}+a_{12} b_{22} \\ a_{21} b_{11}+a_{22} b_{21} & a_{21} b_{12}+a_{22} b_{22} \end{matrix} \right]$$
Exchanging $a$ and $b$ we obtain the expression for the other product $BA$.
Now, if we allow zero, negative and/or rational entries we can probably decompose any matrix in an infinite number of ways.
However, if we restrict ourselves:
$$a_{jk},~b_{jk} \in \mathbb{N}$$
The problem becomes well defined.
Is there an algorithm to decompose an arbitrary square positive integer matrix into a product of several positive integer matrices of the same dimensions?
There is a set of matrices which can'be be decomposed, just like the prime numbers (or irreducible polynomials, for example). The most trivial one is (remember, zero entries are not allowed):
$$\left[ \begin{matrix} 1 & 1 \\ 1 & 1 \end{matrix} \right]$$
There are no natural numbers $a_{11},b_{11},a_{12},b_{21}$, such that:
$$a_{11} b_{11}+a_{12} b_{21}=1$$
The same extends to any dimension $d$. Any 'composite' $d \times d$ matrix will have all entries $ \geq d$. Thus, for square matrices we can name several more 'primes':
$$\left[ \begin{matrix} 2 & 1 \\ 1 & 1 \end{matrix} \right],~~~\left[ \begin{matrix} 1 & 2 \\ 1 & 1 \end{matrix} \right],~~~\left[ \begin{matrix} 1 & 1 \\ 2 & 1 \end{matrix} \right],~~~\left[ \begin{matrix} 1 & 1 \\ 1 & 2 \end{matrix} \right],~~~\left[ \begin{matrix} 2 & 2 \\ 1 & 1 \end{matrix} \right],~~~\left[ \begin{matrix} 1 & 1 \\ 2 & 2 \end{matrix} \right], \dots$$
And in general, any matrix which has at least one entry equal to $1$.
It makes sense, that most entries in 'composite' matrices will be large, since we are multiplying and adding natural numbers. For example:
$$\left[ \begin{matrix} 1 & 2 & 4 \\ 3 & 3 & 1 \\ 3 & 4 & 4 \end{matrix} \right] \left[ \begin{matrix} 2 & 5 & 5 \\ 4 & 5 & 5 \\ 5 & 1 & 4 \end{matrix} \right]=\left[ \begin{matrix} 30 & 19 & 31 \\ 23 & 31 & 34 \\ 42 & 39 & 51 \end{matrix} \right]$$
$$\left[ \begin{matrix} 2 & 5 & 5 \\ 4 & 5 & 5 \\ 5 & 1 & 4 \end{matrix} \right] \left[ \begin{matrix} 1 & 2 & 4 \\ 3 & 3 & 1 \\ 3 & 4 & 4 \end{matrix} \right] =\left[ \begin{matrix} 32 & 39 & 33 \\ 34 & 43 & 41 \\ 20 & 29 & 37 \end{matrix} \right]$$
If no decomposition algorithm for this case exists, is it at least possible to recognize a matrix that can't be decomposed according to the above rules?
It's a strange question... Let $A\in M(N)$ s.t. $A=PQ$ where $P,Q\in M(N)$ are random. I calculate "the" Smith normal decomposition of $A$: $A=UDV$ where $U,V\in GL(\mathbb{Z})$ and $D$ is a diagonal in $M(\mathbb{Z})$. During each Maple test, I consider the matrix $UD=[C_1,\cdots,C_n]$, where $(C_i)_i$ are its columns; curiously,
(P) for every $i$, $C_i\geq 0$ or $C_i\geq 0$. Is it true for any such matrices $A$ ?
EDIT. Answer to @ You're In My Eye . I conjectured that property (P) above and, for every $i,j$, $a_{i,j}\geq n$ characterize the decomposable matrices $A\in M(N)$. Unfortunately, the matrix $A=\begin{pmatrix}10&13\\9&5\end{pmatrix}\in M(N)$ satisfies (P) but is indecomposable.
Remark 1. If $A=UV$ is decomposable, then there are many other decompositions: $A=(UP)(P^TV)$ where $P$ is any permutation.
Remark 2. We can consider the permanent function; if $A=UV$, then $per(A)> per(U)per(V)$ and in particular $per(U)<\dfrac{per(A)}{n!}$. If we look for an eventual decomposition of the $A$ above, then we obtain $\det(U)\in\{\pm 67,\pm 1\}$ and $per(U)\leq 83$.