Почему это SIMD-умножение не быстрее, чем умножение без SIMD?

Предположим, что у нас есть функция, которая умножает два массива по 1000000 удваивается каждый. В C/С++ функция выглядит так:

void mul_c(double* a, double* b)
{
    for (int i = 0; i != 1000000; ++i)
    {
        a[i] = a[i] * b[i];
    }
}

Компилятор создает следующую сборку с -O2:

mul_c(double*, double*):
        xor     eax, eax
.L2:
        movsd   xmm0, QWORD PTR [rdi+rax]
        mulsd   xmm0, QWORD PTR [rsi+rax]
        movsd   QWORD PTR [rdi+rax], xmm0
        add     rax, 8
        cmp     rax, 8000000
        jne     .L2
        rep ret

Из вышеприведенной сборки кажется, что компилятор использует инструкции SIMD, но только умножает на одну двойную каждую итерацию. Поэтому я решил написать ту же функцию в встроенной сборке, где я полностью использую регистр xmm0 и умножаю два удвоения за один раз:

void mul_asm(double* a, double* b)
{
    asm volatile
    (
        ".intel_syntax noprefix             \n\t"
        "xor    rax, rax                    \n\t"
        "0:                                 \n\t"
        "movupd xmm0, xmmword ptr [rdi+rax] \n\t"
        "mulpd  xmm0, xmmword ptr [rsi+rax] \n\t"
        "movupd xmmword ptr [rdi+rax], xmm0 \n\t"
        "add    rax, 16                     \n\t"
        "cmp    rax, 8000000                \n\t"
        "jne    0b                          \n\t"
        ".att_syntax noprefix               \n\t"

        : 
        : "D" (a), "S" (b)
        : "memory", "cc"
    );
}

После измерения времени выполнения отдельно для обеих этих функций кажется, что для обоих из них требуется 1 мс для завершения:

> gcc -O2 main.cpp
> ./a.out < input

mul_c: 1 ms
mul_asm: 1 ms

[a lot of doubles...]

Я ожидал, что реализация SIMD будет по крайней мере вдвое быстрее (0 мс), так как есть только половина команд умножения/памяти.

Итак, мой вопрос: Почему реализация SIMD быстрее, чем обычная реализация C/С++, когда реализация SIMD выполняет только половину количества команд умножения/памяти?

Здесь полная программа:

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

void mul_c(double* a, double* b)
{
    for (int i = 0; i != 1000000; ++i)
    {
        a[i] = a[i] * b[i];
    }
}

void mul_asm(double* a, double* b)
{
    asm volatile
    (
        ".intel_syntax noprefix             \n\t"
        "xor    rax, rax                    \n\t"
        "0:                                 \n\t"
        "movupd xmm0, xmmword ptr [rdi+rax] \n\t"
        "mulpd  xmm0, xmmword ptr [rsi+rax] \n\t"
        "movupd xmmword ptr [rdi+rax], xmm0 \n\t"
        "add    rax, 16                     \n\t"
        "cmp    rax, 8000000                \n\t"
        "jne    0b                          \n\t"
        ".att_syntax noprefix               \n\t"

        : 
        : "D" (a), "S" (b)
        : "memory", "cc"
    );
}

int main()
{
    struct timeval t1;
    struct timeval t2;
    unsigned long long time;

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

    for (int i = 0; i != 1000000; ++i)
    {
        double v;
        scanf("%lf", &v);
        a[i] = v;
        b[i] = v;
        c[i] = v;
    }

    gettimeofday(&t1, NULL);
    mul_c(a, b);
    gettimeofday(&t2, NULL);
    time = 1000 * (t2.tv_sec - t1.tv_sec) + (t2.tv_usec - t1.tv_usec) / 1000;
    printf("mul_c: %llu ms\n", time);

    gettimeofday(&t1, NULL);
    mul_asm(b, c);
    gettimeofday(&t2, NULL);
    time = 1000 * (t2.tv_sec - t1.tv_sec) + (t2.tv_usec - t1.tv_usec) / 1000;
    printf("mul_asm: %llu ms\n\n", time);

    for (int i = 0; i != 1000000; ++i)
    {
        printf("%lf\t\t\t%lf\n", a[i], b[i]);
    }

    return 0;
}

