оптимизация арифметики внутри С для петли

Question

У меня есть две функции с петлями, которые выглядят очень похожими. Количество обрабатываемых данных очень велико, поэтому я стараюсь как можно больше оптимизировать циклы. Время выполнения для второй функции составляет 320 секунд, но первая занимает 460 секунд. Может кто-нибудь, пожалуйста, дайте мне какие-либо предложения, что отличает и как оптимизировать вычисления?

int ii, jj; 
    double c1, c2; 

    for (ii = 0; ii < n; ++ii) { 
     a[jj] += b[ii] * c1; 
     a[++jj] += b[ii] * c2; 
    } 

Второй:

int ii, jj; 
    double c1, c2; 

    for (ii = 0; ii < n; ++ii) { 
     b[ii] += a[jj] * c1; 
     b[ii] += a[++jj] * c2; 
    } 

А вот выход ассемблер для первого цикла:

movl -104(%rbp), %eax 
    movq -64(%rbp), %rcx 
    cmpl (%rcx), %eax 
    jge LBB0_12 
## BB#10:        ## in Loop: Header=BB0_9 Depth=5 
    movslq -88(%rbp), %rax 
    movq -48(%rbp), %rcx 
    movsd (%rcx,%rax,8), %xmm0 ## xmm0 = mem[0],zero 
    mulsd -184(%rbp), %xmm0 
    movslq -108(%rbp), %rax 
    movq -224(%rbp), %rcx  ## 8-byte Reload 
    addsd (%rcx,%rax,8), %xmm0 
    movsd %xmm0, (%rcx,%rax,8) 
    movslq -88(%rbp), %rax 
    movq -48(%rbp), %rdx 
    movsd (%rdx,%rax,8), %xmm0 ## xmm0 = mem[0],zero 
    mulsd -192(%rbp), %xmm0 
    movl -108(%rbp), %esi 
    addl $1, %esi 
    movl %esi, -108(%rbp) 
    movslq %esi, %rax 
    addsd (%rcx,%rax,8), %xmm0 
    movsd %xmm0, (%rcx,%rax,8) 
    movl -88(%rbp), %esi 
    addl $1, %esi 
    movl %esi, -88(%rbp) 

и второй:

movl -104(%rbp), %eax 
    movq -64(%rbp), %rcx 
    cmpl (%rcx), %eax 
    jge LBB0_12 
## BB#10:        ## in Loop: Header=BB0_9 Depth=5 
    movslq -108(%rbp), %rax 
    movq -224(%rbp), %rcx  ## 8-byte Reload 
    movsd (%rcx,%rax,8), %xmm0 ## xmm0 = mem[0],zero 
    mulsd -184(%rbp), %xmm0 
    movslq -88(%rbp), %rax 
    movq -48(%rbp), %rdx 
    addsd (%rdx,%rax,8), %xmm0 
    movsd %xmm0, (%rdx,%rax,8) 
    movl -108(%rbp), %esi 
    addl $1, %esi 
    movl %esi, -108(%rbp) 
    movslq %esi, %rax 
    movsd (%rcx,%rax,8), %xmm0 ## xmm0 = mem[0],zero 
    mulsd -192(%rbp), %xmm0 
    movslq -88(%rbp), %rax 
    movq -48(%rbp), %rdx 
    addsd (%rdx,%rax,8), %xmm0 
    movsd %xmm0, (%rdx,%rax,8) 
    movl -88(%rbp), %esi 
    addl $1, %esi 
    movl %esi, -88(%rbp) 

оригинальная функция намного больше, поэтому здесь я опирались только на куски, отвечающие за эти петли. Остальная часть c-кода и его выход ассемблера абсолютно одинаковы для обеих функций.

Answer 1

Структура этого расчета довольно странная, но его можно значительно оптимизировать. Некоторые проблемы с этим кодом:

• перезагрузка данных от указателя после записи другому указателю, который не известен как псевдоним. Я предполагаю, что они не будут псевдонимом, потому что этот алгоритм будет еще более странным, если бы это было разрешено, но если они действительно должны быть псевдонимом, проигнорируйте это. В общем, структурируйте тело цикла как: сначала загрузите все, выполните вычисления, затем сохраните их обратно. Не смешивайте загрузку и хранение, это делает компилятор более консервативным.
• перезагрузка данных, сохраненных в предыдущей итерации. Компилятор может это немного понять, но это усложняет ситуацию. Не делай этого.
• неявно относится к первому и последнему предметам по-разному. Сначала это выглядит как хорошая однородная петля, но из-за своей странной структуры на самом деле это особый корпус первой и последней вещи.

Итак, давайте сначала исправим второй цикл, который проще. Единственная проблема здесь - это первый магазин до b[ii], который должен действительно совершить ^(tm), потому что он может быть псевдонимом a[jj + 1]. Но это тривиально можно записать так, что эта проблема уходит:

for (ii = 0; ii < n; ++ii) { 
    b[ii] += a[jj] * c1 + a[jj + 1] * c2; 
    jj++; 
} 

Вы можете сказать по выходу сборки, что компилятор является счастливее, и, конечно, бенчмаркинг подтверждает это быстрее.

Старого ASM (только основной цикл, а не дополнительный хлам):

.LBB0_14:        # =>This Inner Loop Header: Depth=1 
    vmulpd ymm4, ymm2, ymmword ptr [r8 - 8] 
    vaddpd ymm4, ymm4, ymmword ptr [rax] 
    vmovupd ymmword ptr [rax], ymm4 
    vmulpd ymm5, ymm3, ymmword ptr [r8] 
    vaddpd ymm4, ymm4, ymm5 
    vmovupd ymmword ptr [rax], ymm4 
    add  r8, 32 
    add  rax, 32 
    add  r11, -4 
    jne  .LBB0_14 

Нового ASM (только основной цикл):

.LBB1_20:        # =>This Inner Loop Header: Depth=1 
    vmulpd ymm4, ymm2, ymmword ptr [rax - 104] 
    vmulpd ymm5, ymm2, ymmword ptr [rax - 72] 
    vmulpd ymm6, ymm2, ymmword ptr [rax - 40] 
    vmulpd ymm7, ymm2, ymmword ptr [rax - 8] 
    vmulpd ymm8, ymm3, ymmword ptr [rax - 96] 
    vmulpd ymm9, ymm3, ymmword ptr [rax - 64] 
    vmulpd ymm10, ymm3, ymmword ptr [rax - 32] 
    vmulpd ymm11, ymm3, ymmword ptr [rax] 
    vaddpd ymm4, ymm4, ymm8 
    vaddpd ymm5, ymm5, ymm9 
    vaddpd ymm6, ymm6, ymm10 
    vaddpd ymm7, ymm7, ymm11 
    vaddpd ymm4, ymm4, ymmword ptr [rcx - 96] 
    vaddpd ymm5, ymm5, ymmword ptr [rcx - 64] 
    vaddpd ymm6, ymm6, ymmword ptr [rcx - 32] 
    vaddpd ymm7, ymm7, ymmword ptr [rcx] 
    vmovupd ymmword ptr [rcx - 96], ymm4 
    vmovupd ymmword ptr [rcx - 64], ymm5 
    vmovupd ymmword ptr [rcx - 32], ymm6 
    vmovupd ymmword ptr [rcx], ymm7 
    sub  rax, -128 
    sub  rcx, -128 
    add  rbx, -16 
    jne  .LBB1_20 

Это также получил раскатал больше (автоматически), но более значительная разница (не то, что разворачивание бесполезно, но сокращение накладных расходов цикла обычно не является таким большим, как правило, его могут обрабатывать только порты, которые не заняты векторными инструкциями) - это сокращение в магазинах, которое принимает его от отношения 2/3 (потенциально узкое место по пропускной способности магазина, где половина магазины бесполезны) до 4/12 (узкое место в чем-то, что действительно должно произойти).

