20

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

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


新程序
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 循环的交换顺序 (swap) 和块大小为 13 (block) 的阻塞算法之间的比较。 替代文字

4

7 回答 7

35

说到加速,如果交换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保持不变。

于 2010-11-29T03:57:52.120 回答
4

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

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

于 2010-11-29T03:47:59.627 回答
4

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

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

如果您根本不想重构您的代码,我怀疑您可以做的不仅仅是您已经做出的更改。其他一切都成为对算法进行微小更改以实现性能提升的权衡。

当然,你还是应该看看编译器生成的asm。这将告诉你更多关于可以做些什么来加速代码。

于 2010-11-29T15:49:02.033 回答
3
  • 如果可以,请使用 SIMD。如果您假设您使用的是能够这样做的平台,那么您绝对必须使用 VMX 寄存器之类的东西,否则您将遭受巨大的性能损失。
  • 不要matrix通过值传递复杂类型 - 使用 const 引用。
  • 不要在每次迭代中调用函数 -dim()在循环之外缓存。
  • 尽管编译器通常会对此进行有效优化,但最好让调用者为您的函数提供矩阵引用来填充,而不是按类型返回矩阵。在某些情况下,这可能会导致昂贵的复制操作。
于 2010-11-29T03:58:18.690 回答
1

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

于 2010-11-29T04:19:42.493 回答
1

这是我对方形浮点矩阵(二维数组)的快速简单乘法算法的实现。它应该比 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++);
            }
        }
    }
}
于 2014-07-18T13:24:14.953 回答
0

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

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

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

于 2010-11-29T03:48:32.557 回答