c++ - 性能:Matlab vs C++ 矩阵 vector 乘法

标签 c++ matlab performance matrix eigen

序言

前段时间,我问了一个关于 Matlab 与 Python 性能的问题(Performance: Matlab vs Python)。我很惊讶 Matlab 比 Python 快,尤其是在 meshgrid 中。在讨论该问题时,有人指出我应该在 Python 中使用包装器来调用我的 C++ 代码,因为我也可以使用 C++ 代码。我在 C++、Matlab 和 Python 中有相同的代码。

在这样做的同时,我再次惊讶地发现,Matlab 在矩阵汇编计算方面比 C++ 更快。我有一个稍大的代码,我正在研究一段矩阵- vector 乘法。较大的代码在多个实例中执行这样的乘法。总体而言,C++ 中的代码比 Matlab 快得多(因为在 Matlab 中调用函数有开销等),但在矩阵 vector 乘法方面,Matlab 似乎优于 C++(代码片段在底部)。

结果

下表显示了组装内核矩阵所需的时间以及将矩阵与 vector 相乘所需的时间的比较。结果针对矩阵大小 NxN 编译,其中 N 从 10,000 到 40,000 不等。这不是那么大。但有趣的是,N 越大,Matlab 的性能就优于 C++。 Matlab 总时间快 3.8 - 5.8 倍。此外,它在矩阵组装计算方面也更快。

 ___________________________________________
|N=10,000   Assembly    Computation  Total  |
|MATLAB     0.3387      0.031        0.3697 |
|C++        1.15        0.24         1.4    |
|Times faster                        3.8    |
 ___________________________________________ 
|N=20,000   Assembly    Computation  Total  |
|MATLAB     1.089       0.0977       1.187  |
|C++        5.1         1.03         6.13   |
|Times faster                        5.2    |
 ___________________________________________
|N=40,000   Assembly    Computation  Total  |
|MATLAB     4.31        0.348        4.655  |
|C++        23.25       3.91         27.16  |
|Times faster                        5.8    |
 -------------------------------------------

问题

在 C++ 中有更快的方法吗?我错过了什么吗?我知道 C++ 正在使用 for 循环,但我的理解是 Matlab 也会在 meshgrid 中做类似的事情。

代码片段

Matlab 代码:

%% GET INPUT DATA FROM DATA FILES ------------------------------------------- %
% Read data from input file
Data       = load('Input/input.txt');
location   = Data(:,1:2);           
charges    = Data(:,3:end);         
N          = length(location);      
m          = size(charges,2);       

%% EXACT MATRIX VECTOR PRODUCT ---------------------------------------------- %
kex1=ex1; 
tic
Q = kex1.kernel_2D(location , location);
fprintf('\n Assembly time: %f ', toc);

tic
potential_exact = Q * charges;
fprintf('\n Computation time: %f \n', toc);

类(使用网格):

classdef ex1
    methods 
        function [kernel] = kernel_2D(obj, x,y) 
            [i1,j1] = meshgrid(y(:,1),x(:,1));
            [i2,j2] = meshgrid(y(:,2),x(:,2));
            kernel = sqrt( (i1 - j1) .^ 2 + (i2 - j2) .^2 );
        end
    end       
end

C++ 代码:

编辑

使用带有以下标志的 make 文件编译:

CC=g++ 
CFLAGS=-c -fopenmp -w -Wall -DNDEBUG -O3 -march=native -ffast-math -ffinite-math-only -I header/ -I /usr/include 
LDFLAGS= -g -fopenmp  
LIB_PATH= 

SOURCESTEXT= src/read_Location_Charges.cpp 
SOURCESF=examples/matvec.cpp
OBJECTSF= $(SOURCESF:.cpp=.o) $(SOURCESTEXT:.cpp=.o)
EXECUTABLEF=./exec/mykernel
mykernel: $(SOURCESF) $(SOURCESTEXT) $(EXECUTABLEF)
$(EXECUTABLEF): $(OBJECTSF)
    $(CC) $(LDFLAGS) $(KERNEL) $(INDEX) $(OBJECTSF) -o $@ $(LIB_PATH)
.cpp.o:
    $(CC) $(CFLAGS) $(KERNEL) $(INDEX) $< -o $@

