性能:Matlab与C++矩阵矢量乘法

前导

不久前，我问了一个关于Matlab与Python性能的问题(性能：Matlab与Python)。我很惊讶Matlab比Python快，尤其是在meshgrid中。在讨论这个问题时，有人告诉我，我应该使用Python中的包装器来调用我的C++代码，因为我也可以使用C++代码。我在C++、Matlab和Python中都有相同的代码。

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

结果

下表显示了组装内核矩阵所需的时间与将矩阵与向量相乘所需时间的比较。针对矩阵大小CCD_ 2编译结果，其中CCD_ 3在10000到40000之间变化。没有那么大。但有趣的是，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;
}