Я также попытался использовать все регистры xmm (0-7) и удалить зависимости команд, чтобы получить лучшие параллельные вычисления:

void mul_asm(double* a, double* b)
{
    asm volatile
    (
        ".intel_syntax noprefix                 \n\t"
        "xor    rax, rax                        \n\t"
        "0:                                     \n\t"
        "movupd xmm0, xmmword ptr [rdi+rax]     \n\t"
        "movupd xmm1, xmmword ptr [rdi+rax+16]  \n\t"
        "movupd xmm2, xmmword ptr [rdi+rax+32]  \n\t"
        "movupd xmm3, xmmword ptr [rdi+rax+48]  \n\t"
        "movupd xmm4, xmmword ptr [rdi+rax+64]  \n\t"
        "movupd xmm5, xmmword ptr [rdi+rax+80]  \n\t"
        "movupd xmm6, xmmword ptr [rdi+rax+96]  \n\t"
        "movupd xmm7, xmmword ptr [rdi+rax+112] \n\t"
        "mulpd  xmm0, xmmword ptr [rsi+rax]     \n\t"
        "mulpd  xmm1, xmmword ptr [rsi+rax+16]  \n\t"
        "mulpd  xmm2, xmmword ptr [rsi+rax+32]  \n\t"
        "mulpd  xmm3, xmmword ptr [rsi+rax+48]  \n\t"
        "mulpd  xmm4, xmmword ptr [rsi+rax+64]  \n\t"
        "mulpd  xmm5, xmmword ptr [rsi+rax+80]  \n\t"
        "mulpd  xmm6, xmmword ptr [rsi+rax+96]  \n\t"
        "mulpd  xmm7, xmmword ptr [rsi+rax+112] \n\t"
        "movupd xmmword ptr [rdi+rax], xmm0     \n\t"
        "movupd xmmword ptr [rdi+rax+16], xmm1  \n\t"
        "movupd xmmword ptr [rdi+rax+32], xmm2  \n\t"
        "movupd xmmword ptr [rdi+rax+48], xmm3  \n\t"
        "movupd xmmword ptr [rdi+rax+64], xmm4  \n\t"
        "movupd xmmword ptr [rdi+rax+80], xmm5  \n\t"
        "movupd xmmword ptr [rdi+rax+96], xmm6  \n\t"
        "movupd xmmword ptr [rdi+rax+112], xmm7 \n\t"
        "add    rax, 128                        \n\t"
        "cmp    rax, 8000000                    \n\t"
        "jne    0b                              \n\t"
        ".att_syntax noprefix                   \n\t"

        : 
        : "D" (a), "S" (b)
        : "memory", "cc"
    );
}

Но он все еще работает на 1 мс, с той же скоростью, что и обычная реализация C/С++.

ОБНОВЛЕНИЯ

Как было предложено ответами/комментариями, я реализовал другой способ измерения времени выполнения:

#include <stdio.h>
#include <stdlib.h>

void mul_c(double* a, double* b)
{
    for (int i = 0; i != 1000000; ++i)
    {
        a[i] = a[i] * b[i];
    }
}

void mul_asm(double* a, double* b)
{
    asm volatile
    (
        ".intel_syntax noprefix             \n\t"
        "xor    rax, rax                    \n\t"
        "0:                                 \n\t"
        "movupd xmm0, xmmword ptr [rdi+rax] \n\t"
        "mulpd  xmm0, xmmword ptr [rsi+rax] \n\t"
        "movupd xmmword ptr [rdi+rax], xmm0 \n\t"
        "add    rax, 16                     \n\t"
        "cmp    rax, 8000000                \n\t"
        "jne    0b                          \n\t"
        ".att_syntax noprefix               \n\t"

        : 
        : "D" (a), "S" (b)
        : "memory", "cc"
    );
}

