Здравствуйте, друзья, меня зовут Ерохин Кирилл, я программист-любитель, а по совместительству популяризатор российского программного и аппаратного обеспечения, и в этом сентябре я провожу второе (теперь ежегодное) соревнование по алгоритмическому программированию на C/C++ для платформы Эльбрус (e2k), для студентов и выпускников со всей России «Кубок СЭРПАС 2026». Сегодня мне нужно дать участникам соревнования пример оптимизации кода под процессор Эльбрус-8СВ, а Хабр мне в этом поможет, ему не впервой.

Оглавление

1. Введение

2. Эльбрус‑8СВ и что нужно знать для оптимизации под него

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

4. Доступный инструментарий

5. Шаг 0. Базовое решение

6. Шаг 1. Быстрый ввод‑вывод и правильная структура данных

7. Шаг 2. Программная конвейеризация (подсказываем компилятору)

8. Шаг 3. Блочная обработка (борьба с пропускной способностью памяти)

9. Шаг 4. Разворот с объединением

10. Шаг 5. Многопоточность

11. Заключение

1. Введение

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

Сама идея автоматизации соревнований и олимпиад по программированию не нова, и в России существует достаточно много крупных онлайн-площадок для организации и проведения подобных мероприятий, однако ни одна из них не имела опыта работы с такой специфической аппаратной платформой и не адаптировала программный стек своей платформы под архитектуру Эльбрус (e2k).

После общения с тремя крупнейшими российскими площадками мне повезло: создатели Codenrock заинтересовались моей самодеятельностью и, изучив историю и результаты первого соревнования, согласились помочь, отказавшись при этом брать плату как за предстоящую адаптацию программного стека своей платформы под Эльбрусы, так и в целом за проведение и техническое сопровождение предстоящего «Кубка СЭРПАС 2026» на ней, за что им отдельная благодарность.

Веб-платформа автоматизированной проверки и учёта решений и участников соревнования снимала проблему роста их количества на корню, смещая бутылочное горлышко проведения соревнования непосредственно к аппаратному обеспечению, а точнее к количеству доступных серверов на процессорах Эльбрус (e2k), с которыми у меня после значительно оживившегося в этом году вторичного рынка материнских плат на Эльбрусах проблем не было. Однако она же (платформа) вносила важный нюанс в сам процесс отладки кода и анализа своих решений (программ) участниками, поскольку в отборочном этапе исчезал прямой доступ к серверам через SSH (в финале всё будет по-прежнему через SSH), а значит, традиционными утилитами (GDB, Valgrind, perf, gprof и т.д.) участникам воспользоваться уже не получится и придётся обходиться ручным инструментированием (да, оказывается именно так называется старая добрая printf-отладка).

К счастью, ещё в прошлом году МЦСТ выпустили долгожданный QEMU для платформы Эльбрус (e2k), что значительно расширило возможности участников соревнования (да и не только их) по отладке своих программ, позволив для проверки корректности работы программы использовать QEMU, а для оценки скорости её работы — printf-отладку через веб-платформу Codenrock.

В качестве же примера, который мы с читателем будем разбирать, я решил выбрать первую задачу отборочного этапа прошлогоднего соревнования на решение СЛАУ вида A·x=b (вторую задачу ребята из НИЦ ЦТ и так без нас разобрали). Единственное отличие от оригинала: величину N, определяющую размер матрицы A и длину b, мы установим равной 2000. При меньших значениях N некоторые оптимизации (например многопоточность) могут банально не успеть показать себя.

Отправной точкой изучения программирования и оптимизации кода под Эльбрус (e2k) является фундаментальный труд сотрудников МЦСТ «Руководство по эффективному программированию на платформе “Эльбрус”», доступный всем абсолютно безвозмездно.

2. Эльбрус‑8СВ и что нужно знать для оптимизации под него

Чтобы нам с тобой, дорогой читатель, понять, почему один код работает быстрее, а другой — нет, нужно разобраться, что из себя представляет наше целевое аппаратное обеспечение, то есть как устроена архитектура Эльбрус и конкретно наш пациент процессор. Эльбрус‑8СВ — это восьмиядерный процессор общего назначения семейства «Эльбрус», реализующий VLIW-архитектуру и систему команд e2k пятой версии, а это значит, что нам с тобой, дорогой читатель, предстоит работать с:

Широким командным словом (архитектура VLIW). В привычных читателю процессорах x86-64 порядок выполнения операций (в том числе параллельных) определяет сам процессор. В Эльбрусе всё иначе: процессор за один такт исполняет широкую команду, в которую заранее, ещё на этапе трансляции, упаковано несколько элементарных операций (арифметических, обращений к памяти, переходов). Решение о том, какие операции выполнять параллельно, принимает компилятор (его планировщик команд на этапе компиляции), а не процессор. Это ключевая мысль: производительность на Эльбрусе прямо зависит от того, насколько хорошо компилятор сумел запаковать команды в широкую команду (командное слово). Если код написан так, что зависимости между операциями мешают их совмещать, широкие команды остаются полупустыми, и ядро процессора простаивает. Отсюда практический вывод: значительная часть оптимизации на Эльбрусе не «хитрости» и «хаки», а помощь компилятору: избавление его от сомнений в независимости данных и предоставление ему длинных, регулярных потоков однотипных операций. Подробно об этом можно прочитать в главах 4.3 и 4.4 Руководства.

На этом месте дорогой читатель может подумать, что компилятор Эльбрусов (в нашем случае это LCC) — это такой ленивый болван, за которого программист должен делать всю работу, но это он просто не разобрался:

• сама суть VLIW‑архитектуры заключена в переносе значительной части задач по управлению ресурсами и планированию инструкций с процессора на компилятор. Компилятор должен заранее спланировать выполнение операций, опираясь преимущественно на данные, доступные на этапе компиляции, и не имея информации о событиях, возникающих во время исполнения программы (промахи кеша, фактические задержки памяти и результаты ветвлений).

• компилятор LCC непрерывно развивается с 1997 года, и отдел разработки языкового компилятора является одним из ключевых в МЦСТ. Альма-матер большинства его сотрудников, регулярно выступающих на профильных мероприятиях, — это ведущие вузы России в области системного программирования (МФТИ, МИФИ и т.д.), а количеству публикаций на тему оптимизирующих компиляторов в научных рецензируемых журналах (например «Труды Института системного программирования РАН») можно только позавидовать. Но важно помнить, что, несмотря на колоссальное количество сил и знаний, вкладываемых в LCC, он не волшебник и не может заранее и безошибочно предсказывать поведение кода сколь угодно большой сложности и объёма.

