矩阵乘法与块
这是我加速矩阵乘法的代码，但它比简单乘法快5％。 我能做些什么来尽可能地提升它？
*正在访问这些表格，例如： C [i，j，n]]用于C [i，j]位置。
void matrixMultFast(float * const C, /* output matrix */ float const * const A, /* first matrix */ float const * const B, /* second matrix */ int const n, /* number of rows/cols */ int const ib, /* size of i block */ int const jb, /* size of j block */ int const kb) /* size of k block */ { int i=0, j=0, jj=0, k=0, kk=0; float sum; for(i=0;i<n;i++) for(j=0;j<n;j++) C[sub2ind(i,j,n)]=0; for(kk=0;kk<n;kk+=kb) { for(jj=0;jj<n;jj+=jb) { for(i=0;i<n;i++) { for(j=jj;j<jj+jb;j++) { sum=C[sub2ind(i,j,n)]; for(k=kk;k<kk+kb;k++) sum += A[sub2ind(i,k,n)]*B[sub2ind(k,j,n)]; C[sub2ind(i,j,n)]=sum; } } } } } // end function 'matrixMultFast4' *它是用C语言编写的，它需要支持C99
您可以做很多很多事情来提高矩阵乘法的效率。
为了研究如何改进基本算法，让我们先来看看我们当前的选项。 当然，天真的实现具有3个循环，其时间复杂度为O(n^3) 。 还有另一种称为Strassen方法的方法，它实现了明显的加速并且具有O(n^2.73)的阶数（但实际上是无用的，因为它没有提供明显的优化方法）。
这是理论上的。 现在考虑矩阵如何存储在内存中。 Row major是标准，但你也找到了专业。 根据方案，转置矩阵可能会因缓存未命中次数减少而提高速度。 理论上的矩阵乘法只是一堆矢量点积和加法。 相同的向量将由多个向量操作，因此将该向量保持在高速缓存中以便更快地访问是有意义的。
现在，随着多核，并行和缓存概念的引入，游戏发生了变化。 如果我们仔细观察一下，我们会看到一个点积只不过是一堆乘法，然后是求和。 这些乘法可以并行完成。 因此，我们现在可以查看数字的并行加载。
现在让我们让事情变得更复杂一些。 在谈论矩阵乘法时，单个浮点和双浮点的大小有所不同。 通常前者是32位，而后者是64位（当然，这取决于CPU）。 每个CPU只有固定数量的寄存器，这意味着您的数字越大，您在CPU中的适应性就越小。 故事的道德是，坚持单浮点，除非你真的需要加倍。
现在我们已经了解了如何调整矩阵乘法的基础知识，不用担心。 您不需要对上面讨论的内容执行任何操作，因为已经有子程序可以执行此操作。 正如评论中所提到的，有GotoBLAS，OpenBLAS，Intel的MKL和Apple的Accelerate框架。 MKL / Accelerate是专有的，但OpenBLAS是一个非常有竞争力的选择。
这是一个很好的小例子，它在我的Macintosh上在几毫秒内乘以2 8k x 8k矩阵：
#include #include #include #include #include int SIZE = 8192; typedef float point_t; point_t* transpose(point_t* A) { point_t* At = (point_t*) calloc(SIZE * SIZE, sizeof(point_t)); vDSP_mtrans(A, 1, At, 1, SIZE, SIZE); return At; } point_t* dot(point_t* A, point_t* B) { point_t* C = (point_t*) calloc(SIZE * SIZE, sizeof(point_t)); int i; int step = (SIZE * SIZE / 4); cblas_sgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans, SIZE/4, SIZE, SIZE, 1.0, &A[0], SIZE, B, SIZE, 0.0, &C[0], SIZE); cblas_sgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans, SIZE/4, SIZE, SIZE, 1.0, &A[step], SIZE, B, SIZE, 0.0, &C[step], SIZE); cblas_sgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans, SIZE/4, SIZE, SIZE, 1.0, &A[step * 2], SIZE, B, SIZE, 0.0, &C[step * 2], SIZE); cblas_sgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans, SIZE/4, SIZE, SIZE, 1.0, &A[step * 3], SIZE, B, SIZE, 0.0, &C[step * 3], SIZE); return C; } void print(point_t* A) { int i, j; for(i = 0; i < SIZE; i++) { for(j = 0; j < SIZE; j++) { printf("%f ", A[i * SIZE + j]); } printf("\n"); } } int main() { for(; SIZE <= 8192; SIZE *= 2) { point_t* A = (point_t*) calloc(SIZE * SIZE, sizeof(point_t)); point_t* B = (point_t*) calloc(SIZE * SIZE, sizeof(point_t)); srand(getpid()); int i, j; for(i = 0; i < SIZE * SIZE; i++) { A[i] = ((point_t)rand() / (double)RAND_MAX); B[i] = ((point_t)rand() / (double)RAND_MAX); } struct timeval t1, t2; double elapsed_time; gettimeofday(&t1, NULL); point_t* C = dot(A, B); gettimeofday(&t2, NULL); elapsed_time = (t2.tv_sec - t1.tv_sec) * 1000.0; // sec to ms elapsed_time += (t2.tv_usec - t1.tv_usec) / 1000.0; // us to ms printf("Time taken for %d size matrix multiplication: %lf\n", SIZE, elapsed_time/1000.0); free(A); free(B); free(C); } return 0; }
此时我还应该提到SSE（流式SIMD扩展），这基本上是你不应该做的事情，除非你使用汇编。 基本上，你正在向你的C代码进行矢量化 ，以使用向量而不是整数。 这意味着您可以操作数据块而不是单个值。 编译器放弃并只是按原样转换代码而不进行自己的优化。 如果做得好，它可以像以前一样加速你的代码 - 你甚至可以触摸O(n^2)的理论底线！ 但很容易滥用SSE，不幸的是大多数人都这样做，最终结果比以前更糟糕。
我希望这能激励你深入挖掘。 矩阵乘法的世界是一个庞大而迷人的世界。 下面，我附上链接以供进一步阅读。
- OpenBLAS
- 更多关于SSE的信息
- 英特尔内在函数