Война с компилятором и собой: об оптимизациях вещественной арифметики на Эльбрусе24.01.2022 11:02

22f472776219b7880fb46b04caa3632d.jpg

Недавно в процессе выполнения учебного задания мне потребовалось реализовать метод конечных разностей для нахождения приближённого решения краевой задачи. По сути, я впервые столкнулся с вычислениями с плавающей точкой и не мог не попробовать запустить свою программу на Эльбрусе, зная о его больших возможностях и заточенности под вычисления такого рода. Хотите удивиться? Отправляйтесь со мной в увлекательное путешествие!

Вкратце о математике

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

Итак, рассматривается дифференциальное уравнение Пуассона в прямоугольной области [0; 4] x [0; 3]:

-\Delta u + q(x, y)u = F(x, y),\spaceгде \space \Delta u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}

Для данного уравнения рассматривается краевая задача, причём на каждом отрезке границы заданы граничные условия третьего типа:

\frac{\partial u}{\partial n}(x, y) + u(x, y) = \psi(x, y),

где n— единичная внешняя нормаль к границе.

Решение u(x, y)будем искать с помощью разностной схемы, которая может быть представлена в виде системы линейных алгебраических уравнений A\omega = B. Вычисления представляют из себя итерационный метод, на каждом шаге kвычисляется новое приближение сеточной функции \omega^{(k + 1)}_{ij} = \omega^{(k)}_{ij} - \tau_{k + 1} r^{(k)}_{ij}, где невязка r^{(k)} = A \omega^{(k)} - B, а итерационный параметр

\tau_{(k + 1)} = \frac{[Ar^{(k)}, r^{(k)}]}{||Ar^{(k)}||^2}.

Для проверки условия остановки вычисляется норма разности между новым и старым приближением. За более подробным описанием отправляю читателя к профильной литературе, например, книге А.А. Самарского [1].

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

for (int i = start_i; i <= stop_i; ++i) {
    for (int j = start_j; j <= stop_j; ++j) {
        int idx = i * (run_config->domain_n + 2) + j;
        double wij  = src[idx];
        double wipj = src[(i + 1) * (run_config->domain_n + 2) + j];
        double wimj = src[(i - 1) * (run_config->domain_n + 2) + j];
        double wijp = src[i       * (run_config->domain_n + 2) + j + 1];
        double wijm = src[i       * (run_config->domain_n + 2) + j - 1];
        double x = (run_config->start_i + i - 1) * run_config->h1;
        double y = (run_config->start_j + j - 1) * run_config->h2;
        double laplacian = (wipj + wimj - 2 * wij) * run_config->sqinv_h1 + (wijp + wijm - 2 * wij) * run_config->sqinv_h2;
        dst[idx] = q(x, y) * wij - laplacian - alpha * F(x, y);
    }
}

Стоит отметить, что постановка учебной задачи подразумевает реализацию параллельной MPI+OpenMP программы, а также (при желании) дополнение программы MPI+CUDA реализацией. Необходимость разбивки области вычислений на домены вносит небольшую специфику в организацию программы, которую я решил здесь сохранить, несмотря на то, что статья посвящена только оптимизации последовательных вычислений и не касается вопросов распараллеливания на несколько процессов.

Дедлайн был вчера: самая обычная программа

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

Полная реализация лежит на github [3]. Основные шаги, проделанные в статье, можно отследить по истории коммитов. Например, вот начальная версия.

Для тестовых замеров будем смотреть на время выполнения 1000 шагов итерационного метода на одном ядре процессора на сетке размером 1000×1000. За эффектом от оптимизаций будем следить на трёх платформах:

  1. «домашний компьютер» — Intel Core i7–9700K @ 4.8 ГГц и компилятор icc (Intel C++ Compiler Classic) 2021.5.0, базовые флаги оптимизации »-O3 -no-prec-div -ansi-alias -xHost -restrict»;

  2. «целевой вычислительный кластер» — IBM POWER8 @ 4.0 ГГц и компилятор IBM XL C/C++ 13.1.6, базовые флаги оптимизации »-O5» (именно на таких процессорах требовалось запускать задачу; с удовольствием бы взял более новую систему, но доступа к свежим процессорам IBM у меня нет);

  3. «главный герой» — Эльбрус-8СВ @ 1.55 ГГц и компилятор lcc 1.25.19, базовые флаги оптимизации »-O4 -ffast -ffast-math -mtune=elbrus-8c2».

Итак, для тестов на Эльбрусе мы взяли наиболее производительный на данный момент процессор (инженерные образцы 6 поколения архитектуры не рассматриваем). У данного процессора есть аппаратная поддержка циклов, векторные регистры, возможности программной конвейеризации циклов, механизм подкачки массивов (APB) и вообще 6 «двухэтажных» АЛУ для вычислений с плавающей точкой, а мы взяли достаточно свежий компилятор и хорошие опции оптимизации. Несмотря на то, что для векторизации требуется соблюдение выравнивания данных, мы надеемся получить хороший результат.

Смотрим на результаты запуска и немного огорчаемся: Эльбрус отстаёт в 5–16 раз от Intel и IBM, впрочем, последние тоже не радуют скоростью. Здесь можно отправиться читать Руководство по эффективному программированию на платформе «Эльбрус» [2], а я просто воспользуюсь опытом, полученным при реализации криптографических примитивов.

Реализация

i7–9700K

POWER8

Эльбрус-8СВ

№0 — базовая

14.1 с

42.9 с

223.6 с

И всё-таки линейкой по пальцам: подстановка «горячих» функций компилятором

Об этом просто невозможно промолчать: нельзя вызывать простые функции в циклах из других единиц трансляции. Надеюсь, дорогой читатель никогда так не делает, если хочет производительности от своей программы. Да, речь сейчас о q (x, y) и F (x, y), если мы решили вычислять их явно на каждой итерации. Не рассматривая случай какой-нибудь универсальности, эти функции должны быть определены строго в том же файле, где вызываются, чтобы компилятор имел возможность подставить их (inlining) в тело цикла при обычной компиляции (без Link Time Optimization, межпроцедурного анализа в масштабах всей программы или режима -fwhole на Эльбрусе).

