Как стать автором
Обновить
484.21
YADRO
Тут про железо и инженерную культуру

Как оптимизировать код на С для x86-процессоров: подсистема кэша и памяти, инструкции AVX-512

Уровень сложностиСредний
Время на прочтение12 мин
Количество просмотров12K

Меня зовут Андрей Бакшаев, я ведущий инженер-программист в YADRO. Моя команда занимается разработкой и оптимизацией математических библиотек под архитектуру x86. До этого я 15 лет работал в Intel. Значительная часть моих задач заключалась в том, чтобы реализовывать некоторые алгоритмы обработки изображений и сигналов в довольно известной математической библиотеке IPP, максимально эффективно используя возможности процессоров. Я также исследовал производительность этих алгоритмов в процессорах на ранней стадии проектирования. 

В статье я поделюсь своим опытом оптимизации низкоуровневого кода на языке C. Рассмотрим подсистему кэша и памяти процессоров и новые инструкции AVX-512. Разберем пример ускорения копирования байтового массива данных и посмотрим, как векторизованный код позволяет сократить время работы широко используемого алгоритма замены байтов по таблице с 619 до 34 мс, то есть примерно в 18 раз.

Подсистема кэша и памяти

Представим, что перед нами стоит типовая задача — скопировать массив байтов. Начнем с обычного кода на С и постараемся постепенно модифицировать его, повышая производительность.

Скалярный и векторизованный код

Реализуем простую функцию copy

// скалярный С
copy(char* src, char* dst, int  len)
{
  for(i=0;i<len;i++){
    dst[i] = src[i];
   }
}

Чем отличается скалярный код от векторизованного? Исполняемые арифметические инструкции в процессоре можно поделить на две категории:

  1. Скалярные — инструкция выполняет действие, например, сложение или умножение, над одним элементом. К этой категории относятся большинство инструкций с регистрами общего назначения в x86.

  2. Векторные — инструкция загружает из памяти или выполняет какие-то арифметические инструкции сразу с N элементами. N зависит от типа данных: char, short, float и так далее. Смысл в том, что обрабатывается целый вектор, то есть несколько элементов одновременно. За счет этого и получается ускорение вычислений.

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

Первое векторное расширение, появившееся в x86, называлось MMX (Multimedia Extension). Регистры были 64 бит, и одновременно обрабатывалось по 64/8=8 байт. Далее появилось расширение SSE (Streaming SIMD Extensions) с векторами по 128 бит (16 байт). На данный момент актуальны расширения AVX2 и AVX512 с векторами по 256 и 512 бит. 

Посмотрим на примере, как выглядит векторизованный код:

// векторизованный avx2
#include "immintrin.h"
copy(char* src, char* dst, int  len)
{
 for(i=0;i<len;i+=32){
  __m256i x0 =_mm256_loadu_si256(src+i);
  _mm256_storeu_si256(dst+i, x0);
 }
}

Здесь:

  • “immintrin.h” — файл с интринсиками языка C.

  • x0 — объявление и использование векторного регистра. 

  • __m256i x0 — переменная, соответствующая 256-битному векторному регистру.  

  • _mm256_loadu_si256 — интринсик загрузки из памяти в 256-битный регистр.

Мы подключаем векторное расширение и немного модифицируем наш скалярный код функции copy. В цикле i++ заменяем на i+=32, поскольку теперь за одну итерацию одной командой loadu мы загружаем 32 байта из памяти. Команда storeu записывает данные в память. Такой код называется AVX2 векторизованным, поскольку он использует 256-битные регистры расширения AVX2. 

Можно написать и AVX512 векторизованный код, который использует 512-битный регистр. Это позволит брать на одной итерации в 2 раза больше данных, по 64 байта:

// векторизованный avx512 
#include "immintrin.h"
copy(char* src, char* dst, int  len)
{
  for(i=0;i<len;i+=64){
   __m512i x0=_mm512_loadu_si512(src+i);
   _mm512_storeu_si512(dst+i, x0);
  }
}

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

Сравнение производительности кода

Мы написали три разных варианта кода для копирования байтового массива данных: скалярный С и векторизованные AVX2, AVX512. Давайте сравним их между собой по производительности. 