void mul_asm2(double* a, double* b)
{
    asm volatile
    (
        ".intel_syntax noprefix                 \n\t"
        "xor    rax, rax                        \n\t"
        "0:                                     \n\t"
        "movupd xmm0, xmmword ptr [rdi+rax]     \n\t"
        "movupd xmm1, xmmword ptr [rdi+rax+16]  \n\t"
        "movupd xmm2, xmmword ptr [rdi+rax+32]  \n\t"
        "movupd xmm3, xmmword ptr [rdi+rax+48]  \n\t"
        "movupd xmm4, xmmword ptr [rdi+rax+64]  \n\t"
        "movupd xmm5, xmmword ptr [rdi+rax+80]  \n\t"
        "movupd xmm6, xmmword ptr [rdi+rax+96]  \n\t"
        "movupd xmm7, xmmword ptr [rdi+rax+112] \n\t"
        "mulpd  xmm0, xmmword ptr [rsi+rax]     \n\t"
        "mulpd  xmm1, xmmword ptr [rsi+rax+16]  \n\t"
        "mulpd  xmm2, xmmword ptr [rsi+rax+32]  \n\t"
        "mulpd  xmm3, xmmword ptr [rsi+rax+48]  \n\t"
        "mulpd  xmm4, xmmword ptr [rsi+rax+64]  \n\t"
        "mulpd  xmm5, xmmword ptr [rsi+rax+80]  \n\t"
        "mulpd  xmm6, xmmword ptr [rsi+rax+96]  \n\t"
        "mulpd  xmm7, xmmword ptr [rsi+rax+112] \n\t"
        "movupd xmmword ptr [rdi+rax], xmm0     \n\t"
        "movupd xmmword ptr [rdi+rax+16], xmm1  \n\t"
        "movupd xmmword ptr [rdi+rax+32], xmm2  \n\t"
        "movupd xmmword ptr [rdi+rax+48], xmm3  \n\t"
        "movupd xmmword ptr [rdi+rax+64], xmm4  \n\t"
        "movupd xmmword ptr [rdi+rax+80], xmm5  \n\t"
        "movupd xmmword ptr [rdi+rax+96], xmm6  \n\t"
        "movupd xmmword ptr [rdi+rax+112], xmm7 \n\t"
        "add    rax, 128                        \n\t"
        "cmp    rax, 8000000                    \n\t"
        "jne    0b                              \n\t"
        ".att_syntax noprefix                   \n\t"

        : 
        : "D" (a), "S" (b)
        : "memory", "cc"
    );
}

unsigned long timestamp()
{
    unsigned long a;

    asm volatile
    (
        ".intel_syntax noprefix \n\t"
        "xor   rax, rax         \n\t"
        "xor   rdx, rdx         \n\t"
        "RDTSCP                 \n\t"
        "shl   rdx, 32          \n\t"
        "or    rax, rdx         \n\t"
        ".att_syntax noprefix   \n\t"

        : "=a" (a)
        : 
        : "memory", "cc"
    );

    return a;
}

int main()
{
    unsigned long t1;
    unsigned long t2;

    double* a;
    double* b;

    a = (double*)malloc(sizeof(double) * 1000000);
    b = (double*)malloc(sizeof(double) * 1000000);

    for (int i = 0; i != 1000000; ++i)
    {
        double v;
        scanf("%lf", &v);
        a[i] = v;
        b[i] = v;
    }

    t1 = timestamp();
    mul_c(a, b);
    //mul_asm(a, b);
    //mul_asm2(a, b);
    t2 = timestamp();
    printf("mul_c: %lu cycles\n\n", t2 - t1);

    for (int i = 0; i != 1000000; ++i)
    {
        printf("%lf\t\t\t%lf\n", a[i], b[i]);
    }

    return 0;
}

Когда я запускаю программу с этим измерением, я получаю этот результат:

mul_c:    ~2163971628 cycles
mul_asm:  ~2532045184 cycles
mul_asm2: ~5230488    cycles <-- what???

Здесь стоит обратить внимание на две вещи: прежде всего, количество циклов варьируется LOT, и я предполагаю, что из-за операционной системы, позволяющей другим процессам работать между ними. Есть ли способ предотвратить или считать только циклы во время моей программы? Кроме того, mul_asm2 производит идентичный вывод по сравнению с двумя другими, но он намного быстрее, как?

Я попробовал программу Z boson в своей системе вместе с моими 2 реализациями и получил следующий результат:

> g++ -O2 -fopenmp main.cpp
> ./a.out
mul         time 1.33, 18.08 GB/s
mul_SSE     time 1.13, 21.24 GB/s
mul_SSE_NT  time 1.51, 15.88 GB/s
mul_SSE_OMP time 0.79, 30.28 GB/s
mul_SSE_v2  time 1.12, 21.49 GB/s
mul_v2      time 1.26, 18.99 GB/s
mul_asm     time 1.12, 21.50 GB/s
mul_asm2    time 1.09, 22.08 GB/s

Ответы

Ответ 1

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

Ваша операция связана с пропускной способностью памяти. Это означает, что процессор тратит большую часть своего времени, ожидая медленного чтения и записи в памяти. Отличное объяснение этого можно найти здесь: Почему векторизация цикла не повышает производительность.

Однако я должен не согласиться с одним утверждением в этом ответе.

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

Фактически, векторизация , разворачивание, и несколько потоков могут значительно увеличить пропускную способность даже в операциях с привязкой к пропускной способности памяти. Причина в том, что трудно получить максимальную пропускную способность памяти. Хорошее объяснение этого можно найти здесь: fooobar.com/questions/9124/....

Остальная часть моего ответа покажет, как векторизация и несколько потоков могут приблизиться к максимальной пропускной способности памяти.

Моя тестовая система: Ubuntu 16.10, Skylake ([email protected]), 32 ГБ оперативной памяти, двухканальный DDR4 @2400 ГГц. Максимальная пропускная способность моей системы составляет 38,4 ГБ/с.

Из приведенного ниже кода я создаю следующие таблицы. Я устанавливаю количество потоков, используя OMP_NUM_THREADS, например. export OMP_NUM_THREADS=4. Эффективность bandwidth/max_bandwidth.

-O2 -march=native -fopenmp
Threads Efficiency
1       59.2%
2       76.6%
4       74.3%
8       70.7%

-O2 -march=native -fopenmp -funroll-loops
1       55.8%
2       76.5%
4       72.1%
8       72.2%

-O3 -march=native -fopenmp
1       63.9%
2       74.6%
4       63.9%
8       63.2%

-O3 -march=native -fopenmp -mprefer-avx128
1       67.8%
2       76.0%
4       63.9%
8       63.2%

-O3 -march=native -fopenmp -mprefer-avx128 -funroll-loops
1       68.8%
2       73.9%
4       69.0%
8       66.8%

После нескольких итераций работы из-за неопределенностей в измерениях я сделал следующие выводы:

  • однопоточные скалярные операции получают более 50% пропускной способности.
  • две потоковые скалярные операции получают максимальную пропускную способность.
  • однопоточные векторные операции быстрее, чем однопоточные скалярные операции.
  • однопоточные SSE-операции быстрее, чем однопоточные операции AVX.
  • развернуть не полезно.
  • разворачивание однопоточных операций выполняется медленнее, чем без разворота.
  • больше потоков, чем ядер (Hyper-Threading) дает более низкую пропускную способность.

Решение, обеспечивающее наилучшую пропускную способность, представляет собой скалярные операции с двумя потоками.

Код, который я использовал для сравнения:

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <omp.h>

#define N 10000000
#define R 100

void mul(double *a, double *b) {
  #pragma omp parallel for
  for (int i = 0; i<N; i++) a[i] *= b[i];
}

int main() {
  double maxbw = 2.4*2*8; // 2.4GHz * 2-channels * 64-bits * 1-byte/8-bits 
  double mem = 3*sizeof(double)*N*R*1E-9; // GB

  double *a = (double*)malloc(sizeof *a * N);
  double *b = (double*)malloc(sizeof *b * N);

  //due to copy-on-write b must be initialized to get the correct bandwidth
  //also, GCC will convert malloc + memset(0) to calloc so use memset(1)
  memset(b, 1, sizeof *b * N);

  double dtime = -omp_get_wtime();
  for(int i=0; i<R; i++) mul(a,b);
  dtime += omp_get_wtime();
  printf("%.2f s, %.1f GB/s, %.1f%%\n", dtime, mem/dtime, 100*mem/dtime/maxbw);

  free(a), free(b);
}