Большим регистровым файлом и аппаратной поддержкой циклов. Регистров у Эльбруса непривычно много: регистровый файл ядра насчитывает 256 регистров, и окно одной процедуры может занимать десятки из них — против шестнадцати архитектурных регистров общего назначения в x86‑64. А для счётных циклов есть и особая аппаратная поддержка: часть регистрового окна может «вращаться», автоматически смещая нумерацию регистров от итерации к итерации. Это позволяет компилятору строить глубокие программные конвейеры без лишних пересылок между регистрами.

Что это значит для нас на практике? Это значит, что вычисления, которым нужно удерживать в регистрах много промежуточных значений (например, обработка нескольких строк матрицы сразу), на Эльбрусе окупаются легче, чем на x86-64, где регистров хватает не всегда. Подробнее про это можно почитать в главах 9.3 и 10.11 Руководства.

Аппаратной подкачкой массивов (APB). В широкой команде чтение памяти может выполняться в четырёх каналах, а запись — только в двух, кроме того, каждое ядро Эльбрус-8СВ имеет специальный механизм аппаратной предварительной подкачки массивов (APB). Если программа читает массив подряд, регулярным шагом, APB заранее подгружает следующие порции данных в собственный буфер, пока ядро занято счётом, скрывая значительную часть задержки памяти. Условие срабатывания APB — именно регулярный, предсказуемый обход памяти. Это ещё одна причина, по которой плоские массивы с последовательным доступом на Эльбрусе выигрывают у «разрозненных» структур данных. Механизм APB подробно разобран в главе 6.6 Руководства.

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

Дана невырожденная целочисленная квадратная матрица A размера 2000×2000 со строгим диагональным преобладанием и вектор b длиной 2000 свободных членов. Требуется найти вектор x (элементы которого кратны 0,25), удовлетворяющий уравнению A·x = b, и вывести его компоненты с шестью знаками после запятой. При этом библиотека EML, как и другие высокооптимизированные математические инструменты (BLAS, LAPACK и т.д.), недоступны.

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

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

• все компоненты точного ответа кратны 0.25 и представимы в двоичном виде без ошибки. Однако промежуточные деления и вычитания округляются, поэтому перестановка операций, векторизация и распараллеливание могут слегка изменить младшие разряды результата. Большое расстояние между соседними допустимыми ответами позволяет сохранить правильный вывод с шестью знаками, но не гарантирует побитового совпадения. Иными словами, нам разрешено переупорядочивать вычисления, но после каждого такого изменения результат необходимо проверять.

Я думаю, ни для кого не открытие, что для прошлогоднего соревнования я специально подбирал задачи, хорошо оптимизируемые под Эльбрусы.

4. Доступный инструментарий

Очевидно, что ПК с процессором Эльбрус у большинства участников соревнования (да и у тебя, уважаемый читатель) нет. К счастью, он нам сейчас и не нужен. Всё, что мы будем делать, прекрасно воспроизводится и запускается с доступными всем участникам соревнования средствами, ранее упомянутыми мной в введении:

Кросс-компилятор + эмулятор QEMU (e2k). Кросс-компилятор LCC позволяет собирать приложения для процессоров архитектуры Эльбрус (e2k), работая на ПК (Linux) с процессорами архитектуры x86-64. Флаги кросс-компилятора идентичны нативному e2k-компилятору LCC: https://dev.mcst.ru/кросс-компилятор-lcc-1-29-16/

QEMU (e2k) — программный эмулятор для запуска программ, предназначенных для аппаратной платформы Эльбрус (e2k). Если проще, то он запускает программы под Linux, собранные для архитектуры e2k, на обычном x86-64 компьютере: https://dev.mcst.ru/эльбрус-qemu-user-версия-1-2/

Эмулятор — это инструмент функциональной проверки: он воспроизводит поведение программы, но не её скорость. Запуск выглядит примерно так:

qemu-e2k -cpu elbrus-v5 ./solver < in2000.txt > out.txt

Флаг -cpu elbrus-v5 задаётся явно: по умолчанию эмулятор моделирует более свежую 6-ю версию системы команд (v6), а Эльбрус‑8СВ реализует 5-ю версию системы команд e2k. При этом проверяемую программу (бинарник) проще собирать статически (флаг -static). Память гостевой программе по умолчанию отводится в объёме 1 ГБ (ключ -R), а нашей задаче при N = 2000 (матрица ≈ 32 МБ) этого вполне хватит.

Стоит помнить, что у текущей версии QEMU (e2k) есть три ограничения:

• потоки исполняются поочерёдно, а не одновременно. Правильность многопоточной программы проверить можно (ответ обязан совпасть с эталонным), ускорение — нельзя; более того, поочерёдное исполнение способно маскировать гонки данных. Для быстрых функциональных прогонов многопоточного варианта стоит задавать OMP_NUM_THREADS=1;

• время исполнения под эмулятором ничего не говорит о времени исполнения на реальном Эльбрусе. Все замеры скорости исполнения кода в эмуляции происходят на x86-64 процессоре, к тому же без настоящей многопоточности. Поэтому любые замеры времени в QEMU бессмысленны;

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

Помимо простого запуска, у эмулятора есть ещё несколько полезных возможностей. Ключ -strace печатает журнал системных вызовов — по нему видно, например, сколько раз и какими порциями программа читает ввод. Ключ -d in_asm,exec ведёт журнал исполнения: на уменьшенном размере задачи по нему можно находить горячие участки кода (журналы велики, для больших N способ непригоден). Счёт вызовов функций (инструментирование gprof) работает; выборка по времени — нет. Отладчик подключается через -g <порт>, но состояние гостевой программы доступно только на чтение: точки останова, пошаговое исполнение, осмотр регистров и памяти — да, изменение — нет. Наконец, ключ -cpu позволяет функционально проверить программу на версиях системы команд от elbrus-v2 до elbrus-v6, не имея соответствующих машин.

