使用 GSL 计算矩阵 Kronecker 积的有效方法答案

【问题标题】：Efficient way to compute Kronecker product of matrices with GSL使用 GSL 计算矩阵 Kronecker 积的有效方法
【发布时间】：2012-12-05 11:44:59
【问题描述】：

我算法的瓶颈是我的函数 Kronecker Product，叫做 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 的高效实现吗？我找不到合适的例程。

【问题讨论】：

你最后解决了吗？我有同样的问题，任何帮助将不胜感激。

标签： c algorithm matrix product gsl

【解决方案1】：

仅从表面上看，我就可以看到您日常生活中可能存在的许多瓶颈：

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

【讨论】：

【解决方案2】：

您可以通过“阻塞”和更有效地利用缓存内存来显着提高性能。

看看这个paper。 Is 具有伪代码，我认为您将能够轻松地将其转换为 C 代码。它还有一种算法可以在给定缓存大小和矩阵参数的情况下计算出最佳块大小。

【讨论】：