更快的方法来做多维矩阵加法？

Question

我有一个大小为 (m * l * 4) 的矩阵 A，m 的大小约为 100,000 且 l=100。列表的大小总是等于 n 并且 n

void MatrixAddition(int l, int n, vector<int>& list, int ***A,int ***C,int cluster)
{
    for(int i=0;i<l;i++)
    {
         for(int j=0;j<4;j++)
              C[cluster][i][j]=0;
    }   

for (int i = 0; i < l; i++)
{
        for(int j=0;j<n;++j)
    {
        for(int k=0;k<4;k++)
            C[cluster][i][k]+=A[list[j]][i][k];
    }
}

}

我使用 gprof 计算整个代码中每个函数花费的时间，我发现 MatrixAddition 函数花费了我 60% 的时间。是否有任何替代方法来编写此函数，以减少我的运行时间。

时间秒秒呼叫毫秒/呼叫毫秒/呼叫名称
52.00 7.85 7.85 20 392.60 405.49 MatrixAddition(int, int, std::vector >&, int***, int***, int)

Answer 1

在第二部分交换循环 i 和循环 j。这将使函数对缓存更加友好。

for(int j=0;j<n;++j)
{
    for (int i = 0; i < l; i++)
    {
        for(int k=0;k<4;k++)
            C[cluster][i][k]+=A[list[j]][i][k];
    }
}

另外，我希望你没有忘记 -O3 标志。

Answer 2

（更新：早期版本的索引错误。这个版本很容易自动矢量化）。

使用 C 多维数组（而不是指向指针的数组），或以i*cols + j 为索引的平面一维数组，因此内存是连续的。这将对硬件预取以充分利用内存带宽的有效性产生巨大影响。具有来自另一个负载的地址的负载确实会降低性能，或者相反，提前知道可预测的地址会很有帮助，因为负载可以在需要它们之前就开始（由于无序执行）。

另外，@user31264 的答案是正确的，您需要交换循环，因此j 上的循环是最外层的。这很好，但仅靠它自己还远远不够。

这也将允许编译器自动矢量化。实际上，我很难让 gcc 很好地自动矢量化。 （但那可能是因为我的索引错误，因为我第一次只查看代码。所以编译器不知道我们在连续内存上循环。）

---

我在Godbolt compiler explorer 上玩过它。

我终于从这个版本中得到了很好的编译器输出，它将 A 和 C 作为平面一维数组并自己进行索引：

void MatrixAddition_contiguous(int rows, int n, const  vector<int>& list,
                               const int *__restrict__ A, int *__restrict__ C, int cluster)
  // still auto-vectorizes without __restrict__, but checks for overlap and
  // runs a scalar loop in that case
{
  const int cols = 4;  // or global constexpr or something
  int *__restrict__ Ccluster = C + ((long)cluster) * rows * cols;

  for(int i=0;i<rows;i++)
    //#pragma omp simd  
    for(int k=0;k<4;k++)
      Ccluster[cols*i + k]=0;

  for(int j=0;j < cols;++j) { // loop over clusters in A in the outer-most loop
    const int *__restrict__ Alistj = A + ((long)list[j]) * rows * cols;
    // #pragma omp simd    // Doesn't work: only auto-vectorizes with -O3
    // probably only -O3 lets gcc see through the k=0..3 loop and treat it like one big loop
    for (int i = 0; i < rows; i++) {
      long row_offset = cols*i;
      //#pragma omp simd  // forces vectorization with 16B vectors, so it hurts AVX2
      for(int k=0;k<4;k++)
        Ccluster[row_offset + k] += Alistj[row_offset + k];
    }
  }
}

手动提升list[j] 绝对有助于编译器意识到存储到C 不会影响将从list[j] 加载的索引。可能不需要手动吊起其他东西。

提升A[list[j]]，而不仅仅是list[j]，是previous approach where I had the indexing wrong 的产物。只要我们尽可能从list[j]提升负载，编译器即使不知道list不与C重叠也能做好。

内部循环，gcc 5.3 针对 x86-64 -O3 -Wall -march=haswell -fopenmp（和 -fverbose-asm）是：

.L26:
    vmovdqu ymm0, YMMWORD PTR [r8+rax]        # MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B], MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B]
    vpaddd  ymm0, ymm0, YMMWORD PTR [rdx+rax]   # vect__71.75, MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B], MEM[base: vectp.73_90, index: ivtmp.91_26, offset: 0B]
    add     r12d, 1   # ivtmp.88,
    vmovdqu YMMWORD PTR [r8+rax], ymm0        # MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B], vect__71.75
    add     rax, 32   # ivtmp.91,
    cmp     r12d, r9d # ivtmp.88, bnd.66
    jb      .L26        #,

所以它一次执行 8 次加法，使用 AVX2 vpaddd，未对齐的加载和未对齐的存储回 C。

由于这是自动-矢量化，它应该可以使用 ARM NEON、PPC Altivec 或任何可以进行打包 32 位加法的任何东西编写好的代码。

我无法让 gcc 用 -ftree-vectorizer-verbose=2 告诉我任何事情，但 clang 的 -Rpass-analysis=loop-vectorize 有点帮助。