在C中有效计算kronecker产品

Question

我对C很陌生,对于我的大多数研究来说,没有太多需要比python更快的东西.然而,事实证明我最近所做的工作需要计算相当大的向量/矩阵,因此可能有一个C + MPI解决方案.

从数学上讲,任务很简单.我有很多维数〜40k的向量,并希望计算这些向量的选定对的Kronecker积,然后对这些kronecker积进行求和.

问题是,如何有效地做到这一点？以下代码结构是否有任何问题,使用for循环或获得效果？

该函数kron如下所述传递载体A和B长度vector_size,并计算它们的Kronecker乘积,它存储在C,一个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非常缺乏经验,但谷歌搜索给我带来了这个查询的一点乐趣.

Answer 1

由于你的循环体完全独立,因此肯定有一种方法可以加速它.在考虑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];
    }
}

顺便说一下,索引int通常不是正确的做法.对于与索引和对象大小有关的所有事情都是size_t正确的typedef.

Answer 2

对于双精度向量（单精度和复数相似），您可以使用 BLAS 例程DGER（秩一更新）或类似的例程一次计算一个乘积，因为它们都在向量上。你要乘多少个向量？请记住，添加一堆向量外积（您可以将其视为克罗内克积）最终会形成矩阵-矩阵乘法，BLASDGEMM可以有效地处理该乘法。不过，如果您确实需要整数运算，您可能需要编写自己的例程。