基於Mex混合編程的分子動力學模擬（MD）

11 人贊了文章

什麼是Mex？

Mex即"MATLAB executable"的簡寫，顧名思義，是使用Matlab運行C, C++ 和 Fortran程序的混合編程語言。

Why Mex?

Mex可以同時繼承Matlab和C++的優點，並可以通過合理的編程設計規避二者的缺點。確切地說是用一方的優點去彌補對方的缺點，最大化程序的運行效率。

Matlab的優點：

• 函數集成度高，直接支持矩陣的運算
• WorkSpace顯示數據直觀明了
• Plot方便快捷

Matlab的缺點：

• 付費，程序本身閉源（可以選用替代品：基於GNU協議的Octave，但是不推薦）
• 低到令人髮指的循環運行效率

C++的優點：

一千個程序員能有一千種理由誇獎C++，而作者能做的也只是在這之間再加上一種：它好厲害啊它好快啊它好全能啊

C++的缺點：

• 數據的輸入嚴重依賴於預設，這樣就造成每次修改Parameter都需要重新編譯；或者通過讀取ASCII文件，總之就是很麻煩。
• 數據的輸出嚴重依賴ASCII文件
• 沒有簡單易行並高度集成的Plot方案，需要依靠其他支持ASCII文件的繪圖軟體。

設計Mex的基本思路

• 所有的非矩陣運算全部交由重新編寫的Mex（C++內核）完成
• 所有的for循環全部交由重新編寫的Mex（C++內核）完成
• 所有的矩陣運算交由Matlab完成

如何安裝、編譯和調用Mex？

在安裝好C++編譯器後在Matlab命令行里輸入

mex -setup

Matlab應該會自動完成接下來的事，出現任何錯誤請善用搜索引擎，解決初始化Mex設置中出現的錯誤並不在本文的討論的範圍內。

在初始化Mex完畢後，將Matlab路徑定位到你的Mex Project目錄，我們就可以在Matlab命令行里輸入

mex *.cpp

將已經使用Mex語法修改完畢的*.cpp文件編譯成Matlab可以調用的二進位文件。

如果成功完成編譯，在Matlab的命令行或者腳本里，我們就可以像調用Matlab自身的函數一樣調用Mex函數。在此之前需要確認Matlab路徑已經定位到你的Mex二進位文件所在的目錄。

使用Mex進行混合編程的實例

下文將通過逐行解釋使用Mex語法修改完畢的control_velocity_berendsen.cpp文件來簡單介紹Mex混合編程的程序邏輯。

本文中所採用的示常式序是作者在Numerical simulation課程上還未發表的Project的一部分，故目前還無法公布完整的源代碼。完整的代碼將在結課後（2017年聖誕節前）公布。

由於採用了樊哲勇前輩在A C code for calculating diffusion coefficient in MD中的部分代碼，在此對前輩和原作者表示萬分的感激！