Онлайн-компилятор. Это веб-интерфейс, который мы вместе с командой платформы Codenrock подготовили для участников соревнования, в рамках которого вставляемый в веб-форму исходный код программы отправляется на настоящий вычислительный сервер с Эльбрус‑8СВ (он не один на самом деле), на нём собирается компилятором LCC (версия 1.29.16, флаги -O3 -ffast -fopenmp) и запускается (не один раз для фильтрации от шума). Обратно возвращаются стандартный вывод или поток ошибок (если таковые были). Никакого иного доступа к серверу нет, а значит, замеры встраиваются в саму программу (printf-отладка). Участки анализируемого кода оборачиваются вызовами (например clock_gettime), а результаты печатаются в поток вывода.

По сути, QEMU (e2k) отвечает на вопрос в отношении работы программы (решения) «правильно ли?», а онлайн-компилятор — на вопрос «как быстро?»

Для удобства я подготовил вот такую таблицу-шпаргалку:

Внимательный читатель наверняка заметил, что, судя по количеству «Да», printf-отладка во всех отношениях выглядит предпочтительнее варианта QEMU + gdb. Однако стоит помнить, что в нашем случае printf-отладка осуществляется через онлайн-компилятор в рамках платформы Codenrock, пусть и без ограничений по количеству попыток, а запуск и отладка программы с использованием QEMU через привычную IDE будет куда удобнее и быстрее для большинства участников.

5. Шаг 0. Базовое решение

Решение любой задачи алгоритмического программирования начинается с анализа задачи и выбора оптимального (внезапно!) алгоритма, а в нашем случае — метода решения СЛАУ. Конечно, сложно найти более хрестоматийную задачу в линейной алгебре, чем решение СЛАУ, и методов, которыми её решают, уже есть не менее дюжины, а то и двух десятков. Принято выделять две основные группы таких методов:

• прямые (точные) методы. Дают решение за конечное число шагов. Примеры: метод Гаусса, метод Крамера, матричный метод, LU-разложение, метод Гаусса-Жордана, метод квадратных корней и так далее.

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

Какой же метод подойдёт нам лучше всего?

Достаточно спросить нейронку проанализировать, что лучше подходит для архитектуры Эльбрус (e2k) независимо от конкретной модели процессора:

• заранее выявляемый параллелизм. Транслятору нужно, чтобы независимые операции были видны в коде — тогда он наполнит ими широкую команду.

• последовательный доступ к памяти (единичный шаг). Тогда включается аппаратный массив предварительной подкачки (APB) и заранее подтягивает данные, скрывая задержку памяти.

• высокая плотность вычислений. Много арифметики на единицу данных, чтобы насытить многочисленные вещественные каналы.

• комбинированные операции. Умножение со сложением/вычитанием за один такт.

• регулярный, предсказуемый поток управления. Ядро исполняет заранее составленное транслятором расписание и не перестраивает его динамически, как современные x86; поэтому «скачущий» код транслятору значительно труднее эффективно спланировать.

• длинные независимые цепочки. Материал для программной конвейеризации, при которой транслятор перекрывает соседние итерации.

• большой регистровый файл. Много значений «в полёте» одновременно.

Теперь сопоставим с основными и самыми популярными методами решения СЛАУ (да, я знаю их гораздо больше, но принципиально для нашей задачи они не отличаются):

Поздравляю, победил метод Гаусса! Это первый или второй курс любого технического вуза или колледжа (или техникума? Как их сейчас называют-то?).

Если вдруг читатель забыл (или даже с ним незнаком), я любезно напомню. Метод состоит из двух этапов. Сначала выполняется прямой ход, в котором с помощью элементарных преобразований строк матрица приводится к треугольному виду, то есть под главной диагональю получаются нули. Затем выполняется обратный ход, в котором, двигаясь снизу вверх по уже треугольной системе, мы последовательно выражаем неизвестные. Основная вычислительная работа сосредоточена в прямом ходе, и его трудоёмкость составляет порядка (2/3)·N³ арифметических операций. Именно этот кубический по N объём и делает задачу интересной для оптимизации: при нашем N = 2000 речь идёт примерно о 2,67 миллиарда умножений-вычитаний или 5,33 миллиарда отдельных операций.

Краеугольный камень метода Гаусса в рамках нашей задачи (мы тут код под Эльбрус оптимизируем, если что) — это тот самый «прямой ход», который реализуется одним циклом. Из каждой нижележащей строки вычитается ведущая строка, помноженная на коэффициент ( A[i][j] -= factor * A[k][j]). Отсюда получаем его свойства, так хорошо подходящие архитектуре Эльбрус (e2k):

• независимые итерации внутреннего цикла, элементы строки обрабатываются независимо друг от друга;

• проход по строке подряд. Память читается и пишется единичным шагом;

• объём работы ⅔×N³ операций над N² данными. Задача ограничена счётом, а не памятью;

• одно умножение и одно сложение на элемент. Эталонная комбинированная операция;

• отсутствие выбора ведущего элемента (следствие строгого диагонального преобладания). Нет перестановок строк, а единственное заметное зависящее от данных ветвление проверяет factor == 0.0;

• простой цикл с известными границами — его легко развернуть и конвейеризовать;

• свобода контролируемого переупорядочивания. Ответы кратны 0.25, но после изменения порядка операций результат нужно проверить.

Для первого нулевого (мы же программисты) шага напишем базовое решение без оптимизаций прямо как в учебнике:

• потоковый ввод‑вывод стандартной библиотеки;

• хранение матрицы в виде «вектора векторов»;

• прямой и обратный ход.

#include <iostream>
#include <iomanip>
#include <vector>

int main() {
    // Чтение размера матрицы
    int n;
    std::cin >> n;

    // Чтение матрицы A размера N×N
    std::vector<std::vector<double>> A(n, std::vector<double>(n));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            std::cin >> A[i][j];

    // Чтение вектора свободных членов b
    std::vector<double> b(n);
    for (int i = 0; i < n; ++i)
        std::cin >> b[i];

    // Прямой ход метода Гаусса: приведение к верхнетреугольному виду
    // Строгое диагональное преобладание гарантирует ненулевые диагональные элементы
    for (int k = 0; k < n; ++k) {
        // Нормировка k-й строки: деление на диагональный элемент
        double pivot = A[k][k];
        for (int j = k; j < n; ++j)
            A[k][j] /= pivot;
        b[k] /= pivot;

        // Исключение переменной x_k из нижележащих уравнений
        for (int i = k + 1; i < n; ++i) {
            double factor = A[i][k];
            if (factor == 0.0)
                continue;
            for (int j = k; j < n; ++j)
                A[i][j] -= factor * A[k][j];
            b[i] -= factor * b[k];
        }
    }

    // Обратный ход: вычисление компонент вектора решения
    std::vector<double> x(n);
    for (int i = n - 1; i >= 0; --i) {
        double sum = b[i];
        for (int j = i + 1; j < n; ++j)
            sum -= A[i][j] * x[j];
        x[i] = sum;
    }

    // Вывод результата с шестью знаками после запятой
    std::cout << std::fixed << std::setprecision(6);
    for (int i = 0; i < n; ++i) {
        std::cout << x[i] << "\n";
    }
    return 0;
}

