(This is my second attempt at a solution. My initial attempt can be found below.)
Your optimization problem can be expressed as
$$
\operatorname{minimize}_{p \in \mathbb R^{n \times m}} \quad F(p) + G(p)
$$
where
$$
F(p) = \sum_{i,j} p_{ij}(\log(p_{ij}) - a_{ij}) + I_\Omega(p)
$$
and $I_\Omega$ is the indicator function of the set
$$
\Omega = \{ p \in \mathbb R^{n \times m} \mid p_{ij} \geq 0 \, \forall \, i,j\}
$$
(so $I_\Omega$ enforces nonnegativity on $p$)
and
$G$ is the indicator function of the set
$$
S = \{ p \in \mathbb R^{n \times m} \mid
\sum_j p_{ij} = 1 \, \forall \, i \quad \text{and}\quad
\sum_{i=1}^n p_{ij}(b_i - 1) = 0 \, \forall j \}.
$$
Evaluating the prox-operator of $G$ requires projecting onto $S$, which is a linear algebra problem with a standard closed-form solution.
We do not have a closed-form expression for the prox-operator of $F$, but $F$ is still a very simple function (for example, it is fully separable). Can we not evaluate the prox-operator of $F$ numerically to high precision very efficiently?
Evaluating the prox-operator of $F$ reduces to evaluating the prox-operator of the function $w:\mathbb R \to \mathbb R \cup \{ \infty \}$ defined by
$$
w(x) = \begin{cases} x \log(x) & \quad \text{if } x \geq 0, \\
\infty & \quad \text{otherwise.}
\end{cases}
$$
(We interpret $x \log(x) = 0$ when $x = 0$.)
To evaluate the prox-operator of $w$ at $\hat x$ (with parameter $t > 0$), we must solve
\begin{align}
\operatorname{minimize}_x & \quad x \log(x) + \frac{1}{2t} (x - \hat x)^2 \\
\text{subject to} & \quad x \geq 0.
\end{align}
If we visualize the graph of $ v(x) = x \log(x) + \frac{1}{2t} (x - \hat x)^2$, we see that $v$ is initially decreasing (as we move away from the origin, to the right) but then eventually increases to $+\infty$. So, $v$ has a minimizer in the interval $(0,\infty)$ which can be found by setting $v'(x) = 0$. This yields:
$$
1 + \log(x) + \frac{1}{t}(x - \hat x) = 0
\iff x + t \log(x) = \hat x - t.
$$
The function $x + t \log(x)$ is strictly increasing, and it ranges from $-\infty$ to $+\infty$ as $x$ ranges from $0$ to $+\infty$.
Thus, there is a single value of $x$ for which
$x + t \log(x) = \hat x - t$.
I think this value of $x$ can be found to high accuracy with Newton's method (and I think it won't take too many iterations). This allows us to evaluate the prox-operator of $F$ efficiently.
Alternatively (and even better), we can express the solution to
$$
\frac{x}{t} + \log(x) = c
$$
in terms of the Lambert W-function.
Note that
\begin{align}
\frac{x}{t} + \log(x) = c & \iff
\log(e^{x/t}) + \log(x) = \log(e^c) \\
&\iff x e^{x/t} = e^c \\
&\iff (x/t) e^{x/t} = \frac{e^c}{t}\\
&\iff x/t = W(e^c/t) \\
&\iff x = t W(e^c/t).
\end{align}
In Matlab, the Lambert W-function can be evaluated using lambertw.
Now that we have seen how to evaluate the prox-operators of both $F$ and $G$ efficiently, we are able to minimize $F(p) + G(p)$ using the Douglas-Rachford method (which solves convex problems of exactly this form).
Here is my initial attempt to solve the problem, which had a certain flaw (or issue, at least) which is discussed at the end.
The Chambolle-Pock algorithm solves convex optimization problems in the canonical form (or "graph form") minimize $f(x) + g(Cx)$, where the convex functions $f$ and $g$ are "simple" (meaning they have easy prox-operators) and $C$ is a matrix or a linear transformation. The Chambolle-Pock algorithm was extended by Condat (in this paper) to solve convex optimization problems of the form
$$
\tag{$\spadesuit$} \operatorname{minimize}_x \quad h(x) + f(x) + g(Cx)
$$
where $f,g$ and $C$ are as above and the convex function $h$ is differentiable with a Lipschitz continuous gradient. Your problem has the form $(\spadesuit)$ where $C$ is the identity mapping and
\begin{align}
h(p) = -\sum_{ij} \, p_{ij}(a_{ij} - \log(p_{ij})) \\
\end{align}
and $f$ is the indicator function of
$$
\{ p \in \mathbb R^{n \times m}_+ \mid \sum_j p_{ij} = 1 \, \forall i \}
$$
and
$g$ is the indicator function of
$$
\{ p \in \mathbb R^{m \times n} \mid \sum_{i=1}^n p_{ij}(b_i - 1) = 0 \, \forall j \}.
$$
NOTE: A problem (or at least an issue) with this approach is that $h$ does not have a Lipschitz continuous gradient and cannot be extended to a differentiable function $\mathbb R^{n \times m}$. That might mean this approach doesn't work.
I'd be interested to see what @MichaelGrant says about this problem.
