c++ - 为什么单纯的C++矩阵乘法比BLAS慢100倍？

Question

我正在研究大型矩阵乘法，并运行以下实验以形成基准测试:

从std normal(0平均值，1 stddev)随机生成两个4096x4096矩阵X，Y。

Z = X * Y

Z的元素总和(以确保它们被访问)并输出。

这是朴素的C++实现:

#include <iostream>
#include <algorithm>

using namespace std;

int main()
{
    constexpr size_t dim = 4096;

    float* x = new float[dim*dim];
    float* y = new float[dim*dim];
    float* z = new float[dim*dim];

    random_device rd;
    mt19937 gen(rd());
    normal_distribution<float> dist(0, 1);

    for (size_t i = 0; i < dim*dim; i++)
    {
        x[i] = dist(gen);
        y[i] = dist(gen);
    }

    for (size_t row = 0; row < dim; row++)
        for (size_t col = 0; col < dim; col++)
        {
            float acc = 0;

            for (size_t k = 0; k < dim; k++)
                acc += x[row*dim + k] * y[k*dim + col];

            z[row*dim + col] = acc;
        }

    float t = 0;

    for (size_t i = 0; i < dim*dim; i++)
        t += z[i];

    cout << t << endl;

    delete x;
    delete y;
    delete z;
}

编译并运行:

$ g++ -std=gnu++11 -O3 test.cpp
$ time ./a.out

这是Octave/matlab实现:

X = stdnormal_rnd(4096, 4096);
Y = stdnormal_rnd(4096, 4096);
Z = X*Y;
sum(sum(Z))

跑:

$ time octave < test.octave

Octave 使用BLAS(我假设 sgemm函数)

硬件是Linux x86-64上的i7 3930X，内存为24 GB。 BLAS似乎正在使用两个核心。也许是超线程对？

我发现在 -O3上用GCC 4.7编译的C++版本执行了9分钟:

real    9m2.126s
user    9m0.302s
sys         0m0.052s

Octave 版本用了6秒:

real    0m5.985s
user    0m10.881s
sys         0m0.144s

我了解到BLAS是经过优化的，并且幼稚的算法完全忽略了缓存等等，但是很严重-90倍？

谁能解释这种差异？ BLAS实现的架构到底是什么？我看到它正在使用Fortran，但是在CPU级别发生了什么？它使用什么算法？如何使用CPU缓存？它调用什么x86-64机器指令？ (它使用的是AVX等高级CPU功能吗？)它从何处获得这种额外的速度？

C++算法的哪些关键优化可以使其与BLAS版本相提并论？

我在gdb下运行了 Octave ，并在计算中途停止了几次。它已经启动了第二个线程，这里是堆栈(所有的停顿看起来都相似):

(gdb) thread 1
#0  0x00007ffff6e17148 in pthread_join () from /lib/x86_64-linux-gnu/libpthread.so.0
#1  0x00007ffff1626721 in ATL_join_tree () from /usr/lib/libblas.so.3
#2  0x00007ffff1626702 in ATL_join_tree () from /usr/lib/libblas.so.3
#3  0x00007ffff15ae357 in ATL_dptgemm () from /usr/lib/libblas.so.3
#4  0x00007ffff1384b59 in atl_f77wrap_dgemm_ () from /usr/lib/libblas.so.3
#5  0x00007ffff193effa in dgemm_ () from /usr/lib/libblas.so.3
#6  0x00007ffff6049727 in xgemm(Matrix const&, Matrix const&, blas_trans_type, blas_trans_type) () from /usr/lib/x86_64-linux-gnu/liboctave.so.1
#7  0x00007ffff6049954 in operator*(Matrix const&, Matrix const&) () from /usr/lib/x86_64-linux-gnu/liboctave.so.1
#8  0x00007ffff7839e4e in ?? () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#9  0x00007ffff765a93a in do_binary_op(octave_value::binary_op, octave_value const&, octave_value const&) () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#10 0x00007ffff76c4190 in tree_binary_expression::rvalue1(int) () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#11 0x00007ffff76c33a5 in tree_simple_assignment::rvalue1(int) () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#12 0x00007ffff76d0864 in tree_evaluator::visit_statement(tree_statement&) () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#13 0x00007ffff76cffae in tree_evaluator::visit_statement_list(tree_statement_list&) () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#14 0x00007ffff757f6d4 in main_loop() () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#15 0x00007ffff7527abf in octave_main () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1