#include "math.h"#include "mex.h"/* include任何你需要的頭文件 * 其中"mex.h"是必須要包含的 */#define K_B 8.625e-5 // Boltzmanns constant in eV K^(-1)/* Mex語法中允許使用define來定義常數 */void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])/* mexFunction為此Mex源文件的主函數，它的變數固定，其中 * int nlhs為輸出矩陣個數，是一個數組，其內容為指針，指向數據類型mxArray * mxArray *plhs[]為輸出矩陣，是一個數組類型的指針，指向輸出矩陣所在的內存地址 * int nrhs 和 const mxArray *prhs[] 與上文所指帶類似 * 不過指帶的是輸入矩陣，並且用const保證輸入矩陣的地址不會被改變 */{ if (nrhs != 5) mexErrMsgTxt("control_velocity: Wrong number of input arguments.
"); if (nlhs > 1) mexErrMsgTxt("control_velocity: Too many output argumnents.
");/* 可省略但強烈建議保留的一步 * 目的是為了保證輸入和輸出矩陣數量能夠滿足C++內核的運算要求 * 如果省略則可能會在運行時造成內存讀寫錯誤，導致Matlab崩潰 */ #define v_IN prhs[0] #define m_IN prhs[1] #define T_0_IN prhs[2] #define ts_IN prhs[3] #define tau_T_IN prhs[4] #define v_OUT plhs[0]/* 可省略但強烈建議保留的一步 * 定義名稱可以更好地指代輸入矩陣和輸出矩陣 * 自上而下依次是 輸入矩陣1、輸入矩陣2、輸入矩陣3、輸入矩陣4、輸入矩陣5、輸出矩陣1、 */ int M_v_I = mxGetM(v_IN); int N_v_I = mxGetN(v_IN); int N_m_I = mxGetN(m_IN);/* mxGetM(v_IN)函數可以獲得v_IN矩陣的行數 * mxGetN(v_IN)函數可以獲得v_IN矩陣的列數 */ v_OUT = mxCreateDoubleMatrix(M_v_I, N_v_I, mxREAL);/* 為輸出矩陣初始化內存空間 */ double *v_I = mxGetPr(v_IN); double *m_I = mxGetPr(m_IN); double *T_0_I = mxGetPr(T_0_IN); double T_0 = T_0_I[0]; double *ts_I = mxGetPr(ts_IN); double ts = ts_I[0]; double *tau_T_I = mxGetPr(tau_T_IN); double tau_T = ts_I[0]; double *v_O = mxGetPr(v_OUT);/* 讀取輸入和輸出矩陣的數據 * 在下文中就可以如同操作C++的數組一樣操作數據 * 其中 T_0、ts 和 tau 均為數值 * 故將其二次定義為double類型的數據 *//* 所有的矩陣都按照一維矩陣進行保存 * 和一般的C++源文件不同，多維數組在Mex語法下會按列讀取對應的一維矩陣 * 而不是和一般的C++多維數組一樣按行讀取 * 這一點尤為重要 */ double temperature = 0.0; for (int i = 0; i < N_v_I; i++) { double v2 = v_I[M_v_I * i] * v_I[M_v_I * i] + v_I[M_v_I * i + 1] * v_I[M_v_I * i + 1] + v_I[M_v_I * i + 2] * v_I[M_v_I * i + 2]; temperature = temperature + m_I[i] * v2; } temperature = temperature / (3.0 * K_B * N_v_I); double scale_factor = sqrt(1 + ts / tau_T * ( T_0 / temperature - 1 )); for (int i = 0; i < M_v_I * N_v_I; i++) { v_O[i] = v_I[i] * scale_factor; }}/* 簡單的數據運算 * 在這裡使用了sqrt()函數 * 從側面驗證了Mex語法對C++原生函數庫的良好兼容 */

在Matlab下完成對control_velocity_berendsen.cpp的編譯後，如果需要調用control_velocity_berendsen則只需要像調用原生函數一樣

v = control_velocity_berendsen(v,m,T_0,ts,tau);

請在務必在運行前再次確認輸入變數和輸出變數的個數和對應性。

運行效率

機器配置：

相同邏輯的各種編程語言的運行時間：

• C++:

• Mex:

• Matlab:很抱歉在本文的撰寫過程當中，作者嘗試過運行Matlab所編寫的程序，目前已經過了五十分鐘並沒有得到任何輸出的結果。

綜合來看Mex比原生C++程序效率高的原因在於數據的輸入與輸出全部由程序內的內存交換而來，省去了生成和讀取ASCII文件的過程。

總的來說，Mex的編程比Matlab和C++都來得費勁。

但是由於在MD模擬中需要不斷地修改Parameter，使得基於Matlab進行數據輸入的Mex在預設參數方面更方便。

再者，在操作者的感官上，所有的數據都在程序內的內存中進行交換，並直接通過WorkSpace展示給操作者，顯得十分直觀。可以直接在Matlab內進行Plot輸出，並不依靠其他軟體。

Mex的函數可以是集成在同一個主Mex函數里，但是在實際操作中我們更願意模塊化Mex函數，使得我們在不同的Project里能夠調用以前寫過的Mex函數，僅僅只需要將Mex編譯後的二進位文件複製粘貼即可。或者直接維護一個屬於自己的Mex函數庫，如果有時間你甚至可以將Matlab的原生函數用Mex重新編寫。

由於Matlab腳本的特性，將腳本程序發送給其他人，他人相當於得到了我們程序的源代碼。在某些情況下我們並不希望公開我們的源代碼，此時我們就可以通過Mex將我們的程序加殼，在對方不進行長時間的反編譯的情況下我們的源代碼是相當安全的。即便如此，基於GNU協議的精神作者還是鼓勵公開源代碼，畢竟知識無價嘛！

這篇科普小短文是作者在失戀後失眠寫的，曾想過要在github上做一個開源的Mex based MD項目，但是由於臨近研究生畢業和準備語言考試沒有太多空閑的時間，僅當作拋磚引玉。

在日後作者會繼續寫幾篇關於Mex編程的文章，如果有時間的話。

（有心人可以看看我在知乎上的回復，大部分坑都沒有填完！）