如何加速C++中的矩阵乘法？

Question

我用这个简单的算法进行矩阵乘法.为了更灵活,我使用了包含动态创建数组的matricies对象.

将此解决方案与我的第一个解决方案与静态数组进行比较,速度慢了4倍.我该怎么做才能加快数据访问速度？我不想改变算法.

 matrix mult_std(matrix a, matrix b) {
 matrix c(a.dim(), false, false);
 for (int i = 0; i < a.dim(); i++)
  for (int j = 0; j < a.dim(); j++) {
   int sum = 0;
   for (int k = 0; k < a.dim(); k++)
    sum += a(i,k) * b(k,j);
   c(i,j) = sum;
  }

 return c;
}

---

编辑

---

我纠正了我的问题!我在下面添加了完整的源代码并尝试了一些建议:

• 交换k和j循环迭代 - >性能改进
• 声明dim()并operator()() 作为inline- >性能改进
• 通过const引用传递参数 - > 性能损失!为什么？所以我不使用它.

现在的表现与现在的表现几乎相同.也许应该有一点改进.

但我有另一个问题:我在函数中出现内存错误mult_strassen(...).为什么？
terminate called after throwing an instance of 'std::bad_alloc'
what(): std::bad_alloc

---

旧程序
main.c http://pastebin.com/qPgDWGpW

c99 main.c -o matrix -O3

---

新程序
matrix.h http://pastebin.com/TYFYCTY7
matrix.cpp http://pastebin.com/wYADLJ8Y
main.cpp http://pastebin.com/48BSqGJr

g++ main.cpp matrix.cpp -o matrix -O3.

---

编辑

---

这是一些结果.标准算法(std),j和k循环(交换)的交换顺序与块大小为13(块)的阻塞算法之间的比较.

Answer 1

说到加速,如果你交换k和j循环迭代的顺序,你的函数将更加缓存友好:

matrix mult_std(matrix a, matrix b) {
 matrix c(a.dim(), false, false);
 for (int i = 0; i < a.dim(); i++)
  for (int k = 0; k < a.dim(); k++)
   for (int j = 0; j < a.dim(); j++)  // swapped order
    c(i,j) += a(i,k) * b(k,j);

 return c;
}

那是因为k最内层循环上的索引会导致b每次迭代都出现缓存缺失.随着j作为最内层的指数,无论是c与b被连续访问,而a原地踏步.

Answer 2

确保成员dim()和operator()()被声明为内联，并且编译器优化已打开。然后使用-funroll-loops（在 gcc 上）等选项。

到底有多大a.dim()？如果矩阵的一行不能容纳几个缓存行，那么最好使用块访问模式而不是一次整行。

Answer 3

首先通过常量引用传递参数：

matrix mult_std(matrix const& a, matrix const& b) {

为了向您提供更多详细信息，我们需要了解所使用的其他方法的详细信息。
要回答为什么原始方法快 4 倍，我们需要查看原始方法。

这个问题无疑是你的，因为这个问题之前已经被解决了一百万次。

此外，在提出此类问题时，始终提供带有适当输入的可编译源代码，以便我们可以实际构建和运行代码并查看发生了什么。

如果没有代码，我们只是猜测。

编辑

修复了原始 C 代码中的主要错误（缓冲区溢出）后

我已经更新了代码以进行公平比较并排运行测试：

 // INCLUDES -------------------------------------------------------------------
 #include <stdlib.h>
 #include <stdio.h>
 #include <sys/time.h>
 #include <time.h>

 // DEFINES -------------------------------------------------------------------
 // The original problem was here. The MAXDIM was 500. But we were using arrays
 // that had a size of 512 in each dimension. This caused a buffer overrun that
 // the dim variable and caused it to be reset to 0. The result of this was causing
 // the multiplication loop to fall out before it had finished (as the loop was
 // controlled by this global variable.
 //
 // Everything now uses the MAXDIM variable directly.
 // This of course gives the C code an advantage as the compiler can optimize the
 // loop explicitly for the fixed size arrays and thus unroll loops more efficiently.
 #define MAXDIM 512
 #define RUNS 10

 // MATRIX FUNCTIONS ----------------------------------------------------------
 class matrix
 {
 public:
 matrix(int dim)
       : dim_(dim)
 {
         data_ = new int[dim_ * dim_];

 }

     inline int dim() const {
                         return dim_;
                 }
                 inline int& operator()(unsigned row, unsigned col) {
                         return data_[dim_*row + col];
                 }

                 inline int operator()(unsigned row, unsigned col) const {
                         return data_[dim_*row + col];
                 }

 private:
     int dim_;
     int* data_;
 };

// ---------------------------------------------------
 void random_matrix(int (&matrix)[MAXDIM][MAXDIM]) {
         for (int r = 0; r < MAXDIM; r++)
                 for (int c = 0; c < MAXDIM; c++)
                         matrix[r][c] = rand() % 100;
 }
 void random_matrix_class(matrix& matrix) {
         for (int r = 0; r < matrix.dim(); r++)
                 for (int c = 0; c < matrix.dim(); c++)
                         matrix(r, c) = rand() % 100;
 }

 template<typename T, typename M>
 float run(T f, M const& a, M const& b, M& c)
 {
         float time = 0;

         for (int i = 0; i < RUNS; i++) {
                 struct timeval start, end;
                 gettimeofday(&start, NULL);
                 f(a,b,c);
                 gettimeofday(&end, NULL);

                 long s = start.tv_sec * 1000 + start.tv_usec / 1000;
                 long e = end.tv_sec * 1000 + end.tv_usec / 1000;

                 time += e - s;
         }
         return time / RUNS;
 }
 // SEQ MULTIPLICATION ----------------------------------------------------------
  int* mult_seq(int const(&a)[MAXDIM][MAXDIM], int const(&b)[MAXDIM][MAXDIM], int (&z)[MAXDIM][MAXDIM]) {
          for (int r = 0; r < MAXDIM; r++) {
                  for (int c = 0; c < MAXDIM; c++) {
                          z[r][c] = 0;
                          for (int i = 0; i < MAXDIM; i++)
                                  z[r][c] += a[r][i] * b[i][c];
                  }
          }
  }
  void mult_std(matrix const& a, matrix const& b, matrix& z) {
          for (int r = 0; r < a.dim(); r++) {
                  for (int c = 0; c < a.dim(); c++) {
                          z(r,c) = 0;
                          for (int i = 0; i < a.dim(); i++)
                                  z(r,c) += a(r,i) * b(i,c);
                  }
          }
  }

  // MAIN ------------------------------------------------------------------------
  using namespace std;
  int main(int argc, char* argv[]) {
          srand(time(NULL));

          int matrix_a[MAXDIM][MAXDIM];
          int matrix_b[MAXDIM][MAXDIM];
          int matrix_c[MAXDIM][MAXDIM];
          random_matrix(matrix_a);
          random_matrix(matrix_b);
          printf("%d ", MAXDIM);
          printf("%f \n", run(mult_seq, matrix_a, matrix_b, matrix_c));

          matrix a(MAXDIM);
          matrix b(MAXDIM);
          matrix c(MAXDIM);
          random_matrix_class(a);
          random_matrix_class(b);
          printf("%d ", MAXDIM);
          printf("%f \n", run(mult_std, a, b, c));

          return 0;
  }

现在的结果：

$ g++ t1.cpp
$ ./a.exe
512 1270.900000
512 3308.800000

$ g++ -O3 t1.cpp
$ ./a.exe
512 284.900000
512 622.000000

由此我们可以看出，在完全优化后，C 代码的速度大约是 C++ 代码的两倍。我在代码中看不到原因。