重叠数组求和、自动向量化和 restrict

9
Arstechnia最近发表了一篇文章《为什么有些编程语言比其他编程语言更快》,比较了Fortran和C,并提到了数组求和。在Fortran中,假定数组不重叠,因此可以进一步优化。在C/C++中,类型相同的指针可能会重叠,因此通常无法使用此优化。但是,在C/C++中,可以使用restrict__restrict关键字告诉编译器不要假定指针重叠。因此我开始研究这个问题与自动向量化相关的内容。
以下代码能在GCC和MSVC中进行向量化。
void dot_int(int *a, int *b, int *c, int n) {
    for(int i=0; i<n; i++) {
        c[i] = a[i] + b[i];
    }
}

我测试了包含和不包含重叠数组的情况,并且得到了正确的结果。然而,如果我手动使用SSE向量化这个循环的方式不能处理重叠的数组。

int i=0;    
for(; i<n-3; i+=4) {
    __m128i a4 = _mm_loadu_si128((__m128i*)&a[i]);
    __m128i b4 = _mm_loadu_si128((__m128i*)&b[i]);
    __m128i c4 = _mm_add_epi32(a4,b4);
    _mm_storeu_si128((__m128i*)c, c4);
}
for(; i<n; i++) {
    c[i] = a[i] + b[i];
}

接下来我尝试使用__restrict。我认为由于编译器可以假设数组不重叠,因此它不会处理重叠的数组,但无论是GCC还是MSVC都可以在使用__restrict时得到正确的结果,即使数组重叠。

void dot_int_restrict(int * __restrict a, int * __restrict b, int * __restrict c, int n) {
    for(int i=0; i<n; i++) {
        c[i] = a[i] + b[i];
    }
}

为什么使用和不使用__restrict的自动向量化代码可以正确处理重叠数组?

以下是我用来测试的完整代码:

#include <stdio.h>
#include <immintrin.h>
void dot_int(int *a, int *b, int *c, int n) {
    for(int i=0; i<n; i++) {
        c[i] = a[i] + b[i];
    }
    for(int i=0; i<8; i++) printf("%d ", c[i]); printf("\n"); 
}

void dot_int_restrict(int * __restrict a, int * __restrict b, int * __restrict c, int n) {
    for(int i=0; i<n; i++) {
        c[i] = a[i] + b[i];
    }
    for(int i=0; i<8; i++) printf("%d ", c[i]); printf("\n"); 
}

void dot_int_SSE(int *a, int *b, int *c, int n) {
    int i=0;    
    for(; i<n-3; i+=4) {
        __m128i a4 = _mm_loadu_si128((__m128i*)&a[i]);
        __m128i b4 = _mm_loadu_si128((__m128i*)&b[i]);
        __m128i c4 = _mm_add_epi32(a4,b4);
        _mm_storeu_si128((__m128i*)c, c4);
    }
    for(; i<n; i++) {
        c[i] = a[i] + b[i];
    }
    for(int i=0; i<8; i++) printf("%d ", c[i]); printf("\n"); 
}

int main() {
    const int n = 100;
    int a[] = {1,1,1,1,1,1,1,1};
    int b1[] = {1,1,1,1,1,1,1,1,1};
    int b2[] = {1,1,1,1,1,1,1,1,1};
    int b3[] = {1,1,1,1,1,1,1,1,1};

    int c[8];
    int *c1 = &b1[1];
    int *c2 = &b2[1];
    int *c3 = &b3[1];

    dot_int(a,b1,c, 8);
    dot_int_SSE(a,b1,c,8);

    dot_int(a,b1,c1, 8);
    dot_int_restrict(a,b2,c2,8);
    dot_int_SSE(a,b3,c3,8);

}

输出结果(来自 MSVC)

2 2 2 2 2 2 2 2 //no overlap default
2 2 2 2 2 2 2 2 //no overlap with manual SSE vector code
2 3 4 5 6 7 8 9 //overlap default
2 3 4 5 6 7 8 9 //overlap with restrict
3 2 2 2 1 1 1 1 //manual SSE vector code

编辑:

这里还有一个更简单的版本,它生成的代码非常简洁。

