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

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

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

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

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

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

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

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

Для проверки условия остановки вычисляется норма разности между новым и старым приближением. За более подробным описанием отправляю читателя к профильной литературе, например, книге А.А. Самарского [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], а я просто воспользуюсь опытом, полученным при реализации криптографических примитивов.

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

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

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

Примечательно, что на Эльбрусе, который имеет фиксированные и относительно большие штрафы за вызов функций, ускорение составляет 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, чтобы не плодить платформозависимых изменений.

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

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

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

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

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

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

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

Разобраться совсем не так сложно, как может показаться на первый взгляд. Работает универсальный приём: просто берём строчку компиляции нужного файла и заменяем флаг -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,sm0x8, %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   alcalcf=1, alct=1! solve_cpu.c : 303   abnabnf=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.

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

.L38838:! solve_cpu.c : 312 ! <0237> {   loop_mode   fmuld,5,sm%dr7, %db[11], %db[15]! solve_cpu.c : 314   movad,0area=0, ind=0, am=0, be=0, %db[0]! solve_cpu.c : 315   movad,1area=0, ind=8, am=1, be=0, %db[1]! solve_cpu.c : 315   movad,3area=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   alcalcf=1, alct=1! solve_cpu.c : 312   abnabnf=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,3area=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 говорится о том, что компилятор умеет сам делать вынос инварианта. Что ж, замеры показывают, что надёжнее это сделать руками:

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

.L19038:! solve_cpu.c : 270 ! <0259> {   loop_mode   movad,2area=0, ind=8, am=1, be=0, %db[1]! solve_cpu.c : 272   movad,3area=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,0area=0, ind=8, am=1, be=0, %db[6]! solve_cpu.c : 272   movad,1area=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   alcalcf=1, alct=1! solve_cpu.c : 270   abnabnf=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   alcalcf=1, alct=1! solve_cpu.c : 272   abnabnf=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,1area=0, ind=0, am=1, be=0, %db[26]! solve_cpu.c : 274   movad,3area=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,1area=1, ind=0, am=1, be=0, %xb[0]! solve_cpu.c : 242 } ! <0248> {   loop_mode   alcalcf=1, alct=1! solve_cpu.c : 234   abnabnf=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,0area=0, ind=0, am=1, be=0, %xb[21]! solve_cpu.c : 253   movaqp,1area=2, ind=0, am=1, be=0, %xb[72]! solve_cpu.c : 237   movaqp,2area=0, ind=0, am=1, be=0, %xb[59]! solve_cpu.c : 252   movaqp,3area=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 пользы от ручного выбора параметра раскрутки не обнаружилось.

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

Остался один очевидный шаг, который можно было бы сделать ещё на этапе алгоритмических оптимизаций, но я сознательно пропустил его в угоду унификации. В функции 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,0area=0, ind=0, am=0, be=0, %xb[43]! solve_cpu.c : 186   movaqp,1area=0, ind=16, am=1, be=0, %xb[50]! solve_cpu.c : 186   movaqp,3area=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,0area=1, ind=0, am=0, be=0, %xb[88]! solve_cpu.c : 175   movaqp,1area=1, ind=16, am=1, be=0, %xb[31]! solve_cpu.c : 175   movaqp,2area=0, ind=0, am=0, be=0, %xb[75]! solve_cpu.c : 185   movaqp,3area=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,0area=3, ind=0, am=1, be=0, %xb[17]! solve_cpu.c : 185   movaqp,1area=4, ind=0, am=1, be=0, %xb[18]! solve_cpu.c : 171   movaqp,2area=4, ind=0, am=1, be=0, %xb[24]! solve_cpu.c : 170   movaqp,3area=1, ind=16, am=0, be=0, %xb[30]! solve_cpu.c : 171 } ! <0057> {   loop_mode   alcalcf=1, alct=1! solve_cpu.c : 168   abnabnf=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,0area=2, ind=16, am=1, be=0, %xb[55]! solve_cpu.c : 170   movaqp,1area=2, ind=0, am=0, be=0, %xb[68]! solve_cpu.c : 170   movaqp,2area=1, ind=0, am=1, be=0, %xb[67]! solve_cpu.c : 171   movaqp,3area=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,1area=0, ind=16, am=0, be=0, %xb[1]! solve_cpu.c : 324   movaqp,2area=0, ind=0, am=0, be=0, %xb[0]! solve_cpu.c : 323   movaqp,3area=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,1area=0, ind=0, am=1, be=0, %xb[7]! solve_cpu.c : 324 } ! <0067> {   loop_mode   alcalcf=1, alct=1! solve_cpu.c : 321   abnabnf=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,2area=0, ind=0, am=0, be=0, %xb[1]! solve_cpu.c : 393   movaqp,3area=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,0area=0, ind=16, am=1, be=0, %xb[6]! solve_cpu.c : 394   movaqp,1area=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   alcalcf=1, alct=1! solve_cpu.c : 391   abnabnf=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,0area=2, ind=0, am=1, be=0, %xb[15]! solve_cpu.c : 211   movaqp,1area=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,0area=0, ind=16, am=1, be=0, %xb[27]! solve_cpu.c : 211   movaqp,1area=3, ind=0, am=1, be=0, %xb[48]! solve_cpu.c : 197   movaqp,2area=3, ind=0, am=1, be=0, %xb[53]! solve_cpu.c : 196   movaqp,3area=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,0area=1, ind=16, am=1, be=0, %xb[73]! solve_cpu.c : 197   movaqp,1area=1, ind=0, am=0, be=0, %xb[79]! solve_cpu.c : 197   movaqp,2area=1, ind=0, am=1, be=0, %xb[85]! solve_cpu.c : 196   movaqp,3area=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,2area=0, ind=0, am=0, be=0, %xb[60]! solve_cpu.c : 201   movaqp,3area=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   alcalcf=1, alct=1! solve_cpu.c : 194   abnabnf=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 в обоих случаях и получить итоговый результат. Ниже сводная таблица со всеми основными результатами.

А что там у -fwhole: компиляция в режиме «вся программа»

Иногда задумываешься: а нужно ли столько трудов по оптимизации? Может быть, есть волшебный ключик, при использовании которого компилятор сам выполнит все оптимизации? На роль такой опции лучше всего претендует межпроцедурный анализ в контексте целиковой программы на этапе линковки и, при возможности, профилировка. В случае xlc достаточно продублировать -O5 во флагах линковщика, для icc и lcc надо подать -ipo и -fwhole, соответственно, при компиляции и линковке. А вот результаты:

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

Заключение

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

С помощью нехитрых приёмов мне удалось получить на Эльбрусе впечатляющее ускорение в 44 раза и приблизить время вычисления на Эльбрусе к результату на современном Intel и IBM POWER8, причём с помощью тех же оптимизаций результат на Intel был улучшен в 4.1 раза, а на IBM — в 11.5 раза (13 раз с учётом межпроцедурного анализа).

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