(gdb) thread 2
#0  0x00007ffff14ba4df in ATL_dJIK56x56x56TN56x56x0_a1_b1 () from /usr/lib/libblas.so.3
(gdb) bt
#0  0x00007ffff14ba4df in ATL_dJIK56x56x56TN56x56x0_a1_b1 () from /usr/lib/libblas.so.3
#1  0x00007ffff15a5fd7 in ATL_dmmIJK2 () from /usr/lib/libblas.so.3
#2  0x00007ffff15a6ae4 in ATL_dmmIJK () from /usr/lib/libblas.so.3
#3  0x00007ffff1518e65 in ATL_dgemm () from /usr/lib/libblas.so.3
#4  0x00007ffff15adf7a in ATL_dptgemm0 () from /usr/lib/libblas.so.3
#5  0x00007ffff6e15e9a in start_thread () from /lib/x86_64-linux-gnu/libpthread.so.0
#6  0x00007ffff6b41cbd in clone () from /lib/x86_64-linux-gnu/libc.so.6
#7  0x0000000000000000 in ?? ()

它正在按预期方式调用BLAS gemm。

第一个线程似乎正在加入第二个线程，因此我不确定这两个线程是否占观察到的200％CPU使用率。

ATL_dgemm libblas.so.3是哪个库，其代码在哪里？

$ ls -al /usr/lib/libblas.so.3
/usr/lib/libblas.so.3 -> /etc/alternatives/libblas.so.3

$ ls -al /etc/alternatives/libblas.so.3
/etc/alternatives/libblas.so.3 -> /usr/lib/atlas-base/atlas/libblas.so.3

$ ls -al /usr/lib/atlas-base/atlas/libblas.so.3
/usr/lib/atlas-base/atlas/libblas.so.3 -> libblas.so.3.0

$ ls -al /usr/lib/atlas-base/atlas/libblas.so.3.0
/usr/lib/atlas-base/atlas/libblas.so.3.0

$ dpkg -S /usr/lib/atlas-base/atlas/libblas.so.3.0
libatlas3-base: /usr/lib/atlas-base/atlas/libblas.so.3.0

$ apt-get source libatlas3-base

是ATLAS 3.8.4

以下是我稍后实现的优化:

使用平铺方法，将X，Y和Z的64x64块预加载到单独的数组中。

更改每个块的计算，以使内部循环如下所示:

for (size_t tcol = 0; tcol < block_width; tcol++)
    bufz[trow][tcol] += B * bufy[tk][tcol];

这允许GCC自动矢量化为SIMD指令，并且还允许指令级并行性(我认为)。

打开 march=corei7-avx。这样可以获得30％的额外速度，但由于我认为BLAS库是预先构建的，因此很容易出错。

这是代码:

#include <iostream>
#include <algorithm>

using namespace std;

constexpr size_t dim = 4096;
constexpr size_t block_width = 64;
constexpr size_t num_blocks = dim / block_width;

double X[dim][dim], Y[dim][dim], Z[dim][dim];

double bufx[block_width][block_width];
double bufy[block_width][block_width];
double bufz[block_width][block_width];

void calc_block()
{
    for (size_t trow = 0; trow < block_width; trow++)
        for (size_t tk = 0; tk < block_width; tk++)
        {
            double B = bufx[trow][tk];

            for (size_t tcol = 0; tcol < block_width; tcol++)
                bufz[trow][tcol] += B * bufy[tk][tcol];
        }
}