Поскольку входной и выходной векторы имеют одинаковую длину, в кэш будет попадать суммарный объем данных, равный удвоенной длине, то есть 2 × длина вектора. Ее и нужно учитывать при изучении эффектов, связанных с кэшем. Посмотрим на график ниже. В нем используется единица измерения CPE (clock per element). Такты, затраченные на копирование вектора, делятся на его длину: CPE = clocks / length. Чем меньше CPE, тем быстрее.

По вертикали — такты на элемент вектора. По горизонтали — длина копируемых данных.
По вертикали — такты на элемент вектора. По горизонтали — длина копируемых данных.

Посмотрим детально, что происходит в L1-, L2-, L3-кэше и памяти:

Данные в L1. Объем L1 в нашем случае — 48 Кб. Данные помещаются сюда до длины ~48 Кб/2 = 24 Кб. AVX512-код обгоняет AVX2 почти в два раза.  

Данные в L2. Объем данных превышает размер L1, но не превышает размер L2 примерно до длины 1,25 Мб/2 = 625 Кб. Здесь производительность сравнялась.

Данные в L3. До длины 4096 Кб данные помещаются в L3. Производительность всех трех версий кода снова на равных.

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

Выравнивание массива данных 

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

Давайте измерим все смещения входного массива относительно такого адреса. Сначала выделим выровненную на 64 байта память с достаточным запасом. А далее будем измерять время копирования в таком цикле, смещаясь на 1 байт.

// измеряем все смещения 0...63 входного массива
src = _aligned_malloc( len+63 , 64);
dst = _aligned_malloc( len    , 64);
for(off=0;off<=64;off++){
   copy(src+off, dst, len);
}

Получаем такой график:

Время копирования элемента на невыровненных данных входного массива 0...63 байта.
Время копирования элемента на невыровненных данных входного массива 0...63 байта.

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

Почему так? Для AVX-512 любое чтение по адресу, не кратному 64-м, затрагивает две кэш-линии (cache split). Это происходит потому, что длина кэш-линии — 64 байта (512 бит) и регистр в нашем AVX512-коде — 64 байта. А читать две кэш-линии дольше, чем одну.

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

Non-temporal инструкции

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

// обычная запись _mm256_storeu_si256
copy(char* src, char* dst, int  len)
{
	for(int i=0;i<len;i+=32){
    	__m256i 
x0=_mm256_loadu_si256(src+i);
        _mm256_storeu_si256(dst+i, x0);
    }
}
// non-temporal запись _mm256_stream_si256 
copy(char* src, char* dst, int  len)
{
	for(int i=0;i<len;i+=32){
 	__m256i 
x0=_mm256_loadu_si256(src+i);
    _mm256_stream_si256(dst+i, x0);
    }
}

Важно отметить, что инструкция _mm256_stream_si256 работает только с выровненной памятью.

Если мы измерим производительность обеих версий кода, то увидим следующий график:

Время записи элементов в зависимости от длины вектора данных. Синий — обычное copy (например, AVX-512), голубой — non-temporal код.
Время записи элементов в зависимости от длины вектора данных. Синий — обычное copy (например, AVX-512), голубой — non-temporal код.

Объем данных у нас больше, чем помещается в память. Non-temporal код значительно лучше работает в ней, но в кэше работает медленнее, чем код «обычной» записи. 

Совет: если вы точно знаете объем данных, к примеру он заведомо гораздо больше всех кэшей в CPU, то можете переключаться между storeu и stream. То есть записывать через кэш или же использовать non-temporal запись.

Рекомендую следующий алгоритм использования:

