我先来抛砖吧，
设对应的未知数是x11,x12,x21,x22,x23,x31,x32,x33，各自与行列的位置对应，条件方程有六个，分别对应三行与三列，如果用穷举法，大概的组合有21^8之多，大约300多亿，非常耗时。不过，幸好这些条件方程是相互联系的，因此，可以用替换法来缩减，每一个条件方程式缩减一个变量，最终可以用5个条件方程缩减5个变量，变成下面的等式：

x22*x32-(4+x22*(9-x31-x32))/x31=5

保留下来需要穷举的变量只有三个，循环量大概21^2*20，约8000多个，这么少的循环，即便是R语言，都表示毫无压力。代码如下：

结果出来后，还需要把NA值和其他变量不符合0-20条件的行去除，
整个过程在R上的运行速度非常快，在我的i3电脑上都只需要大概11毫秒，因此，我也懒得再去优化了。

厉害?

非常有趣

楼主什么垃圾水平，for循环用个不停，一看就是初学者。什么素质，还好意思来我帖子下灌水，公开我代码？