int main()
{
    random_device rd;
    mt19937 gen(rd());
    normal_distribution<double> dist(0, 1);

    for (size_t row = 0; row < dim; row++)
        for (size_t col = 0; col < dim; col++)
        {
            X[row][col] = dist(gen);
            Y[row][col] = dist(gen);
            Z[row][col] = 0;
        }

    for (size_t block_row = 0; block_row < num_blocks; block_row++)
        for (size_t block_col = 0; block_col < num_blocks; block_col++)
        {
            for (size_t trow = 0; trow < block_width; trow++)
                for (size_t tcol = 0; tcol < block_width; tcol++)
                    bufz[trow][tcol] = 0;

            for (size_t block_k = 0; block_k < num_blocks; block_k++)
            {
                for (size_t trow = 0; trow < block_width; trow++)
                    for (size_t tcol = 0; tcol < block_width; tcol++)
                    {
                        bufx[trow][tcol] = X[block_row*block_width + trow][block_k*block_width + tcol];
                        bufy[trow][tcol] = Y[block_k*block_width + trow][block_col*block_width + tcol];
                    }

                calc_block();
            }

            for (size_t trow = 0; trow < block_width; trow++)
                for (size_t tcol = 0; tcol < block_width; tcol++)
                    Z[block_row*block_width + trow][block_col*block_width + tcol] = bufz[trow][tcol];

        }

    double t = 0;

    for (size_t row = 0; row < dim; row++)
        for (size_t col = 0; col < dim; col++)
            t += Z[row][col];

    cout << t << endl;
}

所有操作都在calc_block函数中-花费了90％以上的时间。

新时间是:

real    0m17.370s
user    0m17.213s
sys 0m0.092s

更接近。

calc_block函数的反编译如下:

0000000000401460 <_Z10calc_blockv>:
  401460:   b8 e0 21 60 00          mov    $0x6021e0,%eax
  401465:   41 b8 e0 23 61 00       mov    $0x6123e0,%r8d
  40146b:   31 ff                   xor    %edi,%edi
  40146d:   49 29 c0                sub    %rax,%r8
  401470:   49 8d 34 00             lea    (%r8,%rax,1),%rsi
  401474:   48 89 f9                mov    %rdi,%rcx
  401477:   ba e0 a1 60 00          mov    $0x60a1e0,%edx
  40147c:   48 c1 e1 09             shl    $0x9,%rcx
  401480:   48 81 c1 e0 21 61 00    add    $0x6121e0,%rcx
  401487:   66 0f 1f 84 00 00 00    nopw   0x0(%rax,%rax,1)
  40148e:   00 00 
  401490:   c4 e2 7d 19 01          vbroadcastsd (%rcx),%ymm0
  401495:   48 83 c1 08             add    $0x8,%rcx
  401499:   c5 fd 59 0a             vmulpd (%rdx),%ymm0,%ymm1
  40149d:   c5 f5 58 08             vaddpd (%rax),%ymm1,%ymm1
  4014a1:   c5 fd 29 08             vmovapd %ymm1,(%rax)
  4014a5:   c5 fd 59 4a 20          vmulpd 0x20(%rdx),%ymm0,%ymm1
  4014aa:   c5 f5 58 48 20          vaddpd 0x20(%rax),%ymm1,%ymm1
  4014af:   c5 fd 29 48 20          vmovapd %ymm1,0x20(%rax)
  4014b4:   c5 fd 59 4a 40          vmulpd 0x40(%rdx),%ymm0,%ymm1
  4014b9:   c5 f5 58 48 40          vaddpd 0x40(%rax),%ymm1,%ymm1
  4014be:   c5 fd 29 48 40          vmovapd %ymm1,0x40(%rax)
  4014c3:   c5 fd 59 4a 60          vmulpd 0x60(%rdx),%ymm0,%ymm1
  4014c8:   c5 f5 58 48 60          vaddpd 0x60(%rax),%ymm1,%ymm1
  4014cd:   c5 fd 29 48 60          vmovapd %ymm1,0x60(%rax)
  4014d2:   c5 fd 59 8a 80 00 00    vmulpd 0x80(%rdx),%ymm0,%ymm1
  4014d9:   00 
  4014da:   c5 f5 58 88 80 00 00    vaddpd 0x80(%rax),%ymm1,%ymm1
  4014e1:   00 
  4014e2:   c5 fd 29 88 80 00 00    vmovapd %ymm1,0x80(%rax)
  4014e9:   00 
  4014ea:   c5 fd 59 8a a0 00 00    vmulpd 0xa0(%rdx),%ymm0,%ymm1
  4014f1:   00 
  4014f2:   c5 f5 58 88 a0 00 00    vaddpd 0xa0(%rax),%ymm1,%ymm1
  4014f9:   00 
  4014fa:   c5 fd 29 88 a0 00 00    vmovapd %ymm1,0xa0(%rax)
  401501:   00 
  401502:   c5 fd 59 8a c0 00 00    vmulpd 0xc0(%rdx),%ymm0,%ymm1
  401509:   00 
  40150a:   c5 f5 58 88 c0 00 00    vaddpd 0xc0(%rax),%ymm1,%ymm1
  401511:   00 
  401512:   c5 fd 29 88 c0 00 00    vmovapd %ymm1,0xc0(%rax)
  401519:   00 
  40151a:   c5 fd 59 8a e0 00 00    vmulpd 0xe0(%rdx),%ymm0,%ymm1
  401521:   00 
  401522:   c5 f5 58 88 e0 00 00    vaddpd 0xe0(%rax),%ymm1,%ymm1
  401529:   00 
  40152a:   c5 fd 29 88 e0 00 00    vmovapd %ymm1,0xe0(%rax)
  401531:   00 
  401532:   c5 fd 59 8a 00 01 00    vmulpd 0x100(%rdx),%ymm0,%ymm1
  401539:   00 
  40153a:   c5 f5 58 88 00 01 00    vaddpd 0x100(%rax),%ymm1,%ymm1
  401541:   00 
  401542:   c5 fd 29 88 00 01 00    vmovapd %ymm1,0x100(%rax)
  401549:   00 
  40154a:   c5 fd 59 8a 20 01 00    vmulpd 0x120(%rdx),%ymm0,%ymm1
  401551:   00 
  401552:   c5 f5 58 88 20 01 00    vaddpd 0x120(%rax),%ymm1,%ymm1
  401559:   00 
  40155a:   c5 fd 29 88 20 01 00    vmovapd %ymm1,0x120(%rax)
  401561:   00 
  401562:   c5 fd 59 8a 40 01 00    vmulpd 0x140(%rdx),%ymm0,%ymm1
  401569:   00 
  40156a:   c5 f5 58 88 40 01 00    vaddpd 0x140(%rax),%ymm1,%ymm1
  401571:   00 
  401572:   c5 fd 29 88 40 01 00    vmovapd %ymm1,0x140(%rax)
  401579:   00 
  40157a:   c5 fd 59 8a 60 01 00    vmulpd 0x160(%rdx),%ymm0,%ymm1
  401581:   00 
  401582:   c5 f5 58 88 60 01 00    vaddpd 0x160(%rax),%ymm1,%ymm1
  401589:   00 
  40158a:   c5 fd 29 88 60 01 00    vmovapd %ymm1,0x160(%rax)
  401591:   00 
  401592:   c5 fd 59 8a 80 01 00    vmulpd 0x180(%rdx),%ymm0,%ymm1
  401599:   00 
  40159a:   c5 f5 58 88 80 01 00    vaddpd 0x180(%rax),%ymm1,%ymm1
  4015a1:   00 
  4015a2:   c5 fd 29 88 80 01 00    vmovapd %ymm1,0x180(%rax)
  4015a9:   00 
  4015aa:   c5 fd 59 8a a0 01 00    vmulpd 0x1a0(%rdx),%ymm0,%ymm1
  4015b1:   00 
  4015b2:   c5 f5 58 88 a0 01 00    vaddpd 0x1a0(%rax),%ymm1,%ymm1
  4015b9:   00 
  4015ba:   c5 fd 29 88 a0 01 00    vmovapd %ymm1,0x1a0(%rax)
  4015c1:   00 
  4015c2:   c5 fd 59 8a c0 01 00    vmulpd 0x1c0(%rdx),%ymm0,%ymm1
  4015c9:   00 
  4015ca:   c5 f5 58 88 c0 01 00    vaddpd 0x1c0(%rax),%ymm1,%ymm1
  4015d1:   00 
  4015d2:   c5 fd 29 88 c0 01 00    vmovapd %ymm1,0x1c0(%rax)
  4015d9:   00 
  4015da:   c5 fd 59 82 e0 01 00    vmulpd 0x1e0(%rdx),%ymm0,%ymm0
  4015e1:   00 
  4015e2:   c5 fd 58 80 e0 01 00    vaddpd 0x1e0(%rax),%ymm0,%ymm0
  4015e9:   00 
  4015ea:   48 81 c2 00 02 00 00    add    $0x200,%rdx
  4015f1:   48 39 ce                cmp    %rcx,%rsi
  4015f4:   c5 fd 29 80 e0 01 00    vmovapd %ymm0,0x1e0(%rax)
  4015fb:   00 
  4015fc:   0f 85 8e fe ff ff       jne    401490 <_Z10calc_blockv+0x30>
  401602:   48 83 c7 01             add    $0x1,%rdi
  401606:   48 05 00 02 00 00       add    $0x200,%rax
  40160c:   48 83 ff 40             cmp    $0x40,%rdi
  401610:   0f 85 5a fe ff ff       jne    401470 <_Z10calc_blockv+0x10>
  401616:   c5 f8 77                vzeroupper 
  401619:   c3                      retq   
  40161a:   66 0f 1f 44 00 00       nopw   0x0(%rax,%rax,1)

