См.
A008300Сгенерировать такие матрицы можно, заведя массив сумм в столбцах, изначально равных

, и в рекурсивном цикле по строкам генерировать все сочетания столбцов с ненулевыми значениями сумм по

, уменьшая значения сумм в столбцах из сочетания на 1.
На PARI/GP примерная программка выглядит так (осуществляется только подсчет):
Код:
{ T(n,k) = local(m);
m=vector(n,i,k);
count=0;
genall(m,k);
count
}
{ genall(m,k) = local(m2,S);
S=select(j->m[j],vector(#m,j,j));
if(#S==0,count++;return);
forvec(v=vector(k,i,[1,#S]),
m2=m;
for(j=1,k,
t=S[v[j]];
m2[t]--;
);
genall(m2,k)
,2)
}