Давайте просто посмотрим, как изменится время выполнения программы, если перенести определения функций q (x, y) и F (x, y) из отдельного файла (куда они могли попасть из соображений удобства) поближе к месту их вызова.

Реализация

i7–9700K

POWER8

Эльбрус-8СВ

№0 — базовая

14.1 с

42.9 с

223.6 с

№1 — возможен inline

4.8 c

26.2 c

102.0 c

Примечательно, что на Эльбрусе, который имеет фиксированные и относительно большие штрафы за вызов функций, ускорение составляет 2.2 раза, а на Intel и IBM — 2.9 и 1.5 раза, соответственно. Так что на предсказатель переходов и внеочередное исполнение надейся, а сам не плошай.

Удар ниже пояса: правильный тип для счётчика цикла

Так вышло, что использовать сейчас в качестве счётчика цикла на Эльбрусе любой тип, существенно отличающийся от long — крайне неудачная идея. Связано это с тем, что компилятору необходимо следовать стандарту языка Си, а также с рядом особенностей в реализации компилятора.

В Руководстве нет подробных комментариев, только сказано, что рекомендуется использовать тип long (знаковый). Наверное, идеологически более правильно было бы рекомендовать ssize_t, но тут надо уточнить наличие определения ssize_t на других платформах. Да и с точки зрения компилятора конкретно на Эльбрусе отличий нет: это всё знаковый 64-битный тип.

Что ж, заменяем все границы и счётчики на ssize_t (определение через my_int для удобства экспериментов) и получаем ускорение в 5.4 раза: на 1000 итераций теперь требуется 18.8 секунд вместо 102. В зависимости от ряда факторов, ускорение от подобной замены может быть как почти незаметным, так и более весомым.

Так чем плохи другие типы? Для успешного применения механизма подкачки массивов APB необходимо гарантировать линейность изменения адресов элементов массива в цикле. К сожалению, это сразу отрезает возможность эффективного использования беззнакового типа: в нём допустимо переполнение, а значит на этапе компиляции невозможно гарантировать условие i + 1 > i. Так что привычный многим беззнаковый size_t отпадает. Хотя конкретно в данном примере size_t немного выигрывает, при дальнейших оптимизациях, как правило, проявляется небольшая просадка времени выполнения на беззнаковом типе. Впрочем, разница не так велика и можно продолжать спокойно использовать size_t.

А у меня всё в int, что я делаю не так? Это пример ситуации, когда формально проблемы нет, но есть некоторые сложности в реализации компилятора. Оптимизации применяются на разных стадиях работы компилятора и не всегда успешно согласуются друг с другом. В нашем случае компилятор на более раннем этапе видит цикл и расширяет счётчик до естественного для процессора 64-битного типа. Затем при попытке применить APB внутри цикла происходит вычисление смещения с использованием другой переменной типа int — run_config→domain_n. Получается несогласованность типов переменных, необходимо преобразование, что и приводит к вопиющему снижению производительности. Интересно, что можно вынести из внутреннего цикла чтение domain_n в локальную int переменную и тоже получить ускорение, так как в этом случае расширение переменной до 64 бит происходит в более ранней фазе компиляции.

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

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

Реализация

i7–9700K

POWER8

Эльбрус-8СВ

№1

4.8 c

26.2 c

102.0 c

№1 + локальные int переменные

4.8 c

26.1 c

24.0 с

№1 + size_t счётчики

5.3 с

29.2 с

18.4 с

№2 — ssize_t счётчики

5.3 с

27.9 с

18.8 с

Галопом по Европам: несколько алгоритмических оптимизаций

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

Во-первых, это предвычисления значений функций q (x, y) и F (x, y). На протяжении всего итерационного процесса они не меняются, так что нет смысла тратить время на лишнюю арифметику, тем более, эти функции на практике могут быть более сложными, чем в рассматриваемом примере. Вообще говоря, может оказаться, что считать эти функции будет быстрее, чем загружать значения из памяти, но я не уверен, что такое поведение сохранится после применения дальнейших оптимизаций.

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

В-третьих, это обмен указателей вместо копирования матрицы для перехода к новой итерации. Копирование матрицы вместо обмена указателей увеличивает время работы программы на всех рассматриваемых платформах на 3–8%.

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

Реализация

i7–9700K

POWER8

Эльбрус-8СВ

№2

5.3 с

27.9 с

18.8 с

+ предвычисления q, F

5.98 с

10.72 с

14.23 с

+ меньше проходов по массиву

5.32 с

8.88 с

12.35 с

+ отказ от memcpy — №3

4.93 с

8.42 с

12.03 с

Панки грязи не боятся: читаем на языке ассемблера

Разобраться совсем не так сложно, как может показаться на первый взгляд. Работает универсальный приём: просто берём строчку компиляции нужного файла и заменяем флаг -c на -S, чтобы получить файл с ассемблерным кодом. Не забываем добавить -fverbose-asm, эта опция значительно упрощает чтение на языке ассемблера: для каждой инструкции строится привязка к строке исходного кода, а также включается статический подсчёт тактов (позволяет точно оценить время работы при отсутствии блокировок во время исполнения). Дальше дело за малым: выбираем строчку внутри интересующего нас цикла и находим её в ассемблерном коде по номеру.

Возьмём, например, последний вложенный цикл (из функции update_w_calc_partial_error). Так как компилятор делает ряд оптимизаций по конвейеризации циклов, расщеплению по условиям и создаёт участки компенсирующего кода, то нас прежде всего интересуют вхождения нашей строки в широкие команды (ШК) с меткой loop_mode. Таких оказывается 2, второй выглядит получше, но начнём с первого: 

.L53534:						! solve_cpu.c : 303
	! <0073>
	{
	  loop_mode
	  rbranch	.L63957
	  ldd,0,sm	%dr4, %db[3], %db[22], mas=0x4			! solve_cpu.c : 305
	  fmuld,1,sm	%dr1, %db[26], %db[27]				! solve_cpu.c : 305
	  ldd,2	%dr4, %db[17], %db[36], mas=0x3 ? %pcnt0
	  addd,3,sm	0x8, %db[3], %db[1]				! solve_cpu.c : 303
	  fsubd,4,sm	%db[21], %db[33], %db[38]				! solve_cpu.c : 306
	  ldd,5	%dr3, %db[17], %db[25], mas=0x3 ? %pcnt0
	}