Answer 1

造成代码与BLAS之间性能差异的三个因素(加上有关Strassen算法的注释)。

在您的内部循环中，在k上，您具有y[k*dim + col]。由于内存缓存的安排方式，具有相同k和dim的col的连续值映射到相同的缓存集。缓存的结构方式，每个内存地址都有一个缓存集，在缓存中必须保留其内容。每个高速缓存集都有几行(通常是四行)，这些行中的每一行都可以容纳映射到该特定高速缓存集的任何内存地址。

因为您的内部循环以这种方式遍历y，所以每次使用y中的元素时，它都必须将该元素的内存加载到与上一次迭代相同的集合中。这将强制清除集合中先前的缓存行之一。然后，在col循环的下一个迭代中，y的所有元素都已从缓存中清除，因此必须重新加载它们。

因此，循环每次加载y元素时就需要，它必须从内存中加载，这需要很多CPU周期。

高性能代码通过两种方式避免了这种情况。第一，它将工作分成较小的块。将行和列划分为较小的大小，并使用较短的循环进行处理，这些循环可以使用高速缓存行中的所有元素，并可以在继续进行下一个块之前多次使用每个元素。因此，对x元素和y元素的大多数引用通常来自缓存，通常是在单个处理器周期内。二，在某些情况下，代码会将数据从矩阵的一列(由于几何结构而影响缓存)复制到临时缓冲区的行中(这避免了影响)。这再次允许大多数内存引用从缓存而不是从内存提供。

另一个因素是使用单指令多数据(SIMD)功能。许多现代处理器的指令都在一条指令中加载多个元素(典型的是四个float元素，但现在有八个)，存储多个元素，添加多个元素(例如，对于这四个元素中的每个元素，将其添加到那些元素中的相应一个元素上)四)，乘以多个元素，依此类推。只要您能够安排使用这些指令的工作，只需使用这些指令即可立即使代码快四倍。

这些指令无法在标准C中直接访问。某些优化器现在尝试在可能的情况下使用这些指令，但是这种优化很困难，并且从中获得很多好处并不常见。许多编译器提供了对语言的扩展，这些语言可以访问这些指令。就个人而言，我通常更喜欢使用汇编语言编写以使用SIMD。

另一个因素是在处理器上使用指令级并行执行功能。请注意，在您的内部循环中，acc已更新。下一个迭代不能添加到acc中，直到上一个迭代完成了acc的更新为止。高性能代码将使多个和并行运行(甚至多个SIMD和)。这样的结果是，在执行一个和的加法运算时，将开始另一个和的加法运算。在当今的处理器上，通常一次支持四个或更多正在进行的浮点运算。如所写，您的代码根本无法做到这一点。一些编译器将尝试通过重新排列循环来优化代码，但这要求编译器能够看到特定循环的迭代是彼此独立的，或者可以与另一个循环互换，等等。