if( data_volume < cpu_cache_size )
   	storeu
 else {
   	stream

Сравнение производительности с функцией memcpy

Мы можем также вызвать обычную системную memcpy.  В реализации функции memcpy могут использоваться все перечисленные выше способы, а также специальные rep mov инструкции, которые могут неявно улучшаться в новых моделях процессоров.

Посмотрим на график:

Время записи элементов в зависимости от длины вектора данных.
Время записи элементов в зависимости от длины вектора данных.

График производительности memcpy лежит посередине между векторизованным AVX512-кодом и non-temporal. Следовательно, memcpy тоже хорошо оптимизирована, и можно использовать ее.

Влияние каналов памяти на копирование

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

Если мы возьмем какую-нибудь функцию копирования, то теоретически производительность этой функции должна быть 1/N, где N — количество потоков данных. Однако рано или поздно потоки данных упираются в количество каналов в памяти.  

Такие memory bounded функции (то есть функции с небольшим количеством вычислений по сравнению с загрузками) перестают ускоряться на N потоках гораздо раньше, чем число доступных ядер в процессоре. Например, число ядер на Xeon достигает 24, 48 и более на один процессор, а каналов в памяти гораздо меньше — всего 6 или 8. Для функций вида «мало вычислений, много данных» пропускная способность памяти является ограничителем. 

Покажу для примера график соотношения времени к числу потоков для процессора Skylake:

Ближайшее будущее: High Bandwidth Memory

Мы выяснили, что память до сих пор является ограничителем производительности. Чтобы решить эту проблему, индустрия рассматривает новый стандарт памяти — High Bandwidth Memory (HBM). Сейчас ширина шины памяти составляет 64 байта, а на HBM она достигнет 1024 байт. 

Источник: Forbes.com.
Источник: Forbes.com.

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

  • Wccftech: Intel Sapphire Rapids Xeon Platinum 8472C HBM на 32% быстрее, чем без HBM 8480H.

  • Ixbt.com: новый стандарт высокоскоростной памяти High Bandwidth Memory.

Краткий вывод о системе кэша и памяти

Чтобы оптимизировать код копирования байтового массива данных, можно:

  1. Использовать векторные инструкции.

  2. Выровнять все массивы на 64 байта.

  3. Использовать non-temporal запись.

  4. Использовать функцию memcpy для сравнения.

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

Использование инструкций AVX-512

AVX-512 — это расширение архитектуры x86, которое позволяет использовать векторные инструкции с длиной 512 бит. Давайте посмотрим, чем такие инструкции могут быть полезны при оптимизации.

Табличная замена

Частый прием оптимизации — это замена одного байта другим — например, для:

  • подсчета числа ненулевых бит в байте,

  • поиска первого/последнего ненулевого бита/байта,

  • перестановки бит.

То есть вместо сложных вычислений можно просто встроить таблицу замены в код программы для всех 256 возможных значений байта.

Допустим, у нас есть какие-нибудь вычисления и нам нужно узнать, какой первый ненулевой бит в байте. Алгоритм перестановки младшей и старшей половинок байта выглядит так: 

Байт 56 будет заменен на 65, 34 — на 43.
Байт 56 будет заменен на 65, 34 — на 43.

Как реализовать замену байт? Каждый байт имеет всего 2^8 = 256 возможных значений. Создадим специальную lookup-таблицу (LUT). Она позволит получить те самые 256 значений: 

// создание таблицы
char tbl[256];
make_tbl()
{
   for(i=0;i<256;i++){
	 tbl[i]= 
(i<<4)|(i>>4);
   }
}

И теперь заменим входной массив нашей таблицы:

// замена по таблице
lut_c(uchar* src, uchar* dst, int len)
{
  for(i=0;i<len;i++)
     dst[i] = tbl[ src[i] ];
}

Здесь:

  • dst[i] — это значение в таблице, 

  • src[i] — входной элемент, который будет индексом в таблице. 

Этот код и называется «табличная замена», или lookup table (LUT).

Байтовая перестановка 16 байт

Много лет назад в процессоре Core 2 Duo появилась инструкция _mm_shuffle_epi8, которая позволяет реализовать LUT для 16 значений байт. Вот как выглядит ее применение:

__m128i _mm_shuffle_epi8 (__m128i a, __m128i b)

Мы загружаем в векторный регистр a значения, в b загружаются позиции, из которых нужно взять элементы в a. В с получаем результат:

Если мы используем такую инструкцию, то можем заменить 16 байт. Но этого недостаточно для полноценной байтовой замены.

Байтовая перестановка 128 байт

В наборе VBMI (Vector Bit Manipulation Instructions) появилась новая инструкция байтовой перестановки. Она использует два регистра по 512 бит, то есть два по 64 байта, и сквозную индексацию:

__m512i _mm512_permutex2var_epi8 (__m512i a, __m512i idx, __m512i b)

Здесь:

  • idx — индекс.

  • a, b — регистры. 

Такая инструкция позволяет реализовать LUT для 128 байт — от 0 до 127: 

Например, если индекс 2, то мы загружаем из регистра 0xAA.
Например, если индекс 2, то мы загружаем из регистра 0xAA.

Если нам нужно заменить всего 128 значений, мы можем загрузить посчитанную таблицу и вызвать нашу инструкцию байтовой замены:

// LUT на 128 значений
uchar tbl[128]
...
t0 = _mm512_load(tbl     )
t1 = _mm512_load(tbl+64)
x  = _mm512_permutex2var_epi8 (t0, x, t1)

Такой код позволяет правильно заменить по таблице все байты в регистре x, при условии что все x[i]<128. Старший бит в инструкции не используется. 

Байтовая перестановка 256 байт 

Как тогда заменить всю таблицу длиной 256 байт?

// LUT на все 256 значений
t0 = mm512_load (tbl       )
t1 = mm512_load (tbl + 64)
t2 = mm512_load (tbl +128)
t3 = mm512_load  (tbl +192)
m  = mm512_cmpgt   (x, 127)
a  = mm512_permute (t0, x, t1)
b  = mm512_permute (t2, x, t3)
x  = _mm512_blend (m, a, b);

В t0, t1, t2, t3 мы полностью загружаем всю таблицу, все 256 значений. Далее в битовой маске m сравниваем векторные значения с числом 127 и делаем две перестановки в зависимости от диапазона, где лежит значение. Одна из них будет правильной. Инструкция blend позволяет получить правильное значение в x, используя битовую маску в m

Если перевести на язык математики, то:

m:      m[i]=0,  если x[i] <= 127

          m[i]=1,  если x[i] >  127, i=0..64

t0 и t1:  корректная замена для x[i]  <= 127

t2 и t3:                                   для x[i]  >   127

blend:   x[i]=a[i] если m[i]=0

             x[i]=b[i] если m[i]=1

В x — корректная замена байта на байт для всех 256 возможных значений.

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

Сравнение производительности кода

Вот наш изначальный код табличной замены на С:

// С-код
char tbl[256];
make_tbl(){
 for(i=0;i<256;i++){
   tbl[i]= 
(i<<4)|(i>>4);
 }
}
lut_c(uchar* src, uchar* dst, int len){
  for(i=0;i<len;i++)
     dst[i] = tbl[src[i]];
}

И вот как он стал выглядеть с применением инструкций AVX-512:

// векторизованный avx512 
lut_vec(uchar* src, uchar* dst, int len){
 __m512i x0, x1, x2, x3, t0, t1, t2, t3;
 __mmask64 m;
 t0 = _mm512_loadu_si512(tbl+0*64);
 t1 = _mm512_loadu_si512(tbl+1*64);
 t2 = _mm512_loadu_si512(tbl+2*64);
 t3 = _mm512_loadu_si512(tbl+3*64);
 for(i=0;i<(len&(~63));i+=64){
  x0=_mm512_loadu_si512(src+i);
  x1=_mm512_permutex2var_epi8 (t0, x0, t1);
  x2=_mm512_permutex2var_epi8 (t2, x0, t3);
  m =_mm512_cmpgt_epu8_mask (x0, _mm512_set1_epi8(127));
  x3=_mm512_mask_blend_epi8 (m, x1, x2);
  _mm512_storeu_si512(dst+i,x3);
 }
 for(;i<len;i++)
    dst[i] = tbl[src[i]];
}

Запускаем на ноутбуке (Сore i5-1135G7, 2.4 ГГц, 1024 байта) и получаем следующие результаты: 

  1. C-код: 619 мс.

  2. AVX512-код: 34 мс.

Ускорение:  619 мс / 34 мс = 18 раз.

Масочные регистры в AVX-512 

Все массивы данных обрабатываются по 16 элементов. Но что если нам нужно только 9? Для решения таких задач можно использовать масочные регистры, которые появились в AVX-512:

__m512 _mm512_mask_add_ps (__m512 src, __mmask16 k, __m512 a, __m512 b)

Здесь:

  • __m512 src, __m512 a, __m512 b — это регистры.

  • __mmask16 k — это битовая маска.

Вот схема того, как работают инструкции с масками:

У нас есть регистр a и регистр b. В каждом из них — свои значения. Мы выполняем операцию сложения add. Посмотрим для примера на второй столбец справа. Мы складываем 1.1 и 2.1, получаем в промежуточном итоге операции tmp число 3.2. 

Мы использовали бы значения из tmp, если бы у нас не было масочного регистра k. Но он есть. Битовая маска накладывается на значение промежуточного результата операции. 1 означает, что полученный результат — тот, что нам нужен. Но если в масочном регистре 0, то мы берем не результат сложения, а берем значение из src. В нашем случае правильным результатом будет не 3.2, а 5.0. 

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

Сравнение производительности кода

Возьмем наш AVX512-код из прошлого примера:

//avx512 + «С код остатка» 
 t0 = _mm512_loadu_si512(tbl)
 t1 = _mm512_loadu_si512(tbl+64)
 t2 = _mm512_loadu_si512(tbl+128)
 t3 = _mm512_loadu_si512(tbl+192)
 for(i=0;i<len;i+=64){
  x0=_mm512_loadu_si512(src+i)
  x1=_mm512_permutex2var_epi8 (t0, x0, t1)
  x2=_mm512_permutex2var_epi8 (t2, x0, t3)
  m =_mm512_cmpgt_epu8_mask (x0, 127)
  x3=_mm512_mask_blend_epi8 (m, x1, x2)
  _mm512_storeu_si512(dst+i,x3)
 }
 for(;i<len;i++)
    dst[i] = tbl[src[i]];
}

