I've successfully implemented randomly generating such matrices for R using a modeling based on network flow. I intend to write up those ideas one day, but haven't found the time yet. Reasearching for that, I read in Randomization of Presence–absence Matrices: Comments and New Algorithms by Miklós and Podani:
The Havel-Hakimi theorem (Havel 1955, Hakimi 1962) states that there exists a matrix Xn,m of 0’s and 1’s with row totals a0=(a1, a2,… , an) and column totals b0=(b1, b2,… , bm) such that bi ≥ bi+1 for every 0 < i < m if and only if another matrix Xn−1,m of 0’s and 1’s with row totals a1=(a2, a3,… , an) and column totals b1=(b1−1, b2−1,… ,ba1−1, ba1+1,… , bm) also exists.
I guess that should be the best method to recursively decide your question.
Phrased in my own words: Choose any row, remove it from the list of totals. Call that removed number k. Also subtract one from the k columns with larges sums. You obtain a description of a smaller matrix, and recurse. If at any point you don't have k columns with non-zero sums, then no such matrix can exist. Otherwise you can recursively build a matching matrix using the reverse process: take the matrix returned by the recursive call, then add one more row with k ones, placed in the columns from whose counts you originally subtracted one.
Implementation
bool satisfiable(std::vector<int> a, std::vector<int> b) {
while (!a.empty()) {
std::sort(b.begin(), b.end(), std::greater<int>());
int k = a.back();
a.pop_back();
if (k > b.size()) return false;
if (k == 0) continue;
if (b[k - 1] == 0) return false;
for (int i = 0; i < k; i++)
b[i]--;
}
for (std::vector<int>::iterator i = b.begin(); i != b.end(); i++)
if (*i != 0)
return false;
return true;
}