有效地使用高速缓存可以使性能提高十倍，SIMD可以提供另外四项，而指令级并行性可以提供另外四项，总共达到160，这是完全可行的。

这是根据this Wikipedia page对Strassen算法效果的非常粗略的估算。 Wikipedia页面上说Strassen优于n = 100周围的直接乘法。这表明执行时间的常数因子之比为1003/1002.807≈2.4。显然，这将根据处理器模型，与缓存效果相互作用的矩阵大小等而有很大差异。但是，简单的外推法显示，在n = 4096((4096/100)3-2.807≈2.05)时，斯特拉森的倍数约为直接乘法的两倍。再次，这只是一个大概的估计。

至于以后的优化，请在内部循环中考虑以下代码:

bufz[trow][tcol] += B * bufy[tk][tcol];

一个潜在的问题是bufz通常可以与bufy重叠。由于您对bufz和bufy使用了全局定义，因此编译器可能知道在这种情况下它们不重叠。但是，如果将此代码移到作为参数传递bufz和bufy的子例程中，尤其是如果在单独的源文件中编译该子例程，则编译器不太可能知道bufz和bufy不重叠。在这种情况下，编译器无法对代码进行矢量化或以其他方式重新排序，因为此迭代中的bufz[trow][tcol]可能与另一迭代中的bufy[tk][tcol]相同。

即使编译器可以看到在当前源模块中使用不同的bufz和bufy调用了子例程，如果该例程具有extern链接(默认)，那么编译器也必须允许从外部模块调用该例程，因此，如果bufz和bufy重叠，它必须生成可以正常工作的代码。 (编译器处理此问题的一种方法是生成例程的两个版本，一个版本从外部模块调用，一个版本从当前正在编译的模块调用。是否这样做取决于您的编译器，优化会切换，等)。如果将该例程声明为static，则编译器知道无法从外部模块调用该例程(除非您获取其地址，并且该地址有可能传递到当前模块之外)。

另一个潜在的问题是，即使编译器将这些代码矢量化，它也不一定会为您执行的处理器生成最佳代码。查看生成的汇编代码，看来编译器在重复使用%ymm1。一次又一次地，它将内存中的值乘以%ymm1，将内存中的值乘以%ymm1，并将值从%ymm1存入内存。这有两个问题。

第一，您不希望将这些部分总和频繁存储到内存中。您希望将许多累加累积到一个寄存器中，并且该寄存器只会被很少地写入存储器。要说服编译器执行此操作，可能需要重写代码，以明确地将部分和保留在临时对象中，并在循环完成后将其写入内存。

第二，这些指令名义上是串行相关的。在乘法完成之前，添加操作无法开始，并且在添加完成之前，存储区无法写入内存。 Core i7具有出色的无序执行功能。因此，尽管它具有等待开始执行的加法，但稍后会在指令流中查看乘法并启动它。 (即使该乘法也使用%ymm1，处理器也会动态地重新映射寄存器，以便它使用不同的内部寄存器来执行此乘法。)即使您的代码充满了连续的依赖性，处理器也会尝试在以下位置执行几条指令一旦。但是，许多因素都会干扰这一点。您可以用尽处理器用于重命名的内部寄存器。您使用的内存地址可能会出现错误冲突。 (处理器查看大约十二个内存地址的低位，以查看该地址是否与它试图从较早的指令加载或存储的另一个地址相同。如果位相等，则处理器具有延迟当前的加载或存储，直到它可以验证整个地址是否不同为止。这种延迟可能比当前的加载或存储延迟更多。)因此，最好使指令相互独立。

这是我更喜欢在汇编中编写高性能代码的另一个原因。要在C语言中执行此操作，您必须说服编译器通过执行一些操作(例如编写一些自己的SIMD代码(使用它们的语言扩展名)和手动展开循环(编写多个迭代))来向您提供此类指令。

复制进出缓冲区时，可能会有类似的问题。但是，您报告90％的时间都花在calc_block中，因此我没有仔细研究。

另外，斯特拉森的算法是否解决了其余的差异？