void dot_int(int * __restrict a, int * __restrict b, int * __restrict c, int n) {
    a = (int*)__builtin_assume_aligned (a, 16);
    b = (int*)__builtin_assume_aligned (b, 16);
    c = (int*)__builtin_assume_aligned (c, 16);
    for(int i=0; i<n; i++) {
        c[i] = a[i] + b[i];
    }
}

2
你的测试用例没有涵盖输入数据别名可能出现问题的情况。你需要添加一个这样的测试用例:输入和输出指针别名相同的数组,并且加载存储顺序很重要。 - Paul R
1
@PaulR,我认为我已经这样做了。我将输出添加到了我的问题末尾。输出的c数组与输入的b数组重叠。手动SSE代码失败的原因是因为加载存储顺序很重要。 - Z boson
1
好的 - 抱歉 - 或许我看错了你的代码。 - Paul R
你如何知道编译器正在对代码进行向量化?编译器为这个函数生成了什么样的向量化代码? - Apriori
@Apriori,GCC和MSVC两个编译器都报告说它们将循环向量化了。这就是我知道的。汇编比我预期的要长得多,但初步看来,编译器会生成一个向量化版本和一个非向量化版本,然后根据数组是否重叠(内存地址和范围)来决定采取哪条路径。但这仍然只是猜测,我明天会有更多时间来检查。 - Z boson
@Z boson:有道理。我会对你的发现感兴趣。性能数字也可能展示代码正在执行的优化级别。我自己花了一些时间使用 __restrict 关键字,但从未看到任何区别,至少在 x86/x64 上;我想知道在哪些情况下它产生了正面影响。我的理解是在 PowerPC 上更为重要,因为负载命中存储器的问题更大。当然,我们想看到的优化完全不同。 - Apriori
3个回答

5

我不明白问题出在哪里。在Linux/64位、GCC 4.6、-O3、-mtune=native、-msse4.1(也就是一个非常老的编译器/系统)上进行测试时,这段代码:

void dot_int(int *a, int *b, int *c, int n) {
    for(int i=0; i<n; ++i) {
        c[i] = a[i] + b[i];
    }
}

编译成以下内部循环:

.L4:
    movdqu  (%rdi,%rax), %xmm1
    addl    $1, %r8d
    movdqu  (%rsi,%rax), %xmm0
    paddd   %xmm1, %xmm0
    movdqu  %xmm0, (%rdx,%rax)
    addq    $16, %rax
    cmpl    %r8d, %r10d
    ja      .L4
    cmpl    %r9d, %ecx
    je      .L1

虽然这段代码

void dot_int_restrict(int * __restrict a, int * __restrict b, int * __restrict c, int n) {
    for(int i=0; i<n; ++i) {
        c[i] = a[i] + b[i];
    }
}

编译后的结果如下:

.L15:
    movdqu  (%rbx,%rax), %xmm0
    addl    $1, %r8d
    paddd   0(%rbp,%rax), %xmm0
    movdqu  %xmm0, (%r11,%rax)
    addq    $16, %rax
    cmpl    %r10d, %r8d
    jb      .L15
    addl    %r12d, %r9d
    cmpl    %r12d, %r13d
    je      .L10

显然可以看出减少了一个负载。我猜测它正确估计到在执行求和之前没有必要明确加载内存,因为结果不会覆盖任何内容。

还有更多的优化空间--GCC不知道参数比如128位对齐,因此必须生成一个巨大的前导代码来检查是否存在对齐问题(结果可能会有所不同),以及一个后处理代码来处理额外的未对齐部分(或小于128位)。这实际上发生在上述两个版本中。下面是dot_int生成的完整代码:

