Какой самый быстрый способ вычисления смешанных вещественно-сложных матрично-векторных произведений?

Question

Я сталкиваюсь с проблемой как можно быстрее вычислить произведения матрицы-вектора, где матрица строго вещественная, а вектор сложный. Чтобы быть более точным, реальная матрица является (анти-) симметричной c и имеет разреженную структуру ненулевых плотных блоков.

До сих пор моя стратегия состояла в том, чтобы разделить вектор на его реальный и мнимый разделить и вычислить для каждого плотного блока произведения матрицы-вектора как для вещественных, так и для мнимых частей. Поскольку матрица является (анти-) симметричной c, я намерен одновременно вычислять произведения для блока и его транспонирование, поэтому я могу просто повторно использовать строку кэша, в которой находится матрица. Поэтому для каждого блока я вычислить 4 матрично-векторных произведения для блока и его транспонирование для каждой действительной и мнимой части.

Мой код для вычисления этих 4 произведений для одного блока в конечном итоге выглядит следующим образом:

#define no_alias __restrict__

template <typename VecType, typename MatType>
void trans_mul( const VecType* const no_alias re_in,
                VecType* const no_alias re_out, 
                const VecType* const no_alias im_in,
                VecType* const no_alias im_out, 
                const VecType* const no_alias re_tin,
                VecType* const no_alias re_tout,
                const VecType* const no_alias im_tin,
                VecType* const no_alias im_tout,
                const MatType* no_alias mat, // the matrix block
                const int rows, 
                const int cols)
{
    for(int j = 0; j < cols; ++j) {
        for(int i = 0; i < rows; ++i) {
            const auto m = *mat++;        // this is mat[i, j]
            re_tout[j]  += m * re_tin[i]; // transposed
            im_tout[j]  += m * im_tin[i]; // transposed
            re_out[i]   -= m * re_in[j];
            im_out[i]   -= m * im_in[j];
        }
    }
}

Типичный размер матрицы порядка 10 ^ 2. Я компилирую свой код, используя G CC 9.2.1 с -Ofast -march=native. Из вывода сборки я вижу, что компилятор автоматически векторизуется и использует инструкции SIMD.

Я конкурирую с аналогичным кодом, написанным на Fortran, который все еще работает примерно на 25% быстрее , Конечно, мой код чрезвычайно наивен, но все же я не смог придумать ничего более быстрого, так как агрессивная оптимизация кажется очень эффективной. Я также пытался использовать четыре cblas_dgemv звонка, что, однако, было значительно медленнее, чем мой наивный подход. Есть ли что-то еще, что я могу сделать, или есть какая-нибудь полезная процедура BLAS, которая могла бы удовлетворить мой случай?

Answer 1

Для довольно больших матриц (например,> = 1k) вы можете использовать блокировку регистров для повышения производительности (уменьшает количество загрузки / сохранения памяти по сравнению с арифметическими c операциями). Для маленьких матриц трудно сделать что-то лучше, чем исходный код.

Вот результирующий код с блокировкой регистров:

#define no_alias __restrict__
#define mini(a, b) (((a)<(b)) ? (a) : (b))

template <typename VecType, typename MatType>
void trans_mul_v2( const VecType* const no_alias re_in,
                VecType* const no_alias re_out, 
                const VecType* const no_alias im_in,
                VecType* const no_alias im_out, 
                const VecType* const no_alias re_tin,
                VecType* const no_alias re_tout,
                const VecType* const no_alias im_tin,
                VecType* const no_alias im_tout,
                const MatType* no_alias mat, // the matrix block
                const int rows, 
                const int cols)
{
    // Block size (tuned for Clang/GCC on Intel Skylake processors)
    const int si = 16;
    const int sj = 8;

    for(int bj = 0; bj < cols; bj+=sj) {
        for(int bi = 0; bi < rows; bi+=si) {
            if(bi+si <= rows && bj+sj <= cols)
            {
                // The underlying loops are expected to be unrolled by the compiler
                for(int j = bj; j < bj+sj; ++j) {
                    for(int i = bi; i < bi+si; ++i) {
                        const auto m = mat[j*rows+i]; // Assume a column major ordering
                        re_tout[j]  += m * re_tin[i];
                        im_tout[j]  += m * im_tin[i];
                        re_out[i]   -= m * re_in[j];
                        im_out[i]   -= m * im_in[j];
                    }
                }
            }
            else
            {
                // General case (borders)
                for(int j = bj; j < mini(bj+sj,cols); ++j) {
                    for(int i = bi; i < mini(bi+si,rows); ++i) {
                        const auto m = mat[j*rows+i];
                        re_tout[j]  += m * re_tin[i];
                        im_tout[j]  += m * im_tin[i];
                        re_out[i]   -= m * re_in[j];
                        im_out[i]   -= m * im_in[j];
                    }
                }
            }
        }
    }
}

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

Вот результаты (с G CC 9, использующим типы двойной точности):

With a row=100 and cols=100:
trans_mul_v1: 2.438 us
trans_mul_v2: 2.842 us

With a row=1k and cols=1k:
trans_mul_v1: 452 us
trans_mul_v2: 296 us

With a row=10k and cols=10k:
trans_mul_v1: 71.2 ms
trans_mul_v2: 35.9 ms