Algorithm for a linear optimization problem For the vectors $X=(x_1,\cdots, x_n),~ Y=(y_1,\cdots, y_n)$ and $\alpha=(\alpha_1,\cdots,\alpha_n),~ \beta=(\beta_1,\cdots, \beta_n)\in\mathbb R^n_+$ s.t. $\sum_{k=1}^n\alpha_k~~=~~\sum_{k=1}^n\beta_k~~=~~1$, let $\mathcal P$ be the collection of matrices $p=(p_{i,j})_{1\le i\le n, 1\le j\le n}\in\mathbb R^{n\times n}_+$ satisfying 
\begin{eqnarray}
&&\sum_{j=1}^np_{i,j}~=~\alpha_i,~ \forall 1\le i\le n;~~~~~\sum_{i=1}^np_{i,j}~=~\beta_j,~ \forall 1\le j\le n; \\
&&\sum_{j=1}^np_{i,j}y_j~=~\alpha_ix_i,~ \forall 1\le i\le n.
\end{eqnarray}
Given a matrix $c=(c_{i,j})_{1\le i\le n, 1\le j\le n}$, my question is how to solve numerically the following optimization problem:
\begin{eqnarray}
\sup_{p=(p_{i,j})_{1\le i\le n, 1\le j\le n}\in\mathcal P}~\sum_{i=1}^n\sum_{j=1}^n p_{i,j}c_{i,j}.
\end{eqnarray}
Here assume $\mathcal P\neq \emptyset$, does someone know some related algorithm? Thank you very much!
 A: I am not sure I understand the question. The constraints are a system of linear equations and inequalties in the variables $p_{ij}$ and the objective function is linear. So, this is a box generic linear programming problem.
A: We have the following linear program (LP) in $\mathrm P \in \mathbb R^{n \times n}$
$$\begin{array}{ll} \text{maximize} & \langle \mathrm C, \mathrm P \rangle\\ \text{subject to} & \mathrm P 1_n = \mathrm a\\ & 1_n^{\top} \mathrm P = \mathrm b^{\top}\\ & \mathrm P \mathrm y = \mbox{diag} (\mathrm a) \, \mathrm x\\ & \mathrm P \geq \mathrm O_n\end{array}$$
where $\mathrm a, \mathrm b, \mathrm C, \mathrm x, \mathrm y$ are given. Vectorizing, we obtain the following LP in $\mbox{vec} (\mathrm P) \in \mathbb R^{n^2}$
$$\begin{array}{ll} \text{maximize} & \langle \mbox{vec} (\mathrm C), \mbox{vec} (\mathrm P) \rangle\\ \text{subject to} & (1_n^{\top} \otimes \mathrm I_n) \, \mbox{vec} (\mathrm P) = \mathrm a\\ & (\mathrm I_n \otimes 1_n^{\top}) \, \mbox{vec} (\mathrm P) = \mathrm b\\ & (\mathrm y^{\top} \otimes \mathrm I_n) \, \mbox{vec} (\mathrm P) = \mbox{diag} (\mathrm a) \, \mathrm x\\ & \mbox{vec} (\mathrm P) \geq \mathrm 0_{n^2}\end{array}$$