dot_int:
.LFB626:
        .cfi_startproc
        testl   %ecx, %ecx
        pushq   %rbx
        .cfi_def_cfa_offset 16
        .cfi_offset 3, -16
        jle     .L1
        leaq    16(%rdx), %r11
        movl    %ecx, %r10d
        shrl    $2, %r10d
        leal    0(,%r10,4), %r9d
        testl   %r9d, %r9d
        je      .L6
        leaq    16(%rdi), %rax
        cmpl    $6, %ecx
        seta    %r8b
        cmpq    %rax, %rdx
        seta    %al
        cmpq    %r11, %rdi
        seta    %bl
        orl     %ebx, %eax
        andl    %eax, %r8d
        leaq    16(%rsi), %rax
        cmpq    %rax, %rdx
        seta    %al
        cmpq    %r11, %rsi
        seta    %r11b
        orl     %r11d, %eax
        testb   %al, %r8b
        je      .L6
        xorl    %eax, %eax
        xorl    %r8d, %r8d
        .p2align 4,,10
        .p2align 3
.L4:
        movdqu  (%rdi,%rax), %xmm1
        addl    $1, %r8d
        movdqu  (%rsi,%rax), %xmm0
        paddd   %xmm1, %xmm0
        movdqu  %xmm0, (%rdx,%rax)
        addq    $16, %rax
        cmpl    %r8d, %r10d
        ja      .L4
        cmpl    %r9d, %ecx
        je      .L1
.L3:
        movslq  %r9d, %r8
        xorl    %eax, %eax
        salq    $2, %r8
        addq    %r8, %rdx
        addq    %r8, %rdi
        addq    %r8, %rsi
        .p2align 4,,10
        .p2align 3
.L5:
        movl    (%rdi,%rax,4), %r8d
        addl    (%rsi,%rax,4), %r8d
        movl    %r8d, (%rdx,%rax,4)
        addq    $1, %rax
        leal    (%r9,%rax), %r8d
        cmpl    %r8d, %ecx
        jg      .L5
.L1:
        popq    %rbx
        .cfi_remember_state
        .cfi_def_cfa_offset 8
        ret
.L6:
        .cfi_restore_state
        xorl    %r9d, %r9d
        jmp     .L3
        .cfi_endproc

现在在你的情况下,整数实际上是未对齐的(因为它们在堆栈上),但是如果您可以使它们对齐并告诉GCC,那么您可以改进代码生成:

typedef int intvec __attribute__((vector_size(16)));

void dot_int_restrict_alig(intvec * restrict a, 
                           intvec * restrict b, 
                           intvec * restrict c, 
                           unsigned int n) {
    for(unsigned int i=0; i<n; ++i) {
        c[i] = a[i] + b[i];
    }
}

这将生成没有前导的代码:
dot_int_restrict_alig:
.LFB628:
        .cfi_startproc
        testl   %ecx, %ecx
        je      .L23
        subl    $1, %ecx
        xorl    %eax, %eax
        addq    $1, %rcx
        salq    $4, %rcx
        .p2align 4,,10
        .p2align 3
.L25:
        movdqa  (%rdi,%rax), %xmm0
        paddd   (%rsi,%rax), %xmm0
        movdqa  %xmm0, (%rdx,%rax)
        addq    $16, %rax
        cmpq    %rcx, %rax
        jne     .L25
.L23:
        rep
        ret
        .cfi_endproc

请注意对齐的 128 位加载指令的用法 (movdqaa 表示对齐,与 movdqu,不对齐相反)。

我认为代码还有一个非向量化版本。至少在http://gcc.godbolt.org/上是这样的。我的猜测是编译器会制作两个版本,然后当它进入函数时,决定数组是否重叠,然后选择向量化或非向量化版本。无论如何,明天我会仔细查看你的答案。 - Z boson
你的函数 dot_int_restrict_alig 生成的代码与我手动编写的 SSE 函数几乎完全相同,只是它使用了对齐的加载/存储,并且在重叠数组上核心转储。你发布的汇编代码无法处理重叠数组。必须有更多的代码(你没有发布)来处理重叠数组的情况。 - Z boson
不,我的意思是从其他函数中。 - Z boson
实际上,如果您不进行矢量化编译(只需将 -O3 更改为 -O2),则代码仅为您的代码中的 .L5 部分,因此这显然是非矢量化部分。 - Z boson
我使用__builtin_assume_aligned更新了我的答案结尾,这样可以得到更小的前导部分。 - Z boson
显示剩余4条评论