Представим, что пользователь подал вектор длиной 1087 элементов. Наши итерации в AVX-части — по 64 байта, следовательно: 1087 = 16 × 64 + 63.

Мы выполним основные 16 итераций по 64 байта, а «хвост», то есть последние элементы, мы будем обрабатывать за оставшиеся 63 итерации в цикле на С. 

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

//avx512 + «маски для остатка»
t0 = _mm512_loadu_si512(tbl)
t1 = _mm512_loadu_si512(tbl+64)
t2 = _mm512_loadu_si512(tbl+128)
t3 = _mm512_loadu_si512(tbl+192)
for(i=0;i<len;i+=64){
 x0=_mm512_loadu_si512(src+i)
 x1=_mm512_permutex2var_epi8 (t0, x0, t1)
 x2=_mm512_permutex2var_epi8 (t2, x0, t3)
 m =_mm512_cmpgt_epu8_mask (x0,127)
 x3=_mm512_mask_blend_epi8 (m, x1, x2)
 _mm512_storeu_si512(dst+i,x3)
 }
int tail=len-i;
__mmask64 mtail=(0xFFFFFFFFFFFFFFFFU>>(64 - tail))
if(tail>0){
  x0 = _mm512_maskz_loadu_epi8(mtail, src+i);
  x1 = _mm512_permutex2var_epi8 (t0, x0, t1)
  x2 = _mm512_permutex2var_epi8 (t2, x0, t3)
  m  = _mm512_cmpgt_epu8_mask (x0,127)
  x3 = _mm512_mask_blend_epi8 (m, x1, x2)
  _mm512_mask_storeu_epi8(dst+i, mtail, x3)
}