`

# include"environment.hpp"
using namespace std;
using namespace Eigen;

class ex1 
{
public:
    void kernel_2D(const unsigned long M, double*& x, const unsigned long N,  double*&  y, MatrixXd& kernel)    {   
        kernel  =   MatrixXd::Zero(M,N);
        for(unsigned long i=0;i<M;++i)  {
            for(unsigned long j=0;j<N;++j)  {
                        double X =   (x[0*N+i] - y[0*N+j]) ;
                        double Y =   (x[1*N+i] - y[1*N+j]) ;
                        kernel(i,j) = sqrt((X*X) + (Y*Y));
            }
        }
    }
};

int main()
{
    /* Input ----------------------------------------------------------------------------- */
    unsigned long N = 40000;          unsigned m=1;                   
    double* charges;                  double* location;
    charges =   new double[N * m]();  location =    new double[N * 2]();
    clock_t start;                    clock_t end;
    double exactAssemblyTime;         double exactComputationTime;

    read_Location_Charges ("input/test_input.txt", N, location, m, charges);

    MatrixXd charges_           =   Map<MatrixXd>(charges, N, m);
    MatrixXd Q;
    ex1 Kex1;

    /* Process ------------------------------------------------------------------------ */
    // Matrix assembly
    start = clock();
        Kex1.kernel_2D(N, location, N, location, Q);
    end = clock();
    exactAssemblyTime = double(end-start)/double(CLOCKS_PER_SEC);

    //Computation
    start = clock();
        MatrixXd QH = Q * charges_;
    end = clock();
    exactComputationTime = double(end-start)/double(CLOCKS_PER_SEC);

    cout << endl << "Assembly     time: " << exactAssemblyTime << endl;
    cout << endl << "Computation time: " << exactComputationTime << endl;


    // Clean up
    delete []charges;
    delete []location;
    return 0;
}

最佳答案

正如评论中所说,MatLab 依赖英特尔的 MKL 库来处理矩阵产品,这是此类操作最快的库。尽管如此,仅 Eigen 就应该能够提供类似的性能。为此,请确保使用最新的 Eigen(例如 3.4)和正确的编译标志以启用 AVX/FMA(如果可用)和多线程:

-O3 -DNDEBUG -march=native

由于 charges_ 是一个 vector ,最好使用 VectorXd 来让 Eigen 知道您想要矩阵 vector 乘积而不是矩阵矩阵乘积。

如果你有 Intel 的 MKL,那么你也可以让 Eigen uses it在这种精确操作中获得与 MatLab 完全相同的性能。

关于汇编,最好将两个循环反转以启用矢量化,然后使用 OpenMP 启用多线程(添加 -fopenmp 作为编译器标志)以使最外层循环并行运行,最后可以简化您使用 Eigen 的代码:

void kernel_2D(const unsigned long M, double* x, const unsigned long N,  double*  y, MatrixXd& kernel)    {
    kernel.resize(M,N);
    auto x0 = ArrayXd::Map(x,M);
    auto x1 = ArrayXd::Map(x+M,M);
    auto y0 = ArrayXd::Map(y,N);
    auto y1 = ArrayXd::Map(y+N,N);
    #pragma omp parallel for
    for(unsigned long j=0;j<N;++j)
      kernel.col(j) = sqrt((x0-y0(j)).abs2() + (x1-y1(j)).abs2());
}

使用多线程,您需要测量挂钟时间。在这里(Haswell 有 4 个物理内核以 2.6GHz 运行)对于 N=20000,组装时间下降到 0.36 秒,矩阵 vector 乘积需要 0.24 秒,因此总共需要 0.6 秒,这比 MatLab 快,而我的 CPU 似乎更慢比你的。

关于c++ - 性能:Matlab vs C++ 矩阵 vector 乘法,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/46714235/

相关文章:

c++ - 使用 C++ 类成员函数(不能是静态的)作为 C 回调函数

c++ - 返回另一个成员函数调用的const引用的函数

matlab - 向录制的音频添加回声

java - 在Java中调用空方法会占用资源吗?

performance - Hyperledger Fabric可扩展性

c++ - 如何获取md5sum命令并获取代码中的字符串输出

c++ - 加载.so库C++

matlab - 使符号系数相等

arrays - Matlab - 数组中重复次数最多的值(不仅仅是模式)

angularjs - 如何验证 Angular JS 一次绑定(bind)性能增益?