2
我现在明白了情况。原来MSVC和GCC在使用__restrict时给出不同的结果。MSVC能够正确处理重叠,但GCC不能。因此可以得出结论,MSVC忽略了__restrict关键字而GCC则利用它进一步优化。
GCC的输出结果:
2 2 2 2 2 2 2 2 //no overlap default
2 2 2 2 2 2 2 2 //no overlap with manual SSE vector code
2 3 4 5 6 7 8 9 //overlap without __restrict
2 2 2 2 3 2 2 2 //overlap with __restrict
3 2 2 2 1 1 1 1 //manual SSE vector code

我们可以生成一个纯向量化的函数,它产生的汇编行数几乎与C语言相同(在GCC 4.9.0中使用-O3参数生成所有代码):

void dot_int(int * __ restrict a, int * __restrict b, int * __restrict c) {
    a = (int*)__builtin_assume_aligned (a, 16);
    b = (int*)__builtin_assume_aligned (b, 16);
    c = (int*)__builtin_assume_aligned (c, 16);
    for(int i=0; i<1024; i++) {
        c[i] = a[i] + b[i];
    }
}

生成

dot_int(int*, int*, int*):
    xorl    %eax, %eax
.L2:
    movdqa  (%rdi,%rax), %xmm0
    paddd   (%rsi,%rax), %xmm0
    movaps  %xmm0, (%rdx,%rax)
    addq    $16, %rax
    cmpq    $4096, %rax
    jne .L2
    rep ret

然而,如果我们去掉了允许a与c重叠的__restrict on a,我通过查看汇编代码确定:

void dot_int(int * a, int * __restrict b, int * __restrict c) {
        a = (int*)__builtin_assume_aligned (a, 16);
        b = (int*)__builtin_assume_aligned (b, 16);
        c = (int*)__builtin_assume_aligned (c, 16);
        for(int i=0; i<1024; i++) {
            c[i] = a[i] + b[i];
        }
    }

等同于

__attribute__((optimize("no-tree-vectorize")))
inline void dot_SSE(int * __restrict a, int * __restrict b, int * __restrict c) {
    for(int i=0; i<1024; i+=4) {    
        __m128i a4 = _mm_load_si128((__m128i*)&a[i]);
        __m128i b4 = _mm_load_si128((__m128i*)&a[i]);
        __m128i c4 = _mm_add_epi32(a4,b4);
        _mm_store_si128((__m128i*)&c[i],c4);
    }
}
__attribute__((optimize("no-tree-vectorize")))
void dot_int(int * __restrict a, int * __restrict b, int * __restrict c) {
    a = (int*)__builtin_assume_aligned (a, 16);
    b = (int*)__builtin_assume_aligned (b, 16);
    c = (int*)__builtin_assume_aligned (c, 16);
    int pass = 1;
    if((c+4)<a || (a+4)<c) pass = 0;
    if(pass) {
        for(int i=0; i<1024; i++) {
            c[i] = a[i] + b[i];
        }   
    }
    else {
        dot_SSE(a,b,c);
    }
}

换句话说,如果a和c指针之间的距离小于16个字节(|a-c|<4),则代码会分支到非矢量化形式。这证实了我的猜测,即自动矢量化代码包括矢量化和非矢量化版本以处理重叠。
在重叠的数组上,这将得到正确的结果:2 3 4 5 6 7 8 9。

第一个区间 i1 = a-b 我认为是不必要的,因为如果输入重叠也没关系。 - Z boson

2
如果您在重叠的数组上使用 "restrict",则会出现未定义的行为。这就是“重叠限制”情况下出现的情况。未定义行为意味着任何事情都可能发生。确实发生了。仅仅是巧合,行为与没有使用“restrict”时相同。绝对正确。它完全符合“任何事情都可能发生”的定义。没有什么好抱怨的。

我不同意。当数组不重叠时,我手动完成的解决方案是最优的。如果编译器假定数组不重叠,它应该生成无法获得相同结果的代码。我认为更有可能的是编译器忽略了restrict关键字。 - Z boson
另一个问题是,即使没有限制,为什么要将此函数向量化?当您无法假设数组不重叠时,如何有效地进行向量化?我不认为这可以有效地完成。 - Z boson

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接