.L63965:
	! <0074>
	{
	  loop_mode
	  alc	alcf=1, alct=1			! solve_cpu.c : 303
	  abn	abnf=1, abnt=1			! solve_cpu.c : 303
	  ct	%ctpr1 ? %NOT_LOOP_END				! solve_cpu.c : 303
	  fmuld,1,sm	%dr0, %db[30], %db[39]				! solve_cpu.c : 307
	  ldd,3,sm	%dr3, %db[9], %db[17], mas=0x4 ? %pcnt4			! solve_cpu.c : 306
	  fmul_addd,4,sm	%db[45], %db[37], %db[20], %db[12]			! solve_cpu.c : 307
	  staad,5	%db[42], %aad0[ %aasti1 ]
	  incr,5	%aaincr0	! solve_cpu.c : 306
	}

Видим, что основная часть тела цикла вроде как занимает всего 2 такта, но постойте-ка: инструкция rbranch делает безусловный переход на метку .L63957, где мы видим ещё минимум 11 тактов: 

.L63957:
	! <0790>
	{
	  nop 3
	  fmuld,0,sm	%dr1, %db[36], %db[37]				! solve_cpu.c : 305
	  fmuld,1,sm	%dr0, %db[36], %db[45]				! solve_cpu.c : 307
	}
	! <0794>
	{
	  nop 5
	  fsubd,0,sm	%db[25], %db[37], %db[42]				! solve_cpu.c : 306
	}
	! <0800>
	{
	  ibranch	.L63965				! solve_cpu.c : 307
	}

В этом отдельном куске только арифметика с выдерживанием задержек от инструкций, так что вернёмся к анализу первой части. Замечаем, что загрузки из памяти делаются инструкцией ldd с тегами mas=0×4 и mas=0×3. Ага, замечательно, это у нас применился динамический механизм разрыва зависимостей с использованием аппаратной таблицы DAM (Memory Access Disambiguator). Секундочку…

Соблюдаем социальную дистанцию: разрыв зависимостей

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

О необходимости разрыва зависимостей очень легко забыть, но нельзя недооценивать важность этой информации для эффективной оптимизации. В случае суперскалярных процессоров эффект может быть заметен не всегда, но он тоже присутствует, как будет видно из замеров. Пожалуй, основное проблемное место — это возможности автовекторизации. Действительно, как можно считать одновременно 2 соседние итерации цикла с помощью векторных регистров, если нет уверенности, что результат следующей итерации не зависит от записи на текущей?

В случае Эльбруса проявляется ещё и другая особенность: из-за статического планирования без информации о независимости итераций мы не можем начать исполнять следующую итерацию до окончания текущей. Это сразу лишает нас важнейшего способа быстрого исполнения циклов в условиях статического планирования — программной конвейеризации цикла.

Так что же делать? Начать можно с простого пути: меньше всего компилятор знает об указателях, которые приходят в качестве параметров функции, поэтому мы можем воспользоваться ключом -frestrict-params, говоря компилятору: «считай, что на всех указателях, являющихся параметрами функций, есть модификатор restrict». Только вот наш случай сложнее: часть указателей для удобства упакована в структуру и их разыменование выполняется в 2 шага (сначала читаем указатель из структуры, потом уже по нему обращаемся к массиву). В результате одного restrict на «внешнем» указателе уже недостаточно. Есть более более мощная комбинация флагов -frestrict-all -frestrict-unsafe, но её всерьёз рассматривать не будем.

В lcc 1.25 впервые появилась поддержка прагмы ivdep, предназначенной для игнорирования возможных зависимостей в цикле. Можно попробовать воспользоваться ею, но я предпочитаю вариант с копированием нужных указателей и явным добавлением модификатора restrict. На мой взгляд, такой вариант самый прозрачный и даже немного упрощает код программы. Вот так, например, будет выглядеть начало одной из 3 подправленных функций: 

double update_w_calc_partial_error(const struct RunConfig *run_config, double tau)
{
    double * restrict cur_w = run_config->cur_w;
    double * restrict next_w = run_config->next_w;
    double * restrict residual = run_config->residual;
    // <...>
}

В результате небольшой модификации мы добились ускорения на 33–80% на всех платформах. Подчеркну, что наша оптимизация является общей, а не специфической для Эльбруса. Более того, наиболее выраженный эффект от её применения наблюдается на IBM.

Реализация

i7–9700K

POWER8

Эльбрус-8СВ

№3

4.93 с

8.42 с

12.03 с

№4 — используем restrict

3.69 с

4.69 с

7.48 с

И проверим, как же теперь выглядит исследуемый цикл на языке ассемблера. Отлично, все вхождения цикла стали компактными, раскрученными на 2 и спланированными в 4 такта на тело (то есть 2 такта на итерацию):

Цикл на языке ассемблера

