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

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

至强 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 3 倍 比这里的 gcc 快？）

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

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

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

Answer 1

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

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

Answer 2

Why is icc 3x faster than gcc here?

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 下的 -O3 -march-native、TS = 32 由英特尔酷睿 i5 5200U 提供支持的结果是：

原始版本：

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

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

Why is this code slower when the array size is even?

c

optimization

performance

gcc

icc

尝试 1：循环重组

尝试2：前缀和向量化

版本 A：128 位宽向量

版本 B：256 位宽向量