C ++ - Оптимизированное умножение матриц из бумаги Казушиге Гото работает хуже, чем наивное в флаге O3 - PullRequest
Новая
       41

C ++ - Оптимизированное умножение матриц из бумаги Казушиге Гото работает хуже, чем наивное в флаге O3

0 голосов
/ 18 сентября 2018

Соответствующий документ здесь . Я пытаюсь воспроизвести основную статью Казушиге Гото для быстрого умножения матрицы , распав его на подпрограммы умножения gepp (общая панель-панель) и gebp (общая панель-панель) , что, по-видимому, является самым быстрым строительные блоки для гемм. Я написал код ниже, чтобы проверить его, и с флагом -O3 я увидел, что производительность моего кода на хуже , чем умножение наивной матрицы: 

(~0.5x increase)
Time elapsed : 3.82941 <-- naive
Time elapsed : 6.21072 <-- more complex gebp subroutine

Однако без флага -O3 мы видим, что скорость действительно выше, чем наивная версия: 

(~4x increase)
Time elapsed : 53.4537 <-- naive
Time elapsed : 15.603 <-- more complex gebp subroutine

По предложению @ ztik я попробовал его без флагов -mavx2 -O3 и добавил -O2, который показал аналогичные результаты без каких-либо флагов оптимизации: 

(~4x increase)
Time elapsed : 26.4217 <-- naive
Time elapsed : 6.42583 <-- more complex gebp subroutine

Мне хорошо известно, что -O3 включает десятки сотен флагов оптимизации в gcc, каждый из которых в отдельности увеличивает производительность простого метода больше, чем сложный вариант документа Goto. Однако MKL, ATLAS и т. Д. Являются некоторыми известными вариантами BLAS, которые используют метод Гото (и на несколько порядков быстрее, чем наивный), хотя и с ядрами сборки, а не с кодом C ++.

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

Окружающая среда

Я работаю на MacBook Pro с инструкциями AVX2, с поддержкой Intel:

Intel(R) Core(TM) i7-6660U CPU @ 2.40GHz

Я использую g++-8 для компиляции приведенного ниже кода, и в настоящее время он v8.2.0.

Воспроизводимый код

Воспроизводимый код ниже:

Makefile: 

CC=g++-8
FLAGS=-std=c++17 -ffast-math -mavx2 -O3

run: main
    ./main

main: main.cpp
    $(CC) -o main main.cpp $(FLAGS)

main.cpp, который содержит некоторые утилиты для бенчмаркинга: 

#include <iostream>
#include <cstdio>
#include <ctime>
#include <cmath>

#define TIME_IT(EXPR) \
{ \
    std::clock_t __start; \
    double __duration; \
    __start = std::clock(); \
    EXPR; \
    __duration = (std::clock() - __start) / (double) CLOCKS_PER_SEC; \
    std::cout << "Time elapsed : " << __duration << std::endl; \
}

static constexpr auto X_SIZE = 2048;
static constexpr auto Y_SIZE = 1024;
static constexpr auto Z_SIZE = 2048;
// = 32 floats
static constexpr auto BLK_BYTES = 128;
static constexpr auto BLK_SIZE = BLK_BYTES / 4; 

template <size_t row, size_t mid, size_t col>
void initialize_matrices(float (&a)[row][mid], float (&b)[mid][col], float (&c)[row][col]){
    // Initialize matrices
    for(auto i = 0; i < row; i++){
        for(auto j = 0; j < mid; j++){
            a[i][j] = ((float) (i*Y_SIZE + j)) / (row * mid);
        }
    }
    for(auto i = 0; i < mid; i++){
        for(auto j = 0; j < col; j++){
            b[i][j] = ((float) (i*Z_SIZE + j)) / (mid * col);
        }
    }
    for(auto i = 0; i < row; i++){
        for(auto j = 0; j < col; j++){
            c[i][j] = 0;
        }
    }
}

