当数组大小是偶数时，为什么这段代码会变慢？

Question

警告：实际上这不是由于 2 的幂，而是由于奇偶性。请参阅编辑部分。

---

我发现了一段代码，它表现出非常奇怪的行为。

代码使用二维数组（大小 x 大小）。 当 size 是 2 的幂时，代码会慢 10% 到 40%，最慢的是 size=32。

我使用英特尔编译器获得了这些结果。如果我使用 gcc 5.4 进行编译，2 的力量就会消失，但 代码会慢 3 倍。在不同的 CPU 上对其进行了测试，我认为它应该具有足够的重现性。

代码：

#define N 10000000

unsigned int tab1[TS][TS];

void simul1() {

  for(int i=0; i<TS; i++)
  for(int j=0; j<TS; j++) {

    if(i > 0)
      tab1[i][j] += tab1[i-1][j];

    if(j > 0)
      tab1[i][j] += tab1[i][j-1];
  }
}

int main() {

  for(int i=0; i<TS; i++)
  for(int j=0; j<TS; j++)
    tab1[j][i] = 0;


  for(int i=0; i<N; i++) {
    tab1[0][0] = 1;
    simul1();
  }

  return tab1[10][10];
}

编译：

icc:
        icc main.c -O3 -std=c99 -Wall -DTS=31 -o i31
        icc main.c -O3 -std=c99 -Wall -DTS=32 -o i32
        icc main.c -O3 -std=c99 -Wall -DTS=33 -o i33

gcc:
        gcc main.c -O3 -std=c99 -Wall -DTS=31 -o g31
        gcc main.c -O3 -std=c99 -Wall -DTS=32 -o g32
        gcc main.c -O3 -std=c99 -Wall -DTS=33 -o g33

Xeon E5 上的结果：

time ./icc31
4.549s

time ./icc32
6.557s

time ./icc33
5.188s


time ./gcc31
13.335s

time ./gcc32
13.669s

time ./gcc33
14.399s

我的问题：

• 为什么数组的大小正好是 32 时会出现性能损失？
• （可选：这里为什么 icc 比 gcc 快 3 倍？）

---

编辑：实际上这是由于奇偶性，并且仅适用于大小> = 32。偶数和奇数之间的性能差异是一致的，并且随着大小变大而减小。这是一个更详细的基准：

（y 轴是每秒的元素数，以百万为单位，用 TS² * N / size / 1000000 获得）

我的 CPU 有 32KB L1 缓存和 256 KB L2

Answer 1

为什么这里的 icc 比 gcc 快 3 倍？

GCC 无法向量化内部循环，因为它报告说数据引用之间存在依赖关系。英特尔的编译器足够聪明，可以将内部循环分成两个独立的部分：

for (int j = 1; j < TS; j++)
    tab1[i][j] += tab1[i-1][j];  // this one is vectorized

for (int j = 1; j < TS; j++)
    tab1[i][j] += tab1[i][j-1];

尝试 1：循环重组

通过将simul1 重写为：

您可能会在 GCC 中获得更好的性能：

void simul1(void)
{
    for (int j = 1; j < TS; j++)
        tab1[0][j] += tab1[0][j-1];

    for (int i = 1; i < TS; i++) {
        for (int j = 0; j < TS; j++)
            tab1[i][j] += tab1[i-1][j];
        for (int j = 1; j < TS; j++)
            tab1[i][j] += tab1[i][j-1];
    }
}

我在 GCC 6.3.0 下使用 Intel Core i5 5200U 驱动的 -O3 -march-native、TS = 32 的结果是：

原版：

real    0m21.110s
user    0m21.004s
sys 0m0.000s

修改：

real    0m8.588s
user    0m8.536s
sys 0m0.000s

尝试 2：前缀和的向量化

经过一番研究，我发现有可能通过向量加法和移位来向量化第二个内循环。算法呈现here。

版本 A：128 位宽向量

#include "emmintrin.h"

void simul1(void)
{
    for (int j = 1; j < TS; j++)
        tab1[0][j] += tab1[0][j-1];

    for (int i = 1; i < TS; i++) {
        for (int j = 0; j < TS; j++)
            tab1[i][j] += tab1[i-1][j];

        for (int stride = 0; stride < TS; stride += 4) {
            __m128i v;
            v = _mm_loadu_si128((__m128i*) (tab1[i] + stride));
            v = _mm_add_epi32(v, _mm_slli_si128(v, sizeof(int)));
            v = _mm_add_epi32(v, _mm_slli_si128(v, 2*sizeof(int)));
            _mm_storeu_si128((__m128i*) (tab1[i] + stride), v);
        }

        for (int stride = 4; stride < TS; stride += 4)
            for (int j = 0; j < 4; j++)
                tab1[i][stride+j] += tab1[i][stride-1];
    }
}

结果：

real    0m7.541s
user    0m7.496s
sys 0m0.004s

版本 B：256 位宽向量

这个比较复杂。考虑ints 的八元素向量：

V = (a, b, c, d, e, f, g, h)

我们可以将其视为两个压缩向量：

(a, b, c, d), (e, f, g, h)

首先，算法执行两个独立的求和：

(a, b, c, d), (e, f, g, h)
+
(0, a, b, c), (0, e, f, g)
=
(a, a+b, b+c, c+d), (e, e+f, f+g, g+h)
+
(0, 0, a, a+b), (0, 0, e, e+f)
=
(a, a+b, a+b+c, a+b+c+d), (e, e+f, e+f+g, e+f+g+h)

然后它将第一个向量的最后一个元素传播到第二个向量的每个元素中，因此它最终产生：

(a, a+b, a+b+c, a+b+c+d), (a+b+c+d+e, a+b+c+d+e+f, a+b+c+d+e+f+g, a+b+c+d+e+f+g+h)

我怀疑这些内在函数可以写得更好，所以有一些改进的潜力。

#include "immintrin.h"

void simul1(void)
{
    for (int j = 1; j < TS; j++)
        tab1[0][j] += tab1[0][j-1];

    for (int i = 1; i < TS; i++) {
        for (int j = 0; j < TS; j++)
            tab1[i][j] += tab1[i-1][j];

        for (int stride = 0; stride < TS; stride += 8) {
            __m256i v;
            v = _mm256_loadu_si256((__m256i*) (tab1[i] + stride));
            v = _mm256_add_epi32(v, _mm256_slli_si256(v, sizeof(int)));
            v = _mm256_add_epi32(v, _mm256_slli_si256(v, 2*sizeof(int)));

            __m256i t = _mm256_setzero_si256();
            t = _mm256_insertf128_si256(t, 
                _mm_shuffle_epi32(_mm256_castsi256_si128(v), 0xFF), 1);

            v = _mm256_add_epi32(v, t);
            _mm256_storeu_si256((__m256i*) (tab1[i] + stride), v);
        }

        for (int stride = 8; stride < TS; stride += 8)
            for (int j = 0; j < 8; j++)
                tab1[i][stride+j] += tab1[i][stride-1];
    }
}

结果（Clang 3.8）：

real    0m5.644s
user    0m5.364s
sys 0m0.004s

Answer 2

看起来像是典型的缓存争用案例。您的代码是这样编写的，以便对相邻的矩阵行和列进行操作。当矩阵行与缓存行对齐并且将存储在同一缓存行中时，这可能会很痛苦。

但是没有很多数据。如果一行被踢出快速 L1 缓存，它可能仍然适合相当快的 L2 缓存。对于 GCC 发出的代码，L2 显然足够快，但 L2 跟不上 ICC 的（矢量化）代码。