В настоящее время многоядерные процессоры с гетерогенными архитектурами, в которых сочетаются ядра с различной производительностью, становятся всё более и более распространенными. Если ещё пару лет назад такие архитектуры были в основном распространены в мобильном сегменте (см. ARM BIG.little), то с анонсом в 2022 году компанией Intel процессоров 12-го поколения линейки Intel Core, такие процессоры стали распространяться в сегменте десктопов и рабочих станций. Однако, до сих пор остается открытым вопрос — необходимо ли каким‑то специальным образом учитывать особенности данных архитектур для достижения максимальной многопоточной производительности?

Так, например, компания Intel в своей документации к библиотеке MKL рекомендует производить вычисления только на производительных ядрах. Хотя и замечает, что по мере роста размера задач может быть оправдано использование энергоэффективных ядер. Попробуем и мы разобраться в данном вопросе.

Как уже стало понятно из названия - в качестве технологии на примере которой будем проводить эксперименты выберем OpenMP. Данная технология представляет собой совокупность директив компилятора и библиотечных процедур для языков C, C++, Fortran. “Конёк” OpenMP - параллельное выполнение цикла for. Если итерации цикла могут быть выполнены независимо, то разработчику достаточно добавить строчку #pragma omp parallel for перед циклом, чтобы добиться ускорения за счет использования множества ядер в рамках одной системы. Благодаря этой особенности OpenMP стал особенно популярен в вычислительных задачах математики и физики. Исходя из этого выберем и тестовую задачу - реализуем алгоритм решения СЛАУ методом сопряженных градиентов. Горячим участком данного метода является этап перемножения матрицы на вектор, его и будем распараллеливать описанным выше способом.

void conjugateGradient(const float *A, const float *b, float *x,
                       std::ptrdiff_t n, std::ptrdiff_t maxIterations = 1000) {
    auto r = new (std::align_val_t(64)) float[n];
    auto p = new (std::align_val_t(64)) float[n];
    auto Ap = new (std::align_val_t(64)) float[n];

    // Инициализация: x = 0, r = b - A*x = b, p = r
    std::fill(x, x + n, 0.f);
    std::copy(b, b + n, r);
    std::copy(r, r + n, p);

    float rsold = 0.f;
    for (std::ptrdiff_t i = 0; i < n; ++i) {
        rsold += r[i] * r[i];
    }

    float norm_b = 0.f;
    for (std::ptrdiff_t i = 0; i < n; ++i) {
        rsold += b[i] * b[i];
    }

    for (std::ptrdiff_t iter = 0; iter < maxIterations; ++iter) {

        #pragma omp parallel for
        for (std::ptrdiff_t i = 0; i < n; ++i) {
            float res = 0.f;
            for (std::ptrdiff_t j = 0; j < n; ++j) {
                res += A[i * n + j] * p[j];
            }
            Ap[i] = res;
        }

        float tmp = 0.f;
        for (std::ptrdiff_t i = 0; i < n; ++i) {
            tmp += Ap[i] * p[i];
        }

        float alpha = rsold / tmp;

        // x = x + alpha * p
        // r = r - alpha * Ap
        float rsnew = 0.f;
        for (std::ptrdiff_t i = 0; i < n; ++i) {
            x[i] = x[i] + alpha * p[i];
            r[i] = r[i] - alpha * Ap[i];
            rsnew += r[i] * r[i];
        }

        float beta = rsnew / rsold;

        // p = r + beta * p
        for (std::ptrdiff_t i = 0; i < n; ++i) {
            p[i] = r[i] + p[i] * beta;
        }
        rsold = rsnew;
    }

    delete[] r;
    delete[] p;
    delete[] Ap;
}

И тут самое время упомянуть про необязательный параметр используемой нами прагмы - schedule. Этот параметр отвечает за способ распределения итераций по потоком. Стандартом определен следующий набор базовых способов:

• static - распределение итераций цикла между потоками происходит заранее, равными или почти равными блоками, что минимизирует накладные расходы, но может привести к дисбалансу нагрузки;

• dynamic - итерации цикла динамически распределяются между потоками небольшими блоками (по умолчанию — по 1), что улучшает балансировку нагрузки, но увеличивает накладные расходы на синхронизацию;

• guided - итерации сначала распределяются крупными блоками, которые постепенно уменьшаются, что сочетает преимущества static и dynamic;

• auto - распределение итераций определяется компилятором или системой;

• runtime - тип распределения (static, dynamic или guided) задаётся во время выполнения через переменную окружения OMP_SCHEDULE, что позволяет настраивать поведение без перекомпиляции.

Так как этот параметр не является обязательным, в случае его отсутствия компилятор сам решает какое распределение использовать. Как показывает практика, большинство реализаций в таком случае использует статическое распределение итераций.

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

Эксперименты

Будем сравнивать три версии реализованной процедуры: 

• «Static P+E» — статическое распределение итераций, работают все ядра (поведение по умолчанию);

• «Static P» — статическое распределение, работают только производительные ядра. Будем использовать переменную окружению OMP_PLACES, чтобы контролировать на каких именно ядрах исполняются потоки.

• «Dynamic P+E» — динамическое распределение, работают все ядра. 

Для полноты сравнения также будем варьировать размер матрицы от 500×500 до 16000×16000 элементов типа float. 