template <size_t row, size_t mid, size_t col>
void matmul1(float (&a)[row][mid], float (&b)[mid][col], float (&c)[row][col]){
    for(auto i = 0; i < row; i++){
        for(auto j = 0; j < col; j++){
            float sum = 0;
            for(auto k = 0; k < mid; k++){
                sum += a[i][k] * b[k][j]; 
            }
            c[i][j] = sum;
        }
    } 
}

template <size_t col>
inline void gebp(float (&Ab)[BLK_SIZE][BLK_SIZE], float (&Bp)[BLK_SIZE][col], float (&Cp)[BLK_SIZE][col]){
    // We can optimize this subroutine but that'd be overkill for now.  
    for(auto j = 0; j < col; j++){
        for(auto i = 0; i < BLK_SIZE; i++){
            float sum = 0;
            for(auto k = 0; k < BLK_SIZE; k++){
                sum += Ab[i][k] * Bp[k][j];
            }
            Cp[i][j] += sum;
        }
    }
}

template <size_t row, size_t col>
inline void packb(float (&a)[row][col], float (&b)[BLK_SIZE][BLK_SIZE], size_t m, size_t n){
    // size_t m, n in this case means the m,n-th block to pack.
    auto start_row = m * BLK_SIZE;
    auto start_col = n * BLK_SIZE;
    for(size_t i = 0; i < BLK_SIZE; i++){
        for(size_t j = 0; j < BLK_SIZE; j++){
            b[i][j] = a[start_row + i][start_col + j];
        }
    }
}

template <size_t row, size_t mid, size_t col>
inline void gemm(float (&a)[row][mid], float (&b)[mid][col], float (&c)[row][col]){
    // Divide up the matrix into panels:
    // Suppose row / BLK_SIZE = M
    //         col / BLK_SIZE = N
    //         mid / BLK_SIZE = K
    //
    // For the rest of the function, we assume A, B, C as a, b, c variables,
    // and {var}p = panel of var
    //     {var}b = block of var
    //
    // (TODO: We assume it's perfectly divisible for now)

    // Layout: A = (M,K), B = (K,N), C = (M,N)
    auto M = row / BLK_SIZE;
    auto N = col / BLK_SIZE;
    auto K = mid / BLK_SIZE;
    for(auto p = 0; p < K; p++){
        // Reassign B[p*BLK_SIZE : (p+1)*BLK_SIZE][:] into Bp
        float (&Bp)[BLK_SIZE][col] = *(float (*)[BLK_SIZE][col]) &b[p * BLK_SIZE];

        for(auto i = 0; i < M; i++){
            // Pack A[i*BLK_SIZE : (i+1)*BLK_SIZE][p*BLK_SIZE : (p+1)*BLK_SIZE] into Ab
            float Ab[BLK_SIZE][BLK_SIZE];

            // Reassign C[i*BLK_SIZE : (i+1)*BLK_SIZE][:] into Cp
            float (&Cp)[BLK_SIZE][col] = *(float (*)[BLK_SIZE][col]) &c[i * BLK_SIZE];

            packb(a, Ab, i, p); 
            // The result of Ab and Bp should be in Cp
            gebp(Ab, Bp, Cp);
        } 
    }
}

template <size_t row, size_t col>
bool allclose(float (&a)[row][col], float (&b)[row][col], float threshold = 1e-5, bool verbose = true){
    bool is_equal = true;
    for(auto i = 0; i < row; i++){
        for(auto j = 0; j < col; j++){
            bool current_element = std::abs(a[i][j] - b[i][j]) < threshold;
            if(verbose && !current_element){
                std::cerr << "Element at [" << i << "][" << j << "] is incorrect : " 
                    << a[i][j] << " vs. " << b[i][j] << "." << std::endl;
            }
            is_equal = is_equal && current_element;
        }
    } 
    return is_equal;
}

