c - 使用 GSL 计算矩阵的克罗内克积的有效方法

Question

我的算法的瓶颈是我的函数 Kronecker Product called KPro:

gsl_matrix *KPro(gsl_matrix *a, gsl_matrix *b) {
    int i, j, k, l;
    int m, p, n, q;
    m = a->size1;
    p = a->size2;
    n = b->size1;
    q = b->size2;

    gsl_matrix *c = gsl_matrix_alloc(m*n, p*q);
    double da, db;

     for (i = 0; i < m; i++)    {
          for (j = 0; j < p; j++)   {
              da = gsl_matrix_get (a, i, j);
              for (k = 0; k < n; k++)   {
                  for (l = 0; l < q; l++)   {
                      db = gsl_matrix_get (b, k, l);
                      gsl_matrix_set (c, n*i+k, q*j+l, da * db);                
                  }
              }
          }
      }

    return c;
}

您知道使用 GSL 的高效实现吗？我找不到合适的套路。

Answer 1

仅从表面上看，我可以看到您的日常工作中有很多可能的瓶颈:

1. 重用矩阵 c 而不是每次都重新分配它，即将函数堆栈变量提升为类成员或静态提升为文件。分配一次并分配到最大可能的问题大小。
2. 调用所有这些 gsl_matrix_get 和 gsl_matrix_set 肯定会阻止编译器对您的代码进行自动向量化，请考虑使用基于模板的矩阵实现以及重载或内联运算符和直接内存访问。
3. 想想您正在使用的矩阵排序:它是行主要的吗？还是专栏专业？缓存未命中比你在那里做的任何其他事情都要昂贵。您想利用空间局部性和重用，通过以最内层循环(计算发生的地方)访问已预取的相邻矩阵元素的方式重新排序循环来实现。
4. 进行对齐的内存分配，使矢量化更容易、更高效。
5. 考虑使用循环展开和阻塞