Время исполнения: 24 237 мс (просьба преждевременно не падать в обморок)

Ускорение: ×1,0 (база отсчёта).

6. Шаг 1. Быстрый ввод‑вывод и правильная структура данных

Наблюдательный читатель уже обратил своё внимание на неожиданно большое время выполнения базового решения (программы) в предыдущем шаге и сейчас мы с вами разберёмся в чём тут причина.

Замерим время выполнения основных операций, функций и блоков кода, используя printf-отладку. Добавим заголовочные файлы chrono для замеров времени и stdio для printf(), а также вспомогательные элементы для удобства замеров:

#include <chrono>
#include <cstdio>

using clk = std::chrono::steady_clock;

static double ms(clk::time_point a, clk::time_point b) {
    return std::chrono::duration<double, std::milli>(b - a).count();
}

Основными операциями (блоками) программы у нас являются:

1. Чтение (ввод) N

2. Чтение (ввод) матрицы A

3. Чтение (ввод) вектора b

4. Прямой ход

5. Обратный ход

6. Вывод результата (вектора)

Используя операции вида clk::time_point t = clk::now(); и нашу вспомогательную функцию ms() для вычисления разницы между двумя метками времени, узнаём скорость выполнения основных операций:

Первое, что должно удивить читателя, не знакомого с разработкой под платформу Эльбрус (e2k), — это то, что самым узким (горячим) местом оказалось чтение (ввод) матрицы A, а вовсе не операции вычисления «прямого хода». Причина этого заключается в том, что потоковый ввод стандартной библиотеки std::cin устроен достаточно тяжело: на каждое число он проверяет форматные флаги, состояние потока, локаль и прочее, и делает это через цепочку вызовов, которую у компилятора LCC не получается хорошо оптимизировать. Для процессоров Эльбрус такой поток нерегулярной, ветвящейся работы — прямая противоположность тому, что они ожидают получить на обработку.

Данный пример при всей своей простоте ярко демонстрирует, почему проекты (код), не требующие явного портирования под «Эльбрус» (e2k) и без проблем компилируемые LCC, всё же требуют анализа и оптимизации. Это прямое следствие того, что последние 40 лет для персональных компьютеров (25 лет для серверов) архитектура x86 (включая x86-64) является доминирующей, и все приложения, библиотеки и алгоритмы создаются с фокусом на x86, делая весь окружающий нас IT-мир x86-центричным.

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

При этом, несмотря на то что внутри каждой строки элементы расположены последовательно (что, казалось бы, хорошо), при переходах между строками, особенно при обходе матрицы по столбцам, процессору приходится обращаться к менее регулярной последовательности адресов. Механизм асинхронной подкачки данных в буфер предподкачки массива (APB) рассчитан на регулярные последовательности адресов, поэтому использование двумерных векторов (std::vector<std::vector<T>>) может снижать его эффективность.

Первым делом разберёмся с медленным вводом (чтением). Чтобы уменьшить накладные расходы, входные данные читаются не по одному символу или числу, а сразу большими блоками, в нашем случае будем читать по 64 КБ (выбираем эмпирически), а прочитанный блок будем сохранять в обычном массиве:

static char inputBuffer[1 << 16]; // 64 КБ (2^16 байт)
static int bufLen = 0; // количество данных в буфере
static int bufPos = 0; // позиция следующего символа

// Функция возвращает следующий символ входных данных
static int readChar() {
    // Если все символы текущего блока уже обработаны, читаем следующий блок
    if (bufPos == bufLen) {
        bufLen = static_cast<int>(
            std::fread(inputBuffer, 1, sizeof(inputBuffer), stdin)
        );
        bufPos = 0;

        // Нулевой результат означает конец файла или ошибку чтения
        if (bufLen <= 0) {
            return -1;
        }
    }

    // Возвращаем очередной символ и переходим к следующему.
    return inputBuffer[bufPos++];
}

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

1. читает очередной символ из массива inputBuffer;

2. увеличивает bufPos;

3. возвращает прочитанный символ.

Когда bufPos достигает bufLen, текущий блок заканчивается. Тогда функция std::fread запрашивает следующую порцию данных и снова заполняет буфер. Таким образом, сравнительно дорогая операция чтения выполняется редко, примерно один раз на каждые 64 КБ входных данных, а не для каждого отдельного символа или числа.

Сама функция readChar() ещё не преобразует текст в числа. Она только последовательно выдаёт символы. Такой способ ввода удобен для Эльбруса не из-за какой-либо специальной команды процессора, а благодаря более простой и регулярной обработке:

• данные читаются крупными блоками;

• символы в буфере просматриваются последовательно;

• уменьшаются накладные расходы форматированного ввода;

• основной цикл работает с обычным непрерывным массивом;

• компилятору проще оптимизировать простой цикл разбора символов.

Далее нам нужно, используя readChar(), читать конкретные числа. Матрица A по условию задачи целочисленная, и этим грех не воспользоваться: для неё достаточно простейшего целочисленного разбора без всякой работы с дробной частью:

// Функция чтения целого числа со знаком
static long readInteger() {
    int c = readChar();
    while (c != '-' && (c < '0' || c > '9')) c = readChar(); // пропуск пробелов

    bool negative = (c == '-');
    if (negative) c = readChar();

    long value = 0;
    while (c >= '0' && c <= '9') {
        value = value * 10 + (c - '0');
        c = readChar();
    }
    return negative ? -value : value;
}

Схема прозрачна: пропустить разделители, запомнить знак, накопить число по цифрам (классическое value = value \* 10 + цифра). Ни настроек локали, ни универсального разбора форматов, только то, что действительно нужно для решения нашей конкретной задачи.

Теперь переходим к вектору b. Его элементы могут содержать дробную часть, например «-284.500000 15.250000 7.750000». Поэтому функции readInteger() уже недостаточно. Новая функция должна распознать четыре части записи числа:

1. Пропустить пробелы и другие разделители;

2. Определить знак числа;

3. Прочитать целую часть;

4. При наличии точки прочитать дробную часть.

// Функция читает число с возможной дробной частью.
static double readNumber() {
    int c = readChar();

    // Пропускаем пробелы, переводы строк и другие разделители.
    // Останавливаемся на знаке минус или первой цифре числа.
    while (c != '-' && (c < '0' || c > '9')) {
        c = readChar();
    }

    // Запоминаем знак числа.
    bool negative = (c == '-');
    if (negative) {
        c = readChar();
    }

    // Читаем целую часть.
    double value = 0.0;
    while (c >= '0' && c <= '9') {
        value = value * 10.0 + (c - '0');
        c = readChar();
    }

    // Если после целой части встретилась точка,
    // читаем цифры дробной части.
    if (c == '.') {
        c = readChar();
        double fraction = 0.0;
        double scale = 1.0;
        while (c >= '0' && c <= '9') {
            fraction = fraction * 10.0 + (c - '0');
            scale *= 10.0;
            c = readChar();
        }
        value += fraction / scale;
    }

    // Применяем сохранённый знак.
    return negative ? -value : value;
}

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

1. читает очередной символ;

2. проверяет, является ли он цифрой;

3. преобразует цифру вычитанием ‘0’;

4. добавляет её к накапливаемому значению.

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

После подготовки функций быстрого ввода переходим к оптимизации чтения в «сплющенную» матрицу (элементы которой логически записаны в одномерном массиве). Сначала читаем размер квадратной матрицы:

int n = (int)readInteger();

Далее главное структурное изменение. Вместо двумерного вектора используем один непрерывный массив длиной N×N:

// чтение матрицы A в «сплющенную» матрицу элемент A[i][j] -> A[i*n + j]
std::vector<double> A((size_t)n * n);

for (int i = 0; i < n; ++i) {
    double* row = &A[(size_t)i * n]; // строка i лежит подряд
    for (int j = 0; j < n; ++j)
        row[j] = (double)readInteger();
}

При этом мы по-прежнему используем std::vector, со всеми его удобствами, но один, а не тысяча вложенных. Указатель row наводится на начало строки i, и дальше строка заполняется по единичному шагу. Дальше читаем вектор свободных членов:

// чтение вектора свободных членов b
std::vector<double> b(n);

for (int i = 0; i < n; ++i)
    b[i] = readNumber();

Теперь займёмся прямым ходом. Сам алгоритм по сравнению с базовой версией изменять не нужно. Меняем только способ хранения и обращения к элементам матрицы, храня всё в одном std::vector<double>, а начало нужной строки будем определять с помощью указателя:

// Прямой ход
// Матрица A хранится в сплющенном массиве.
for (int k = 0; k < n; ++k) {
    // Указатель на начало ведущей строки k.
    double* rowK =
        A.data() + static_cast<std::size_t>(k) * n;

    // Ведущий элемент A[k][k].
    double pivot = rowK[k];

    // Делим ведущую строку на ведущий элемент.
    for (int j = k; j < n; ++j) {
        rowK[j] /= pivot;
    }
    b[k] /= pivot;

    // Обнуляем элементы под ведущим элементом.
    for (int i = k + 1; i < n; ++i) {
        // Указатель на начало обрабатываемой строки i.
        double* rowI =
            A.data() + static_cast<std::size_t>(i) * n;

        // Элемент A[i][k], который нужно обнулить.
        double factor = rowI[k];

        // Если элемент уже равен нулю, строку можно не изменять.
        if (factor == 0.0) {
            continue;
        }

        // Вычитаем из строки i ведущую строку,
        // умноженную на factor.
        for (int j = k; j < n; ++j) {
            rowI[j] -= factor * rowK[j];
        }

        // Такую же операцию выполняем с правой частью системы.
        b[i] -= factor * b[k];
    }
}

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

Далее ведущий элемент rowK[k] сохраняется в переменной pivot. Затем на него делятся элементы ведущей строки, начиная с позиции k, и соответствующий элемент b[k]. В результате ведущий элемент становится равен единице, а уравнение остаётся равносильным исходному. Элементы левее k не обрабатываются, поскольку они уже были преобразованы на предыдущих шагах.

После нормализации ведущей строки обрабатываются все строки ниже неё. Для каждой строки один раз вычисляется указатель rowI, а элемент rowI[k] сохраняется в переменной factor. Затем из строки rowI вычитается ведущая строка rowK, умноженная на factor, и так же изменяется b[i]. В результате элемент под главной диагональю становится равен нулю. Если factor уже равен нулю, обработка строки пропускается, поскольку вычитание нулевого множителя ничего не изменит.

Для обратного хода делаем всё аналогично. Строка берётся указателем на непрерывный участок:

// Обратный ход: вычисление компонент вектора решения
std::vector<double> x(n);

for (int i = n - 1; i >= 0; --i) {
    const double* rowI = &A[(size_t)i * n];
    double sum = b[i];
    for (int j = i + 1; j < n; ++j)
        sum -= rowI[j] * x[j];
    x[i] = sum;
}

Вывод реализуем через printf():

for (int i = 0; i < n; ++i) {
    double v = x[i];
    std::printf("%.6f\n", v);
}

Тут нужно озвучить один нюанс, особенно важный для участников соревнования. Мы с ребятами из Codenrock столкнулись с ним ещё во время отладки платформы для проверки решений участников. Дело в том, что нулевые значения после вычислений — не точный ноль, а крошечный остаток вроде -3e-13 и printf() напечатал бы его как «-0.000000». Поэтому всё, что округляется до нуля в шести знаках, нужно выводить как «0.000000» без знака. Однако, чтобы не изводить участников соревнования, мы решили делать такую проверку и соответствующую замену автоматически на стороне сервера. У нас нет цели за█████ть их до смерти.

Теперь повторяем замеры скорости (времени) и изучаем результат:

Вот это другое дело: ускорение чтения матрицы в 55,8 раза, с 7 921,7 до 141,9 мс, вот такое оно программирование под e2k в x86-мире. Прямой ход тоже удалось неплохо ускорить в 6,6 раза, с 16 286,9 до 2 445,8 мс. Теперь он читает память последовательно, лучше использует APB и значительно меньше страдает от задержек памяти.

