性能：Matlab的VS C++矩阵向量乘法

Question

序言

前段时间我问一个关于VS的Python（Performance: Matlab vs Python）Matlab的性能问题。我很惊讶Matlab比Python更快，特别是在meshgrid。在讨论这个问题时，有人指出我应该在Python中使用包装来调用我的C++代码，因为我也可以使用C++代码。我在C++，Matlab和Python中拥有相同的代码。

虽然这样做，我再次惊讶地发现matlab比矩阵程序集和计算中的C++更快。我有一个稍大的代码，我正在调查一段矩阵向量乘法。较大的代码在多个实例中执行这种乘法。总的来说，C++中的代码比Matlab快得多（因为在Matlab中调用的函数有一个开销等），但是Matlab似乎在矩阵向量乘法（底部的代码片段）中优于C++。

结果

下表显示了所花费的时间来组装内核矩阵，它需要与矢量乘以矩阵的时间的比较。结果汇编为矩阵大小NxN，其中N从10,000到40,000变化。这不是那么大。但有趣的是，Matlab的性能优于C++，N的性能就越好。 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); 

类（使用meshgrid）：

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

Answer 1

如MatLab的依赖英特尔的MKL库矩阵产品的评论说，这是这种类型的操作最快速的图书馆。尽管如此，Eigen本身应该能够提供类似的性能。为此，请确保使用最新的Eigen（例如3。4），以及适当的编译标志启用AVX/FMA（如果可用）和多线程：

-O3 -DNDEBUG -march=native 

由于charges_是一个载体，更好的使用VectorXd来征知道你想要一个矩阵向量的产品，而不是一个矩阵定矩阵一。

如果您拥有英特尔的MKL，那么您也可以让Eigen uses it获得与MatLab完全相同的性能，用于此精确操作。

关于大会，更好地反两个循环，使量化，然后要使OpenMP多线程（添加-fopenmp为编译器标志），使并行的最外层循环运行，最后你可以使用简化代码征：

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.36s，矩阵向量乘积为0.24s，因此总共为0.6s，这比MatLab快，而我的CPU似乎比较慢比你的。