I am trying to solve a tricky (for me) 2D binary LP problem. Suppose I have a rectangular matrix B (for benefits), that is I by J and all elements are known and non-stochastic. I want to pick rows and columns to maximize the sum the elements in the intersections of the selected rows and columns, subject to some constraints.

Making all this a bit more formal:

B is matrix with I=927 rows and J=137 columns. It's currently a Stata dataset.
Let binary decision variable $r_i$ indicate whether row $i$ is selected
Let binary decision variable $c_j$ indicate whether column $j$ is selected.
Let a binary decision variable $x_{i,j}$ indicate whether entry $(i,j)$ is selected

The problem is to maximize $\sum_{i,j} B_{i,j} x_{i,j}$ subject to linear constraints:

\begin{align}
x_{i,j} &\le r_i &&\text{for all $i,j$}\\
x_{i,j} &\le c_j &&\text{for all $i,j$}\\
\sum_i r_i & \le 167 \\
20 \le \sum_j c_j & \le 30
\end{align}

If it makes the math any easier, I want 20-30 columns and about 5000 elements selected all together.

Can something like this be done in Mata/Stata?