代码使用2D数组(大小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
解决方法
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由Intel Core i5 5200U供电：
原始版本：
real 0m21.110s user 0m21.004s sys 0m0.000s
改性：
real 0m8.588s user 0m8.536s sys 0m0.000s
尝试2：前缀和的矢量化
在一些reaseach之后,我发现有可能通过向量加法和移位来向量化第二个内环.算法呈现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,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位宽矢量
这个更复杂.考虑整数的八元素向量：
V = (a,b,c,d,e,f,g,h)
我们可以将它视为两个打包向量：
(a,d),(e,h)
首先,算法执行两个独立的求和：
(a,h) + (0,a,c),(0,g) = (a,a+b,b+c,c+d),e+f,f+g,g+h) + (0,a+b),e+f) = (a,a+b+c,a+b+c+d),e+f+g,e+f+g+h)
然后它将第一个向量的最后一个元素传播到第二个向量的每个元素中,因此它最终产生：
(a,(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,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
