在C中有效计算kronecker产品
我对C很陌生，对于我的大多数研究来说，没有太多需要比python更快的东西。 然而，事实certificate我最近所做的工作需要计算相当大的向量/矩阵，因此可能需要C + MPI解决方案。
从数学上讲，任务很简单。 我有很多维数〜40k的向量，并希望计算这些向量的选定对的Kronecker积，然后对这些kronecker积进行求和。
问题是，如何有效地做到这一点？ 以下代码结构是否有任何问题，使用for循环或获得效果？
下面描述的函数vector_size传递长度为vector_size向量A和B ，并计算它们存储在C的kronecker乘积，即vector_size*vector_size矩阵。
void kron(int *A, int *B, int *C, int vector_size) { int i,j; for(i = 0; i < vector_size; i++) { for (j = 0; j < vector_size; j++) { C[i*vector_size+j] = A[i] * B[j]; } } return; }
这对我来说似乎很好，当然（如果我没有做出一些愚蠢的语法错误）产生正确的结果，但我有一种潜在的怀疑，即嵌入式循环不是最佳的。 如果我还有另一种方法，请告诉我。 建议欢迎。
我感谢你的耐心和任何建议。 再一次，我对C非常缺乏经验，但谷歌搜索给我带来了这个查询的一点乐趣。
对于双精度向量（单精度和复数相似），您可以使用BLAS例程DGER （rank-one update）或类似的方法一次一个地进行产品，因为它们都在向量上。 你乘以多少个向量？ 请记住，添加一堆矢量外部产品（您可以将Kronecker产品视为）将最终作为矩阵矩阵乘法，BLAS的DGEMM可以有效处理。 但是，如果真的需要整数运算，则可能需要编写自己的例程。
由于你的循环体完全独立，因此肯定有一种方法可以加速它。 在考虑MPI之前，最简单的就是利用几个核心。 OpenMP应该做得很好。
#pragma omp parallel for for(int i = 0; i < vector_size; i++) { for (int j = 0; j < vector_size; j++) { C[i][j] = A[i] * B[j]; } }
现在许多编译器都支持这一点。
您也可以尝试将一些常见的表达式拖出内部循环，但是像gcc，icc或clang这样的优秀编译器应该自己完成这些：
#pragma omp parallel for for(int i = 0; i < vector_size; ++i) { int const x = A[i]; int * vec = &C[i][0]; for (int j = 0; j < vector_size; ++j) { vec[j] = x * B[j]; } }
BTW，用int索引通常不是正确的做法。 size_t是与索引和对象大小有关的所有内容的正确typedef 。
如果您的编译器支持C99（并且您从未传递与A和B相同的向量），请考虑在支持C99的模式下编译并将函数签名更改为：
void kron(int * restrict A, int * restrict B, int * restrict C, int vector_size);
restrict关键字向编译器承诺A ， B和C指向的数组不是别名（重叠）。 在编写代码时，编译器必须在内循环的每次执行时重新加载A[i] ，因为它必须是保守的，并假设您的C[]存储可以修改A[]值。 在restrict ，编译器可以假设这不会发生。
找到解决方案 （感谢@Jeremiah Willcock）： GSL的BLAS绑定看起来很漂亮。 如果我们逐步选择向量A和B对并将它们添加到某个“运行总计”向量/矩阵C ，则上述kron函数的以下修改版本
void kronadd(int *A, int *B, int *C, int vector_size, int alpha) { int i,j; for(i = 0; i < vector_size; i++) { for (j = 0; j < vector_size; j++) { C[i*vector_size+j] = alpha * A[i] * B[j]; } } return; }
从function上讲，它与BLAS DGER函数（可以作为gsl_blas_dger访问）精确对应。 初始kron函数是DGER， alpha = 0 ， C是正确维度的未初始化（归零）矩阵/向量。
事实certificate，最终可能更容易为这些库简单地使用python绑定。 但是，我想我在尝试解决这些问题时已经学到了很多东西。 在其他回复中有一些更有用的建议，如果你有同样的问题需要处理，请检查出来。 感谢大家！
这是数值计算界常见的问题，最好的办法是使用像Matlab （或其自由软件克隆之一 ）这样经过良好调试的软件包。
你甚至可能找到一个python绑定它，所以你可以摆脱C.
以上所有（可能）比在python中严格编写的代码更快。 如果你需要更快的速度，我会建议一些事情：
- 研究使用Fortran而不是C. Fortran编译器倾向于更好地优化数值计算（如果你使用gcc则会有一个例外，因为它的C和Fortran编译器都使用相同的后端）。
- 考虑并行化您的算法。 我知道Fortran的变种有并行循环语句。 我认为周围有一些C插件可以做同样的事情。 如果您使用的是PC（和单精度），您还可以考虑使用video卡的GPU，这本质上是一个非常便宜的arrays处理器。
另一个易于实现的优化是，如果您知道数组的内部维度可以被n整除，那么将n个赋值语句添加到循环体中，减少必要迭代次数，并对循环进行相应更改数数。
这个策略可以通过在外部循环周围使用switch语句来推广，其中数组大小的大小可以被2,3,4和5整除，或者是最常见的。 这可以带来相当大的性能，并且与建议1和3兼容，以进一步优化/并行化。 一个好的编译器甚至可以为你做这样的事情（也就是循环展开）。
另一个优化是利用指针算法来避免数组索引。 像这样的东西应该做的伎俩：
int i, j; for(i = 0; i < vector_size; i++) { int d = *A++; int *e = B; for (j = 0; j < vector_size; j++) { *C++ = *e++ * d; } }
这也避免了通过将其缓存在局部变量中多次访问A [i]的值，这可能会给你一个小的速度提升。 （请注意，此版本不可并行，因为它会改变指针的值，但仍可以使用循环展开。）
为了解决你的问题，我认为你应该尝试使用Eigen 3，它是一个使用所有矩阵函数的C ++库！
如果你有时间，去看看它的文档！ =）
祝好运 ！