Месяц назад мне в телеграм написал человек и предложил доступ к системе с процессором Эльбрус-8СВ.
И, конечно же, я согласился. Так как мне интересно.
Не каждый день неизвестные люди в Интернете предлагают доступ к удалённым хостам.

Меня зовут Леонид Лагунов и я математик-программист.
Я решил написать про своё знакомство с процессором Эльбрус.

Содержание:

• Предшествующий опыт работы с VLIW

• Начинаю изучать Эльбрус-8СВ

• Постановка задачи

• Пишем программу

  • 0. Простейшая реализация

  • 1. Включение режима APB (Array Prefetch Buffer)

  • 2. Раскрутка цикла (x2)

  • 3. Раскрутка цикла (x3)

  • 4. Раскрутка цикла (x4)

  • 5. Использование SIMD64 (movah, pshufb, paddw)

  • 6. SIMD64 раскрутка цикла (x2)

  • 7. SIMD64 раскрутка цикла (x3)

  • 8. SIMD64 раскрутка цикла (x3) + 16-битные "half"-счётчики

  • 9. Использование SIMD128 (movaw, qpshufb, qpaddd)

  • 10. SIMD128 раскрутка цикла (x2*2)

  • 11. SIMD128 раскрутка цикла (x2) + 16-битные "half"-счётчики

  • 12. Использование SIMD64 (2) (movad, psadbw, addd)

  • 13. SIMD64 (2) раскрутка цикла (x2)

  • 14. Использование SIMD128 (2) (movaqp, qpsadbw, qpaddd)

  • 15. SIMD128 (2) раскрутка цикла (x2)

  • 16. SIMD128 (2) раскрутка цикла (x2*2)

  • 17. Использование EML (eml_Vector_Sum_8U)

• Итоговые результаты

• Исследование просадки скорости

• Динамическое изменение частоты процессора

• Зачем вращаются регистры в цикле

• Заключение

Предшествующий опыт работы с VLIW

Ранее мне приходилось работать с системой TI KeyStone II (4 ядра ARM Cortex-A15 и 8 ядер TMS320C66). Для этой системы мной были реализованы алгоритмы FFT-1D (для 2¹⁸ точек) и FFT-2D (для 2⁹ точек), с помощью которых решалась задача поиска картинки в картинке. Управляющий код работал на 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
        }
        ...

Поизучал ассемблерный код. Появились вопросы.

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

Тех.поддержка дала следующие советы:

• использовать int64_t в счётчике цикла

• использовать #pragma loop count(100) для активации режима APB

• использовать интринсики pshufb и paddw для векторных вычислений

• использовать раскрутку цикла для заполнения всех вычислительных юнитов

• использовать #pragma unroll(1) или ключ -fno-loop-unroll (отмена автоматической раскрутки циклов) для упрощения ассемблерного вывода

• использовать EML — там уже всё есть

Постановка задачи

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

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

Дано: указатель на массив и количество элементов в массиве.
Требуется: найти сумму элементов массива.
Элементами массива пусть будут целочисленные байты (uint8_t).

int t0 = clock();
/***  здесь измеряемый код  ***/
int t1 = clock();
int usec = t1 - t0;

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 итераций. В предыдущих вариантах внутренний цикл состоял из 2²⁴ итераций и там активация 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 итераций. В предыдущих вариантах внутренний цикл состоял из 2²⁴ итераций и там активация 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.

Итоговые результаты

Кто хочет посмотреть весь код, вот исходный файл 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.

Выражаю благодарность сообществу СЭРПАС и лично его руководителю Кириллу Ерохину за предоставление доступа к оборудованию.