前言

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

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

结果

下表显示了组装内核矩阵所需的时间与将矩阵与向量相乘所需的时间的比较.结果编译为矩阵大小N×N,其中N在10,000到40,000之间变化.哪个不是那么大.但有趣的是,Mat获得的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

for(unsigned long j=0;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(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;

}