Оптимизация ядра диагонального умножения матрицы на вектор (? Diamv)

В качестве (совершенно необязательного) задания для вводного курса программирования на C ++ я пытаюсь реализовать диагональное умножение матрицы на вектор (?diamv) ядро, т.е. математически
$$ mathbf {y} leftarrow alpha mathbf {y} + beta mathbf {M} mathbf {x} $$
для диагонально кластерной матрицы $ mathbf {M} $, плотные векторы $ mathbf {x} $ и $ mathbf {y} $, и скаляры $ alpha $ и $ beta $. Я считаю, что могу обоснованно мотивировать следующие предположения:

  1. Процессоры, выполняющие вычислительные потоки, способны выполнять расширение набора инструкций SSE4.2 (но не обязательно AVX2),
  2. Схема доступа матрицы $ mathbf {M} $ не влияет на вычисления, и поэтому не нужно учитывать временную локальность кеша между вызовами ядра,
  3. Матрица $ mathbf {M} $ не помещается в кеш, очень диагонально сгруппирован с диагональным шаблоном, который известен во время компиляции, и квадрат,
  4. Матрица $ mathbf {M} $ не содержит регулярно встречающихся последовательностей по диагоналям, которые допускали бы сжатие по оси,
  5. Для структуры матрицы не существует функции переупорядочивания. $ mathbf {M} $ это привело бы к продукту без кеширования с более низкой стоимостью, чем идеальный алгоритм, оптимизированный для многоуровневой памяти,
  6. Исходные данные выровнены по адекватной границе,
  7. OpenMP, выбранный из-за его популярности, доступен для обеспечения параллелизма с общей памятью. Никакого параллелизма с распределенной памятью не требуется, поскольку предполагается, что алгоритм декомпозиции предметной области, например DP-FETI, будет декомпозировать обработку до уровня узла из-за типичного размера проблемы.

Проведя обзор литературы, я пришел к следующим выводам относительно ее дизайна и реализации (это краткое изложение, с возрастающей степенью детализации, с обширным обзором литературы, доступным по запросу для экономии места):

  1. «Для достижения высокой производительности параллельная реализация умножения разреженных матриц на вектор должна поддерживать масштабируемость», согласно Уайту и Садаяппану, 1997.
  2. Схема хранения диагональной матрицы,
    $$ DeclareMathOperator { vec} {vec} DeclareMathOperator { val} {val} vec left ( val {(i, j)} Equiv a_ {i, i + j} right) $$
    куда $ vec $ – это оператор векторизации матрицы, который получает вектор путем наложения столбцов матрицы операнда друг на друга. Сохраняя матрицу в этом формате, я считаю, что расположение кеша является максимально оптимальным для обеспечения параллелизма по строкам. Разделение шахматной доски сводится к построчному для диагональных матриц. Кроме того, это позволяет повторно использовать исходный вектор, что необходимо, если матрица не используется повторно, пока она еще находится в кеше (Frison 2016).
  3. Я считаю, что вышесказанное должно всегда держите, до того как векторизация вообще рассматривается? Нерегулярные заполненные области матрицы, то есть верхний левый и нижний правый, можно обрабатывать отдельно без дополнительных затрат в асимптотическом смысле (поскольку матрица сгруппирована по диагонали и очень велика).
  4. Поскольку доступ к этой матрице является линейным, программная предварительная выборка не требуется. Я все равно включил его для проверки кода в том месте, которое считал наиболее логичным.

Следующий фрагмент представляет мои лучшие усилия с учетом вышеупомянутого:

#include <algorithm>
#include <stdint.h>
#include <type_traits>

#include <xmmintrin.h>
#include <emmintrin.h>

#include <omp.h>

#include "tensors.hpp"


#define CEIL_INT_DIV(num, denom)        1 + ((denom - 1) / num)

