性能:Matlab的VS C++矩阵向量乘法
前段时间我问一个关于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 [email protected] $(LIB_PATH)
.cpp.o:
$(CC) $(CFLAGS) $(KERNEL) $(INDEX) $< -o [email protected]
`
# 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_
是一个载体,更好的使用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似乎比较慢比你的。
谢谢你的回答。添加'march = native'不会执行任何操作。关于'VectorXd'的建议,使其速度更快,但仍然比Matlab慢。尽管我尝试了你的建议代码片段,但是我收到了错误。拳头是'auto'我把它改成'ArrayXd'。但我在编译时遇到的主要错误是对omp_get_num_threads的未定义引用,我在源文件中添加了#include omp.h,并添加了'-fopenmp'标志。仍然没有你能帮我实施吗?谢谢 –
没关系,经过一些更多的Google搜索之后,我通过添加'LDFLAGS = -g -fopenmp'来完成工作。我已将部分构建文件放入问题中。它现在正在工作,但我没有看到任何时间上的变化。他们同样不幸。还有其他建议吗? –
不,您不能用'ArrayXd'替换auto,因为这会执行深层复制(昂贵)。你应该在你的'CFLAGS'中加入'-std = C++ 11'来启用C++ 11支持。然后,正如我所说的,您需要测量挂墙时间(而不是CPU时间),例如使用C++ 11 [std :: chrono](http://www.cplusplus.com/reference/chrono/)。 – ggael
你应该把你用来编译你的C++代码的标志 – yakoudbz
你可以清楚地改进你初始化矩阵的方式。首先,不要打电话::零,你正在浪费时间初始化一切。其次,尝试查看矩阵是以行 - 主还是列 - 主的顺序存储的。如果它是列主要的,内部循环应该在每一行迭代! – yakoudbz
由于'm'是一个,使用'VectorXd'可能会更快。 – m7913d