Старое решение с временной ошибкой

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

Одно внутреннее решение для встроенного сборочного подхода просто:

void mul_SSE(double*  a, double*  b) {
  for (int i = 0; i<N/2; i++) 
      _mm_store_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
}

Позвольте мне определить некоторый тестовый код

#include <x86intrin.h>
#include <string.h>
#include <stdio.h>
#include <x86intrin.h>
#include <omp.h>

#define N 1000000
#define R 1000

typedef __attribute__(( aligned(32)))  double aligned_double;
void  (*fp)(aligned_double *a, aligned_double *b);

void mul(aligned_double* __restrict a, aligned_double* __restrict b) {
  for (int i = 0; i<N; i++) a[i] *= b[i];
}

void mul_SSE(double*  a, double*  b) {
  for (int i = 0; i<N/2; i++) _mm_store_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
}

void mul_SSE_NT(double*  a, double*  b) {
  for (int i = 0; i<N/2; i++) _mm_stream_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
}

void mul_SSE_OMP(double*  a, double*  b) {
  #pragma omp parallel for
  for (int i = 0; i<N; i++) a[i] *= b[i];
}

void test(aligned_double *a, aligned_double *b, const char *name) {
  double dtime;
  const double mem = 3*sizeof(double)*N*R/1024/1024/1024;
  const double maxbw = 34.1;
  dtime = -omp_get_wtime();
  for(int i=0; i<R; i++) fp(a,b);
  dtime += omp_get_wtime();
  printf("%s \t time %.2f s, %.1f GB/s, efficency %.1f%%\n", name, dtime, mem/dtime, 100*mem/dtime/maxbw);
}

int main() {
  double *a = (double*)_mm_malloc(sizeof *a * N, 32);
  double *b = (double*)_mm_malloc(sizeof *b * N, 32);

  //b must be initialized to get the correct bandwidth!!!
  memset(a, 1, sizeof *a * N);
  memset(b, 1, sizeof *a * N);

  fp = mul,         test(a,b, "mul        ");
  fp = mul_SSE,     test(a,b, "mul_SSE    ");
  fp = mul_SSE_NT,  test(a,b, "mul_SSE_NT ");
  fp = mul_SSE_OMP, test(a,b, "mul_SSE_OMP");

  _mm_free(a), _mm_free(b);
}

Теперь первый тест

g++ -O2 -fopenmp test.cpp
./a.out
mul              time 1.67 s, 13.1 GB/s, efficiency 38.5%
mul_SSE          time 1.00 s, 21.9 GB/s, efficiency 64.3%
mul_SSE_NT       time 1.05 s, 20.9 GB/s, efficiency 61.4%
mul_SSE_OMP      time 0.74 s, 29.7 GB/s, efficiency 87.0%

Итак, с -O2, который не векторизовать циклы, мы видим, что внутренняя версия SSE намного быстрее, чем простое решение C mul. efficiency = bandwith_measured/max_bandwidth, где max составляет 34,1 ГБ/с для моей системы.

Второй тест

g++ -O3 -fopenmp test.cpp
./a.out
mul              time 1.05 s, 20.9 GB/s, efficiency 61.2%
mul_SSE          time 0.99 s, 22.3 GB/s, efficiency 65.3%
mul_SSE_NT       time 1.01 s, 21.7 GB/s, efficiency 63.7%
mul_SSE_OMP      time 0.68 s, 32.5 GB/s, efficiency 95.2%

С -O3 векторизация цикла, и внутренняя функция не предлагает практически никаких преимуществ.

Третий тест

g++ -O3 -fopenmp -funroll-loops test.cpp
./a.out
mul              time 0.85 s, 25.9 GB/s, efficency 76.1%
mul_SSE          time 0.84 s, 26.2 GB/s, efficency 76.7%
mul_SSE_NT       time 1.06 s, 20.8 GB/s, efficency 61.0%
mul_SSE_OMP      time 0.76 s, 29.0 GB/s, efficency 85.0%

С помощью -funroll-loops GCC разворачивает контуры восемь раз, и мы видим значительное улучшение, за исключением нерезидентного решения магазина, а не реального преимущества для решения OpenMP.