.L38838:						! solve_cpu.c : 312
	! <0237>
	{
	  loop_mode
	  fmuld,5,sm	%dr7, %db[11], %db[15]				! solve_cpu.c : 314
	  movad,0	area=0, ind=0, am=0, be=0, %db[0]	! solve_cpu.c : 315
	  movad,1	area=0, ind=8, am=1, be=0, %db[1]	! solve_cpu.c : 315
	  movad,3	area=0, ind=0, am=0, be=0, %db[8]	! solve_cpu.c : 314
	}
	! <0238>
	{
	  loop_mode
	  fmuld,3,sm	%dr0, %db[11], %db[20]				! solve_cpu.c : 316
	  fmuld,4,sm	%dr7, %db[12], %db[16]				! solve_cpu.c : 314
	  fmuld,5,sm	%dr0, %db[12], %db[21]				! solve_cpu.c : 316
	}
	! <0239>
	{
	  loop_mode
	  fmuld,3,sm	%db[22], %db[17], %db[11]				! solve_cpu.c : 316
	  fmuld,4,sm	%db[23], %db[18], %db[24]				! solve_cpu.c : 316
	  fsubd,5,sm	%db[7], %db[17], %db[12]				! solve_cpu.c : 315
	}
	! <0240>
	{
	  loop_mode
	  alc	alcf=1, alct=1			! solve_cpu.c : 312
	  abn	abnf=1, abnt=1			! solve_cpu.c : 312
	  ct	%ctpr1 ? %NOT_LOOP_END				! solve_cpu.c : 312
	  fsubd,1,sm	%db[6], %db[18], %db[17]				! solve_cpu.c : 315
	  staad,2	%db[19], %aad2[ %aasti1 ]		! solve_cpu.c : 315
	  faddd,3,sm	%dr3, %db[13], %dr3				! solve_cpu.c : 316
	  faddd,4,sm	%dr1, %db[26], %dr1				! solve_cpu.c : 316
	  staad,5	%db[14], %aad2[ %aasti1 + _f32s,_lts0 0x8 ]
	  incr,5	%aaincr2	! solve_cpu.c : 315
	  movad,3	area=0, ind=8, am=1, be=0, %db[7]	! solve_cpu.c : 314
	}

От работы кони дохнут: вынос инварианта из цикла

Обратимся теперь к циклу из функции calc_tau_part: 

for (my_int j = 2; j <= run_config->domain_n - 1; ++j) {
    my_int idx = i * (run_config->domain_n + 2) + j;
    num += rhox * a_residual[idx] * residual[idx];
    div += rhox * a_residual[idx] * a_residual[idx];
}

Легко видеть, что здесь есть 2 умножения на величину rhox, значение которой не меняется в цикле. Её можно вынести за пределы цикла, так как мы изначально допускаем нарушение порядка арифметических операций. Конечно, это всего лишь одно умножение, но важна сама идея, так как в других примерах легко можно встретить вычисление в цикле, скажем, синуса от постоянной величины. Примечательно, что в документации к IBM XLC говорится о том, что компилятор умеет сам делать вынос инварианта. Что ж, замеры показывают, что надёжнее это сделать руками:

Реализация

i7–9700K

POWER8

Эльбрус-8СВ

№4

3.69 с

4.69 с

7.48 с

№5 — вынос инварианта

3.67 с

4.37 с

7.08 с

Эффект на Эльбрусе хорошо иллюстрируется: теперь цикл спланирован в 1 такт без раскрутки вместо 4 тактов при раскрутке на 2.

Для любителей ассемблера

Было: 

.L19038:						! solve_cpu.c : 270
	! <0259>
	{
	  loop_mode
	  movad,2	area=0, ind=8, am=1, be=0, %db[1]	! solve_cpu.c : 272
	  movad,3	area=0, ind=0, am=0, be=0, %db[0]	! solve_cpu.c : 272
	}
	! <0260>
	{
	  loop_mode
	  fmuld,4,sm	%dr5, %db[3], %db[15]				! solve_cpu.c : 272
	  fmuld,5,sm	%dr5, %db[2], %db[12]				! solve_cpu.c : 272
	  movad,0	area=0, ind=8, am=1, be=0, %db[6]	! solve_cpu.c : 272
	  movad,1	area=0, ind=0, am=0, be=0, %db[9]	! solve_cpu.c : 272
	}
	! <0261>
	{
	  loop_mode
	  fmuld,3,sm	%db[17], %db[10], %db[19]				! solve_cpu.c : 272
	  fmuld,4,sm	%db[17], %db[5], %db[16]				! solve_cpu.c : 273
	  fmuld,5,sm	%db[14], %db[13], %db[20]				! solve_cpu.c : 272
	}
	! <0262>
	{
	  loop_mode
	  alc	alcf=1, alct=1			! solve_cpu.c : 270
	  abn	abnf=1, abnt=1			! solve_cpu.c : 270
	  ct	%ctpr1 ? %NOT_LOOP_END				! solve_cpu.c : 270
	  fmuld,1,sm	%db[14], %db[4], %db[5]				! solve_cpu.c : 273
	  faddd,2,sm	%dr0, %db[7], %dr0				! solve_cpu.c : 273
	  faddd,3,sm	%dr1, %db[21], %dr1				! solve_cpu.c : 272
	  faddd,4,sm	%dr3, %db[18], %dr3				! solve_cpu.c : 273
	  faddd,5,sm	%dr2, %db[22], %dr2				! solve_cpu.c : 272
	}

Стало: 

.L18685:						! solve_cpu.c : 272
	! <0083>
	{
	  loop_mode
	  alc	alcf=1, alct=1			! solve_cpu.c : 272
	  abn	abnf=1, abnt=1			! solve_cpu.c : 272
	  ct	%ctpr1 ? %NOT_LOOP_END				! solve_cpu.c : 272
	  fmuld,1,sm	%db[37], %db[36], %db[39]				! solve_cpu.c : 274
	  faddd,2,sm	%db[25], %db[47], %db[17]				! solve_cpu.c : 274
	  fmuld,4,sm	%db[37], %db[37], %db[38]				! solve_cpu.c : 275
	  faddd,5,sm	%db[24], %db[46], %db[16]				! solve_cpu.c : 275
	  movad,1	area=0, ind=0, am=1, be=0, %db[26]	! solve_cpu.c : 274
	  movad,3	area=0, ind=0, am=1, be=0, %db[27]	! solve_cpu.c : 274
	}

Ну, кошечка, ну ещё капельку: добиваемся векторизации

Существенное ограничение на Эльбрус-8СВ при работы с векторными регистрами — данные должны быть выровнены (по 16-байтной границе). Без этого компилятор не сможет сгенерировать код с использованием 128-битных регистров. А вот как получить векторный код при соблюдении необходимых требований — «а вот это науке ещё не известно!» Точнее, рецепт примерно стандартный: в текущей версии компилятору нужна информация о количестве итераций цикла, чтобы «решиться» на векторизацию. Эта информация может быть получена из подсказки с помощью прагмы loop count (N) или из профиля исполнения программы. Тем не менее, я пока не выработал для себя удобного способа работы, так что покажу способ для продвинутых: как с помощью intrinsic-функций гарантированно получить векторный код.

