2017-10-12 109 views
6

序言性能: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; 
} 
+5

你應該把你用來編譯你的C++代碼的標誌 – yakoudbz

+5

你可以清楚地改進你初始化矩陣的方式。首先,不要打電話::零,你正在浪費時間初始化一切。其次,嘗試查看矩陣是以行 - 主還是列 - 主的順序存儲的。如果它是列主要的,內部循環應該在每一行迭代! – yakoudbz

+0

由於'm'是一個,使用'VectorXd'可能會更快。 – m7913d

回答

9

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

+0

謝謝你的回答。添加'march = native'不會執行任何操作。關於'VectorXd'的建議,使其速度更快,但仍然比Matlab慢。儘管我嘗試了你的建議代碼片段,但是我收到了錯誤。拳頭是'auto'我把它改成'ArrayXd'。但我在編譯時遇到的主要錯誤是對omp_get_num_threads的未定義引用,我在源文件中添加了#include omp.h,並添加了'-fopenmp'標誌。仍然沒有你能幫我實施嗎?謝謝 –

+0

沒關係,經過一些更多的Google搜索之後,我通過添加'LDFLAGS = -g -fopenmp'來完成工作。我已將部分構建文件放入問題中。它現在正在工作,但我沒有看到任何時間上的變化。他們同樣不幸。還有其他建議嗎? –

+0

不,您不能用'ArrayXd'替換auto,因爲這會執行深層複製(昂貴)。你應該在你的'CFLAGS'中加入'-std = C++ 11'來啓用C++ 11支持。然後,正如我所說的,您需要測量掛牆時間(而不是CPU時間),例如使用C++ 11 [std :: chrono](http://www.cplusplus.com/reference/chrono/)。 – ggael