Месяц назад мне в телеграм написал человек и предложил доступ к системе с процессором Эльбрус-8СВ.
И, конечно же, я согласился. Так как мне интересно.
Не каждый день неизвестные люди в Интернете предлагают доступ к удалённым хостам.
Меня зовут Леонид Лагунов и я математик-программист.
Я решил написать про своё знакомство с процессором Эльбрус.
Содержание:
Предшествующий опыт работы с VLIW
Ранее мне приходилось работать с системой TI KeyStone II (4 ядра ARM Cortex-A15 и 8 ядер TMS320C66). Для этой системы мной были реализованы алгоритмы FFT-1D (для 218 точек) и FFT-2D (для 29 точек), с помощью которых решалась задача поиска картинки в картинке. Управляющий код работал на ARM (написан на Си), а рабочий код для TMS320C66 писался на ассемблере. В итоге был опубликован материал, с которым можно ознакомиться вот тут.
Процессор TMS320C66 является VLIW-процессором. В привычных архитектурах (x86-64, ARM, MIPS и т.д.) инструкции выполняются последовательно (с точки зрения программиста). Для ускорения работы применяется аппаратная конвейеризация и внеочередное исполнение инструкций. В VLIW-процессоре конвейеризацией инструкций управляет сама программа, т.е. в коде явно указано в какой такт какую инструкцию следует запустить. В общем случае одновременно выполняются несколько инструкций (запущены независимо и завершатся независимо друг от друга).
Чтобы написать программу на ассемблере для VLIW-процессора, требуется подробно знать, что же каждая инструкция делает в каждый такт своего исполнения. Для TMS320C66 данная информация находится в официальной документации (SPRUGH7).
Начинаю изучать Эльбрус-8СВ
Эльбрус тоже является VLIW-процессором. Поэтому у меня появилась мысль решить на нём эту же задачу (написать реализацию алгоритма FFT на ассемблере).
Полной документации на инструкции процессора мне найти не удалось. В руководстве на сайте МЦСТ показано, что можно написать код на Си, скомпилировать из него ассемблерный вывод (lcc -S) и проанализировать этот вывод.
Используя такой подход, можно лишь догадываться о длительности исполнения инструкций. Поначалу я решил, что такая информация о каждой инструкции отсутствует и у меня не будет возможности писать эффективный код на ассемблере. Однако, позже я нашёл в руководстве раздел с этой информацией (глава 12. Планирование последовательности операций).
Как я понял, в Эльбрус «длительность исполнения инструкции» не является фиксированной величиной и зависит от того, куда пойдут результаты этой инструкции. Мне показалось, что этот факт резко усложняет написание кода человеком (в TMS320C66 было проще).
В итоге я решил начать с реализации какого-нибудь простого алгоритма на Си, чтобы изучать его ассемблерный вывод.
Написал вот такое:
double avg(uint8_t *data, size_t count)
{
int sum = 0;
for(size_t i = 0; i < count; ++i)
sum += data[i];
return (double)sum / count;
}Скомпилировал с ключами -O3 -S и получил код на ассемблере.
Код на ассемблере
...
{
setwd wsz = 0x9, nfx = 0x1, dbl = 0x0
setbn rsz = 0x4, rbs = 0x4, rcur = 0x0
setbp psz = 0x3
disp %ctpr1, .L124
rwd,0 _f64,_lts1 0xff0000000000, %lsr
}
{
cmpbdb,0,sm 0x3, %r1, %pred1
cmpbdb,1,sm 0x1, %r1, %pred3
ldb,2,sm %r0, %r4, %b[8]
ldb,3,sm %r0, 0x2, %b[4]
cmpbdb,4,sm 0x2, %r1, %pred2
ldb,5,sm %r0, 0x1, %b[6]
}
{
nop 2
ldb,0,sm %r0, 0x3, %b[2]
addd,1,sm 0x4, 0x0, %b[5]
addd,2,sm 0x4, 0x1, %b[3]
}
.L124:
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
abp abpf=1, abpt=1
ct %ctpr1 ? %pred3
addd,2,sm %b[3], 0x1, %b[1]
ldb,3,sm %r0, %b[5], %b[0]
cmpbdb,4,sm %b[5], %r1, %pred0
adds,5,sm %r3, %b[8], %r3
}
...Поизучал ассемблерный код. Появились вопросы.
Анализ ассемблерного кода
Фигурными скобками выделены инструкции, которые запускаются одновременно.
%r0 = data (первый аргумент функции)%r1 = count (второй аргумент функции)%r3 = sum
Видно, что делается prolog (до метки .L124).
Виден цикл (метка .L124).
Видно, что имена регистров %pred крутятся с помощью abp.
Вероятно, номера %b крутятся с помощью abn.
При этом номера %b прокручиваются по 2 раза за такт (как это получается?).
Зачем вообще здесь вращаются регистры?
В TMS320C66 я обходился без вращения регистров.
Далее был запрос в тех.поддержку МЦСТ с просьбой помочь разобраться с тем, что здесь накомпилировалось, и указать, как можно ускорить данный код.
Тех.поддержка дала следующие советы:
использовать
int64_tв счётчике циклаиспользовать
#pragma loop count(100)для активации режима APBиспользовать интринсики
pshufbиpaddwдля векторных вычисленийиспользовать раскрутку цикла для заполнения всех вычислительных юнитов
использовать
#pragma unroll(1)или ключ-fno-loop-unroll(отмена автоматической раскрутки циклов) для упрощения ассемблерного выводаиспользовать EML — там уже всё есть
Постановка задачи
Так как я пытаюсь изучить ассемблер Эльбруса, буду писать всё самостоятельно без EML, постепенно улучшая результаты.
По ходу изучения было решено заменить вычисление среднего арифметического на просто сумму элементов массива. Поэтому далее всё будет описываться так, будто бы задача изначально была такая:
Дано: указатель на массив и количество элементов в массиве.
Требуется: найти сумму элементов массива.
Элементами массива пусть будут целочисленные байты (uint8_t).
Требования ко входным данным
В полном решении задачи массив следует разбить на 3 участка: начало, середина и хвост. Середина выбирается исходя из свойств алгоритма и обрабатывается наиболее эффективно, а начало и хвост являются обрезками и обрабатываются менее эффективно до и после середины соответственно.
Свойства алгоритма могут выдвигать требования к:
— кратности количества элементов массива
— выровненности (align) указателя на массив
В представленных ниже реализациях отсутствует обработка начала и хвоста массива. Считается, что входные данные целиком являются серединой, т.е. удовлетворяют требованиям алгоритма.
Наименьшее общее кратное по требованиям к:
— кратности количества элементов массива равно 192
— выровненности (align) указателя на массив равно 16 байт (128 бит)
Это значит, что для тестирования алгоритмов на одних и тех же данных:
— количество элементов должно быть кратно 192
— указатель на массив должен быть выровнен (align) на 16 байт (128 бит)
Флаги компиляции
При компиляции использовались следующие ключи lcc:
-Wall -O3 -ffast-math
Как измерялось время
Замеры времени делались двумя способами (ощутимой разницы не замечено):
1) с помощью функции clock():
int t0 = clock();
/*** здесь измеряемый код ***/
int t1 = clock();
int usec = t1 - t0;2) с помощью функции clock_gettime():
struct timespec t0, t1;
clock_gettime(CLOCK_REALTIME, &t0);
/*** здесь измеряемый код ***/
clock_gettime(CLOCK_REALTIME, &t1);
int usec = (t1.tv_sec - t0.tv_sec)*1000000 + (t1.tv_nsec - t0.tv_nsec)/1000;Также использовалось чтение счётчика тактов процессора:
uint64_t get_clock_count()
{
uint64_t dst;
#pragma asm_inline
asm ("rrd %%clkr, %0" : "=r" (dst));
return dst;
}Итоговые результаты сведены в таблицу, которая находится после описания всех шагов.
После таблицы лежат ссылки на код.
Пишем программу
Далее везде перед основным циклом используется #pragma unroll(1), чтобы компилятор не раскручивал этот цикл. Раскрутку будем делать самостоятельно, иначе можно запутаться.
Сначала рассмотрим простейшую реализацию подсчёта суммы элементов массива.
Сразу применим совет «использовать int64_t в счётчике цикла».
0. Простейшая реализация
uint64_t base(uint8_t *data, uint32_t count)
{
uint64_t sum = 0;
#pragma unroll(1)
for(int64_t i = 0; i < count; ++i)
sum += data[i];
return sum;
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
ldb,3,sm %r4, %b[3], %b[0]
addd,4,sm %b[3], 0x1, %b[1]
addd,5,sm %r3, %b[6], %r3
}Основной цикл на ассемблере состоит из двух инструкций (остальные нам не интересны): - ldb (запросить чтение байта из памяти в регистр) - addd (сложить два 64-битных регистра в третий)
Ещё один addd используется для инкремента индекса массива (int64_t i).
Все инструкции, заключённые в фигурные скобки, запускаются одновременно в один такт.
В текущий такт работы цикла инструкция addd использует данные, запрошенные инструкцией ldb за несколько тактов до текущего.
Чтобы это работало, перед циклом происходит предварительный запрос нескольких первых байтов (не видно в данном коде).
Замеры времени: [count = 7080000] 9837 usec [count = 600000000] 1785000 usec
Далее посмотрим, как повлияет на время работы вставка #pragma loop count(100).
1. Включение режима APB (Array Prefetch Buffer)
uint64_t base_apb(uint8_t *data, uint32_t count)
{
uint64_t sum = 0;
#pragma unroll(1)
#pragma loop count(100) // enable apb (ld -> mova)
for(int64_t i = 0; i < count; ++i)
sum += data[i];
return sum;
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
addd,5,sm %r3, %b[10], %r3 ? %pcnt0
movab,1 area=0, ind=0, am=1, be=0, %b[0]
}Основной цикл на ассемблере состоит из двух инструкций: - movab (считать байт из APB в регистр) - addd (сложить два 64-битных регистра в третий)
Здесь произошла замена инструкции ldb на movab.
Пропал инкремент индекса массива (данные читаются последовательно из APB).
Также перед циклом происходит активация APB (не видно в данном коде).
Замеры времени: [count = 7080000] 6078 usec (100%) [count = 600000000] 551000 usec (100%)
Видим ускорение по сравнению с предыдущим вариантом из-за отправки запросов на чтение памяти заранее.
Данный вариант берём за контрольный (100%), в дальнейшем будем сравниваться с ним.
Все дальнейшие реализации сразу используют int64_t в счётч��ке цикла и вставку #pragma loop count(100). Есть ощущение, что эти действия являются стандартным началом любого длинного цикла обработки данных и можно не проверять, как будет без них.
Теперь поиграем с раскруткой цикла, чтобы посмотреть, что это даст.
2. Раскрутка цикла (x2)
// count % 2 == 0
uint64_t base_x2(uint8_t *data, uint32_t count)
{
uint64_t sum[2] = {0,0};
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < count; i += 2)
{
sum[0] += data[i + 0];
sum[1] += data[i + 1];
}
return sum[0] + sum[1];
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
addd,4,sm %r4, %b[11], %r4 ? %pcnt0
addd,5,sm %r3, %b[10], %r3 ? %pcnt0
movab,0 area=0, ind=0, am=0, be=0, %b[1]
movab,1 area=0, ind=1, am=1, be=0, %b[0]
}В основном цикле стало по две инструкции movab и addd.
Теперь выполняются 2 сложения за один такт.
Замеры времени: [count = 7080000] 3075 usec (1/2) [count = 600000000] 281000 usec (1/2)
Ускорение в 2 раза относительно контрольного варианта.
Ещё раскрутим цикл.
3. Раскрутка цикла (x3)
// count % 3 == 0
uint64_t base_x3(uint8_t *data, uint32_t count)
{
uint64_t sum[3] = {0,0,0};
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < count; i += 3)
{
sum[0] += data[i + 0];
sum[1] += data[i + 1];
sum[2] += data[i + 2];
}
return sum[0] + sum[1] + sum[2];
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
addd,3,sm %r4, %b[10], %r4 ? %pcnt0
addd,4,sm %r3, %b[20], %r3 ? %pcnt0
addd,5,sm %r5, %b[11], %r5 ? %pcnt0
movab,0 area=0, ind=0, am=0, be=0, %b[0]
movab,1 area=0, ind=1, am=1, be=0, %b[10]
movab,3 area=0, ind=0, am=1, be=0, %b[1]
}В основном цикле стало по три инструкции movab и addd.
Теперь выполняются 3 сложения за один такт.
Замеры времени: [count = 7080000] 2029 usec (1/3) [count = 600000000] 188000 usec (1/3)
Ускорение в 3 раза относительно контрольного варианта.
И ещё немного раскрутки.
4. Раскрутка цикла (x4)
// count % 4 == 0
uint64_t base_x4(uint8_t *data, uint32_t count)
{
uint64_t sum[4] = {0,0,0,0};
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < count; i += 4)
{
sum[0] += data[i + 0];
sum[1] += data[i + 1];
sum[2] += data[i + 2];
sum[3] += data[i + 3];
}
return sum[0] + sum[1] + sum[2] + sum[3];
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
addd,2,sm %r6, %b[20], %r6 ? %pcnt0
addd,3,sm %r4, %b[10], %r4 ? %pcnt0
addd,4,sm %r5, %b[11], %r5 ? %pcnt0
addd,5,sm %r3, %b[21], %r3 ? %pcnt0
movab,0 area=0, ind=0, am=0, be=0, %b[10]
movab,1 area=0, ind=1, am=1, be=0, %b[0]
movab,2 area=0, ind=0, am=0, be=0, %b[1]
movab,3 area=0, ind=1, am=1, be=0, %b[11]
}В основном цикле стало по четыре инструкции movab и addd.
Теперь выполняются 4 сложения за один такт.
Замеры времени: [count = 7080000] 1522 usec (1/4) [count = 600000000] 142000 usec (1/4)
Ускорение в 4 раза относительно контрольного варианта.
Дальнейшая раскрутка цикла не даёт ускорения, потому что не может быть больше четырёх movab в одном такте.
При дальнейшей раскрутке тело цикла разбивается на несколько тактов.
Теперь попробуем использовать интринсики pshufb и paddw для векторных вычислений.
5. Использование SIMD64 (movah, pshufb, paddw)
// count % 2 == 0
// alignof(data) == 2
uint64_t SIMD64_1(uint8_t *data, uint32_t count)
{
uint64_t sum = 0;
uint16_t *p16 = (uint16_t *)data;
while(count > 0)
{
uint64_t sum2 = 0;
int max_block_size = 2 * (1<<24);
int block_size = count < max_block_size ? count : max_block_size;
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < block_size/2; ++i)
{
uint64_t mask_zzz1zzz0 = 0x8080800180808000;
uint64_t data2 = __builtin_e2k_pshufb(0, *p16++, mask_zzz1zzz0);
sum2 = __builtin_e2k_paddw(sum2, data2);
}
sum += (uint32_t)sum2 + (sum2 >> 32);
count -= block_size;
}
return sum;
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
paddw,0,sm %r1, %g16, %r1 ? %pcnt0
pshufb,1,sm 0x0, %b[10], %r3, %g16
movah,1 area=0, ind=0, am=1, be=0, %b[0]
}Основной цикл на ассемблере состоит из трёх инструкций: - movah (считать "half" из APB в регистр) - pshufb (переложить байты из одного 64-битного регистра в другой в соответствии с маской) - paddw (сложить два 64-битных регистра-вектора-"word" в третий)
"half" = 2 байта (16 бит) "word" = 4 байта (32 бита)
В данном варианте используются 64-битные SIMD инструкции.
Здесь учтено возможное переполнение "word"-половинок.
Выполняются 2 сложения за один такт.
Замеры времени: [count = 7080000] 3040 usec (1/2) [count = 600000000] 281000 usec (1/2)
Ускорение в 2 раза относительно контрольного варианта.
Попробуем раскрутить цикл.
6. SIMD64 раскрутка цикла (x2)
// count % 4 == 0
// alignof(data) == 2
uint64_t SIMD64_1_x2(uint8_t *data, uint32_t count)
{
uint64_t sum = 0;
uint16_t *p16 = (uint16_t *)data;
while(count > 0)
{
uint64_t sum2[2] = {0,0};
int max_block_size = 2 * 2 * (1<<24);
int block_size = count < max_block_size ? count : max_block_size;
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < block_size/2; i += 2)
{
uint64_t mask_zzz1zzz0 = 0x8080800180808000;
uint64_t data2_0 = __builtin_e2k_pshufb(0, *p16++, mask_zzz1zzz0);
uint64_t data2_1 = __builtin_e2k_pshufb(0, *p16++, mask_zzz1zzz0);
sum2[0] = __builtin_e2k_paddw(sum2[0], data2_0);
sum2[1] = __builtin_e2k_paddw(sum2[1], data2_1);
}
sum += (uint32_t)sum2[0] + (sum2[0] >> 32);
sum += (uint32_t)sum2[1] + (sum2[1] >> 32);
count -= block_size;
}
return sum;
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
paddw,0,sm %r3, %g17, %r3 ? %pcnt0
pshufb,1,sm 0x0, %b[11], %r4, %g17
paddw,3,sm %r1, %g16, %r1 ? %pcnt0
pshufb,4,sm 0x0, %b[10], %r4, %g16
movah,0 area=0, ind=0, am=0, be=0, %b[1]
movah,1 area=0, ind=2, am=1, be=0, %b[0]
}В основном цикле стало по две инструкции movah, pshufb и paddw.
Теперь выполняются 4 сложения за один такт.
Замеры времени: [count = 7080000] 1520 usec (1/4) [count = 600000000] 141000 usec (1/4)
Ускорение в 4 раза относительно контрольного варианта.
Дальнейшая раскрутка цикла не даёт ускорения, потому что не может быть больше двух paddw в одном такте.
При дальнейшей раскрутке тело цикла разбивается на несколько тактов.
Однако, мы можем заменить paddw на просто addd. Разница будет, если в addd произойдёт перенос из младшей половины регистра в старшую. Если не допускать такого переноса, то разницы не будет.
Отсюда следующий пункт.
7. SIMD64 раскрутка цикла (x3)
// count % 6 == 0
// alignof(data) == 2
uint64_t SIMD64_1_x3(uint8_t *data, uint32_t count)
{
uint64_t sum = 0;
uint16_t *p16 = (uint16_t *)data;
while(count > 0)
{
uint64_t sum2[3] = {0,0,0};
int max_block_size = 3 * 2 * (1<<24);
int block_size = count < max_block_size ? count : max_block_size;
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < block_size/2; i += 3)
{
uint64_t mask_zzz1zzz0 = 0x8080800180808000;
sum2[0] += __builtin_e2k_pshufb(0, *p16++, mask_zzz1zzz0);
sum2[1] += __builtin_e2k_pshufb(0, *p16++, mask_zzz1zzz0);
sum2[2] += __builtin_e2k_pshufb(0, *p16++, mask_zzz1zzz0);
}
sum += (uint32_t)sum2[0] + (sum2[0] >> 32);
sum += (uint32_t)sum2[1] + (sum2[1] >> 32);
sum += (uint32_t)sum2[2] + (sum2[2] >> 32);
count -= block_size;
}
return sum;
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
pshufb,0,sm 0x0, %b[11], %r5, %b[20]
pshufb,1,sm 0x0, %b[10], %r5, %b[19]
addd,2,sm %r3, %b[28], %r3 ? %pcnt0
pshufb,3,sm 0x0, %b[20], %r5, %b[11]
addd,4,sm %r1, %b[27], %r1 ? %pcnt0
addd,5,sm %r4, %b[19], %r4 ? %pcnt0
movah,0 area=0, ind=0, am=0, be=0, %b[1]
movah,1 area=0, ind=2, am=1, be=0, %b[0]
movah,3 area=0, ind=0, am=1, be=0, %b[10]
}В основном цикле стало по три инструкции movah, pshufb и addd.
Теперь выполняются 6 сложений за один такт.
Замеры времени: [count = 7080000] 1014 usec (1/6) [count = 600000000] 95000 usec (1/6)
Ускорение в 6 раз относительно контрольного варианта.
Дальнейшая раскрутка цикла не даёт ускорения, потому что использованы все 6 вычислительных юнитов.
При дальнейшей раскрутке тело цикла разбивается на несколько тактов.
Можно попробовать использовать вектора из четырёх "half"-счётчиков, вместо двух "word"-счётчиков. В этом случае следует учесть возможное переполнение "half"-счётчиков.
Отсюда следующий пункт.
8. SIMD64 раскрутка цикла (x3) + 16-битные "half"-счётчики
// count % 12 == 0
// alignof(data) == 4
uint64_t SIMD64_1_x3_16bit(uint8_t *data, uint32_t count)
{
uint64_t sum = 0;
uint32_t *p32 = (uint32_t *)data;
while(count > 0)
{
uint64_t sum4[3] = {0,0,0};
int max_block_size = 3 * 4 * (1<<8);
int block_size = count < max_block_size ? count : max_block_size;
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < block_size/4; i += 3)
{
uint64_t mask_z3z2z1z0 = 0x8003800280018000;
sum4[0] += __builtin_e2k_pshufb(0, *p32++, mask_z3z2z1z0);
sum4[1] += __builtin_e2k_pshufb(0, *p32++, mask_z3z2z1z0);
sum4[2] += __builtin_e2k_pshufb(0, *p32++, mask_z3z2z1z0);
}
sum += (uint16_t)sum4[0] + (uint16_t)(sum4[0] >> 16) + (uint16_t)(sum4[0] >> 32) + (sum4[0] >> 48);
sum += (uint16_t)sum4[1] + (uint16_t)(sum4[1] >> 16) + (uint16_t)(sum4[1] >> 32) + (sum4[1] >> 48);
sum += (uint16_t)sum4[2] + (uint16_t)(sum4[2] >> 16) + (uint16_t)(sum4[2] >> 32) + (sum4[2] >> 48);
count -= block_size;
}
return sum;
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
pshufb,0,sm 0x0, %b[20], %r5, %b[19]
pshufb,1,sm 0x0, %b[11], %r5, %b[20]
addd,2,sm %r3, %b[27], %r3 ? %pcnt0
pshufb,3,sm 0x0, %b[10], %r5, %b[11]
addd,4,sm %r1, %b[28], %r1 ? %pcnt0
addd,5,sm %r4, %b[19], %r4 ? %pcnt0
movaw,0 area=0, ind=0, am=0, be=0, %b[10]
movaw,1 area=0, ind=4, am=1, be=0, %b[1]
movaw,3 area=0, ind=0, am=1, be=0, %b[0]
}Здесь произошла замена инструкции movah на movaw (считать "word" из APB в регистр).
Теперь выполняются 12 сложений за один такт.
Замеры времени: [count = 7080000] 693 usec [count = 600000000] 100000 usec
Мы ожидаем увидеть ускорение в 12 раз относительно контрольного варианта.
Однако, видим ускорение всего в 6 раз (почти не отличается от предыдущего варианта).
По сравнению с предыдущим вариантом видно ускорение при малом объёме данных и замедление при большом объёме данных.
Ускорение при малом объёме данных, скорее всего, связано с тем, что данные лежат в кэше.
Все предыдущие замеры также были подвержены этому факту.
При большом объёме данных они честно читаются из памяти. Отсутствие ускорения можно объяснить большими накладными расходами на активацию APB. Каждый раз перед внутренним циклом нужно активировать APB, а цикл состоит всего из 256 итераций. В предыдущих вариантах внутренний цикл состоял из 224 итераций и там активация APB не была заметна.
Можно подумать: «нельзя ли один раз активировать APB, а не перед каждым внутренним циклом?» Ответ такой: если бы мы могли так сделать, то это был бы один большой цикл. И в этом одном цикле нужно было бы находить место для инструкций скидывания "half"-счётчиков. Но все вычислительные юниты уже заняты.
Теперь попробуем заменить инструкции SIMD64 на аналогичные SIMD128.
9. Использование SIMD128 (movaw, qpshufb, qpaddd)
// count % 4 == 0
// alignof(data) == 4
uint64_t SIMD128_1(uint8_t *data, uint32_t count)
{
uint64_t sum = 0;
uint32_t *p32 = (uint32_t *)data;
while(count > 0)
{
__v2di sum4 = {0};
int max_block_size = 4 * (1<<24);
int block_size = count < max_block_size ? count : max_block_size;
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < block_size/4; ++i)
{
__v2di mask_zzz3zzz2zzz1zzz0 = {0x8080800180808000, 0x8080800380808002};
sum4 += __builtin_e2k_qpshufb((__v2di){*p32}, (__v2di){*p32}, mask_zzz3zzz2zzz1zzz0);
++p32;
}
union union_128 un = {.i128 = sum4};
sum += (uint64_t)un.u32[0] + un.u32[1] + un.u32[2] + un.u32[3];
count -= block_size;
}
return sum;
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
qpaddd,0,sm %r1, %g16, %r1 ? %pcnt0
qpshufb,1,sm %b[9], %b[9], %r3, %g16
qppackdl,4,sm 0x0, %b[10], %b[1]
movaw,1 area=0, ind=0, am=1, be=0, %b[0]
}Основной цикл на ассемблере состоит из трёх инструкций: - movaw (считать "word" из APB в регистр) - qpshufb (переложить байты из одного 128-битного регистра в другой в соответствии с маской) - qpaddd (сложить два 128-битных регистра-вектора-"dword" в третий)
"word" = 4 байта (32 бита)"dword" = 8 байт (64 бита)
Также видна инструкция qppackdl (формирование 128-битного регистра из двух 64-битных половинок). Эта инструкци�� соответствует созданию переменной типа __v2di из (uint64_t)0 и (uint64_t)(*p32). Мне не удалось придумать, как написать код на Си так, чтобы компилятор её не генерировал.
Единственное найденное решение (использовал далее) — читать сразу по 16 байт (128 бит).
В данном варианте используются 128-битные SIMD инструкции.
Здесь учтено возможное переполнение "word"-четвертинок.
Выполняются 4 сложения за один такт.
Замеры времени: [count = 7080000] 1520 usec (1/4) [count = 600000000] 141000 usec (1/4)
Ускорение в 4 раза относительно контрольного варианта.
Попробуем раскрутить цикл.
10. SIMD128 раскрутка цикла (x2*2)
// count % 16 == 0
// alignof(data) == 16
uint64_t SIMD128_1_x2x2(uint8_t *data, uint32_t count)
{
uint64_t sum = 0;
__v2di *p128 = (__v2di *)data;
while(count > 0)
{
//__v2di sum4[4] = {{0}, {0}, {0}, {0}};
__v2di sum4_0 = {0};
__v2di sum4_1 = {0};
__v2di sum4_2 = {0};
__v2di sum4_3 = {0};
int max_block_size = 16 * (1<<24);
int block_size = count < max_block_size ? count : max_block_size;
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < block_size/16; ++i)
{
__v2di mask_zzz3zzz2zzz1zzz0 = {0x8080800180808000, 0x8080800380808002};
__v2di mask_zzz7zzz6zzz5zzz4 = {0x8080800580808004, 0x8080800780808006};
__v2di mask_zzzBzzzAzzz9zzz8 = {0x8080800980808008, 0x8080800B8080800A};
__v2di mask_zzzFzzzEzzzDzzzC = {0x8080800D8080800C, 0x8080800F8080800E};
sum4_0/*sum4[0]*/ += __builtin_e2k_qpshufb(*p128, *p128, mask_zzz3zzz2zzz1zzz0);
sum4_1/*sum4[1]*/ += __builtin_e2k_qpshufb(*p128, *p128, mask_zzz7zzz6zzz5zzz4);
sum4_2/*sum4[2]*/ += __builtin_e2k_qpshufb(*p128, *p128, mask_zzzBzzzAzzz9zzz8);
sum4_3/*sum4[3]*/ += __builtin_e2k_qpshufb(*p128, *p128, mask_zzzFzzzEzzzDzzzC);
++p128;
}
union union_128 un0 = {.i128 = sum4_0/*sum4[0]*/};
union union_128 un1 = {.i128 = sum4_1/*sum4[1]*/};
union union_128 un2 = {.i128 = sum4_2/*sum4[2]*/};
union union_128 un3 = {.i128 = sum4_3/*sum4[3]*/};
sum += (uint64_t)un0.u32[0] + un0.u32[1] + un0.u32[2] + un0.u32[3];
sum += (uint64_t)un1.u32[0] + un1.u32[1] + un1.u32[2] + un1.u32[3];
sum += (uint64_t)un2.u32[0] + un2.u32[1] + un2.u32[2] + un2.u32[3];
sum += (uint64_t)un3.u32[0] + un3.u32[1] + un3.u32[2] + un3.u32[3];
count -= block_size;
}
return sum;
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
qpaddd,0,sm %r1, %b[5], %r1 ? %pcnt0
qpshufb,1,sm %b[6], %b[6], %r10, %b[3]
qpaddd,3,sm %r5, %b[3], %r5 ? %pcnt0
qpshufb,4,sm %b[6], %b[6], %r9, %b[1]
movaqp,1 area=0, ind=0, am=1, be=0, %b[0]
}
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
qpaddd,0,sm %r4, %b[3], %r4 ? %pcnt0
qpshufb,1,sm %b[4], %b[4], %r7, %b[3]
qpaddd,3,sm %r3, %b[1], %r3 ? %pcnt0
qpshufb,4,sm %b[4], %b[4], %r6, %b[1]
}"quad" = 16 байт (128 бит)
Здесь произошла замена инструкции movaw на movaqp (считать "quad" из APB в регистр).
Это сделано для борьбы с генерацией инструкции qppackdl (см. предыдущий пункт).
Для обработки 16 байт пришлось сделать раскрутку дважды, т.е. раскрутить цикл в 4 раза.
В основном цикле стало по четыре инструкции qpshufb и qpaddd.
Теперь выполняются 16 сложений за 2 такта (8 сложений за один такт).
Замеры времени: [count = 7080000] 760 usec (1/8) [count = 600000000] 67000 usec (1/8)
Ускорение в 8 раз относительно контрольного варианта.
Дальнейшая раскрутка цикла не даёт ускорения, потому что не может быть больше двух qpaddd в одном такте.
При дальнейшей раскрутке тело цикла разбивается на несколько тактов.
Можно попробовать использовать вектора из восьми "half"-счётчиков, вместо четырёх "word"-счётчиков. В этом случае следует учесть возможное переполнение "half"-счётчиков.
Отсюда следующий пункт.
11. SIMD128 раскрутка цикла (x2) + 16-битные "half"-счётчики
// count % 16 == 0
// alignof(data) == 16
uint64_t SIMD128_1_x2_16bit(uint8_t *data, uint32_t count)
{
uint64_t sum = 0;
__v2di *p128 = (__v2di *)data;
while(count > 0)
{
//__v2di sum8[2] = {{0}, {0}};
__v2di sum8_0 = {0};
__v2di sum8_1 = {0};
int max_block_size = 16 * (1<<8);
int block_size = count < max_block_size ? count : max_block_size;
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < block_size/16; ++i)
{
__v2di mask_z7z6z5z4z3z2z1z0 = {0x8003800280018000, 0x8007800680058004};
__v2di mask_zFzEzDzCzBzAz9z8 = {0x800B800A80098008, 0x800F800E800D800C};
sum8_0/*sum8[0]*/ += __builtin_e2k_qpshufb(*p128, *p128, mask_z7z6z5z4z3z2z1z0);
sum8_1/*sum8[1]*/ += __builtin_e2k_qpshufb(*p128, *p128, mask_zFzEzDzCzBzAz9z8);
++p128;
}
union union_128 un0 = {.i128 = sum8_0/*sum8[0]*/};
union union_128 un1 = {.i128 = sum8_1/*sum8[1]*/};
sum += (uint64_t)un0.u16[0] + un0.u16[1] + un0.u16[2] + un0.u16[3] + un0.u16[4] + un0.u16[5] + un0.u16[6] + un0.u16[7];
sum += (uint64_t)un1.u16[0] + un1.u16[1] + un1.u16[2] + un1.u16[3] + un1.u16[4] + un1.u16[5] + un1.u16[6] + un1.u16[7];
count -= block_size;
}
return sum;
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
qpaddd,0,sm %r3, %g17, %r3 ? %pcnt0
qpshufb,1,sm %b[10], %b[10], %r5, %g17
qpaddd,3,sm %r1, %g16, %r1 ? %pcnt0
qpshufb,4,sm %b[10], %b[10], %r4, %g16
movaqp,1 area=0, ind=0, am=1, be=0, %b[0]
}Теперь выполняются 16 сложений за один такт.
Замеры времени: [count = 7080000] 499 usec [count = 600000000] 63000 usec
Здесь такая же ситуация, как в пункте 8.
Мы ожидаем увидеть ускорение в 16 раз относительно контрольного варианта.
Однако, видим ускорение всего в 9 раз (почти не отличается от предыдущего варианта).
По сравнению с предыдущим вариантом видно ускорение при малом объёме данных и слабое ускорение при большом объёме данных.
Ускорение при малом объёме данных, скорее всего, связано с тем, что данные лежат в кэше.
При большом объёме данных они честно читаются из памяти. Слабое ускорение можно объяснить большими накладными расходами на активацию APB. Каждый раз перед внутренним циклом нужно активировать APB, а цикл состоит всего из 256 итераций. В предыдущих вариантах внутренний цикл состоял из 224 итераций и там активация APB не была заметна.
Можно подумать: «нельзя ли один раз активировать APB, а не перед каждым внутренним циклом?» Ответ такой: если бы мы могли так сделать, то это был бы один большой цикл. И в этом одном цикле нужно было бы находить место для инструкций скидывания "half"-счётчиков.
В отличие от пункта 8 здесь не все вычислительные юниты заняты (юниты 2 и 5 свободны).
Можно ли с их помощью организовать скидывание "half"-счётчиков, я не изучал.
Далее была найдена «волшебная» SIMD-инструкция psadbw.
Эта инструкция как будто специально создана для решения этой задачи.
Попробуем использовать её.
12. Использование SIMD64 (2) (movad, psadbw, addd)
// count % 8 == 0
// alignof(data) == 8
uint64_t SIMD64_2(uint8_t *data, uint32_t count)
{
uint64_t sum = 0;
uint64_t *p64 = (uint64_t *)data;
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < count/8; ++i)
{
sum += __builtin_e2k_psadbw(0, p64[i]);
}
return sum;
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
psadbw,4,sm 0x0, %b[10], %b[1]
addd,5,sm %r3, %b[9], %r3 ? %pcnt0
movad,1 area=0, ind=0, am=1, be=0, %b[0]
}Основной цикл на ассемблере состоит из трёх инструкций: - movad (считать "dword" из APB в регистр) - psadbw (сложить модули разностей байтов из двух 64-битных регистров и поместить результат в третий регистр) - addd (сложить два 64-битных регистра в третий)
"dword" = 8 байт (64 бит)
В данном варианте используются 64-битные SIMD инструкции.
Если один из аргументов инструкции psadbw равен нулю, то фактически выполняется сложение модулей всех 8 байтов второго аргумента.
Выполняются 8 сложений за один такт.
Замеры времени: [count = 7080000] 761 usec (1/8) [count = 600000000] 71000 usec (1/8)
Ускорение в 8 раз относительно контрольного варианта.
Посмотрим на раскрутку цикла.
13. SIMD64 (2) раскрутка цикла (x2)
// count % 16 == 0
// alignof(data) == 8
uint64_t SIMD64_2_x2(uint8_t *data, uint32_t count)
{
uint64_t sum[2] = {0,0};
uint64_t *p64 = (uint64_t *)data;
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < count/8; i += 2)
{
sum[0] += __builtin_e2k_psadbw(0, p64[i + 0]);
sum[1] += __builtin_e2k_psadbw(0, p64[i + 1]);
}
return sum[0] + sum[1];
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
psadbw,1,sm 0x0, %b[11], %b[11]
addd,3,sm %r4, %b[19], %r4 ? %pcnt0
psadbw,4,sm 0x0, %b[10], %b[10]
addd,5,sm %r3, %b[18], %r3 ? %pcnt0
movad,0 area=0, ind=0, am=0, be=0, %b[1]
movad,1 area=0, ind=8, am=1, be=0, %b[0]
}В основном цикле стало по две инструкции movad, psadbw и addd.
Теперь выполняются 16 сложений за один такт.
Замеры времени: [count = 7080000] 380 usec (1/16) [count = 600000000] 36000 usec (1/16)
Ускорение в 16 раз относительно контрольного варианта.
Дальнейшая раскрутка цикла не даёт ускорения, потому что не может быть больше двух psadbw в одном такте.
При дальнейшей раскрутке тело цикла разбивается на несколько тактов.
Теперь попробуем использовать «волшебную» SIMD-инструкцию в версии SIMD128.
14. Использование SIMD128 (2) (movaqp, qpsadbw, qpaddd)
// count % 16 == 0
// alignof(data) == 16
uint64_t SIMD128_2(uint8_t *data, uint32_t count)
{
__v2di sum2 = {0};
__v2di *p128 = (__v2di *)data;
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < count/16; ++i)
{
sum2 += __builtin_e2k_qpsadbw((__v2di){0}, p128[i]);
}
union union_128 un = {.i128 = sum2};
return un.u64[0] + un.u64[1];
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
qpaddd,3,sm %b[14], %b[15], %b[12] ? %pcnt0
qpsadbw,4,sm %r3, %b[11], %b[11]
movaqp,1 area=0, ind=0, am=1, be=0, %b[1]
}Основной цикл на ассемблере состоит из трёх инструкций: - movaqp (считать "quad" из APB в регистр) - qpsadbw (применить psadbw к старшей и младшей половинам участвующих регистров) - qpaddd (сложить два 128-битных регистра-вектора-"dword" в третий)
"dword" = 8 байт ( 64 бит)"quad" = 16 байт (128 бит)
В данном варианте используются 128-битные SIMD инструкции.
Выполняются 16 сложений за один такт.
Замеры времени: [count = 7080000] 380 usec (1/16) [count = 600000000] 36000 usec (1/16)
Ускорение в 16 раз относительно контрольного варианта.
И, конечно, раскрутим цикл.
Стоит заметить, что наивная раскрутка путём копирования строчек кода здесь не дала ускорения. После некоторых усилий получилось добиться ускорения, если разбить обрабатываемый массив на две половины.
Анализ ассемблерного вывода показывает, что:
— при наивной раскрутке цикла компилятор генерирует чтение из одного «канала» APB
— при разбиении массива на две половины — генерирует чтение из двух «каналов» APB
15. SIMD128 (2) раскрутка цикла (x2)
// count % 32 == 0
// alignof(data) == 16
uint64_t SIMD128_2_x2(uint8_t *data, uint32_t count)
{
//__v2di sum2[2] = {{0}, {0}};
__v2di sum2_0 = {0};
__v2di sum2_1 = {0};
__v2di *p128_0 = (__v2di *) data;
__v2di *p128_1 = (__v2di *)&data[count/2];
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < count/2/16; ++i)
{
sum2_0/*sum2[0]*/ += __builtin_e2k_qpsadbw((__v2di){0}, p128_0[i]);
sum2_1/*sum2[1]*/ += __builtin_e2k_qpsadbw((__v2di){0}, p128_1[i]);
}
union union_128 un0 = {.i128 = sum2_0/*sum2[0]*/};
union union_128 un1 = {.i128 = sum2_1/*sum2[1]*/};
return un0.u64[0] + un0.u64[1] + un1.u64[0] + un1.u64[1];
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
qpaddd,0,sm %b[15], %b[29], %b[13] ? %pcnt0
qpsadbw,1,sm %r3, %b[25], %b[25]
qpaddd,3,sm %b[14], %b[28], %b[12] ? %pcnt0
qpsadbw,4,sm %r3, %b[24], %b[24]
movaqp,1 area=0, ind=0, am=1, be=0, %b[15]
movaqp,3 area=0, ind=0, am=1, be=0, %b[14]
}В основном цикле стало по две инструкции movaqp, qpsadbw и qpaddd.
Теперь выполняются 32 сложения за один такт.
Замеры времени: [count = 7080000] 190 usec (1/32) [count = 600000000] 29000 usec
Мы ожидаем увидеть ускорение в 32 раза относительно контрольного варианта.
Однако, видим ускорение всего в 19 раз (в 1,25 раз быстрее предыдущего варианта).
Точ��ее, видим ускорение в 32 раза при малом объёме данных и в 19 раз при большом.
Ускорение при малом объёме данных, скорее всего, связано с тем, что данные лежат в кэше.
При большом объёме данных они честно читаются из памяти.
Даже для получения этого результата пришлось разделить чтение из памяти на две половины. Без такого разделения компилятор генерирует чтение из одного «канала» APB и ускорения не происходит совсем.
Предположительно, мы упёрлись в пропускную способность системы памяти.
Здесь я должен был написать:
Дальнейшая раскрутка цикла не даёт ускорения, потому что не может быть больше двух qpsadbw в одном такте.
При дальнейшей раскрутке тело цикла разбивается на несколько тактов.
Всё это правда, однако, ускорение всё же есть.
Оно происходит, потому что при дальнейшей раскрутке удаётся сделать чтение памяти по порядку всеми «каналами» APB.
Смотрим как это выглядит.
16. SIMD128 (2) раскрутка цикла (x2*2)
// count % 64 == 0
// alignof(data) == 16
uint64_t SIMD128_2_x2x2(uint8_t *data, uint32_t count)
{
__v2di sum2_0 = {0};
__v2di sum2_1 = {0};
__v2di sum2_2 = {0};
__v2di sum2_3 = {0};
__v2di *p128 = (__v2di *)data;
#pragma unroll(1)
#pragma loop count(100)
for(int64_t i = 0; i < count/16; i += 4)
{
sum2_0 += __builtin_e2k_qpsadbw((__v2di){0}, p128[i + 0]);
sum2_1 += __builtin_e2k_qpsadbw((__v2di){0}, p128[i + 1]);
sum2_2 += __builtin_e2k_qpsadbw((__v2di){0}, p128[i + 2]);
sum2_3 += __builtin_e2k_qpsadbw((__v2di){0}, p128[i + 3]);
}
union union_128 un0 = {.i128 = sum2_0};
union union_128 un1 = {.i128 = sum2_1};
union union_128 un2 = {.i128 = sum2_2};
union union_128 un3 = {.i128 = sum2_3};
return un0.u64[0] + un0.u64[1] + un1.u64[0] + un1.u64[1] + un2.u64[0] + un2.u64[1] + un3.u64[0] + un3.u64[1];
}Код на ассемблере Эльбрус (основной цикл):
{
loop_mode
qpaddd,0,sm %b[23], %g17, %b[21] ? %pcnt0
qpsadbw,1,sm %r3, %b[15], %g17
qpaddd,3,sm %b[22], %g16, %b[20] ? %pcnt0
qpsadbw,4,sm %r3, %b[14], %g16
movaqp,2 area=0, ind=0, am=0, be=0, %b[23]
movaqp,3 area=0, ind=16, am=1, be=0, %b[22]
}
{
loop_mode
alc alcf=1, alct=1
abn abnf=1, abnt=1
ct %ctpr1 ? %NOT_LOOP_END
qpaddd,0,sm %b[9], %b[31], %b[7] ? %pcnt0
qpsadbw,1,sm %r3, %b[27], %b[27]
qpaddd,3,sm %b[8], %b[30], %b[6] ? %pcnt0
qpsadbw,4,sm %r3, %b[26], %b[26]
movaqp,0 area=0, ind=0, am=0, be=0, %b[9]
movaqp,1 area=0, ind=16, am=1, be=0, %b[8]
}В основном цикле стало по четыре инструкции movaqp, qpsadbw и qpaddd.
Теперь выполняются 64 сложения за 2 такта (32 сложения за один такт).
Замеры времени: [count = 7080000] 190 usec (1/32) [count = 600000000] 26500 usec
Ускорение в 32 раза при малом объёме данных и в 21 раз при большом.
Ускорение при малом объёме данных, скорее всего, связано с тем, что данные лежат в кэше.
При большом объёме данных они честно читаются из памяти.
Ускорение относительно предыдущего варианта происходит из-за чтения памяти по порядку вместо чтения из двух половин.
Однако, это всё ещё не ожидаемое ускорение в 32 раза.
Предположительно, мы упёрлись в пропускную способность системы памяти.
Ну и для полноты картины посмотрим, как работает EML.
17. Использование EML (eml_Vector_Sum_8U)
uint64_t EML(uint8_t *data, uint32_t count)
{
eml_64f sum = 0;
eml_Status status = eml_Vector_Sum_8U(data, count, &sum);
if (status != EML_OK)
puts("EML error");
return sum;
}Для сравнения была использована функция eml_Vector_Sum_8U из библиотеки EML.
Замеры времени: [count = 7080000] 191 usec [count = 600000000] 29000 usec
Скорость работы такая же, как в пункте 15.
Итоговые результаты
№ | Алгоритм | Сложений за такт | Два замера времени (микросекунды) | Реальное ускорение |
0 | base | 1 | 9837 / 1785000 | — |
1 | base_apb | 1 | 6078 / 551000 | 100% |
2 | base_x2 | 2 | 3075 / 281000 | x2 |
3 | base_x3 | 3 | 2029 / 188000 | x3 |
4 | base_x4 | 4 | 1522 / 142000 | x4 |
5 | SIMD64_1 | 2 | 3040 / 281000 | x2 |
6 | SIMD64_1_x2 | 4 | 1520 / 141000 | x4 |
7 | SIMD64_1_x3 | 6 | 1014 / 95000 | x6 |
8 | SIMD64_1_x3_16bit | 12 | 693 / 100000 | x8.8 / x5.5 |
9 | SIMD128_1 | 4 | 1520 / 141000 | x4 |
10 | SIMD128_1_x2x2 | 8 | 760 / 67000 | x8 |
11 | SIMD128_1_x2_16bit | 16 | 499 / 63000 | x12 / x8.7 |
12 | SIMD64_2 | 8 | 761 / 71000 | x8 |
13 | SIMD64_2_x2 | 16 | 380 / 36000 | x16 |
14 | SIMD128_2 | 16 | 380 / 36000 | x16 |
15 | SIMD128_2_x2 | 32 | 190 / 29000 | x32 / x19 |
16 | SIMD128_2_x2x2 | 32 | 190 / 26500 | x32 / x21 |
17 | EML | — | 191 / 29000 | x32 / x19 |
Кто хочет посмотреть весь код, вот исходный файл sum.c и ассемблерный вывод sum.s.
Исследование просадки скорости
Было замечено, что при малых размерах данных можно достичь ускорения в 32 раза по сравнению с контрольным вариантом, но при больших размерах данных такое ускорение не достигается.
Попробуем исследовать, где находится граница между малыми и большими размерами данных.
Предположительно, эта граница связана с кэшем процессора, поэтому позапускаем разные размеры и посмотрим на скорость работы.
Размеры кэшей процессора Эльбрус-8СВ:
L1D = 64 КБ (на каждое ядро)L2 = 512 КБ (на каждое ядро)L3 = 16 МБ (общий на все ядра)
Начнём с L1.
Выглядит всё стабильно.
Переходим к проверке L2.
То же самое, проблем не видно.
Смотрим на L3.
Видим, что просадки начинаются примерно после выхода за пределы L3.
Также небольшое падение скорости видно в пунктах 13 и 14 (раньше не обращал внимание).
После прохождения границы в 24 МБ (1,5 размера L3) падение прекращается.
Данное поведение похоже на работу честного 2-way кэша (2x8 МБ).
Читаем (Данные 24 МБ) = (начало 8 МБ) + (середина 8 МБ) + (хвост 8МБ):
1) читаем начало и середину — кэш заполняется
2) читаем хвост — из кэша выталкиваются самые старые 8 МБ (т.е. начало)
К концу первого прохода в кэше лежат середина и хвост массива.
При следующих проходах:
1) читаем начало — выталкивается середина
2) читаем середину — выталкивается хвост
3) читаем хвост — выталкивается начало
При 1-way кэше полные тормоза достигались бы на двукратном размере данных.
При n-way кэше полные тормоза достигаются при размере данных (1 + 1/n) от размера кэша.
Это лишь мои рассуждения. Не стоит считать меня экспертом в этом вопросе.
Если у кого-то есть идеи, как это можно объяснить по-другому, пишите в комментариях.
Динамическое изменение частоты процессора
В выводе написанной программы не просто так появился столбец Такты/мкс.
Дело в том, что у этого процессора частота меняется от 1000 МГц до 1200 МГц.
Впервые я отметил этот факт, когда посмотрел вывод lscpu:
Впоследствии, во время написания первых версий программы, я обратил внимание на разную скорость работы в зависимости от того, как заполнялся массив. При заполнении константным значением скорость была выше, чем при заполнении случайными числами. Этот эффект проявлялся только на малых размерах массива.
Когда заполнялось константами, программа завершалась ещё до того, как частота вырастала до максимума. Заполнение случайными числами происходит долго, частота успевает подняться до начала измерения времени.
Поэтому в финальных версиях я добавил цикл вокруг вызова исследуемой функции.
Для примера покажу, как может выглядеть вывод, если такой цикл убрать.
Зачем вращаются регистры в цикле
Когда я только начал изучать ассемблерный вывод компилятора, я заметил, что используется вращение регистров в цикле. У меня сразу возник вопрос, зачем это нужно. Ведь до этого я писал на TMS320C66, а там не требовалось ничего вращать.
После обсуждения этого вопроса с @shcher у меня появилось следующее понимание ситуации.
У процессора Эльбрус есть следующее отличие от TMS320C66: при запуске инструкции происходит блокировка (для последующих инструкций) на чтение из регистров, в которые эта инструкция будет писать.
Если попробовать сделать что-то такое (псевдокод):
load R1
add R0, R1, R0то перед инструкцией add произойдет остановка конвейера до момента, когда данные придут в R1.
На первый взгляд это кажется логичным.
Однако, при таком подходе становится неудобно писать конвейеризованные циклы.
Например, посмотрим на такой псевдокод:
load R1
sum:
{
loop sum
add R0, R1, R0
load R1
}
add R0, R1, R0Здесь каждая итерация будет занимать более 1 такта, потому что add, запущенный на текущей итерации, ждёт данных от load, запущенной на предыдущей итерации.
Получается последовательный код, отсутствует конвейеризация.
Для организации конвейера придётся сделать чтение в несколько разных регистров:
load R1
load R2
load R3
sum:
{
add R0, R1, R0
load R1
}
{
add R0, R2, R0
load R2
}
{
loop sum
add R0, R3, R0
load R3
}
add R0, R1, R0
add R0, R2, R0
add R0, R3, R0Теперь происходит 3 сложения за 3 такта (если данные приходят в регистр не дольше, чем через 2 такта после load).
Вращение регистров помогает упростить этот код. Получается примерно так:
load R1
load R2
load R3
sum:
{
loop sum
add R0, R1, R0
load R1
abn
}
add R0, R1, R0
add R0, R2, R0
add R0, R3, R0Здесь инструкция abn вращает имена регистров R1, R2 и R3.
Это эквивалентно перекладыванию данных между ними по кругу:
R1 <-- R2 R2 <-- R3 R3 <-- R1
----------
В TMS320C66 конвейеризованные циклы пишутся проще, потому что у каждой инструкции есть фиксированная длительность исполнения. Например, зафиксировано, что после запуска инструкции load данные придут в R1 в конце 5 такта (первым тактом считаем такт запуска load). Если мы считаем R1 раньше, чем через 5 тактов, то просто прочитается старое значение.
Поэтому можно написать вот такой псевдокод и всё будет работать:
load R1
load R1
load R1
load R1
load R1
sum:
{
loop sum
add R0, R1, R0
load R1
}
add R0, R1, R0
add R0, R1, R0
add R0, R1, R0
add R0, R1, R0
add R0, R1, R0Минусом данного подхода является невозможность запуска программы на железе с отличающимися длительностями инструкций. Если мы хотим ускорить какую-нибудь инструкцию и сохранить обратную совместимость, придётся создать новую "быструю" инструкцию с сохранением старой "медленной".
Такой процесс невозможно продолжать бесконечно.
Вариант, реализованный в Эльбрусе, позволяет программе работать на железе с отличающимися длительностями инструкций. На неродном железе работать будет, быть может, не оптимально. Но будет работать.
----------
Кроме обратной совместимости, вращение регистров может понадобиться в такой задаче:
out[i] = in[i] + in[i-5]
Если решать эту задачу без вращения, получится примерно такой псевдокод:
load R5 // in[i-5]
load R4 // in[i-4]
load R3 // in[i-3]
load R2 // in[i-2]
load R1 // in[i-1]
load R0 // in[i]
{
R10 <-- R0 + R5 // out[i] = in[i] + in[i-5]
load R0
R1 <-- R0
R2 <-- R1
R3 <-- R2
R4 <-- R3
R5 <-- R4
}
sum:
{
loop sum
store R10
R10 <-- R0 + R5 // out[i] = in[i] + in[i-5]
load R0
R1 <-- R0
R2 <-- R1
R3 <-- R2
R4 <-- R3
R5 <-- R4
}
{
store R10
R10 <-- R0 + R5 // out[i] = in[i] + in[i-5]
}
store R10Стрелка обозначает перекладывание содержимого из регистра в регистр какой-нибудь инструкцией. Нам требуется это перекладывание, чтобы не потерять предыдущее значение регистра. Предыдущее значение затирается новым каждую итерацию (каждый такт).
Вращение регистров упрощает псевдокод до такого:
load R5 // in[i-5]
load R4 // in[i-4]
load R3 // in[i-3]
load R2 // in[i-2]
load R1 // in[i-1]
load R0 // in[i]
{
R10 <-- R0 + R5 // out[i] = in[i] + in[i-5]
load R0
abn
}
sum:
{
loop sum
store R10
R10 <-- R0 + R5 // out[i] = in[i] + in[i-5]
load R0
abn
}
{
store R10
R10 <-- R0 + R5 // out[i] = in[i] + in[i-5]
}
store R10Здесь инструкция abn вращает имена регистров R0, R1, R2, R3, R4, R5.
Может быть, в этом описании есть неточности, но идея должна быть понятна.
Заключение
Вместо написания кода на ассемблере можно думать на ассемблере, а писать на Си.
И постоянно сверяться, правильно ли компилятор понял, что я хотел получить.
Я могу писать так, пока не разберусь в ассемблере Эльбруса достаточно хорошо.
После чего можно будет писать сразу на ассемблере.
Компилятор оптимизирует всё подряд и перемешивает код без нужды.
Такой код сложнее понимать, чем написанный человеком вручную.
Наверно, теперь можно подумать о реализации FFT.
Выражаю благодарность сообществу СЭРПАС и лично его руководителю Кириллу Ерохину за предоставление доступа к оборудованию.