Основные усилия мы направим на цикл в функции calc_aw_b. Его особенность в том, что на каждой итерации нужны 3 подряд идущих элемента массива. Соответственно, выровненный доступ к памяти здесь просто так невозможен. Для обхода этой особенности на итерации с номером j будем формировать регистр, содержащий сдвинутое на 1 значение (как бы src[j + 1] элемент), с помощью ассемблерной инструкции qppermb из считанных выровненных 128-битных значений src[j] и src[j + 2]. То есть из src[j] мы берём старшую половину, а из src[j + 2] — младшую.

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

for (my_int j = start_j; j < stop_j; j += 2) {
    my_int idx = i * (run_config->domain_n + 2) + j;
    __v2di wij = *((__v2di *) &src[idx]);
    __v2di wipj = *((__v2di *) &src[(i + 1) * (run_config->domain_n + 2) + j]);
    __v2di wimj = *((__v2di *) &src[(i - 1) * (run_config->domain_n + 2) + j]);
    
    src_0 = wij;
    src_1 = *((__v2di *) &src[idx + 2]);
    __v2di wijp = __builtin_e2k_qppermb(src_1, src_0, mix_doubles);
    
    __v2di t0 = __builtin_e2k_qpfaddd(wipj, wimj);
    __v2di t1 = __builtin_e2k_qpfmuld(const_2, wij);
    __v2di t2 = __builtin_e2k_qpfaddd(wijp, wijm);
    __v2di t3 = __builtin_e2k_qpfsubd(t0, t1);
    __v2di t4 = __builtin_e2k_qpfsubd(t2, t1);
    __v2di t5 = __builtin_e2k_qpfmuld(t3, sqinv_h1);
    __v2di t6 = __builtin_e2k_qpfmuld(t4, sqinv_h2);
    __v2di laplacian = __builtin_e2k_qpfaddd(t5, t6);
    __v2di t7 = __builtin_e2k_qpfmuld(wij, *((__v2di *) &q_mat[idx]));
    __v2di t8 = __builtin_e2k_qpfmuld(v2alpha, *((__v2di *) &b_mat[idx]));
    *((__v2di *) &dst[idx]) = __builtin_e2k_qpfsubd(__builtin_e2k_qpfsubd(t7, laplacian), t8);
    wijm = wijp;
}

Выглядит уже довольно громоздко, а это ещё не видна обработка головы и хвоста цикла. Зато компилируется в красивую конструкцию с итоговым темпом обработки 1 такт на элемент матрицы.

.L15436:						! solve_cpu.c : 234
	! <0247>
	{
	  loop_mode
	  qpfadd_rsubd,0,sm	%xb[38], %xb[40], %xb[95], %xb[40]			! solve_cpu.c : 248
	  qpfadd_rsubd,1,sm	%xb[82], %xb[92], %xb[97], %xb[1]			! solve_cpu.c : 247
	  qpfmuld,2,sm	%xb[48], %xr5, %xb[56]				! solve_cpu.c : 250
	  qpfmul_addd,4,sm	%xb[13], %xr6, %xb[62], %xb[48]			! solve_cpu.c : 251
	  qpfmul_rsubd,5,sm	%xb[32], %xb[89], %xb[56], %xb[13]			! solve_cpu.c : 254
	  movaqp,1	area=1, ind=0, am=1, be=0, %xb[0]	! solve_cpu.c : 242
	}
	! <0248>
	{
	  loop_mode
	  alc	alcf=1, alct=1			! solve_cpu.c : 234
	  abn	abnf=1, abnt=1			! solve_cpu.c : 234
	  ct	%ctpr1 ? %NOT_LOOP_END				! solve_cpu.c : 234
	  qppermb,1,sm	%xb[6], %xb[8], %xr0, %xb[36] ? %pcnt21			! solve_cpu.c : 242
	  qpfmuld,2,sm	%xr2, %xb[4], %xb[89]				! solve_cpu.c : 245
	  qpfmul_subd,4,sm	%xr7, %xb[59], %xb[21], %xb[62]			! solve_cpu.c : 254
	  staaqp,5	%xb[72], %aad5[ %aasti1 ]
	  incr,5	%aaincr0	! solve_cpu.c : 254
	  movaqp,0	area=0, ind=0, am=1, be=0, %xb[21]	! solve_cpu.c : 253
	  movaqp,1	area=2, ind=0, am=1, be=0, %xb[72]	! solve_cpu.c : 237
	  movaqp,2	area=0, ind=0, am=1, be=0, %xb[59]	! solve_cpu.c : 252
	  movaqp,3	area=1, ind=0, am=1, be=0, %xb[82]	! solve_cpu.c : 238
	}

Здесь хорошо видно, что компилятор успешно объединяет последовательные инструкции в сдвоенные — это позволяет более плотно планировать их исполнение (без «второго этажа» 6 АЛУ не хватило бы, чтобы спланировать цикл в 2 такта на тело). И работает немного быстрее: удаётся улучшить результат на 14% и достичь 6.19 секунд (в таблицах этот шаг будет обозначен как »№6 — векторизация»). Заметим, что на рассматриваемых процессорах Intel и IBM нет требований к выравниванию данных, а потому компиляторы ещё на предыдущих этапах смогли справиться с векторизацией всех циклов.

Ломаем законы физики: линейная обработка против объёма вычислений

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

Так зачем мы всё это делаем во внешнем цикле, если по факту двумерный массив представлен в памяти линейным участком? Заменяем двойной цикл (с количеством итераций не больше M для внешнего и N для вложенного) на одинарный с общим количеством итераций MN + 2M — 2 и поднимаем его в начало функции. Граничные точки будут перезаписаны ниже при обработке границ. Да, мы в этом цикле делаем больше вычислений, чем было, и даже что-то считаем для точек, которые нам вообще не нужны (используются при межпроцессных синхронизациях), но избавление от вложенного цикла даёт выигрыш. Напрашивается аналогия с работой кэша: при переходе к линейному доступу к памяти можно делать больше операций, но выигрыш всё равно останется.

После перехода к одному большому циклу на Эльбрусе можно проверить поведение программы при выборе параметра unroll для этого цикла. Обычно требуется посмотреть всего несколько значений (до 4 или 8), в данном случае оптимальным оказалась раскрутка на 3 итерации. На Intel и IBM пользы от ручного выбора параметра раскрутки не обнаружилось.

Реализация

i7–9700K

POWER8

Эльбрус-8СВ

№6 — векторизация

3.67 с

4.37 с

6.19 с

№7 — коллапсирование гнезда циклов

3.60 с

4.25 с

5.46 с

Последний штрих: процессор быстрее памяти

Остался один очевидный шаг, который можно было бы сделать ещё на этапе алгоритмических оптимизаций, но я сознательно пропустил его в угоду унификации. В функции calc_aw_b параметр alpha принимает значения 0 или 1 и отвечает за необходимость вычитания правой части при вычислении AX - \alpha B. Эти 2 случая можно было бы сразу разделить, но я предлагаю взглянуть на проблему под другим углом.

Основные циклы на языке ассемблера сейчас выглядят так: 

.L13023:						! solve_cpu.c : 168
	! <0052>
	{
	  loop_mode
	  qpfadd_rsubd,1,sm	%xb[61], %xb[36], %xb[95], %xb[15]			! solve_cpu.c : 180
	  qpfmsd,3,sm	%xb[15], %xb[87], %xb[91], %xb[103]			! solve_cpu.c : 187
	  qpfnmad,4,sm	%xr11, %xb[50], %xb[108], %xb[104]			! solve_cpu.c : 187
	  qpfmuld,5,sm	%xb[65], %xr0, %xb[105]				! solve_cpu.c : 183
	}
	! <0053>
	{
	  loop_mode
	  qpfadd_rsubd,1,sm	%xb[110], %xb[6], %xb[68], %xb[29]			! solve_cpu.c : 181
	  qpfmuld,2,sm	%xr10, %xb[35], %xb[16]				! solve_cpu.c : 178
	  qpfmsd,3,sm	%xb[43], %xb[29], %xb[109], %xb[106]			! solve_cpu.c : 187
	  qppermb,4,sm	%xb[5], %xb[35], %xr12, %xb[4] ? %pcnt5			! solve_cpu.c : 175
	}
	! <0054>
	{
	  loop_mode
	  qpfadd_rsubd,1,sm	%xb[67], %xb[110], %xb[93], %xb[61]			! solve_cpu.c : 181
	  qpfmuld,2,sm	%xr10, %xb[5], %xb[66]				! solve_cpu.c : 178
	  qpfmad,3,sm	%xb[66], %xr9, %xb[111], %xb[107]			! solve_cpu.c : 184
	  qppermb,4,sm	%xb[33], %xb[90], %xr12, %xb[65]			! solve_cpu.c : 175
	  movaqp,0	area=0, ind=0, am=0, be=0, %xb[43]	! solve_cpu.c : 186
	  movaqp,1	area=0, ind=16, am=1, be=0, %xb[50]	! solve_cpu.c : 186
	  movaqp,3	area=2, ind=0, am=1, be=0, %xb[36]	! solve_cpu.c : 186
	}
	! <0055>
	{
	  loop_mode
	  qpfmuld,2,sm	%xb[31], %xr0, %xg16				! solve_cpu.c : 183
	  qpfmad,3,sm	%xb[76], %xr9, %xg16, %xb[87]			! solve_cpu.c : 184
	  qppermb,4,sm	%xb[90], %xb[5], %xr12, %xb[108]			! solve_cpu.c : 175
	  qpfmuld,5,sm	%xb[75], %xr0, %xb[109]				! solve_cpu.c : 183
	  movaqp,0	area=1, ind=0, am=0, be=0, %xb[88]	! solve_cpu.c : 175
	  movaqp,1	area=1, ind=16, am=1, be=0, %xb[31]	! solve_cpu.c : 175
	  movaqp,2	area=0, ind=0, am=0, be=0, %xb[75]	! solve_cpu.c : 185
	  movaqp,3	area=0, ind=16, am=1, be=0, %xb[76]	! solve_cpu.c : 185
	}
	! <0056>
	{
	  loop_mode
	  qpfadd_rsubd,1,sm	%xb[30], %xb[24], %xb[18], %xb[62]			! solve_cpu.c : 180
	  qpfmuld,2,sm	%xr10, %xb[90], %xb[91]				! solve_cpu.c : 178
	  qpfmad,3,sm	%xb[17], %xr9, %xb[105], %xb[95]			! solve_cpu.c : 184
	  qpfnmad,4,sm	%xr11, %xb[62], %xb[100], %xb[99]			! solve_cpu.c : 187
	  staaqp,5	%xb[102], %aad5[ %aasti1 ]		! solve_cpu.c : 187
	  movaqp,0	area=3, ind=0, am=1, be=0, %xb[17]	! solve_cpu.c : 185
	  movaqp,1	area=4, ind=0, am=1, be=0, %xb[18]	! solve_cpu.c : 171
	  movaqp,2	area=4, ind=0, am=1, be=0, %xb[24]	! solve_cpu.c : 170
	  movaqp,3	area=1, ind=16, am=0, be=0, %xb[30]	! solve_cpu.c : 171
	}
	! <0057>
	{
	  loop_mode
	  alc	alcf=1, alct=1			! solve_cpu.c : 168
	  abn	abnf=1, abnt=1			! solve_cpu.c : 168
	  ct	%ctpr1 ? %NOT_LOOP_END				! solve_cpu.c : 168
	  qpfadd_rsubd,0,sm	%xb[72], %xb[71], %xb[68], %xb[72]			! solve_cpu.c : 180
	  qpfadd_rsubd,1,sm	%xb[4], %xb[67], %xb[16], %xb[71]			! solve_cpu.c : 181
	  staaqp,2	%xb[101], %aad5[ %aasti1 + _f32s,_lts1 0x10 ]	! solve_cpu.c : 187
	  qpfmsd,3,sm	%xb[98], %xb[86], %xb[97], %xb[98]			! solve_cpu.c : 187
	  qpfnmad,4,sm	%xr11, %xb[55], %xb[103], %xb[100]			! solve_cpu.c : 187
	  staaqp,5	%xb[104], %aad5[ %aasti1 + _f32s,_lts0 0x20 ]
	  incr,5	%aaincr2	! solve_cpu.c : 187
	  movaqp,0	area=2, ind=16, am=1, be=0, %xb[55]	! solve_cpu.c : 170
	  movaqp,1	area=2, ind=0, am=0, be=0, %xb[68]	! solve_cpu.c : 170
	  movaqp,2	area=1, ind=0, am=1, be=0, %xb[67]	! solve_cpu.c : 171
	  movaqp,3	area=3, ind=0, am=1, be=0, %xb[1]	! solve_cpu.c : 175
	}


