How to model and solve this problem? I have a matrix $P \in M_n(\mathbb N)$, where   

$$ P = 
\begin{bmatrix}
0 &amp;  P_{12}  &amp; \ldots &amp; P_{1n}\\
P_{21}  &amp;  0 &amp; \ldots &amp; P_{2n}\\
\vdots &amp; \vdots &amp; \ddots &amp; \vdots\\
P_{n1}  &amp;   P_{n2}    &amp; \ldots &amp; 0
\end{bmatrix}$$

with $P_{ii} = 0$ for all $i \in \{1,2,\dots,n\}$. I need to find matrices $A, B \in M_n(\mathbb N)$ that satisfy $A + B = P$ and that satisfy the following constraints

$$\forall i\in [\![ 1,n]\!] \sum_{k=1}^{n} A_{ik} = \sum_{k=1}^{n} A_{ki}$$

such that $\sum_{i=1}^{n} \sum_{k=1}^{n} A_{ik}$ is maximized. 

I need to implement an algorithm to solve this problem. On my input data I have approximately: $ n = 16 $  and  $ \sum_{i=1}^{n}\sum_{k=1}^{n} P_{ik} = 60000 $  so the brute-force approach is out of question.
I don't know what could be a good approximation algorithm so my current approach is to reduce the problem to a binary integer programming problem and then apply Branch-and-Cut but I have serious doubt on it effectiveness for this specific problem.

Finding the optimal solution in polynomial time would be perfect (not sure if it's possible), but I can satisfy myself with a good approximation algorithm. Not having a strong background in CS I'm a bit confused, help would be greatly appreciated !
Choose a modelling language, like JuMP, and solve it as a MILP.

An attempt I made with CPLEX solves a random instance ( $n=16$, $max(P)=60000$ ) to optimality almost instantly. 

I've implemented the following

$$ max(z)$$
$$ z = \sum_{i=1}^n\sum_{k=1}^n A(i,k) $$

$$   \sum_{k=1}^n A(i,k) =\sum_{k=1}^n A(k,i)\quad \forall i \in [1,n]$$
$$ A(i,k) + B(i,k) = P(i,k)\quad \forall (i,k) \in [1,n] \times [1,n]$$

$$A\in[0,max(P)-1]^{n\times n}\subset\mathbb{Z}^{n\times n}$$
$$B\in[0,max(P)-1]^{n\times n}\subset\mathbb{Z}^{n\times n}$$