#if defined(__INTEL_COMPILER)
#define AGNOSTIC_UNROLL(N)              unroll (N)
#elif defined(__CLANG__)
#define AGNOSTIC_UNROLL(N)              clang loop unroll_count(N)
#elif defined(__GNUG__)
#define AGNOSTIC_UNROLL(N)              unroll N
#else
#warning "Compiler not supported"
#endif

/* Computer-specific optimization parameters */
#define PREFETCH                        true
#define OMP_SIZE                        16
#define BLK_I                           8
#define SSE_REG_SIZE                    128
#define SSE_ALIGNMENT                   16
#define SSE_UNROLL_COEF                 3


namespace ranges = std::ranges;


/* Calculate the largest absolute value ..., TODO more elegant? */
template <typename T1, typename T2>
auto static inline largest_abs_val(T1 x, T2 y) {
    return std::abs(x) > std::abs(y) ? std::abs(x) : std::abs(y);
}


/* Define intrinsics agnostically; compiler errors thrown automatically */
namespace mm {
    /* _mm_load_px - [...] */
    inline auto load_px(float const *__p) { return _mm_load_ps(__p); };
    inline auto load_px(double const *__dp) { return _mm_load_pd(__dp); };

    /* _mm_store_px - [...] */
    inline auto store_px(float *__p, __m128 __a) { return _mm_store_ps(__p, __a); };
    inline auto store_px(double *__dp, __m128d __a) { return _mm_store_pd(__dp, __a); };

    /* _mm_set1_px - [...] */
    inline auto set_px1(float __w) { return _mm_set1_ps(__w);};
    inline auto set_px1(double __w) { return _mm_set1_pd(__w); };

    /* _mm_mul_px - [...] */
    inline auto mul_px(__m128 __a, __m128 __b) { return _mm_mul_ps(__a, __b);};
    inline auto mul_px(__m128d __a, __m128d __b) { return _mm_mul_pd(__a, __b); };
}


namespace tensors {
    template <typename T1, typename T2>
    int diamv(matrix<T1> const &M, 
              vector<T1> const &x,
              vector<T1> &y,
              vector<T2> const &d,
              T1 alpha, T1 beta) noexcept {
        /* Initializations */
        /* - Compute the size of an SSE vector */
        constexpr size_t sse_size =  SSE_REG_SIZE / (8*sizeof(T1));
        /* - Validation of arguments */
        static_assert((BLK_I >= sse_size && BLK_I % sse_size == 0), "Cache blocking is invalid");
        /* - Reinterpretation of the data as aligned */
        auto M_ = reinterpret_cast<T1 *>(__builtin_assume_aligned(M.data(), SSE_ALIGNMENT));
        auto x_ = reinterpret_cast<T1 *>(__builtin_assume_aligned(x.data(), SSE_ALIGNMENT));
        auto y_ = reinterpret_cast<T1 *>(__builtin_assume_aligned(y.data(), SSE_ALIGNMENT));
        auto d_ = reinterpret_cast<T2 *>(__builtin_assume_aligned(d.data(), SSE_ALIGNMENT));
        /* - Number of diagonals */
        auto n_diags = d.size();
        /* - Number of zeroes for padding TODO more elegant? */
        auto n_padding_zeroes = largest_abs_val(ranges::min(d), ranges::max(d));
        /* - No. of rows lower padding needs to be extended with */
        auto n_padding_ext = (y.size() - 2*n_padding_zeroes) % sse_size;
        /* - Broadcast α and β into vectors outside of the kernel loop */
        auto alpha_ = mm::set_px1(alpha);
        auto beta_ = mm::set_px1(beta);

        /* Compute y := αy + βMx in two steps */
        /* - Pre-compute the bounding areas of the two non-vectorizable and single vect. areas */
        size_t conds_begin[] = {0, M.size() - (n_padding_ext+n_padding_zeroes)*n_diags};
        size_t conds_end[] = {n_padding_zeroes*n_diags, M.size()};
        /* - Non-vectorizable areas (top-left and bottom-right resp.) */
#pragma AGNOSTIC_UNROLL(2)
        for (size_t NONVEC_LOOP=0; NONVEC_LOOP<2; NONVEC_LOOP++) {
            for (size_t index_M=conds_begin[NONVEC_LOOP]; index_M<conds_end[NONVEC_LOOP]; index_M++) {
                auto index_y = index_M / n_diags;
                auto index_x = d[index_M % n_diags] + index_y;
                if (index_x >= 0)
                    y_[index_y] = (alpha * y_[index_y]) + (beta * M_[index_M] * x_[index_x]);
            }
        }
        /* - Vectorized area - (parallel) iteration over the x parallelization blocks */
#pragma omp parallel for shared (M_, x_, y_) schedule(static)
        for (size_t j_blk=conds_end[0]+1; j_blk<conds_begin[1]; j_blk+=BLK_I*n_diags) {
            /* Iteration over the x cache blocks */
            for (size_t j_bare = 0; j_bare < CEIL_INT_DIV(sse_size, BLK_I); j_bare++) {
                size_t j = j_blk + (j_bare*n_diags*sse_size);
                /* Perform y = ... for this block, potentially with unrolling */
                /* *** microkernel goes here *** */
#if PREFETCH
                /* __mm_prefetch() */
#endif
            }
        }

        return 0;
    };
}
 