float a[X_SIZE][Y_SIZE];
float b[Y_SIZE][Z_SIZE];
float c1[X_SIZE][Z_SIZE];
float c2[X_SIZE][Z_SIZE];

int main(){
    initialize_matrices(a, b, c1);
    TIME_IT(matmul1(a, b, c1))

    // We must guarrantee c is all zeros at first.
    initialize_matrices(a, b, c2);
    TIME_IT(gemm(a, b, c2))

    std::cout << allclose(c1, c2, 1e-1, true) << std::endl;
    return 0;
}

С диагностическими флажками 

g++-8 -o main main.cpp -std=c++17 -ffast-math -mavx2 -O3 -ftree-vectorize -fopt-info-vec-missed
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:1083:16: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/char_traits.h:320:25: note: not vectorized: not enough data-refs in basic block.
main.cpp:41:23: note: not vectorized: complicated access pattern.
main.cpp:41:23: note: bad data access.
main.cpp:42:27: note: Unknown misalignment, naturally aligned
main.cpp:36:23: note: not vectorized: complicated access pattern.
main.cpp:36:23: note: bad data access.
main.cpp:37:27: note: Unknown misalignment, naturally aligned
main.cpp:38:34: note: not vectorized: not enough data-refs in basic block.
main.cpp:38:13: note: not vectorized: no vectype for stmt: MEM[(float *)vectp_a.16_42] = vect__6.15_7;
 scalar_type: vector(8) float
main.cpp:38:13: note: not vectorized: not enough data-refs in basic block.
main.cpp:36:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:43:34: note: not vectorized: not enough data-refs in basic block.
main.cpp:43:13: note: not vectorized: no vectype for stmt: MEM[(float *)vectp_b.10_15] = vect__12.9_38;
 scalar_type: vector(8) float
