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

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

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

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

Итак, рассматривается дифференциальное уравнение Пуассона в прямоугольной области [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 шагов итерационного метода на одном ядре процессора на сетке размером 1000x1000. За эффектом от оптимизаций будем следить на трёх платформах:

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,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.

И проверим, как же теперь выглядит исследуемый цикл на языке ассемблера. Отлично, все вхождения цикла стали компактными, раскрученными на 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 говорится о том, что компилятор умеет сам делать вынос инварианта. Что ж, замеры показывают, что надёжнее это сделать руками:

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

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

Остался один очевидный шаг, который можно было бы сделать ещё на этапе алгоритмических оптимизаций, но я сознательно пропустил его в угоду унификации. В функции 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 в обоих случаях и получить итоговый результат. Ниже сводная таблица со всеми основными результатами.

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

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

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

Заключение

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

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

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