【问题标题】:How to speed up matrix multiplication in C++?如何加快 C++ 中的矩阵乘法?
【发布时间】:2011-05-17 02:03:09
【问题描述】:

我正在使用这个简单的算法执行矩阵乘法。为了更加灵活,我将对象用于包含动态创建的数组的矩阵。

将此解决方案与我的第一个使用静态数组的解决方案进行比较,速度要慢 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;
}


编辑
我更正了我的问题!我在下面添加了完整的源代码并尝试了您的一些建议:
  • 交换 kj 循环迭代 -> 性能提升
  • dim()operator()() 声明为inline -> 性能提升
  • 通过 const 引用传递参数 -> 性能损失! 为什么?所以我不使用它。

现在的表现几乎与旧 porgram 中的表现相同。也许应该有更多的改进。

但我还有另一个问题:函数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


新计划
矩阵.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 循环的交换顺序 (swap) 和块大小为 13 (block) 的阻塞算法之间的比较。

【问题讨论】:

  • 你的意思是写一个只适用于方阵的矩阵乘法吗?只要内部尺寸相等,就可以定义乘法。
  • 您通过引用传递 a 和 b,对吗?你不是为了调用这个函数而复制两个矩阵吗?
  • 你也可以使用Eigen,这个特别好调。 (不要让 LGPL 许可证吓到您 - 它是一个仅包含标头的库,并且 LGPL 的“病毒”条款不成立。请参阅the FAQ。)
  • @Inverse:这并非普遍适用。你见过class matrix 的复制构造函数代码吗?如果不是,那你只是在胡乱猜测。
  • @multiholle:我发现了主要问题。您的 C 代码有一个巨大的缓冲区溢出。这是重置 dim 全局变量,导致您的乘法提前退出。

标签: c++ arrays benchmarking matrix-multiplication


【解决方案1】:

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

a.dim() 到底有多大?如果矩阵的一行不适合仅仅几个缓存行,那么最好使用块访问模式而不是一次整行。

【讨论】:

    【解决方案2】:

    我在这里进行了大胆的猜测,但是如果您动态分配矩阵会产生如此巨大的差异,那么问题可能出在碎片化上。同样,我不知道底层矩阵是如何实现的。

    您为什么不手动为矩阵分配内存,确保它是连续的,并自己构建指针结构?

    另外,dim() 方法是否有任何额外的复杂性?我也会将其声明为内联。

    【讨论】:

      【解决方案3】:

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

      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 作为最内层索引,cb 都被连续访问,而a 保持不变。

      【讨论】:

        【解决方案4】:
        • 如果可以,请使用 SIMD。如果您假设您使用的平台能够执行大量矢量数学运算,那么您绝对必须使用 VMX 寄存器之类的东西,否则您将遭受巨大的性能损失。
        • 不要通过值传递像 matrix 这样的复杂类型 - 使用 const 引用。
        • 不要在每次迭代中调用函数 - 在循环之外缓存 dim()
        • 虽然编译器通常对此进行有效优化,但最好让调用者为您的函数提供矩阵引用来填充,而不是按类型返回矩阵。在某些情况下,这可能会导致昂贵的复制操作。

        【讨论】:

        • VMX?我们知道他在 PowerPC 上运行吗?
        • 因此我的评论是,“假设您使用的平台有能力这样做”。我的愿景有点模糊——这些天我几乎只在 PPC 上编程,我在这里的一些答案看起来很奇怪。
        • 啊,明白了。我只是不确定,因为您首先提到了 SIMD(一般而言),这是有道理的,然后您将其范围缩小到仅 VMX。
        • 是的,这绝对是我的一个错误。我稍微编辑了答案以使其更有意义:)感谢您指出这一点,我知道我有时住在山洞里!
        • 顺便说一句,不确定您是否知道(但其他阅读您答案的人可能不知道):VMX 的 x86“等效”称为 SSE
        【解决方案5】:

        通过 const 引用传递参数以开始:

        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++ 代码的两倍。我在代码中看不到原因。

        【讨论】:

        • 正如 OP 在已编辑问题中所说, const 引用结果变慢了。它总是值得尝试,但它与许多不同的编译器优化交互,以至于它远非一个安全的选择。有时更快,有时更慢。
        • @jalf:总体上我也同意你的 cmets,但在这种特定情况下,提交的代码有问题。
        • @逆:请仔细阅读文章。我完全同意这篇文章,并且使用按值复制确实非常有用,尤其是在编译器执行 RVO 和 NRVO 时。但是这与这里无关,因为我们需要返回一个新数组(我们无法优化它,结果是一个新数组)。
        • @multiholle:在 C++ 中,使用类几乎没有任何开销。
        【解决方案6】:

        你说你不想修改算法,但这究竟是什么意思?

        展开循环算作“修改算法”吗?使用 SSE/VMX 无论您的 CPU 上可用的 SIMD 指令如何?使用某种形式的blocking 来改善缓存局部性怎么样?

        如果您根本不想重构代码,我怀疑除了已经做出的更改之外,您还可以做更多的事情。其他一切都变成了对算法进行细微更改以实现性能提升的权衡。

        当然,你还是应该看看编译器生成的asm。这将告诉您更多关于可以采取哪些措施来加快代码速度。

        【讨论】:

        • 使用阻塞加速算法,太棒了!
        【解决方案7】:

        这是我对方形浮点矩阵(2D 数组)的快速简单乘法算法的实现。它应该比 chrisaycock 代码快一点,因为它节省了一些增量。

        static void fastMatrixMultiply(const int dim, float* dest, const float* srcA, const float* srcB)
        {
            memset( dest, 0x0, dim * dim * sizeof(float) );
        
            for( int i = 0; i < dim; i++ ) {
                for( int k = 0; k < dim; k++ ) 
                {
                    const float* a = srcA + i * dim + k;
                    const float* b = srcB + k * dim;
                    float* c = dest + i * dim;
        
                    float* cMax = c + dim;
                    while( c < cMax ) 
                    {   
                        *c++ += (*a) * (*b++);
                    }
                }
            }
        }
        

        【讨论】:

          猜你喜欢
          • 2011-05-31
          • 1970-01-01
          • 2015-04-08
          • 2017-01-25
          • 2011-11-30
          • 1970-01-01
          • 1970-01-01
          • 1970-01-01
          • 1970-01-01
          相关资源
          最近更新 更多