Теперь для этого первого цикла, когда вы берете первый и последний итерации, это просто добавив два масштабируется b «S для каждого a, а затем мы помещаем первый и последний итерации назад в отдельно:

a[0] += b[0] * c1; 
for (ii = 1; ii < n; ++ii) { 
    a[ii] += b[ii - 1] * c2 + b[ii] * c1; 
} 
a[n] += b[n - 1] * c2; 

это берет его из этого (заметим, что это даже не векторизации):

.LBB0_3:        # =>This Inner Loop Header: Depth=1 
    vmulsd xmm3, xmm0, qword ptr [rsi + 8*rax] 
    vaddsd xmm2, xmm2, xmm3 
    vmovsd qword ptr [rdi + 8*rax], xmm2 
    vmulsd xmm2, xmm1, qword ptr [rsi + 8*rax] 
    vaddsd xmm2, xmm2, qword ptr [rdi + 8*rax + 8] 
    vmovsd qword ptr [rdi + 8*rax + 8], xmm2 
    vmulsd xmm3, xmm0, qword ptr [rsi + 8*rax + 8] 
    vaddsd xmm2, xmm2, xmm3 
    vmovsd qword ptr [rdi + 8*rax + 8], xmm2 
    vmulsd xmm2, xmm1, qword ptr [rsi + 8*rax + 8] 
    vaddsd xmm2, xmm2, qword ptr [rdi + 8*rax + 16] 
    vmovsd qword ptr [rdi + 8*rax + 16], xmm2 
    lea  rax, [rax + 2] 
    cmp  ecx, eax 
    jne  .LBB0_3 

Для этого:

.LBB1_6:        # =>This Inner Loop Header: Depth=1 
    vmulpd ymm4, ymm2, ymmword ptr [rbx - 104] 
    vmulpd ymm5, ymm2, ymmword ptr [rbx - 72] 
    vmulpd ymm6, ymm2, ymmword ptr [rbx - 40] 
    vmulpd ymm7, ymm2, ymmword ptr [rbx - 8] 
    vmulpd ymm8, ymm3, ymmword ptr [rbx - 96] 
    vmulpd ymm9, ymm3, ymmword ptr [rbx - 64] 
    vmulpd ymm10, ymm3, ymmword ptr [rbx - 32] 
    vmulpd ymm11, ymm3, ymmword ptr [rbx] 
    vaddpd ymm4, ymm4, ymm8 
    vaddpd ymm5, ymm5, ymm9 
    vaddpd ymm6, ymm6, ymm10 
    vaddpd ymm7, ymm7, ymm11 
    vaddpd ymm4, ymm4, ymmword ptr [rcx - 96] 
    vaddpd ymm5, ymm5, ymmword ptr [rcx - 64] 
    vaddpd ymm6, ymm6, ymmword ptr [rcx - 32] 
    vaddpd ymm7, ymm7, ymmword ptr [rcx] 
    vmovupd ymmword ptr [rcx - 96], ymm4 
    vmovupd ymmword ptr [rcx - 64], ymm5 
    vmovupd ymmword ptr [rcx - 32], ymm6 
    vmovupd ymmword ptr [rcx], ymm7 
    sub  rbx, -128 
    sub  rcx, -128 
    add  r11, -16 
    jne  .LBB1_6 

Приятный и векторизованный на этот раз, и многое другое хранение и погрузка продолжаются.

Оба изменения в сочетании сделали его примерно в два раза быстрее на моем ПК, но, конечно, YMMV.

Я все еще, что этот код странный. Обратите внимание, как мы изменяем a[n] в последней итерации первого цикла, а затем используем его в первой итерации второго цикла, в то время как другой a - это просто вид сбоку и смотреть. Странно. Может быть, это действительно так, но, откровенно говоря, это похоже на ошибку.