Перед разворачиванием цикла сборка для mul wiht -O3 равна

    xor     eax, eax
.L2:
    movupd  xmm0, XMMWORD PTR [rsi+rax]
    mulpd   xmm0, XMMWORD PTR [rdi+rax]
    movaps  XMMWORD PTR [rdi+rax], xmm0
    add     rax, 16
    cmp     rax, 8000000
    jne     .L2
    rep ret

С -O3 -funroll-loops сборка для mul:

   xor     eax, eax
.L2:
    movupd  xmm0, XMMWORD PTR [rsi+rax]
    movupd  xmm1, XMMWORD PTR [rsi+16+rax]
    mulpd   xmm0, XMMWORD PTR [rdi+rax]
    movupd  xmm2, XMMWORD PTR [rsi+32+rax]
    mulpd   xmm1, XMMWORD PTR [rdi+16+rax]
    movupd  xmm3, XMMWORD PTR [rsi+48+rax]
    mulpd   xmm2, XMMWORD PTR [rdi+32+rax]
    movupd  xmm4, XMMWORD PTR [rsi+64+rax]
    mulpd   xmm3, XMMWORD PTR [rdi+48+rax]
    movupd  xmm5, XMMWORD PTR [rsi+80+rax]
    mulpd   xmm4, XMMWORD PTR [rdi+64+rax]
    movupd  xmm6, XMMWORD PTR [rsi+96+rax]
    mulpd   xmm5, XMMWORD PTR [rdi+80+rax]
    movupd  xmm7, XMMWORD PTR [rsi+112+rax]
    mulpd   xmm6, XMMWORD PTR [rdi+96+rax]
    movaps  XMMWORD PTR [rdi+rax], xmm0
    mulpd   xmm7, XMMWORD PTR [rdi+112+rax]
    movaps  XMMWORD PTR [rdi+16+rax], xmm1
    movaps  XMMWORD PTR [rdi+32+rax], xmm2
    movaps  XMMWORD PTR [rdi+48+rax], xmm3
    movaps  XMMWORD PTR [rdi+64+rax], xmm4
    movaps  XMMWORD PTR [rdi+80+rax], xmm5
    movaps  XMMWORD PTR [rdi+96+rax], xmm6
    movaps  XMMWORD PTR [rdi+112+rax], xmm7
    sub     rax, -128
    cmp     rax, 8000000
    jne     .L2
    rep ret

Четвертый тест

g++ -O3 -fopenmp -mavx test.cpp
./a.out
mul              time 0.87 s, 25.3 GB/s, efficiency 74.3%
mul_SSE          time 0.88 s, 24.9 GB/s, efficiency 73.0%
mul_SSE_NT       time 1.07 s, 20.6 GB/s, efficiency 60.5%
mul_SSE_OMP      time 0.76 s, 29.0 GB/s, efficiency 85.2%

Теперь неотрицательная функция является самой быстрой (исключая версию OpenMP).

Таким образом, в этом случае нет причин использовать встроенные или встроенные сборки, потому что мы можем получить лучшую производительность с соответствующими параметрами компилятора (например, -O3, -funroll-loops, -mavx).

Система тестирования: Ubuntu 16.10, Skylake ([email protected]), оперативная память 32 ГБ. Максимальная пропускная способность памяти (34,1 ГБ/с) https://ark.intel.com/products/88967/Intel-Core-i7-6700HQ-Processor-6M-Cache-up-to-3_50-GHz

Вот еще одно решение, заслуживающее рассмотрения. Инструкция cmp не нужна, если мы рассчитываем от -N до нуля и получаем доступ к массивам как N+i. GCC должен был зафиксировать это давным-давно. Он устраняет одну инструкцию (хотя из-за макро-op слияния cmp и jmp часто считаются одним микрооператором).