Некоторые важные примечания:

  1. tensors.hpp – это простая библиотека только для заголовков, которую я написал для случая, чтобы действовать как единый слой абстракции для тензоров различного порядка (с CRTP), имеющих разные схемы хранения. Он также содержит псевдонимы, например, векторов и плотных матриц.

  2. Я считаю, что для микроядра есть две возможности

    а. Итерировать линейно по векторизованной матрице в каждом блоке кеша; это составило бы построчную итерацию по матрице $ mathbf {M} $ внутри каждого блока кэша и, следовательно, скалярный продукт. Насколько мне известно, скалярные произведения неэффективны в плотных матрично-векторных произведениях как из-за зависимостей данных, так и из-за того, как внутренние функции разлагаются на μops.

    б. Итерации по строкам в блоках кеша в векторизованной матрице, что составляет итерацию по диагоналям в матрице $ mathbf {M} $ в каждом блоке кеша. Из-за того, как матрица $ mathbf {M} $ хранится, то есть в его векторизованной форме, это повлечет за собой расходы на широковещательную передачу чисел с плавающей запятой (что, насколько мне известно, является сложным вопросом), но позволит строкам внутри блоков выполняться параллельно.

    Боюсь, что упустил другие, лучшие варианты. Это основная причина для открытия этого вопроса. Я полностью застрял. Кроме того, я считаю, что различия в том, насколько хорошо используются повторно векторы источника / назначения, слишком близки, чтобы их можно было назвать. Кто-нибудь знает, как я могу подойти к более глубокому пониманию этого вопроса?

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

  4. Есть ли способ сделать “уродливый” невекторизуемый код более элегантным?

Проверяя вышесказанное, я чувствую себя полным любителем; все отзывы очень ценятся. Заранее спасибо.

1 ответ
1

Функции в mm пространство имен, похоже, не имеет никакого значения по сравнению с прямым вызовом мнемоники. Я бы избавился от них, чтобы код был простым. Кроме того, знаете ли вы правила языка, касающиеся использования ведущих символов подчеркивания в именах наизусть? Если нет, я бы посоветовал вообще избегать их использования, в противном случае я бы предостерег от этого в любом случае, потому что не все знают эти неясные правила (ссылаясь на аргументы, начинающиеся с двойного подчеркивания).

Вы приложили значительные усилия, чтобы сделать это быстро. Вы проверили это? Вы знаете, действительно ли ваш код быстрее, чем наивное поэлементное произведение двух векторов? (Предполагая, что матрица действительно диагональна, как вы сначала говорите, а не сгруппирована, как вы позже скажете). Запуск и обмен данными между потоками, которые распределяют рабочую нагрузку, являются накладными расходами, которые, вероятно, не окупятся для небольших матриц. Вы спросите, насколько мала? Зависит от вашего процессора. Тест, тест, тест! В дополнение к сравнению с наивной реализацией (опять же, предполагая, что матрица на самом деле диагональная), я бы сравнил с существующей библиотекой, например Собственный или же Армадилло чтобы увидеть, где находится ваш код в общей схеме вещей.

