Скорость умножения матриц BLAS / Eigen и Numpy - PullRequest
Новая
       21

Скорость умножения матриц BLAS / Eigen и Numpy

0 голосов
/ 25 апреля 2020

Я тестирую матричное умножение для разных библиотек, поскольку я думаю переписать некоторый код на Cython на нативный c ++. Однако простые тесты, по-видимому, подразумевают, что numpy быстрее, чем BLAS или собственные значения для простых умножений матриц.

Я написал следующие файлы: 

#!test_blas.cpp
#include <random>
#include <cstdio>
#include <stdlib.h>
#include <iostream>
#include <cblas.h>

int main ( int argc, char* argv[] ) {

    // Random numbers
    std::mt19937_64 rnd;
    std::uniform_real_distribution<double> doubleDist(0, 1);

    // Create arrays that represent the matrices A,B,C
    const int n = 2000;
    double*  A = new double[n*n];
    double*  B = new double[n*n];
    double*  C = new double[n*n];

    // Fill A and B with random numbers
    for(uint i =0; i <n; i++){
        for(uint j=0; j<n; j++){
            A[i*n+j] = doubleDist(rnd);
            B[i*n+j] = doubleDist(rnd);
        }
    }

    // Calculate A*B=C
    clock_t start = clock();
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1.0, A, n, B, n, 0.0, C, n);
    clock_t end = clock();
    double time = double(end - start)/ CLOCKS_PER_SEC;
    std::cout<< "Time taken : " << time << std::endl; 
    // Clean up
    delete[] A;
    delete[] B;
    delete[] C;

    return 0;
}


#!test_eigen.cpp
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{

    int n_a_rows = 2000;
    int n_a_cols = 2000;
    int n_b_rows = n_a_cols;
    int n_b_cols = 2000;

    MatrixXd a(n_a_rows, n_a_cols);

    for (int i = 0; i < n_a_rows; ++ i)
        for (int j = 0; j < n_a_cols; ++ j)
            a (i, j) = n_a_cols * i + j;

    MatrixXd b (n_b_rows, n_b_cols);
    for (int i = 0; i < n_b_rows; ++ i)
        for (int j = 0; j < n_b_cols; ++ j)
            b (i, j) = n_b_cols * i + j;

    MatrixXd d (n_a_rows, n_b_cols);

    clock_t begin = clock();

    d = a * b;
    clock_t end = clock();
    double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
    std::cout << "Time taken : " << elapsed_secs << std::endl;

}

#!test_numpy.py

import numpy as np
import time
N = 2000

a = np.random.rand(N, N)
b = np.random.rand(N, N)
start = time.time()
c = a.dot(b)
print(f"Time taken : {time.time() - start}")

Наконец, я создал следующий тестовый файл 

#!test_matrix.sh
c++ -O2 -march=native -std=c++11 -I /usr/include/eigen3 test_eigen.cpp -o eigen
c++ -O2 -march=native -std=c++11 test_blas.cpp  -o blas -lcblas

echo "testing BLAS"
./blas
echo "testing Eigen"
./eigen
echo "testing numpy"
python test_numpy.py

, который дает вывод 

testing BLAS
Time taken : 1.63807
testing Eigen
Time taken : 0.795115
testing numpy
Time taken : 0.28397703170776367

Теперь мой вопрос: почему numpy является самым быстрым из этих тестов? Я что-то упустил в отношении оптимизации?

Возможно, numpy использует многопоточность для вычисления матричного произведения. Однако добавление флага компилятора -fopenmp приводит к ухудшению производительности для eigen и BLAS.

Я использую g ++ версии 9.0.3-1. Numpy - это версия 1.18.1, использующая python 3.8.2. Заранее спасибо.

...