void mul_SSE_v2(double*  a, double*  b) {
  for (ptrdiff_t i = -N; i<0; i+=2)
    _mm_store_pd(&a[N + i], _mm_mul_pd(_mm_load_pd(&a[N + i]),_mm_load_pd(&b[N + i])));

Сборка с -O3

mul_SSE_v2(double*, double*):
    mov     rax, -1000000
.L9:
    movapd  xmm0, XMMWORD PTR [rdi+8000000+rax*8]
    mulpd   xmm0, XMMWORD PTR [rsi+8000000+rax*8]
    movaps  XMMWORD PTR [rdi+8000000+rax*8], xmm0
    add     rax, 2
    jne     .L9
    rep ret
}

Эта оптимизация может быть полезна только при использовании массивов, например. кэш L1, то есть не считывающий из основной памяти.

Наконец-то я нашел способ получить решение простой C, чтобы не сгенерировать инструкцию cmp.

void mul_v2(aligned_double* __restrict a, aligned_double* __restrict b) {
  for (int i = -N; i<0; i++) a[i] *= b[i];
}

И затем вызовите функцию из отдельного объектного файла, такого как mul_v2(&a[N],&b[N]), поэтому это, пожалуй, лучшее решение. Однако, если вы вызываете функцию из того же объектного файла (единицы перевода), как тот, который он определил в GCC, снова генерирует команду cmp.

Кроме того,

void mul_v3(aligned_double* __restrict a, aligned_double* __restrict b) {
  for (int i = -N; i<0; i++) a[N+i] *= b[N+i];
}

все еще генерирует инструкцию cmp и генерирует ту же самую сборку, что и функция mul.

Функция mul_SSE_NT глупа. Он использует не временные хранилища, которые полезны только при записи в память, но поскольку функция читает и записывает в один и тот же адрес, невременные хранилища не только бесполезны, они дают более низкие результаты.

Предыдущие версии этого ответа получили неправильную пропускную способность. Причина была в том, что массивы не были инициализированы.

Ответ 2

Ваш код asm действительно в порядке. Что не так, как вы его измеряете. Как я указал в комментариях, вы должны:

a) использовать больше итераций - 1 миллион ничего для современного процессора

b) использовать HPT для измерения

c) использовать RDTSC или RDTSCP для подсчета реальных часов процессора

Кроме того, почему вы боитесь -O3 выбрать? Не забудьте создать код для вашей платформы, поэтому используйте -march = native. Если ваш процессор поддерживает компилятор AVX или AVX2, у вас будет возможность создать еще лучший код.

Следующее - дайте компилятору некоторые подсказки об псевдониме и присваивании, если вы знаете код.

Вот моя версия вашего mul_c - да, это GCC, но вы показали, что использовали GCC

void mul_c(double* restrict a, double* restrict b)
{
   a = __builtin_assume_aligned (a, 16);
   b = __builtin_assume_aligned (b, 16);

    for (int i = 0; i != 1000000; ++i)
    {
        a[i] = a[i] * b[i];
    }
}

Он будет производить:

mul_c(double*, double*):
        xor     eax, eax
.L2:
        movapd  xmm0, XMMWORD PTR [rdi+rax]
        mulpd   xmm0, XMMWORD PTR [rsi+rax]
        movaps  XMMWORD PTR [rdi+rax], xmm0
        add     rax, 16
        cmp     rax, 8000000
        jne     .L2
        rep ret

Если у вас есть AVX2 и убедитесь, что данные выровнены по 32 байта, он станет

mul_c(double*, double*):
        xor     eax, eax
.L2:
        vmovapd ymm0, YMMWORD PTR [rdi+rax]
        vmulpd  ymm0, ymm0, YMMWORD PTR [rsi+rax]
        vmovapd YMMWORD PTR [rdi+rax], ymm0
        add     rax, 32
        cmp     rax, 8000000
        jne     .L2
        vzeroupper
        ret

Так что нет необходимости в asm, если компилятор может сделать это для вас;)

Ответ 3

Я хочу добавить еще одну точку зрения на проблему. Инструкции SIMD дают большую производительность, если ограничений на привязку к памяти нет. Но в текущем примере слишком много операций по загрузке и хранению памяти и слишком мало вычислений ЦП. Таким образом, процессор успевает обрабатывать входящие данные без использования SIMD. Если вы используете данные другого типа (например, 32-битный float) или более сложный алгоритм, пропускная способность памяти не будет ограничивать производительность ЦП, а использование SIMD даст больше преимуществ.