是否有一种就地乘以平方矩阵的算法？
问题内容：
用于乘以4x4矩阵的朴素算法如下所示：
void matrix_mul(double out[4][4], double lhs[4][4], double rhs[4][4]) {
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
out[i][j] = 0.0;
for (int k = 0; k < 4; ++k) {
out[i][j] += lhs[i][k] * rhs[k][j];
}
}
}
}
显然，如果out == lhs或out == rhs（此处==表示引用相等），此算法将给出假结果。是否有一个版本允许其中一种或两种情况不只是复制矩阵？如有必要，我很高兴为每种情况提供不同的功能。
我找到了这篇论文，但它讨论了Strassen-
Winograd算法，该算法对我的小型矩阵而言过于严格。这个问题的答案似乎表明，如果out == lhs && out == rhs（即，我们试图对矩阵求平方），那么就无法就位，但是即使没有令人信服的证据也无法证明。
问题答案：
我对这个答案并不感到兴奋（我主要是为了使“显然不能完成的”人群沉默而发布它），但我怀疑使用真正的就地算法可以做得更好（
O（1）个用于存储两个nxn矩阵的额外存储字）。让我们将两个矩阵相乘A和B。假定A和B没有别名。
如果A是上三角的，那么乘法问题看起来像这样。
[a11 a12 a13 a14] [b11 b12 b13 b14]
[ 0 a22 a23 a24] [b21 b22 b23 b24]
[ 0 0 a33 a34] [b31 b32 b33 b34]
[ 0 0 0 a44] [b41 b42 b43 b44]
我们可以将乘积计算为B，如下所示。将B的第一行乘以a11。将a12B的第二行与第一行相乘。将a13B的第三行与第一行相乘。将a14B的第四行与第一行相乘。
现在，我们用正确的乘积覆盖了B的第一行。幸运的是，我们不再需要它。将B的第二行乘以a22。将a23B的第三行与第二行相乘。（你明白了。）
同样，如果A是下三角单元，则乘法问题看起来像这样。
[ 1 0 0 0 ] [b11 b12 b13 b14]
[a21 1 0 0 ] [b21 b22 b23 b24]
[a31 a32 1 0 ] [b31 b32 b33 b34]
[a41 a42 a43 1 ] [b41 b42 b43 b44]
将a43时间添加到B的第三行到第四行。将a42B的第二行与第四行相乘。将a41B的第一行与第四行相乘。将a32B的第二行与第三行相乘。（你明白了。）
完整的算法是先对A进行LU分解，再将UB乘以B，再将LB乘以B，然后对LU进行适当的分解（我不确定是否有人这样做，但是似乎很容易将A逆转脚步）。实际上有大约一百万个理由不执行此操作，其中两个原因是A可能无法进行LU分解，并且通常不会使用浮点算术精确地重构A。
