为啥 blas 比 numpy 慢

Posted 2023-02-23

【中文标题】为啥 blas 比 numpy 慢

【英文标题】：why blas is slower than numpy为什么 blas 比 numpy 慢

【发布时间】：2018-01-25 05:46:09

【问题描述】：

感谢 Mats Petersson 的帮助。他的 C++ 的运行时间终于看起来正常了！但我有两个新问题。

为什么 Mats Petersson 的代码比我的代码快两倍？

Mats Petersson 的 C++ 代码是：

#include <iostream>
#include <openblas/cblas.h>
#include <array>
#include <iterator>
#include <random>
#include <ctime>
using namespace std;
const blasint m = 100, k = 100, n = 100;
// Mats Petersson's declaration
array<array<double, k>, m> AA[500]; 
array<array<double, n>, k> BB[500]; 
array<array<double, n>, m> CC[500]; 
// My declaration
array<array<double, k>, m> AA1; 
array<array<double, n>, k> BB1; 
array<array<double, n>, m> CC1; 

int main(void) 
    CBLAS_ORDER Order = CblasRowMajor;
    CBLAS_TRANSPOSE TransA = CblasNoTrans, TransB = CblasNoTrans;

    const float alpha = 1;
    const float beta = 0;
    const int lda = k;
    const int ldb = n;
    const int ldc = n;
    default_random_engine r_engine(time(0));
    uniform_real_distribution<double> uniform(0, 1);

    double dur = 0;
    clock_t start,end;
    double total = 0;
    // Mats Petersson's initialization and computation
    for(int i = 0; i < 500; i++) 
        for (array<array<double, k>, m>::iterator iter = AA[i].begin(); iter != AA[i].end(); ++iter) 
            for (double &number : (*iter))
                number = uniform(r_engine);
        
        for (array<array<double, n>, k>::iterator iter = BB[i].begin(); iter != BB[i].end(); ++iter) 
            for (double &number : (*iter))
                number = uniform(r_engine);
        
    
    start = clock();
    for(int i = 0; i < 500; ++i)
        cblas_dgemm(Order, TransA, TransB, m, n, k, alpha, &AA[i][0][0], lda, &BB[i][0][0], ldb, beta, &CC[i][0][0], ldc);
    
    end = clock();
    dur += (double)(end - start);
    cout<<endl<<"Mats Petersson spends "<<(dur/CLOCKS_PER_SEC)<<" seconds to compute it"<<endl<<endl;

    // It turns me!  
    dur = 0;
    for(int i = 0; i < 500; i++)
        for(array<array<double, k>, m>::iterator iter = AA1.begin(); iter != AA1.end(); ++iter)
            for(double& number : (*iter))
                number = uniform(r_engine);
        
        for(array<array<double, n>, k>::iterator iter = BB1.begin(); iter != BB1.end(); ++iter)
            for(double& number : (*iter))
                number = uniform(r_engine);
        
        start = clock();
        cblas_dgemm(Order, TransA, TransB, m, n, k, alpha, &AA1[0][0], lda, &BB1[0][0], ldb, beta, &CC1[0][0], ldc);
        end = clock();
        dur += (double)(end - start);
    

    cout<<endl<<"I spend "<<(dur/CLOCKS_PER_SEC)<<" seconds to compute it"<<endl<<endl;  

结果如下：

Mats Petersson spends 0.215056 seconds to compute it

I spend 0.459066 seconds to compute it

那么，为什么他的代码比我的代码快两倍？

Python 还是更快？

numpy 代码是

import numpy as np
import time
a = 
b = 
c = 
for i in range(500):
    a[i] = np.matrix(np.random.rand(100, 100))
    b[i] = np.matrix(np.random.rand(100, 100))
    c[i] = np.matrix(np.random.rand(100, 100))
start = time.time()
for i in range(500):
    c[i] = a[i]*b[i]
print(time.time() - start)

结果是：

还是看不懂！

【问题讨论】：

也许 numpy 生成随机数更快？如果你想比较乘法的速度，也许你应该去掉随机数生成部分。

试试g++ -O2 -march=native -std=c++11 main.cpp -lopenblas?

你确定 Python 代码中的 a*b 实际上是计算它而不是懒惰吗？ （在 Python3 中很懒惰 range）/ C++ 代码末尾缺少 。 / 您在 C++ 中的 fill 命令正在用一个值填充每一行，这与 Python 行为不同。不过应该不是问题。

@user202729 numpy.ndarray 操作几乎都是急切的。

AA、BB 和 CC 是否连续？ numpy.ndarray 对象是。

【参考方案1】：

因此，我无法使用以下代码重现原始结果：

#include <iostream>
#include <openblas/cblas.h>
#include <array>
#include <iterator>
#include <random>
#include <ctime>
using namespace std;

const blasint m = 100, k = 100, n = 100;
array<array<double, k>, m> AA[500];
array<array<double, n>, k> BB[500];
array<array<double, n>, m> CC[500];

int main(void) 
    CBLAS_ORDER Order = CblasRowMajor;
    CBLAS_TRANSPOSE TransA = CblasNoTrans, TransB = CblasNoTrans;


    const float alpha = 1;
    const float beta = 0;
    const int lda = k; 
    const int ldb = n; 
    const int ldc = n; 
    default_random_engine r_engine(time(0));
    uniform_real_distribution<double> uniform(0, 1);

    double dur = 0;
    clock_t start,end;
    double total = 0;

    for(int i = 0; i < 500; i++)
        for(array<array<double, k>, m>::iterator iter = AA[i].begin(); iter != AA[i].end(); ++iter)
            for(double& number : (*iter))
                number = uniform(r_engine);
        
        for(array<array<double, n>, k>::iterator iter = BB[i].begin(); iter != BB[i].end(); ++iter)
            for(double& number : (*iter))
                number = uniform(r_engine);
        
    

    start = clock();
    for(int i = 0; i < 500; i++)
    
        cblas_dgemm(Order, TransA, TransB, m, n, k, alpha, &AA[i][0][0], lda, &BB[i][0][0], ldb, beta, 
            &CC[i][0][0], ldc);
    total += CC[i][i/5][i/5];
    
    end = clock();
    dur = (double)(end - start);

    cout<<endl<<"It spends "<<(dur/CLOCKS_PER_SEC)<<" seconds to compute it"<<endl<<endl;
    cout << "total =" << total << endl;

还有这段代码：

import numpy as np
import time
a = 
b = 
c = 
for i in range(500):
    a[i] = np.matrix(np.random.rand(100, 100))
    b[i] = np.matrix(np.random.rand(100, 100))
    c[i] = np.matrix(np.random.rand(100, 100))
start = time.time()
for i in range(500):
    c[i] = a[i]*b[i]
print(time.time() - start)

我们知道循环（几乎）做同样的事情。我的结果如下：

python 2.7：0.676353931427 python 3.4：0.6782681941986084 clang++ -O2: 0.117377 g++ -O2: 0.117685

使数组全局化可确保我们不会炸毁堆栈。我还将 rengine1 更改为 rengine，因为它无法按原样编译。

然后我确保两个示例都计算了 500 个不同的数组值。

有趣的是，g++ 的总时间比 clang++ 的总时间短得多——但这是时间测量之外的循环，实际的矩阵乘法是相同的，给或取千分之一秒。 python 的总执行时间介于 clang 和 g++ 之间。

【讨论】：

我很想知道为什么 clang 的总执行时间比 GCC 大。 clang++ 是否使用了 libc++ 和 g++ stdlibc++？

我认为clang只是没有像gcc那样优化初始化代码，但我没有研究它。太容易被卷入那些事情中，而不能及时工作...... :)