Уважаемый читатель может заметить, что такое значительное ускорение чтения (ввода) матрицы A и вектора b — это прекрасно, но у нас тут соревнование по решению и оптимизации алгоритмических задач под Эльбрус (e2k), а не чемпионат по скоростному чтению матриц, и будет прав. Однако этот пример служит очень важным напоминанием о том, что:

• в реальных задачах входные данные редко (мягко говоря) будут попадать к вам в удобном для обработки (вычислений) виде (в том числе и на x86);

• в мире, где последние, как минимум, 20 лет всё строится вокруг x86, лёгкой прогулки для e2k ждать не стоит.

Время исполнения: 2 596,3 мс

Ускорение (общее): ×9,34

7. Шаг 2. Программная конвейеризация (подсказываем компилятору)

Прежде всего давайте посмотрим, где после шага 1 у нас сконцентрировались узкие (горячие) места:

Узкое место теперь одно: прямой ход занимает 94,2 % общего времени. Ввод после шага 1 сжался до 5,5 %, всё остальное — треть процента, и оптимизировать есть смысл только прямой ход.

Аналитический счёт операций. Главный кубический член прямого хода содержит около 2,67 млрд умножений-вычитаний, каждое из которых считается за две операции, то есть (2/3)·N³ ≈ 5,33 млрд операций. Делим на измеренное время и получаем эффективную скорость:

5,33·10⁹ операций / 2,446 с ≈ 2,2 ГФлоп/с

При этом пиковая скорость одного ядра Эльбрус-8СВ на вещественных числах двойной точности — 288 ГФлоп/с, на одно ядро — 36 ГФлоп/с, а для скалярного кода это как минимум вдвое меньше. То есть сейчас процессор недозагружен (мягко говоря). Почему? Давайте разбираться.

У прямого хода метода Гаусса границы всех трёх циклов определяются только величиной N, хотя проверка factor == 0.0 внутри цикла всё же зависит от данных. Сколько раз выполнится тело циклов — вопрос чисто арифметический, ответ на него есть ещё до первого запуска циклов. В этом и есть главное отличие и преимущество в нашем случае метода Гаусса над итерационными методами (Якоби, сопряжённых градиентов и т.д.). Если процессор (программный конвейер) недозагружен при том, что структура обрабатываемых данных «плоская», а сам алгоритм почти не зависит от них, то, вероятнее всего, компилятор, не сумев убедиться в независимости итераций (наложение rowI и rowK), не может эффективно заполнить (упаковать) операции в широкие команды.

По сути, итерации внутри циклов независимы, элемент j никак не связан с элементом j+1. Но компилятор этого не знает, для него rowI и rowK — это два указателя в один и тот же массив A, и формально запись в rowI[j] могла бы изменить память, которую вот-вот прочитает rowK[j+1] (такое перекрытие указателей называют наложением). Пока наложение не исключено, компилятор перестраховывается и не совмещает и не переставляет операции соседних итераций, а именно на этом стоит вся производительность широких команд в Эльбрусах.

Но в отличие от компилятора мы с тобой, дорогой читатель, знаем, что итерации внутри циклов независимы и наложения происходить не может. Поэтому нам нужно явно указать на это компилятору. Для этого мы используем два средства, которые описаны в Руководстве:

• квалификатор типа указателя __restrict, говорящий компилятору, что в данном фрагменте память, доступная через этот указатель, не перекрывается с памятью, доступной через другие указатели (раздел о разрыве зависимостей по памяти, глава 6 Руководства).

• директива #pragma ivdep, указывающая, что зависимостей между итерациями через память нет. Получив её, компилятор может свободнее совмещать операции для разных j в широких командах и строить программный конвейер, если это оказывается выгодно (глава 6.3 Руководства).

Добавляем их в код:

for (int k = 0; k < n; ++k) {
    double* __restrict rowK = &A[(size_t)k * n]; // <-- __restrict
    double pivot = rowK[k];
    for (int j = k; j < n; ++j)
        rowK[j] /= pivot;
    b[k] /= pivot;

    for (int i = k + 1; i < n; ++i) {
        double* __restrict rowI = &A[(size_t)i * n]; // __restrict
        double factor = rowI[k];
        if (factor == 0.0) continue;

#pragma ivdep // <-- #pragma ivdep
        for (int j = k; j < n; ++j)
            rowI[j] -= factor * rowK[j];
        b[i] -= factor * b[k];
    }
}

И делаем контрольный замер:

Прямой ход удалось ускорить в 1,4 раза с 2 445,8 мс до 1 743,5 мс.

Время исполнения: 1 903,3 мс

Ускорение (общее): ×12,73

8. Шаг 3. Блочная обработка (борьба с пропускной способностью памяти)

Возвращаемся к анализу узких (горящих) мест после оптимизаций шага 2:

Самой нагруженной частью с точки зрения затраченного времени по-прежнему остаётся прямой ход метода Гаусса — 1 743,5 мс (91,6 %), и это при том, что мы заметно улучшили эффективность заполнения операций в широкие команды, наладив программный конвейер. А если мы упираемся не в потолок вычислительной мощности процессора, то следующее место, куда стоит смотреть, — это память.

Для этого даже не надо делать сложный подсчёт байтов, которые прямой ход гоняет через память, чтобы понять, где ядро простаивает и чего оно ждёт (хотя кто я такой, дорогой читатель, чтобы запрещать тебе такое). У нас уже есть ресурсоёмкий элемент в лице матрицы A, размер которого мы можем легко корректировать, тем самым выявляя узкие места в работе с памятью. Отследим, как меняется производительность прямого хода (ГФлоп/с) при разных N, намеренно проследив за тем, как при этом ведёт себя кэш ядра:

Пока матрица A помещается в кэш (16 Мб), скорость шага 1 растёт, но стоит ей превысить 16 Мб L3, производительность проседает на треть, с 3,27 до 2,18 ГФлоп/с. Код не менялся, операций столько же, изменилось лишь то, откуда ядро берёт данные.

Такая просадка указывает на то, что ядро процессора простаивает в ожидании, пока очередная строка матрицы проделает долгий путь из оперативной памяти в регистры. Оптимизации, применённые нами на шаге 2, немного выправляют ситуацию (конвейер прячет задержку), но причину не убрали. Матрица A по-прежнему целиком прокачивается через память на каждом шаге исключения, а нужно, чтобы данные, однажды попавшие в кэш, использовались многократно, прежде чем их вытеснят.

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

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

А что, если сгруппировать эти вычитания друг с другом во времени? Возьмём группу из нескольких соседних шагов и пройдём по нижележащим строкам один раз, выполняя для каждой строки сразу все причитающиеся ей вычитания этой группы, пока строка загружена и горяча.

Ведущие строки группы при этом нужны постоянно, а значит, группа должна быть такой, чтобы её ведущие строки целиком жили в кэше: 64 строки по 16 КБ — около мегабайта, помещаются с большим запасом. Каждая строка матрицы будет проходить из памяти один раз на группу, а не один раз на шаг, и 2000 проходов превращаются в 2000/64 ≈ 31. Такой подход называют блочной обработкой или блочным LU-разложением (не путать с дюжиной других LU-разложений). Давайте перенесём всё в код.

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

const int NB = 64; // ширина группы (блока)

for (int kk = 0; kk < n; kk += NB) {
    // если kk+64 ещё в пределах матрицы, то граница kb равна kk+64, иначе n
    const int kb = (kk + NB < n) ? kk + NB : n;

    // ... три фазы обработки блока [kk, kb) ...
}

Группа охватывает шаги исключения с номерами от kk до kb (не включая). Выражение для kb аккуратно обрезает последнюю группу, если N не делится на 64 нацело, последний блок просто получается короче, и никаких особых случаев дальше не понадобится.

Подготовка блока столбцов. Здесь решается тонкость, которую мы отложили при выводе идеи: ведущие строки группы зависят друг от друга. Вторая ведущая строка станет собой только после того, как из неё вычтут первую; третья — после первых двух. Поэтому внутри узкой полосы столбцов [kk, kb) — её называют панелью — сначала выполняется самое обычное исключение, как в шаге 2:

// Подготовка блока столбцов
for (int k = kk; k < kb; ++k) {
    double* __restrict rowK = &A[(size_t)k * n];
    const double inv = 1.0 / rowK[k];

    for (int i = k + 1; i < n; ++i)
        A[(size_t)i * n + k] *= inv; // множитель l_ik на месте

    for (int i = k + 1; i < n; ++i) {
        double* __restrict rowI = &A[(size_t)i * n];
        const double f = rowI[k];

#pragma ivdep
        for (int j = k + 1; j < kb; ++j)
            rowI[j] -= f * rowK[j];
    }
}

Три отличия от оригинального прямого хода:

• внутренний цикл идёт только до kb, а не до конца строки. Всё, что правее панели, мы сознательно откладываем. Этим займёмся чуть позже, когда множители всей группы будут готовы;

• множители сохраняются. Величина l = a[i][k] / a[k][k] записывается прямо на место обнуляемого элемента, под диагональю. В обычном методе Гаусса её не используют, но нам она ещё понадобится для вычитания из «хвостов»;

• деление заменено умножением на обратное. Одно 1.0 / rowK[k] вместо тысяч делений в столбце (деление в разы дороже умножения).

Хвосты строк у блока столбцов. Строки блока столбцов с номерами kk+1 … kb−1 уже получили вычитания в пределах панели, но их правые части, столбцы от kb до конца, ещё нет. Доделываем, заодно обновляя вектор b теми же множителями:

// Хвосты строк у блока столбцов.
for (int k = kk + 1; k < kb; ++k) {
    double* __restrict rowK = &A[(size_t)k * n];

    for (int p = kk; p < k; ++p) {
        const double* __restrict rowP = &A[(size_t)p * n];
        const double f = rowK[p];

#pragma ivdep
        for (int j = kb; j < n; ++j)
            rowK[j] -= f * rowP[j];

        b[k] -= f * b[p];
    }
}

Порядок обхода здесь важен. Строка k обрабатывается только после того, как строки kk … k−1 уже завершены, поэтому цикл по k идёт по возрастанию, и каждая строка вычитает из себя только готовые строки, стоящие выше в той же панели. Множитель f = rowK[p] — это как раз сохранённое в предыдущем блоке «Подготовка блока столбцов» число l из-под диагонали.

Главный проход (одна строка + сразу вся группа). Каждая строка ниже блока столбцов загружается один раз и получает все 64 причитающихся ей вычитания подряд, пока она горячая:

// Главный проход (одна строка и сразу вся группа)
// По одному проходу на строку, сразу все NB множителей блока
// b обновляется той же схемой.
for (int i = kb; i < n; ++i) {
    double* __restrict rowI = &A[(size_t)i * n];

    for (int p = kk; p < kb; ++p) {
        const double* __restrict rowP = &A[(size_t)p * n];
        const double f = rowI[p];

#pragma ivdep
        for (int j = kb; j < n; ++j)
            rowI[j] -= f * rowP[j];

        b[i] -= f * b[p];
    }
}

Обратите внимание на движение данных. Внешний цикл по i идёт по строкам матрицы и каждая проходит из памяти один раз на всю группу. Внутренний цикл по p перебирает 64 ведущие строки блока столбцов, которые уже в кэше и остаются там. Между двумя обращениями к rowP через кэш проходит не вся матрица, как раньше, а одна строка.

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

for (int i = n - 1; i >= 0; --i) {
    const double* rowI = &A[(size_t)i * n];
    double sum = b[i];
    for (int j = i + 1; j < n; ++j)
        sum -= rowI[j] * x[j];
    x[i] = sum / rowI[i]; // ← новое: делим на диагональ
}

Делаем замеры:

Прямой ход удалось ускорить в 1,3 раза с 1 743,5 мс до 1 342,5 мс. Повторим изменения производительности прямого хода при разных N, немного расширив их для отслеживания динамики:

Время исполнения: 1 514,5 мс

Ускорение (общее): ×16,0

9. Шаг 4. Разворот с объединением

Смотрим, где теперь у нас сконцентрировались узкие (горячие) места:

Прямой ход по-прежнему занимает 88,6 % времени вычислений, при этом бутылочное горлышко в виде ограничения пропускной способности ОЗУ (данные переиспользуются из кэша) мы устранили на предыдущем шаге.

Но достигнутая производительность после оптимизаций предыдущего шага хоть и не прекращает рост, но заметно замедляется после N=6000 и явно не имеет тенденции достичь даже 5 ГФлоп/с, что по-прежнему далеко от возможностей ядра процессора. При этом единственным действительно узким (горячим) местом в коде остаётся одна строка:

rowI[j] -= f * rowP[j];

В ней всего две операции: умножение и вычитание. А вот обращения к данным три: чтение rowP[j], чтение rowI[j] и запись rowI[j] обратно. Полтора обращения к данным на каждую операцию при том, что за такт ядро может исполнить до шести вещественных операций, но операции чтения размещаются только в четырёх каналах, а записи — только в двух (глава 4.4 Руководства). То есть АЛУ, вероятнее всего, простаивают, а возможности размещения операций чтения и записи становятся нашим предполагаемым бутылочным горлышком.

Давайте посмотрим, что с этим можно сделать. Чтение ведущей строки rowP[j] не убрать: чтобы произвести вычитание, нужно её знать, по одному чтению на каждое вычитание, и это минимум. А вот пара «чтение rowI[j] и запись rowI[j]» повторяется при каждом вычитании. Проследим судьбу одного элемента за группу — приём, знакомый нам по шагу 3, только этажом ниже, между кэшем и регистрами: элемент вычитает из себя составляющие 64 ведущих строк и ради этого 128 раз пересылается между кэшем и регистрами, хотя итог всех пересылок — одно-единственное число: исходное значение минус сумма составляющих. Но ведь сумму можно накапливать в регистре по частям. Вычтем за две пересылки (пару чтение-запись) вклад сразу двух ведущих строк, и этим одна пара «чтение-запись» обслужит два вычитания: обращений станет 1 на операцию вместо 1,5.

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

for (int p = kk; p + 1 < kb; p += 2) {
    const double* __restrict rowP0 = &A[(size_t)(p + 0) * n];
    const double* __restrict rowP1 = &A[(size_t)(p + 1) * n];
    const double f0 = rowI[p + 0];
    const double f1 = rowI[p + 1];

#pragma ivdep
    for (int j = kb; j < n; ++j)
        rowI[j] -= f0 * rowP0[j] + f1 * rowP1[j];

    b[i] -= f0 * b[p + 0] + f1 * b[p + 1];
}

Замеряем:

Падение в 1,40 раза. Теперь осталось понять оптимальное количество ведущих строк. Подберём замерами (код намеренно не показываю за ненадобностью: что для 4, что для 8 он практически идентичен):

При четырёх ведущих строках получаем оптимальный вариант. Удвоение до восьми добавляет лишь 1,8 %, обращения перестают быть самым узким местом. Берём по четыре: выигрыш почти максимальный, код вдвое проще, да и ширина блока 64 делится на 4 нацело. Фактически в регистрах копится сумма вычитаемых составляющих, четыре произведения f×rowP[j] складываются, не покидая регистров, и лишь готовый итог уходит в ячейку. Вот такой вот разворот с объединением.

Делаем замеры:

Прямой ход удалось ускорить в 1,53 раза, с 1 342,5 до 878,8 мс.

Время исполнения: 1 051,5 мс

Ускорение (общее): ×23,05

10. Шаг 5. Многопоточность

Снова смотрим, где теперь у нас сконцентрировались узкие (горячие) места:

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

Стоит заметить, что не каждый алгоритм хорошо распараллеливается, однако это не наш случай. В главном проходе каждая нижележащая строка получает свои вычитания независимо от остальных: поток, обрабатывающий строку i, пишет только в неё и в элемент b[i], а ведущие строки блока все потоки лишь читают. Аналогично устроены оба цикла по строкам в подготовке блока столбцов. Две части прямого хода распараллеливать не следует. Хвосты строк блока столбцов малы (63 короткие строки на блок): накладные расходы на раздачу работы превысили бы выигрыш. Обратный ход содержит зависимость по данным: каждая неизвестная вычисляется через уже найденные и занимает 5,6 мс, поэтому простое распределение его цикла между потоками невозможно и не нужно по величине.

OpenMP. Для распараллеливания вычислений будем использовать поддерживаемый компилятором и платформой Эльбрус (e2k) программный интерфейс (стандарт) OpenMP. Число потоков задаётся переменной окружения OMP_NUM_THREADS. Мы намеренно будем использовать семь потоков, оставив одно ядро системе.

Подготовка блока столбцов:

for (int k = kk; k < kb; ++k) {
    double* __restrict rowK = &A[(size_t)k * n];
    const double inv = 1.0 / rowK[k];

#pragma omp parallel for schedule(static) // строки — по потокам
    for (int i = k + 1; i < n; ++i)
        A[(size_t)i * n + k] *= inv;

#pragma omp parallel for schedule(static) // и здесь
    for (int i = k + 1; i < n; ++i) {
        double* __restrict rowI = &A[(size_t)i * n];
        const double f = rowI[k];

#pragma ivdep
        for (int j = k + 1; j < kb; ++j)
            rowI[j] -= f * rowK[j];
    }
}

Директива #pragma omp parallel for перед циклом сообщает: итерации независимы, распредели их между потоками, а schedule(static) задаёт статическую раздачу — заранее равными долями, без динамического диспетчера, для однородной по трудоёмкости работы это самый дешёвый вариант.

Главный проход (внутри — развёрнутое по четырём ведущим строкам ядро шага 4, без изменений):

#pragma omp parallel for schedule(static) // строки — по потокам
for (int i = kb; i < n; ++i) {
    double* __restrict rowI = &A[(size_t)i * n];
    ... // тело шага 4
}

При этом мы можем посмотреть оптимальное количество потоков для вычислений:

Разница между семью и восемью потоками — 0,6 %. Восьмой поток почти ничего не добавляет, поэтому и выбор семи потоков оказывается изначально удачным (бывают же совпадения). Делаем полный прогон и строим таблицу распределения времени при семи потоках:

Время исполнения: 451,3 мс

Ускорение (общее): ×53,7

11. Заключение

Безусловно, это не все возможные оптимизации, которые можно выжать из кода, но дальше может встать вопрос их реальной необходимости с учётом сложности и объёма кода, который они добавят, а мы с вами остановимся на достигнутом. Общее ускорение ×53,7, а если считать только вклад вычислительных приёмов (от шага 1, когда ввод уже был приведён в порядок), то ×5,75. В конце концов, моё дело было сделать разбор оптимизации кода под платформу Эльбрус (e2k) для участников соревнования. С чем мы с тобой, дорогой читатель, успешно, как мне кажется, справились.

А если тебе от 18 до 24 лет включительно, то обязательно присоединяйся к соревнованию: cup.serpas.ru

И нет, проживание на Эльбрусе не оплачивается.

---

Не забудьте подписаться на наше сообщество энтузиастов в Telegram, VK и MAX.