.L18027:						! solve_cpu.c : 321
	! <0064>
	{
	  loop_mode
	}
	! <0065>
	{
	  loop_mode
	  movaqp,1	area=0, ind=16, am=0, be=0, %xb[1]	! solve_cpu.c : 324
	  movaqp,2	area=0, ind=0, am=0, be=0, %xb[0]	! solve_cpu.c : 323
	  movaqp,3	area=0, ind=16, am=1, be=0, %xb[4]	! solve_cpu.c : 323
	}
	! <0066>
	{
	  loop_mode
	  qpfmuld,3,sm	%xb[2], %xb[2], %xb[12]				! solve_cpu.c : 326
	  qpfmuld,4,sm	%xb[6], %xb[3], %xb[11]				! solve_cpu.c : 325
	  qpfmuld,5,sm	%xb[6], %xb[6], %xb[8]				! solve_cpu.c : 326
	  movaqp,1	area=0, ind=0, am=1, be=0, %xb[7]	! solve_cpu.c : 324
	}
	! <0067>
	{
	  loop_mode
	  alc	alcf=1, alct=1			! solve_cpu.c : 321
	  abn	abnf=1, abnt=1			! solve_cpu.c : 321
	  ct	%ctpr1 ? %NOT_LOOP_END				! solve_cpu.c : 321
	  qpfmuld,1,sm	%xb[2], %xb[9], %xb[3]				! solve_cpu.c : 325
	  qpfaddd,2,sm	%xr0, %xb[5], %xr0				! solve_cpu.c : 325
	  qpfaddd,3,sm	%xr2, %xb[14], %xr2				! solve_cpu.c : 326
	  qpfaddd,4,sm	%xr3, %xb[13], %xr3				! solve_cpu.c : 325
	  qpfaddd,5,sm	%xr1, %xb[10], %xr1				! solve_cpu.c : 326
	}


.L33648:						! solve_cpu.c : 391
	! <0067>
	{
	  loop_mode
	  movaqp,2	area=0, ind=0, am=0, be=0, %xb[1]	! solve_cpu.c : 393
	  movaqp,3	area=0, ind=16, am=1, be=0, %xb[0]	! solve_cpu.c : 393
	}
	! <0068>
	{
	  loop_mode
	  qpfmuld,4,sm	%xr5, %xb[3], %xb[15]				! solve_cpu.c : 393
	  qpfmuld,5,sm	%xr5, %xb[2], %xb[12]				! solve_cpu.c : 393
	  movaqp,0	area=0, ind=16, am=1, be=0, %xb[6]	! solve_cpu.c : 394
	  movaqp,1	area=0, ind=0, am=0, be=0, %xb[7]	! solve_cpu.c : 394
	}
	! <0069>
	{
	  loop_mode
	  qpfmuld,3,sm	%xb[17], %xb[17], %xb[3]				! solve_cpu.c : 395
	  qpfmuld,4,sm	%xb[14], %xb[14], %xb[2]				! solve_cpu.c : 395
	  qpfsubd,5,sm	%xb[11], %xb[17], %xb[16]				! solve_cpu.c : 394
	}
	! <0070>
	{
	  loop_mode
	  alc	alcf=1, alct=1			! solve_cpu.c : 391
	  abn	abnf=1, abnt=1			! solve_cpu.c : 391
	  ct	%ctpr1 ? %NOT_LOOP_END				! solve_cpu.c : 391
	  qpfsubd,1,sm	%xb[10], %xb[14], %xb[11]				! solve_cpu.c : 394
	  staaqp,2	%xb[13], %aad2[ %aasti1 + _f32s,_lts0 0x10 ]	! solve_cpu.c : 394
	  qpfaddd,3,sm	%xr3, %xb[5], %xr3				! solve_cpu.c : 395
	  qpfaddd,4,sm	%xr4, %xb[4], %xr4				! solve_cpu.c : 395
	  staaqp,5	%xb[18], %aad2[ %aasti1 ]
	  incr,5	%aaincr2	! solve_cpu.c : 394
	}

Видно, что каждый цикл спланирован так, что на обработку одного элемента матрицы в каждом цикле должен тратиться 1 такт. Если же посмотреть на результаты профилировки с помощью утилиты perf, то оказывается, что цикл в функции calc_aw_b исполняется примерно в 1.5 раза дольше, чем другие 2 (функция calc_aw_b вызывается в 2 раза чаще, чем остальные): 

Samples: 7K of event 'task-clock:u', Event count (approx.): 7103000000
  Children      Self  Command     Shared Object           Symbol
+   59.75%    59.75%  prog_e2kv5  prog_e2kv5              [.] calc_aw_b
+   20.37%    20.37%  prog_e2kv5  prog_e2kv5              [.] update_w_calc_partial_error
+   18.16%    18.16%  prog_e2kv5  prog_e2kv5              [.] calc_tau_part

Этот эффект объясняется тем, что данные не успевают подгружаться из оперативной памяти с достаточным для процессора темпом. Для эксперимента разделим случаи умножения на 1 и на 0 и в случае умножения на 0 будем не загружать значение матрицы B, а брать какое-нибудь из уже использованных. Планирование цикла принципиально не поменялось, это те же 6 тактов, хоть в нём и на 3 инструкции movaqp меньше.

Аналогичное планирование: 

.L13226:						! solve_cpu.c : 194
	! <0053>
	{
	  loop_mode
	  qpfadd_rsubd,1,sm	%xb[69], %xb[79], %xb[48], %xb[28]			! solve_cpu.c : 206
	  qpfmsd,3,sm	%xb[15], %xb[28], %xb[94], %xb[94]			! solve_cpu.c : 213
	  qpfnmad,4,sm	%xr12, %xb[14], %xb[101], %xb[97]			! solve_cpu.c : 213
	  qpfmuld,5,sm	%xb[90], %xr0, %xb[98]				! solve_cpu.c : 209
	  movaqp,0	area=2, ind=0, am=1, be=0, %xb[15]	! solve_cpu.c : 211
	  movaqp,1	area=0, ind=0, am=0, be=0, %xb[16]	! solve_cpu.c : 211
	}
	! <0054>
	{
	  loop_mode
	  qpfadd_rsubd,1,sm	%xb[76], %xb[6], %xb[93], %xb[72]			! solve_cpu.c : 207
	  qpfmuld,2,sm	%xr9, %xb[45], %xb[69]				! solve_cpu.c : 204
	  qpfmsd,3,sm	%xb[53], %xb[27], %xb[102], %xb[99]			! solve_cpu.c : 213
	  qppermb,4,sm	%xb[5], %xb[45], %xr11, %xb[4] ? %pcnt5			! solve_cpu.c : 201
	  movaqp,0	area=0, ind=16, am=1, be=0, %xb[27]	! solve_cpu.c : 211
	  movaqp,1	area=3, ind=0, am=1, be=0, %xb[48]	! solve_cpu.c : 197
	  movaqp,2	area=3, ind=0, am=1, be=0, %xb[53]	! solve_cpu.c : 196
	  movaqp,3	area=1, ind=16, am=0, be=0, %xb[63]	! solve_cpu.c : 196
	}
	! <0055>
	{
	  loop_mode
	  qpfadd_rsubd,1,sm	%xb[36], %xb[76], %xb[46], %xb[86]			! solve_cpu.c : 207
	  qpfmuld,2,sm	%xr9, %xb[5], %xb[91]				! solve_cpu.c : 204
	  qpfmad,3,sm	%xb[34], %xr10, %xb[103], %xb[100]			! solve_cpu.c : 210
	  qppermb,4,sm	%xb[43], %xb[62], %xr11, %xb[34]			! solve_cpu.c : 201
	  movaqp,0	area=1, ind=16, am=1, be=0, %xb[73]	! solve_cpu.c : 197
	  movaqp,1	area=1, ind=0, am=0, be=0, %xb[79]	! solve_cpu.c : 197
	  movaqp,2	area=1, ind=0, am=1, be=0, %xb[85]	! solve_cpu.c : 196
	  movaqp,3	area=2, ind=0, am=1, be=0, %xb[1]	! solve_cpu.c : 201
	}
	! <0056>
	{
	  loop_mode
	  qpfmuld,2,sm	%xb[74], %xr0, %xg16				! solve_cpu.c : 209
	  qpfmad,3,sm	%xb[60], %xr10, %xg16, %xb[90]			! solve_cpu.c : 210
	  qppermb,4,sm	%xb[62], %xb[5], %xr11, %xb[74]			! solve_cpu.c : 201
	  qpfmuld,5,sm	%xb[41], %xr0, %xb[101]				! solve_cpu.c : 209
	  movaqp,2	area=0, ind=0, am=0, be=0, %xb[60]	! solve_cpu.c : 201
	  movaqp,3	area=0, ind=16, am=1, be=0, %xb[41]	! solve_cpu.c : 201
	}
	! <0057>
	{
	  loop_mode
	  qpfadd_rsubd,1,sm	%xb[59], %xb[54], %xb[71], %xb[30]			! solve_cpu.c : 206
	  qpfmuld,2,sm	%xr9, %xb[62], %xb[44]				! solve_cpu.c : 204
	  qpfmad,3,sm	%xb[30], %xr10, %xb[98], %xb[54]			! solve_cpu.c : 210
	  qpfnmad,4,sm	%xr12, %xb[44], %xb[95], %xb[59]			! solve_cpu.c : 213
	  staaqp,5	%xb[96], %aad4[ %aasti1 ]		! solve_cpu.c : 213
	}
	! <0058>
	{
	  loop_mode
	  alc	alcf=1, alct=1			! solve_cpu.c : 194
	  abn	abnf=1, abnt=1			! solve_cpu.c : 194
	  ct	%ctpr1 ? %NOT_LOOP_END				! solve_cpu.c : 194
	  qpfadd_rsubd,0,sm	%xb[89], %xb[83], %xb[93], %xb[56]			! solve_cpu.c : 206
	  qpfadd_rsubd,1,sm	%xb[4], %xb[36], %xb[69], %xb[37]			! solve_cpu.c : 207
	  staaqp,2	%xb[61], %aad4[ %aasti1 + _f32s,_lts1 0x10 ]	! solve_cpu.c : 213
	  qpfmsd,3,sm	%xb[70], %xb[37], %xb[56], %xb[93]			! solve_cpu.c : 213
	  qpfnmad,4,sm	%xr12, %xb[84], %xb[94], %xb[94]			! solve_cpu.c : 213
	  staaqp,5	%xb[97], %aad4[ %aasti1 + _f32s,_lts0 0x20 ]
	  incr,5	%aaincr2	! solve_cpu.c : 213
	}

А вот скорость вычислений выросла на 8%, что подтверждает наличие зависимости времени вычисления от скорости работы памяти. Остаётся честно избавиться от умножения на alpha в обоих случаях и получить итоговый результат. Ниже сводная таблица со всеми основными результатами.

Реализация

i7–9700K

POWER8

Эльбрус-8СВ

№0 — базовая

14.1 с

42.9 с

223.6 с

№1 — возможен inline

4.8 c

26.2 c

102.0 c

№2 — ssize_t счётчики

5.3 с

27.9 с

18.8 с

№3 — алгоритмические оптимизации

4.93 с

8.42 с

12.03 с

№4 — используем restrict

3.69 с

4.69 с

7.48 с

№5 — вынос инварианта

3.67 с

4.37 с

7.0

© Habrahabr.ru