C++ - Optimized matrix multiply from Kazushige Goto's paper performs worse than naive in O3 flag

user_3781180 0
Viewed 833

The relevant paper is here. I am trying to reproduce the seminal paper by Kazushige Goto for fast matrix multiplication by decaying it down to subroutines of gepp (general panel-panel) and gebp(general block-panel) multiplication, which apparently is the fastest building blocks for gemm. I wrote the code below to test it and with -O3 flag, I saw the performance of my code is actually worse than the naive matrix multiplication:

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

However, without the -O3 flag we see that the speed is indeed faster than the naive version:

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

With @ztik's suggestion, I tried it without the -mavx2 -O3 flags, and added -O2, which showed similar results to without any optimization flags:

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

I am well aware that -O3 enables tens of hundreds of optimization flags in gcc, each of which will individually increase the performance of the naive method more than the complicated Goto paper variant. However, MKL, ATLAS, etc are some famous BLAS variants which use Goto's method(and are orders of magnitudes faster than naive), albeit with assembly kernels rather than C++ code.

Is this expected, and if so, how can I get the performance improvement(other than having to write self-unrolled assembly)? Did I just write some terrible code, with some kind of bug that doesn't abuse caching properly?

Environment

I am running on a macbook pro with AVX2 instructions, with intel generation:

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

I use g++-8 to compile the below code, and it's currently v8.2.0.

Reproducible code

The reproducible code is below:

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, which contains some utilities for benchmarking:

#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;
}

With diagnosis flags

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
matrix-multiplication blas gcc c++
2 Answers

Hand-wavy answer: Playing around with block size

I tried playing with gebp's block size, and it seems like for 2048x1024 x 1024 x 2048 matrices, a single float seems to be the best for performance under -O2 and no vectorized instructions. This is a bit strange, since I thought bigger blocks could reside inside the tlb/L1/L2 caches.

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

However, here's something very interesting, with vectorized instructions:

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

... and it gets worse from here. It seems like the gcc auto-vectorizer seems to have aggressively optimized for special cases of 4 and 8 for AVX and AVX2 instructions. This is just a speculation, and not a rigorous finding. In extreme cases, the speedup actually goes from 4x to 8x on vectorized instructions!

On a 2048x2048 square matrix multiplication under vectorized environment:

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

This is a whopping >10x speed improvement. Compared to ATLAS, in my benchmarks (using xtensor as the interface), my implementation is only <2x slower.

user_37811800

Just an assumption that doesn't look well in a comment:

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];
    }
}

can be optimized slightly into

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];
    }
}

its likely already done by compiler but who knows without checking assemby output?

Also isn't memcpy good here? It looks like a copy.

user_14700920
Related