I'm expanding my comment to an answer.
You can always do it the brute force way.
Find some permutation matrix $P_\pi$ such that none of the entries $A[i,\pi(i)]$ are $0$.
Subtract $(\min_i A[i,\pi(i)])\, P_\pi$ from $A$.
Repeat steps (1) and (2) until you're left with the $0$ matrix.
As the OP comments, to prove that it terminates, you can note that each time you are left with a multiple of a double stochastic matrix, and step (2) always adds at least one more zero entry.
We still need to show that step (1) is always possible. This follows from Hall's marriage theorem. This theorem states that such a permutation $\pi$ exists unless there is a set $R$ of rows and a set $C$ of columns such that $|R|+|C| < n$ and $A[i,j] = 0$ unless $i \in R$ or $j \in C$. In other words, all the non-zero entries are either in a row in $C$ or a column in $R$.
We will show by contradiction that such a matrix cannot be doubly stochastic. Look at the the sum of all elements $A[i,j]$ with $i \in R$ and $j \in \bar{C}$. Since each of the column sums is $1$,
$$\sum_{i\in R, j\in \bar{C}} A[i,j] = \sum_{j\in \bar{C}} A[i,j] = |\bar{C}|.$$ Also, each of the row sums is $1$, so
$$\sum_{i\in R, j\in \bar{C}} A[i,j] \leq \sum_{i\in R} A[i,j] = |R|.$$
But $|R| < |\bar{C}|$, so
$$\sum_{i\in R, j\in \bar{C}} A[i,j]\leq |R| < |\bar{C}|= \sum_{i\in R, j\in \bar{C}} A[i,j],$$
a contradiction.
Note that this gives a proof of Birkhoff's theorem. To turn it into a real algorithm, you still need an algorithm for step (1) that finds a permutation in Hall's marriage theorem. There are lots of algorithms to do this; I don't know which is the simplest one. (It's a special case of the bipartite maximum matching problem, which in turn is a special case of the assignment problem, so any bipartite maximum matching algorithm, or any algorithm solving the assignment problem, will work.).
Best Answer
The matrix $A$ is almost circulant. If the third and fourth rows were swapped, then it would be. So let's consider an alternative matrix: $B=PA$ where $P$ is the permutation matrix $$P=\begin{bmatrix} 1&0&0&0\\ 0&1&0&0\\ 0&0&0&1\\ 0&0&1&0 \end{bmatrix}.$$ Then $B$ is circulant.
Since it is circulant, it can be diagonalized with the Discrete Fourier Transform (DFT). The DFT values of $(0.4, 0.3, 0.2, 0.1)$ are $$d = \begin{bmatrix} 1.0000 + 0.0000i \\ 0.2000 - 0.2000i \\ 0.2000 + 0.0000i \\ 0.2000 + 0.2000i \end{bmatrix}$$ Therefore, $B=F^\ast D F$, where $D=\text{diag}(d)$, where $F$ is the DFT matrix.
Note that $B^2=B\cdot B = F^\ast D F \cdot F^\ast D F$. But, since $F$ is unitary, $F F^\ast=I$. Therefore, $B^2=F^\ast D^2 F$. Hopefully you can now see that $B^n = F^\ast D^n F$. (This is a common trick.)
Therefore, $$\lim_{N\rightarrow\infty} B^N = \lim_{N\rightarrow\infty} F^\ast D^N F = F^\ast \text{diag}(1,0,0,0)F.$$
Taking the inverse DFT of $(1,0,0,0)$ yields $(0.25,0.25,0.25,0.25)$. And making a circulant matrix from this vector yields your answer.
(Another way to think about this is to ask, what circulant matrix would diagonalize with eigenvalues $(1,0,0,0)$; and that matrix is the specified matrix.)
So how does this relate to our original $A$?
I didn't figure it out. Luckily, there's another person in this thread that appears to have! https://math.stackexchange.com/a/3114768/24205