Запустим оба решения на ноутбуке с процессором Сore i5-1135G7, 2.4 ГГц и сравним результаты.

AVX512 + «C код остатка»:

len=1024: 34 мс

len=1087: 52 мс

AVX512 + «маски для остатка»:

len=1024: 34 мс

len=1087: 36 мс

Ускорение: 52 мс / 36 мс = 1.44 раза.

Время работы функции для вектора длиной 1087 элементов составляло 52 миллисекунды, а благодаря масочным регистрам оно сократилось до 36 миллисекунд. То есть без применения битовой маски мы могли потерять 40% производительности.

Краткий вывод об инструкциях AVX-512

Использование инструкций AVX-512 позволяет значительно ускорить выполнение кода:

  1. Можно сделать байтовую перестановку на все 256 бит в 18 раз быстрее, чем с обычным кодом на C.

  2. Масочные регистры позволяют быстрее обрабатывать число элементов, не кратное 16.

Конечно, инструкции AVX-512 подходят не только для ускорения байтовых перестановок. Инструкций великое множество, и каждая может быть полезна оптимизатору в разных ситуациях. 

Полезные ссылки

Напоследок поделюсь тремя ссылками, которые могут быть полезны оптимизатору. Все материалы — на английском:

  1. Как и когда использовать инструкции AVX-512.

  2. Онлайн-справочник со списком интринсиков

  3. Подробное описание микроархитектур процессоров.

Теги:
Хабы:
Всего голосов 42: ↑39 и ↓3+48
Комментарии23

Публикации

Информация

Сайт
yadro.com
Дата регистрации
Дата основания
Численность
5 001–10 000 человек
Местоположение
Россия
Представитель
Ульяна Соловьева