Я понимаю, что это упражнение, но я просто отмечу, что в противном случае вам следует полностью использовать для этого одну из существующих библиотек.

Я не уверен, действительно ли вам помогают прагмы развертывания цикла, компилятор должен уже выполнять развертывание в нужной степени. Похоже, что эта прагма в лучшем случае ничего не делает или в худшем случае вызывает развертывание, когда вы этого не хотите (например, при отладке или если вы хотите, чтобы размер двоичного файла был меньше). Имея это в виду, это добавляет дополнительную сложность коду. Я бы протестировал с использованием -O2 и -O3 и без них и увидел разницу.

Я часто вижу это с начинающими программистами на C / C ++, они хотят делать все причудливые вещи, чтобы сделать это быстро, вручную развертывать, извлекать общие подвыражения и т. Д. Чего они обычно не знают, так это того, что компилятор уже делает все это и гораздо больше, и он делает это лучше, чем они, или, по крайней мере, не хуже почти во всех случаях. Компиляторы, которые у нас есть сегодня, – это не старый компилятор вашего дедушки, они потрясающе генерируют эффективный код до такой степени, что лучшее, что вы можете сделать, – это четко выразить свой алгоритм и намерения и позволить компилятору генерировать код.

Показательный пример, который быстрее:

int factorial(int n){
    if(1==n)
        return n;
    return n * factorial(n-1);
}

или же

int factorial (int n){
    int ans=1;
    while (n>1){
        ans*=n;
        n--;
    }
    return ans;
}

?

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

Знаете ли вы, что компилятор может генерировать SSE и другие векторные инструкции, если вы скажете ему, что это разрешено? См. Например
этот сайт. Я уверен, что вы можете найти эквивалент в clang и msvc.

Все это для того, чтобы сказать … Вам действительно нужна сложность, которую вы добавляете? Я не знаю, потому что вы не показали мне никаких тестов. 🙂

Что касается предварительной выборки, я никогда не видел случая, чтобы ручная вставка инструкции предварительной выборки оказывала сколько-нибудь заметное влияние, предварительные выборки ЦП обычно очень хороши.

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

Надеюсь, это поможет!

(Написано на телефон, извините за опечатки и автокоррекцию)

  • Большое спасибо! 1) Смысл пространства имен mm заключается в перегрузке встроенных функций, например _mm_load_ps и _mm_load_pd 2) __arg следует соглашениям в файлах заголовков 3) извинения, матрица кластеризована по диагонали, а не строго по диагонали 4) матрицы относительно большие, в элементах> 300k 5) Я понимаю, что вы имеете в виду под словом optim. компиляторы, будут изучать, Hager & Wellein заявляют, что такие случаи трудно оптимизировать автоматически / с подсказками 6) не удалось найти реализацию diamv с открытым исходным кодом, только Intel mkl 7) Я проведу тестирование, как только закончу (автоматическую) векторизацию. Еще раз спасибо!

    – user169291

  • 2) соглашение зарезервировано для внутреннего использования компилятором, вы не должны следовать ему, потому что вы не компилятор, и вы можете столкнуться с конфликтами имен с внутренними макросами компилятора и тем, что не может даже не отображаться как ошибка компилятора. Это то, что я имел в виду, говоря о языковых правилах, касающихся использования ведущих подчеркиваний. 5) За последние годы произошли большие улучшения, не знаю, сколько лет вашему эталону. Вот почему вы должны тестировать 🙂 6) для начала вы можете сравнить с разреженными матрицами в Eigen. 🙂

    – Эмили Л.

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

    – Эмили Л.

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *