Почему векторизация цикла не улучшает производительность

Я изучаю эффект векторизации на производительность программы. В связи с этим я написал следующий код:

#include <stdio.h>
#include <sys/time.h>
#include <stdlib.h>

#define LEN 10000000

int main(){

    struct timeval stTime, endTime;

    double* a = (double*)malloc(LEN*sizeof(*a));
    double* b = (double*)malloc(LEN*sizeof(*b));
    double* c = (double*)malloc(LEN*sizeof(*c));

    int k;
    for(k = 0; k < LEN; k++){
        a[k] = rand();
        b[k] = rand();
    }

    gettimeofday(&stTime, NULL);

    for(k = 0; k < LEN; k++)
        c[k] = a[k] * b[k];

    gettimeofday(&endTime, NULL);

    FILE* fh = fopen("dump", "w");
    for(k = 0; k < LEN; k++)
        fprintf(fh, "c[%d] = %f\t", k, c[k]);
    fclose(fh);

    double timeE = (double)(endTime.tv_usec + endTime.tv_sec*1000000 - stTime.tv_usec - stTime.tv_sec*1000000);

    printf("Time elapsed: %f\n", timeE);

    return 0;
}

В этом коде я просто инициализирую и умножая два вектора. Результаты сохраняются в векторе c. В основном меня интересует эффект векторизации следующего цикла:

for(k = 0; k < LEN; k++)
    c[k] = a[k] * b[k];

Я скомпилирую код, используя следующие две команды:

1) icc -O2 TestSMID.c -o TestSMID -no-vec -no-simd
2) icc -O2 TestSMID.c -o TestSMID -vec-report2

Я ожидаю увидеть улучшение производительности, поскольку вторая команда успешно векторизовает цикл. Тем не менее, мои исследования показывают, что нет улучшения производительности, когда цикл векторизован.

Возможно, я что-то пропустил, так как я не очень хорошо знаком с этой темой. Поэтому, пожалуйста, дайте мне знать, если что-то не так с моим кодом.

Заранее благодарим за помощь.

PS: Я использую Mac OSX, поэтому нет необходимости выровнять данные, поскольку все выделенные ячейки памяти выравниваются по 16 байт.

Изменить: Я хотел бы поблагодарить всех вас за ваши комментарии и ответы. Я подумал о ответе, предложенном @Mysticial, и есть некоторые дополнительные моменты, которые следует упомянуть здесь. Во-первых, как отметил В. Винска, c[k]=a[k]*b[k] не принимает только один цикл. В дополнение к приращению индекса цикла и сравнению для обеспечения того, чтобы k было меньше, чем LEN, для выполнения операции необходимо выполнить другие действия. Посмотрев на код сборки, сгенерированный компилятором, можно увидеть, что простое умножение требует гораздо больше одного цикла. Обозначенная версия выглядит так:

L_B1.9:                         # Preds L_B1.8
        movq      %r13, %rax                                    #25.5
        andq      $15, %rax                                     #25.5
        testl     %eax, %eax                                    #25.5
        je        L_B1.12       # Prob 50%                      #25.5
                                # LOE rbx r12 r13 r14 r15 eax
L_B1.10:                        # Preds L_B1.9
        testb     $7, %al                                       #25.5
        jne       L_B1.32       # Prob 10%                      #25.5
                                # LOE rbx r12 r13 r14 r15
L_B1.11:                        # Preds L_B1.10
        movsd     (%r14), %xmm0                                 #26.16
        movl      $1, %eax                                      #25.5
        mulsd     (%r15), %xmm0                                 #26.23
        movsd     %xmm0, (%r13)                                 #26.9
                                # LOE rbx r12 r13 r14 r15 eax
L_B1.12:                        # Preds L_B1.11 L_B1.9
        movl      %eax, %edx                                    #25.5
        movl      %eax, %eax                                    #26.23
        negl      %edx                                          #25.5
        andl      $1, %edx                                      #25.5
        negl      %edx                                          #25.5
        addl      $10000000, %edx                               #25.5
        lea       (%r15,%rax,8), %rcx                           #26.23
        testq     $15, %rcx                                     #25.5
        je        L_B1.16       # Prob 60%                      #25.5
                                # LOE rdx rbx r12 r13 r14 r15 eax
L_B1.13:                        # Preds L_B1.12
        movl      %eax, %eax                                    #25.5
                                # LOE rax rdx rbx r12 r13 r14 r15
L_B1.14:                        # Preds L_B1.14 L_B1.13
        movups    (%r15,%rax,8), %xmm0                          #26.23
        movsd     (%r14,%rax,8), %xmm1                          #26.16
        movhpd    8(%r14,%rax,8), %xmm1                         #26.16
        mulpd     %xmm0, %xmm1                                  #26.23
        movntpd   %xmm1, (%r13,%rax,8)                          #26.9
        addq      $2, %rax                                      #25.5
        cmpq      %rdx, %rax                                    #25.5
        jb        L_B1.14       # Prob 99%                      #25.5
        jmp       L_B1.20       # Prob 100%                     #25.5
                                # LOE rax rdx rbx r12 r13 r14 r15
L_B1.16:                        # Preds L_B1.12
        movl      %eax, %eax                                    #25.5
                                # LOE rax rdx rbx r12 r13 r14 r15
L_B1.17:                        # Preds L_B1.17 L_B1.16
        movsd     (%r14,%rax,8), %xmm0                          #26.16
        movhpd    8(%r14,%rax,8), %xmm0                         #26.16
        mulpd     (%r15,%rax,8), %xmm0                          #26.23
        movntpd   %xmm0, (%r13,%rax,8)                          #26.9
        addq      $2, %rax                                      #25.5
        cmpq      %rdx, %rax                                    #25.5
        jb        L_B1.17       # Prob 99%                      #25.5
                                # LOE rax rdx rbx r12 r13 r14 r15
L_B1.18:                        # Preds L_B1.17
        mfence                                                  #25.5
                                # LOE rdx rbx r12 r13 r14 r15
L_B1.19:                        # Preds L_B1.18
        mfence                                                  #25.5
                                # LOE rdx rbx r12 r13 r14 r15
L_B1.20:                        # Preds L_B1.14 L_B1.19 L_B1.32
        cmpq      $10000000, %rdx                               #25.5
        jae       L_B1.24       # Prob 0%                       #25.5
                                # LOE rdx rbx r12 r13 r14 r15
L_B1.22:                        # Preds L_B1.20 L_B1.22
        movsd     (%r14,%rdx,8), %xmm0                          #26.16
        mulsd     (%r15,%rdx,8), %xmm0                          #26.23
        movsd     %xmm0, (%r13,%rdx,8)                          #26.9
        incq      %rdx                                          #25.5
        cmpq      $10000000, %rdx                               #25.5
        jb        L_B1.22       # Prob 99%                      #25.5
                                # LOE rdx rbx r12 r13 r14 r15
L_B1.24:                        # Preds L_B1.22 L_B1.20

И неглавная версия:

L_B1.9:                         # Preds L_B1.8
        xorl      %eax, %eax                                    #25.5
                                # LOE rbx r12 r13 r14 r15 eax
L_B1.10:                        # Preds L_B1.10 L_B1.9
        lea       (%rax,%rax), %edx                             #26.9
        incl      %eax                                          #25.5
        cmpl      $5000000, %eax                                #25.5
        movsd     (%r15,%rdx,8), %xmm0                          #26.16
        movsd     8(%r15,%rdx,8), %xmm1                         #26.16
        mulsd     (%r13,%rdx,8), %xmm0                          #26.23
        mulsd     8(%r13,%rdx,8), %xmm1                         #26.23
        movsd     %xmm0, (%rbx,%rdx,8)                          #26.9
        movsd     %xmm1, 8(%rbx,%rdx,8)                         #26.9
        jb        L_B1.10       # Prob 99%                      #25.5
                                # LOE rbx r12 r13 r14 r15 eax

Кроме того, процессор не загружает только 24 байта. В каждом доступе к памяти загружается полная строка (64 байта). Что еще более важно, так как память, необходимая для a, b и c, смежна, prefetcher определенно поможет много и загрузит следующие блоки заранее. Сказав это, я думаю, что ширина памяти, рассчитанная @Mysticial, слишком пессимистична.

Кроме того, использование SIMD для улучшения производительности программы для очень простого добавления упоминается в Руководстве по векторизации Intel. Таким образом, кажется, что мы должны получить некоторое улучшение производительности для этого очень простого цикла.

Edit2: Еще раз спасибо за ваши комментарии. Кроме того, спасибо @Mysticial пример кода, я наконец увидел эффект SIMD на повышение производительности. Проблема, как упоминалось в Mystical, была пропускная способность памяти. Выбрав небольшой размер для a, b и c, которые вписываются в кеш L1, можно видеть, что SIMD может значительно улучшить производительность. Вот результаты, которые я получил:

icc -O2 -o TestSMIDNoVec -no-vec TestSMID2.c: 17.34 sec

icc -O2 -o TestSMIDVecNoUnroll -vec-report2 TestSMID2.c: 9.33 sec

И разворачивание цикла улучшает производительность еще больше:

icc -O2 -o TestSMIDVecUnroll -vec-report2 TestSMID2.c -unroll=8: 8.6sec

Кроме того, я должен упомянуть, что для моего процессора требуется только один цикл для завершения итерации при компиляции с помощью -O2.

PS: Мой компьютер - это ядро Macbook Pro i5 @2.5 ГГц (двухъядерное)

Ответ 1

Этот оригинальный ответ был верным еще в 2013 году. По состоянию на 2017 год все изменилось настолько, что и вопрос, и ответ устарели.

См. конец этого ответа для обновления 2017 года.

Оригинальный ответ (2013):

Потому что у вас узкое место по пропускной способности памяти.

Хотя векторизация и другие микрооптимизации могут улучшить скорость вычислений, они не могут увеличить скорость вашей памяти.

В вашем примере:

for(k = 0; k < LEN; k++)
    c[k] = a[k] * b[k];

Вы делаете один проход по всей памяти, делая очень мало работы. Это максимизирует пропускную способность вашей памяти.

Поэтому, независимо от того, как он оптимизирован, (векторизованный, развернутый и т.д.), он не будет намного быстрее.

Типичный настольный компьютер 2013 года имеет порядок 10 ГБ/с пропускной способности памяти *.
Ваша петля касается 24 байта/итерации.

Без векторизации современный процессор x64, вероятно, может сделать около 1 итерации цикл *.

Предположим, что вы работаете на частоте 4 ГГц:

(4 * 10^9) * 24 bytes/iteration = 96 GB/s

Это почти 10x вашей пропускной способности памяти - без векторизации.

* Неудивительно, что несколько человек сомневались в цифрах, которые я приводил выше, так как я не дал никаких цитат. Хорошо, что они были на вершине моей головы из опыта. Итак, вот некоторые тесты, чтобы доказать это.

Итерация цикла может работать так же быстро, как 1 цикл/итерация:

Мы можем избавиться от узкого места памяти, если мы уменьшим LEN, чтобы он соответствовал кешу.
(Я тестировал это на С++, поскольку это было проще, но это не имеет значения.)

#include <iostream>
#include <time.h>
using std::cout;
using std::endl;

int main(){
    const int LEN = 256;

    double *a = (double*)malloc(LEN*sizeof(*a));
    double *b = (double*)malloc(LEN*sizeof(*a));
    double *c = (double*)malloc(LEN*sizeof(*a));

    int k;
    for(k = 0; k < LEN; k++){
        a[k] = rand();
        b[k] = rand();
    }

    clock_t time0 = clock();

    for (int i = 0; i < 100000000; i++){
        for(k = 0; k < LEN; k++)
            c[k] = a[k] * b[k];
    }

    clock_t time1 = clock();
    cout << (double)(time1 - time0) / CLOCKS_PER_SEC << endl;
}

Процессор: Intel Core i7 2600K @4.2 ГГц
Компилятор: Visual Studio 2012
Время: 6.55 секунд

В этом тесте я провел 25 600 000 000 итераций всего 6.55 секунд.

6.55 * 4.2 GHz= 27 510 000 000 циклов
27,510,000,000 / 25,600,000,000= 1.074 цикла/итерация

Теперь, если вам интересно, как это можно сделать:

2 загрузки
1 магазин
1 умножить
счетчик приращений
сравнить + ветвь

всего за один цикл...

Это потому, что современные процессоры и компиляторы удивительны.

В то время как каждая из этих операций имеет латентность (особенно умножение), процессор может одновременно выполнять несколько итераций. Моя тестовая машина - это процессор Sandy Bridge, способный выдерживать нагрузки 2x128b, 1x128b store и 1x256b vector FP умножать каждый отдельный цикл. И потенциально другой один или два векторных или целочисленных ops, если нагрузки являются операндами источника памяти для микро-fused uops. (2 загрузки + 1 пропускная способность магазина только при использовании 256-битных AVX-загрузок/хранилищ, в противном случае всего два операционных блока памяти за цикл (не более одного хранилища)).

Глядя на сборку (которую я опускаю для краткости), кажется, что компилятор развернул цикл, тем самым уменьшив накладные расходы. Но ему не удалось его векторизовать.

Ширина памяти составляет порядка 10 ГБ/с:

Самый простой способ проверить это - через memset():

#include <iostream>
#include <time.h>
using std::cout;
using std::endl;

int main(){
    const int LEN = 1 << 30;    //  1GB

    char *a = (char*)calloc(LEN,1);

    clock_t time0 = clock();

    for (int i = 0; i < 100; i++){
        memset(a,0xff,LEN);
    }

    clock_t time1 = clock();
    cout << (double)(time1 - time0) / CLOCKS_PER_SEC << endl;
}

Процессор: Intel Core i7 2600K @4.2 ГГц
Компилятор: Visual Studio 2012
Время: 5.811 секунд

Поэтому для записи на 100 ГБ памяти требуется моя машина 5.811 секунд. Это примерно 17,2 ГБ/с.

И мой процессор находится на более высоком конце. Процессоры Nehalem и Core 2 имеют меньшую пропускную способность памяти.

Обновить март 2017:

По состоянию на 2017 год все усложнилось.

Благодаря DDR4 и четырехканальной памяти уже невозможно, чтобы один поток насыщал пропускную способность памяти. Но проблема пропускной способности не обязательно исчезает. Несмотря на то, что пропускная способность повысилась, процессорные ядра также улучшились - и их больше.

Математически:

Каждое ядро имеет ограничение пропускной способности X.
Основная память имеет ограничение полосы пропускания Y.
В старых системах X > Y.
В современных высокопроизводительных системах X < Y. Но X * (# of cores) > Y.

Назад в 2013 году: Sandy Bridge @4 ГГц + двухканальный DDR3 @1333 МГц

Без векторизации (8-байтная загрузка/сохранение): X = 32 GB/s и Y = ~17 GB/s
Векторизованный SSE * (16-байтная загрузка/сохранение): X = 64 GB/s и Y = ~17 GB/s

Теперь в 2017 году: Haswell-E @4 ГГц + четырехканальный DDR4 @2400 МГц

Без векторизации (8-байтная загрузка/сохранение): X = 32 GB/s и Y = ~70 GB/s
Векторизованный AVX * (32-байтовая загрузка/сохранение): X = 64 GB/s и Y = ~70 GB/s

_{(Для Sandy Bridge и Haswell архитектурные ограничения в кеше ограничивают пропускную способность до 16 байтов/циклов независимо от ширины SIMD.)}

Итак, в настоящее время один поток не всегда сможет насытить пропускную способность памяти. И вам нужно будет векторизовать для достижения этого предела X. Но вы по-прежнему будете использовать ограничение полосы пропускания основной памяти Y с двумя или более потоками.

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

Ответ 2

EDIT: Изменен ответ много. Кроме того, пожалуйста, не обращайте внимания на большинство из того, что я написал ранее, о том, что Мистический ответ не совсем корректен. Хотя, я до сих пор не согласен с тем, что это узкое место в памяти, так как, несмотря на очень широкий спектр тестов, я не видел никаких признаков того, что исходный код связан скоростью памяти. Тем временем он продолжал показывать явные признаки привязки к CPU.

Могут быть много причин. И поскольку причина [s] может быть очень аппаратно-зависимой, я решил, что не должен спекулировать на основе догадок. Просто собираюсь описать эти вещи, с которыми я столкнулся во время более поздних тестов, где я использовал гораздо более точный и надежный метод измерения времени процессора и цикл за цикл в 1000 раз. Я считаю, что эта информация может помочь. Но, пожалуйста, возьмите его с солью, поскольку это зависит от оборудования.

При использовании инструкций из семейства SSE векторизованный код, который я получил, был более чем на 10% быстрее по сравнению с не-векторизованным кодом.
Векторизованный код с использованием SSE-семейства и векторизованного кода с использованием AVX работает более или менее с одинаковой производительностью.
При использовании инструкций AVX не-векторизованный код работал быстрее - на 25% или быстрее, чем все, что я пробовал.
Результаты масштабируются линейно с тактовой частотой процессора во всех случаях.
Результаты не сильно повлияли на часы памяти.
Результаты были в значительной степени затронуты задержкой памяти - гораздо больше, чем часы памяти, но не так сильно, как часы процессора повлияли на результаты.

WRT Мистический пример запуска почти 1 итерации за такт - я не ожидал, что планировщик CPU будет таким эффективным и будет принимать 1 итерацию каждые 1,5-2 такта. Но, к моему удивлению, это не так; Я был не прав, сожалею об этом. Мой собственный процессор работал еще эффективнее - 1.048 циклов/итераций. Поэтому я могу подтвердить, что эта часть Мистического ответа определенно правильна.

Ответ 3

Как уже говорилось, Mystential, ограничения пропускной способности основной памяти являются узким местом для больших буферов. Путь вокруг этого заключается в том, чтобы перепроектировать вашу обработку для работы в кусках, которые вписываются в кеш. (Вместо того, чтобы умножить целые 200 Мбайт удвоений, умножьте всего лишь 128kiB, затем сделайте что-то с этим. Таким образом, код, который использует вывод умножения, найдет его в кэше L2. Обычно L2 256kiB и является приватным для каждого ядра ЦП, о последних разработках Intel.)

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

Если вы это сделаете, убедитесь, что компилятор еще не создает потоковые хранилища (movnt...). Эти записи обходят кэши, чтобы не загрязнять их данными, которые не подходят. Следующее чтение этих данных должно касаться основной памяти.

Ответ 4

На всякий случай, когда [] b [] и c [] сражаются за кеш L2::

#include <string.h> /* for memcpy */

 ...

 gettimeofday(&stTime, NULL);

    for(k = 0; k < LEN; k += 4) {
        double a4[4], b4[4], c4[4];
        memcpy(a4,a+k, sizeof a4);
        memcpy(b4,b+k, sizeof b4);
        c4[0] = a4[0] * b4[0];
        c4[1] = a4[1] * b4[1];
        c4[2] = a4[2] * b4[2];
        c4[3] = a4[3] * b4[3];
        memcpy(c+k,c4, sizeof c4);
        }

    gettimeofday(&endTime, NULL);

Сокращает время работы от 98429.000000 до 67213.000000; разворачивание цикла в 8 раз сводит его к 57157.000000 здесь.