main.cpp:43:13: note: not vectorized: not enough data-refs in basic block.
main.cpp:41:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:51:1: note: not vectorized: not enough data-refs in basic block.
main.cpp:121:23: note: not vectorized: multiple nested loops.
main.cpp:121:23: note: bad loop form.
main.cpp:125:27: note: not vectorized: multiple nested loops.
main.cpp:125:27: note: bad loop form.
main.cpp:69:23: note: not vectorized: multiple nested loops.
main.cpp:69:23: note: bad loop form.
main.cpp:70:27: note: step unknown.
main.cpp:70:27: note: not consecutive access *Cp_14[i_63][j_42] = _30;
main.cpp:70:27: note: not vectorized: complicated access pattern.
main.cpp:70:27: note: bad data access.
main.cpp:72:31: note: step unknown.
main.cpp:72:31: note: misalign = 0 bytes of ref Ab[i_63][k_64]
main.cpp:72:31: note: Unknown alignment for access: *Bp_12[k_64][j_42]
main.cpp:72:31: note: vector alignment may not be reachable
main.cpp:72:31: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:72:31: note: not ssa-name.
main.cpp:72:31: note: use not simple.
main.cpp:72:31: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:72:31: note: not ssa-name.
main.cpp:72:31: note: use not simple.
main.cpp:72:31: note: no array mode for V8SF[2048]
main.cpp:72:31: note: single-element interleaving not supported for not adjacent vector loads
main.cpp:72:31: note: not falling back to elementwise accesses
main.cpp:72:31: note: not vectorized: relevant stmt not supported: _24 = *Bp_12[k_64][j_42];
main.cpp:72:31: note: bad operation or unsupported loop bound.
main.cpp:72:31: note: step unknown.
main.cpp:72:31: note: misalign = 0 bytes of ref Ab[i_63][k_64]
main.cpp:72:31: note: Unknown alignment for access: *Bp_12[k_64][j_42]
main.cpp:72:31: note: vector alignment may not be reachable
main.cpp:72:31: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:72:31: note: not ssa-name.
main.cpp:72:31: note: use not simple.
main.cpp:72:31: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:72:31: note: not ssa-name.
main.cpp:72:31: note: use not simple.
main.cpp:72:31: note: no array mode for V4SF[2048]
main.cpp:72:31: note: single-element interleaving not supported for not adjacent vector loads
main.cpp:72:31: note: not falling back to elementwise accesses
main.cpp:72:31: note: not vectorized: relevant stmt not supported: _24 = *Bp_12[k_64][j_42];
main.cpp:72:31: note: bad operation or unsupported loop bound.
main.cpp:97:25: note: not vectorized: loop contains function calls or data references that cannot be analyzed
main.cpp:96:10: note: not vectorized: not enough data-refs in basic block.
main.cpp:95:10: note: not vectorized: not enough data-refs in basic block.
main.cpp:95:10: note: not vectorized: not enough data-refs in basic block.
main.cpp:97:25: note: not vectorized: not enough data-refs in basic block.
main.cpp:72:31: note: not consecutive access _24 = *Bp_12[k_64][j_42];
main.cpp:72:31: note: not consecutive access _23 = Ab[i_63][k_64];
main.cpp:72:31: note: not vectorized: no grouped stores in basic block.
main.cpp:70:27: note: not consecutive access _29 = *Cp_14[i_63][j_42];
main.cpp:70:27: note: not consecutive access *Cp_14[i_63][j_42] = _30;
main.cpp:70:27: note: not vectorized: no grouped stores in basic block.
main.cpp:69:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:125:27: note: not vectorized: not enough data-refs in basic block.
main.cpp:121:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:137:1: note: not vectorized: not enough data-refs in basic block.
main.cpp:140:6: note: not vectorized: control flow in loop.
main.cpp:140:6: note: bad loop form.
main.cpp:144:49: note: not vectorized: control flow in loop.
main.cpp:144:49: note: bad loop form.
main.cpp:145:13: note: not consecutive access _3 = *a_20(D)[i_17][j_131];
main.cpp:145:13: note: not consecutive access _4 = *b_21(D)[i_17][j_131];
main.cpp:145:13: note: not vectorized: no grouped stores in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:228:43: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:228:43: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/basic_ios.h:49:7: note: not consecutive access _118 = MEM[(const struct basic_ios *)_58]._M_ctype;
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/basic_ios.h:49:7: note: not consecutive access _55 = MEM[(struct basic_ostream *)_36]._vptr.basic_ostream;
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/basic_ios.h:49:7: note: not consecutive access _56 = MEM[(long int *)_55 + -24B];
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/basic_ios.h:49:7: note: not vectorized: no grouped stores in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:874:2: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:875:51: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:877:27: note: not consecutive access _127 = MEM[(const struct ctype *)_118].D.35451._vptr.facet;
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:877:27: note: not consecutive access _128 = MEM[(int (*) () *)_127 + 48B];
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:877:27: note: not vectorized: no grouped stores in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:877:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:143:27: note: not vectorized: not enough data-refs in basic block.
main.cpp:142:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:152:12: note: not vectorized: not enough data-refs in basic block.
main.cpp:55:23: note: not vectorized: multiple nested loops.
main.cpp:55:23: note: bad loop form.
main.cpp:56:27: note: step unknown.
main.cpp:56:27: note: inner step doesn't divide the vector alignment.
main.cpp:56:27: note: Unknown alignment for access: a[i_57][k_59]
main.cpp:56:27: note: misalign = 0 bytes of ref b[k_59][j_58]
main.cpp:56:27: note: misalign = 0 bytes of ref c1[i_57][j_58]
main.cpp:56:27: note: Unknown misalignment, naturally aligned
main.cpp:56:27: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:56:27: note: not ssa-name.
main.cpp:56:27: note: use not simple.
main.cpp:56:27: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:56:27: note: not ssa-name.
main.cpp:56:27: note: use not simple.
main.cpp:56:27: note: can't use a fully-masked loop because no conditional operation is available.
main.cpp:162:5: note: not vectorized: not enough data-refs in basic block.
main.cpp:162:5: note: not vectorized: not enough data-refs in basic block.
main.cpp:58:31: note: not vectorized: no vectype for stmt: vect__37.134_72 = MEM[(float *)vectp_b.132_70];
 scalar_type: vector(8) float
