Недавно в процессе выполнения учебного задания мне потребовалось реализовать метод конечных разностей для нахождения приближённого решения краевой задачи. По сути, я впервые столкнулся с вычислениями с плавающей точкой и не мог не попробовать запустить свою программу на Эльбрусе, зная о его больших возможностях и заточенности под вычисления такого рода. Хотите удивиться? Отправляйтесь со мной в увлекательное путешествие!
Вкратце о математике
Для понимания происходящего математическая постановка задачи не является принципиально важной, так что её подробный разбор я позволю себе пропустить и ограничусь коротким описанием.
Итак, рассматривается дифференциальное уравнение Пуассона в прямоугольной области [0; 4] x [0; 3]:
Для данного уравнения рассматривается краевая задача, причём на каждом отрезке границы заданы граничные условия третьего типа:
где 
— единичная внешняя нормаль к границе.
Решение 
будем искать с помощью разностной схемы, которая может быть представлена в виде конечно-разностной операторной задачи
. Вычисления представляют из себя итерационный метод, на каждом шаге 
вычисляется новое приближение сеточной функции 
, где невязка 
, а итерационный параметр
![\tau_{(k + 1)} = \frac{[Ar^{(k)}, r^{(k)}]}{||Ar^{(k)}||^2}.](https://habrastorage.org/getpro/habr/upload_files/81f/3e2/474/81f3e24743e00341f1e9873b152803b0.svg)
Для проверки условия остановки вычисляется норма разности между новым и старым приближением. За более подробным описанием отправляю читателя к профильной литературе, например, книге А.А. Самарского [1].
Формулы аппроксимации (по просьбам трудящихся)
Исходное уравнение во внутренних точках прямоугольной области аппроксимируется разностным уравнением
где 
а разностный оператор Лапласа
Это аппроксимация второго порядка точности. Если сделать для граничных точек аналогично (наивно), то получится аппроксимация первого порядка точности (что испортит точность всей схемы), поэтому применяется более хитрая схема. Рассмотрим на примере правой границы:
С помощью ряда Тейлора выразим производную:
Переходим к сеточной функции, подставляем в граничное условие (вторая производная выражается через исходное уравнение с оператором Лапласа) и домножаем всё выражение на 
получаем:
Здесь 
для правой границы. Аналогичные выражения получаются для остальных границ (можно увидеть в явном виде в исходном коде на github).
В результате основные вычисления на каждом шаге представляют из себя несколько запусков пары вложенных циклов с классической характерной структурой: несколько чтений данных из матриц, арифметические операции и запись в выходную матрицу и/или накопление результата в промежуточной переменной. С производительностью исполнения этих циклов мы и будем разбираться. Вот так выглядит один из них:
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 шагов итерационного метода на одном ядре процессора на сетке размером 1000x1000. За эффектом от оптимизаций будем следить на трёх платформах:
«домашний компьютер» — Intel Core i7-9700K @ 4.8 ГГц и компилятор icc (Intel C++ Compiler Classic) 2021.5.0, базовые флаги оптимизации «-O3 -no-prec-div -ansi-alias -xHost -restrict»;
«целевой вычислительный кластер» — IBM POWER8 @ 4.0 ГГц и компилятор IBM XL C/C++ 13.1.6, базовые флаги оптимизации «-O5» (именно на таких процессорах требовалось запускать задачу; с удовольствием бы взял более новую систему, но доступа к свежим процессорам IBM у меня нет);
«главный герой» — Эльбрус-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 на других платформах (позже мне подсказали, что лучше всего использовать ptrdiff_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=0x4 и mas=0x3. Ага, замечательно, это у нас применился динамический механизм разрыва зависимостей с использованием аппаратной таблицы 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 и отвечает за необходимость вычитания правой части при вычислении 
. Эти 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.08 с |
№6 — ручная векторизация | 3.67 с | 4.37 с | 6.19 с |
№7 — коллапсирование гнезда циклов | 3.60 с | 4.25 с | 5.46 с |
№8 — уменьшаем нагрузку на память | 3.43 с | 3.80 с | 5.04 с |
А что там у -fwhole: компиляция в режиме «вся программа»
Иногда задумываешься: а нужно ли столько трудов по оптимизации? Может быть, есть волшебный ключик, при использовании которого компилятор сам выполнит все оптимизации? На роль такой опции лучше всего претендует межпроцедурный анализ в контексте целиковой программы на этапе линковки и, при возможности, профилировка. В случае xlc достаточно продублировать -O5 во флагах линковщика, для icc и lcc надо подать -ipo и -fwhole, соответственно, при компиляции и линковке. А вот результаты:
Реализация | i7-9700K | POWER8 | Эльбрус-8СВ |
№0 — базовая | 5.0 с | 19.6 с | 94.5 с |
№1 — возможен inline | 5.0 c | 19.8 c | 95.1 c |
№2 — ssize_t счётчики | 5.0 с | 19.8 с | 18.7 с |
№3 — алгоритмические оптимизации | 4.90 с | 7.05 с | 11.01 с |
№4 — используем restrict | 3.70 с | 4.91 с | 7.41 с |
№5 — вынос инварианта | 3.68 с | 4.64 с | 7.03 с |
№6 — ручная векторизация | 3.68 с | 4.64 с | 6.33 с |
№7 — коллапсирование гнезда циклов | 3.65 с | 3.36 с | 5.44 с |
№8 — уменьшаем нагрузку на память | 3.43 с | 3.33 с | 5.04 с |
В целом, всё видно из таблицы, но я подчеркну: межпроцедурный анализ даёт хорошие результаты, но не заменяет оптимизации, которые может применить разработчик.
Заключение
На примере относительно простой учебной программы мне удалось обнаружить основные сложности и продемонстрировать техники, с которыми приходится сталкиваться при оптимизации под Эльбрус. Конечно, есть ещё ряд особенностей, которые не вошли в эту статью, но их вклад я считаю менее заметным и важным на первых порах оптимизации.
С помощью нехитрых приёмов мне удалось получить на Эльбрусе впечатляющее ускорение в 44 раза и приблизить время вычисления на Эльбрусе к результату на современном Intel и IBM POWER8, причём с помощью тех же оптимизаций результат на Intel был улучшен в 4.1 раза, а на IBM — в 11.5 раза (13 раз с учётом межпроцедурного анализа).
Нельзя не отметить, что оптимизация под Эльбрус не отличается принципиально от оптимизации под другие процессоры общего назначения и следует, в основном, тем же принципам. При этом совершенно не надо бояться VLIW-ассемблера: при наличии минимальных навыков работы с любым компилятором не составляет труда разобраться в сгенерированном коде и понять, из-за чего результат не соответствует ожиданиям.