Запускать наше приложение будем на компьютере с процессором Intel Core i9-14900K, который относится к поколению Raptor Lake Refresh и состоит из 8 производительных P‑ядер на микроархитектуре Raptor Cove и 16 энергоэффективных E‑ядер на микроархитектуре Gracemont. Операционная система — GNU/Linux, дистрибутив — Ubuntu 24.04 LTE, ядро версии 6.14.0, компилятор — gcc 13.3.

Получаем следующие результаты. На малых размерах матрицы (до 2000×2000) наблюдаем лидерство подхода с использованием исключительно P‑ядер, что объясняется большими накладными расходами на динамическую балансировку. По мере роста размера матрицы доля накладных расходов уменьшается и подход с использованием всех доступных ядер вкупе с динамической балансировкой вырывается вперед. Разница во времени исполнения может достигать 30% в сравнении с двумя другими подходами. 

Хоть данная формула и игнорирует многие аспекты реального железа (например, напрочь игнорируется пропускная способность памяти), но тем не менее позволяет выполнить некоторую оценку «сверху». Отметим, что при некоторых размерах матрицы мы достигаем теоретической оценки — успех!

Можно ли получить результат лучше?

У простой версии с динамическим распределением итераций есть одна проблема — неэффективная работа с кэшем. Как уже говорилось ранее, по умолчанию единицей планирования является одна итерация, а значит, каждый поток читает только относительно небольшой непрерывный кусок памяти, который относится к ровно одной строке матрицы. И что намного хуже — записывает по одному элементу в результирующий вектор, что увеличивает вероятность так называемый false sharing — ситуации, когда несколько потоков выполняют запись в одну кэш‑строку и процессору приходится тратить дополнительное время на их синхронизацию.

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

И теперь можем построить наш итоговый график. Диапазон размеров матрицы, в котором мы достигаем теоретического ускорения существенно увеличился, более того при некоторых размерах мы его превосходим. С чем это может быть связано? По всей видимости, здесь мы сталкиваемся с тем же эффектом, когда в традиционных архитектурах наблюдаем сверхлинейное ускорение из‑за увеличения размера доступного кэша при увеличении числа используемых ядер. Но всё равно глядя на график ускорения остается несколько вопросов.

Первый из них — почему производительность версии «Dynamic P+E» так быстро падает после некоторого размера? Для того чтобы ответить на этот вопрос воспользуемся утилитой perf. Будем анализировать отношение счетчиков cache‑misses к cache‑references (по сути — процент кэш‑промахов) и LLC‑load‑misses к LLC‑loads (к сожалению, это уже не процент промахов к кэшу последнего уровня, а просто отношение числа промахов к числу загрузок).

По всей видимости здесь действует два разнонаправленных фактора: с одной стороны версия использующая ядра обоих типов располагает суммарно большим объемом кэша, но большее число одновременно работающих потоков начинают уже мешать друг к другу при доступе к общему ресурсу — кэшу последнего уровня. Как итог мы видим более эффективное использование кэша версией «Static P», за исключением диапазона от 3000 до 5000 тысяч элементов. Напомню, именно в этом диапазоне ранее мы наблюдали ускорение выше теоретической оценки.

По мере дальнейшего увеличения размера матрицы число кэш‑промахов в обоих версиях выравнивается и приближается к 100%. Ограничивающим фактором в работе приложения становится пропускная способность памяти, становятся практически равны и их производительности, но они по‑прежнему работают эффективнее версии «Static P+E» приблизительно на 10-15%.

Почему же «практически», а не равную? Если приблизить показатели эффективности на больших размерах задачи, то можно заметить, что версии с динамическим распределением итераций вне зависимости от того используют они энергоэффективные ядра или нет, оказываются на 3% быстрее версии «Static P». 

По всей видимости ещё одними источниками дисбаланса в данном процессоре, помимо ядер разных типов, является особенности работы алгоритма динамического изменения тактовой частоты (aka Turbo Boost). Косвенно это может подтвердить утилита s-tiu. Можем видеть, что частоты первого ядра в среднем оказываются выше, а само ядро больше простаивает при работе версии “Static P”. При динамическом же распределении итераций, можем наблюдать более ровную загрузку ядер.

Скриншоты s-tui

Итоги

Для достижения максимальной производительности приложений для современных центральных процессоров компании Intel необходимо принимать во внимание их ключевую особенность — наличие производительных и энергоэффективных ядер. В зависимости от вашего алгоритма и объема данных с которыми он оперирует, возможно, вам будет достаточно только отказаться от использования энергоэффективных ядер с помощью переменной окружения OMP_PLACES, а возможно вам придется изменять исходный код приложения для смены политики планирования исполнения потоков. В противном случае, вы рискуете недосчитаться вплоть до 30% производительности.

Послесловие

Идея статьи появилась в ходе моего преподавания курса Основ параллельного программирования на факультете Информационных технологий НГУ. В рамках курса из года в год студентам приходится выполнять одни и те же лабораторные работы на знакомство с технологиями OpenMP и MPI. Тем интереснее было обнаружить, что стандартный ответ на вопрос «Какой подход к распределению итераций цикла между потоками стоит выбрать в такой‑то задаче?» заиграл новыми красками, когда у студентов стали появляться первые ноутбуки с гетерогенной архитектурой CPU. Причем у разных студентов как на зло получались разные ответы на данный вопрос — незначительные отличия в реализациях и объеме данных вносили свой вклад, это и стало толчком к тому, чтобы провести собственное исследование.