main.cpp:58:31: note: not consecutive access _36 = a[i_57][k_59];
main.cpp:58:31: note: not vectorized: no grouped stores in basic block.
main.cpp:61:13: note: not vectorized: no vectype for stmt: MEM[(float *)vectp_c1.138_80] = vect_sum_40.136_76;
 scalar_type: vector(8) float
main.cpp:61:13: note: not vectorized: not enough data-refs in basic block.
main.cpp:55:23: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:562:44: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:561:18: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:562:44: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:561:18: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:175:29: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:113:13: note: not vectorized: not enough data-refs in basic block.
main.cpp:170:1: note: not vectorized: not enough data-refs in basic block.
./main
Time elapsed : 4.03692
Time elapsed : 6.33257

Ответы [ 2 ]

0 голосов
/ 18 сентября 2018

Просто предположение, которое не очень хорошо выглядит в комментарии: 

auto start_row = m * BLK_SIZE;
auto start_col = n * BLK_SIZE;
for(size_t i = 0; i < BLK_SIZE; i++){
    for(size_t j = 0; j < BLK_SIZE; j++){
        b[i][j] = a[start_row + i][start_col + j];
    }
}

может быть немного оптимизировано в 

const auto start_row = m * BLK_SIZE;
const auto start_col = n * BLK_SIZE;
const auto end_row   = start_row + BLK_SIZE;
const auto end_col   = start_col + BLK_SIZE;
for(size_t i = start_row; i < end_row; i++){
    for(size_t j = start_col; j < end_col; j++){
        b[i][j] = a[i][j];
    }
}

, вероятно, это уже сделано компилятором, но кто знает, не проверяявыход сборки?

И разве memcpy здесь не хорош?Это похоже на копию.

0 голосов
/ 18 сентября 2018

Hand-wavy answer: Игра с размером блока

Я попытался поиграть с размером блока gebp, и похоже, что для матриц 2048x1024 x 1024 x 2048 один float кажется лучшим для производительности при-O2 и без векторизованных инструкций.Это немного странно, так как я думал, что большие блоки могут находиться внутри кешей tlb / L1 / L2. 

1 float:

Time elapsed : 2.33989

8 floats:

Time elapsed : 4.34408

16 floats:

Time elapsed : 4.522

32 floats:

Time elapsed : 6.42583

64 floats:

Time elapsed : 7.35193

128 floats:

Time elapsed : 8.18502

256 floats:

Time elapsed : 8.1686

Однако вот что очень интересное, с векторизованными инструкциями: 

1 float:

Time elapsed : 1.3284

2 floats:

Time elapsed : 1.62332

4 floats:

Time elapsed : 0.421444

8 floats:

Time elapsed : 0.527425

16 floats:

Time elapsed : 4.03877

... и здесь становится все хуже.Похоже, что авто-векторизатор gcc, по-видимому, агрессивно оптимизирован для особых случаев 4 и 8 для инструкций AVX и AVX2. Это всего лишь предположение, а не строгий вывод. В крайних случаях ускорение фактически увеличивается от 4x до 8x для векторизованных инструкций!

Для умножения квадратной матрицы 2048x2048 в векторизованной среде: 

Time elapsed : 11.9983 <-- naive
Time elapsed : 0.822692 <-- vectorized
Time elapsed : 0.472855 <-- ATLAS blas GEMM

Этоколоссальный> 10-кратное улучшение скорости.По сравнению с ATLAS, в моих тестах (с использованием xtensor в качестве интерфейса) моя реализация всего лишь в 2 раза медленнее.

...