性能：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似乎比較慢比你的。