Matlab以MEX方式呼叫C原始碼
如果我有一個用C語言寫的函式,實現了一個功能,如一個簡單的函式:
double add(double x, double y)
{
return x + y;
}
現在我想要在Matlab中使用它,比如輸入:
>> a = add(1.1, 2.2)
3.3000
要得出以上的結果,那應該怎樣做呢?
解決方法之一是要通過使用MEX檔案,MEX檔案使得呼叫C函式和呼叫Matlab的內建函式一樣方便。MEX檔案是由原C程式碼加上MEX檔案專用的介面函式後編譯而成的。可以這樣理解,MEX檔案實現了一種介面,它把在Matlab中呼叫函式時輸入的自變數通過特定的介面調入了C函式,得出的結果再通過該介面調回Matlab。該特定介面的操作,包含在mexFunction這個函式中,由使用者具體設定。
所以現在我們要寫一個包含add和mexFunction的C檔案,Matlab呼叫函式,把函式中的自變數(如上例中的1.1和2.2)傳給 mexFunction的一個引數,mexFunction把該值傳給add,把得出的結果傳回給mexFunction的另一個引數,Matlab通過該引數來給出在Matlab語句中呼叫函式時的輸出值(如上例中的a)。
值得注意的是,mex檔案是與平臺有關的,以我的理解,mex檔案就是另類的動態連結庫。在matlab6.5中使用mex -v 選項,你可以看到最後mex階段有類似如下的資訊:
--> "del _lib94902.obj"
--> "del "test.exp""
--> "del "test.lib""
也就是說,雖然在matlab6.5生成的是dll檔案,但是中間確實有過lib檔案生成。
比如該C檔案已寫好,名為add.c。那麼在Matlab中,輸入:
>> mex add.c
就能把add.c編譯為MEX檔案(編譯器的設定使用指令mex -setup),在Windows中,MEX檔案型別為mexw32,即現在我們得出add.mexw32檔案。現在,我們就可以像呼叫M函式那樣呼叫 MEX檔案,如上面說到的例子。所以,通過MEX檔案,使用C函式就和使用M函式是一樣的了。
我們現在來說mexFunction怎樣寫。
mexFunction的定義為:
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
}
可以看到,mexFunction是沒返回值的,它不是通過返回值把結果傳回Matlab的,而是通過對引數plhs的賦值。mexFunction的四個引數皆是說明Matlab呼叫MEX檔案時的具體資訊,如這樣呼叫函式時:
>> b = 1.1; c = 2.2;
>> a = add(b, c)
mexFunction四個引數的意思為:
nlhs = 1,說明呼叫語句左手面(lhs-left hand side)有一個變數,即a。
nrhs = 2,說明呼叫語句右手面(rhs-right hand side)有兩個自變數,即b和c。
plhs是一個數組,其內容為指標,該指標指向資料型別mxArray。因為現在左手面只有一個變數,即該陣列只有一個指標,plhs[0]指向的結果會賦值給a。
prhs和plhs類似,因為右手面有兩個自變數,即該陣列有兩個指標,prhs[0]指向了b,prhs[1]指向了c。要注意prhs是const的指標陣列,即不能改變其指向內容。
因為Matlab最基本的單元為array,無論是什麼型別也好,如有double array、 cell array、 struct array……所以a,b,c都是array,b = 1.1便是一個1x1的double array。而在C語言中,Matlab的array使用mxArray型別來表示。所以就不難明白為什麼plhs和prhs都是指向mxArray型別的指標陣列。
完整的add.c如下:
#include "mex.h" // 使用MEX檔案必須包含的標頭檔案
// 執行具體工作的C函式
double add(double x, double y)
{
return x + y;
}
// MEX檔案介面函式
void mexFunction(int nlhs,mxArray *plhs[], int nrhs,const mxArray *prhs[])
{
double *a;
double b, c;
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
a = mxGetPr(plhs[0]);
b = *(mxGetPr(prhs[0]));
c = *(mxGetPr(prhs[1]));
*a = add(b, c);
}
mexFunction的內容是什麼意思呢?我們知道,如果這樣呼叫函式時:
>> output = add(1.1, 2.2);
在未涉及具體的計算時,output的值是未知的,是未賦值的。所以在具體的程式中,我們建立一個1x1的實double矩陣(使用 mxCreateDoubleMatrix函式,其返回指向剛建立的mxArray的指標),然後令plhs[0]指向它。接著令指標a指向plhs [0]所指向的mxArray的第一個元素(使用mxGetPr函式,返回指向mxArray的首元素的指標)。同樣地,我們把prhs[0]和prhs [1]所指向的元素(即1.1和2.2)取出來賦給b和c。於是我們可以把b和c作自變數傳給函式add,得出給果賦給指標a所指向的mxArray中的元素。因為a是指向plhs[0]所指向的mxArray的元素,所以最後作輸出時,plhs[0]所指向的mxArray賦值給output,則
output便是已計算好的結果了。
上面說的一大堆指向這指向那,什麼mxArray,初學者肯定都會被弄到頭暈眼花了。很抱歉,要搞清楚這些亂糟糟的關係,只有多看多練。
實際上mexFunction是沒有這麼簡單的,我們要對使用者的輸入自變數的個數和型別進行測試,以確保輸入正確。如在add函式的例子中,使用者輸入char array便是一種錯誤了。
從上面的講述中我們總結出,MEX檔案實現了一種介面,把C語言中的計算結果適當地返回給Matlab罷了。當我們已經有用C編寫的大型程式時,大可不必在 Matlab裡重寫,只寫個介面,做成MEX檔案就成了。另外,在Matlab程式中的部份計算瓶頸(如迴圈),可通過MEX檔案用C語言實現,以提高計算速度。
以上是對mex檔案的初步認識,下面詳細介紹如何用c語言編寫mex檔案:
1 為什麼要用C語言編寫MEX檔案
MATLAB是矩陣語言,是為向量和矩陣操作設計的,一般來說,如果運算可以用向量或矩陣實現,其運算速度是非常快的。但若運算中涉及到大量的迴圈處理,MATLAB的速度的令人難以忍受的。解決方法之一為,當必須使用for迴圈時,把它寫為MEX檔案,這樣不必在每次執行迴圈中的語句時MATLAB都對它們進行解釋。
2 編譯器的安裝與配置
要使用MATLAB編譯器,使用者計算機上應用事先安裝與MATLAB適配的以下任何一種ANSI C/C++編譯器:
5.0、6.0版的MicroSoft Visual C++(MSVC)
5.0、5.2、5.3、5.4、5.5版的Borland C++
LCC(由MATLAB自帶,只能用來產生MEX檔案)
下面是安裝與配置MATLAB編譯器應用程式MEX的設定的步驟:
(1)在MATLAB命令視窗中執行mex –setup,出現下列提示:
Please choose your compiler for building external interface (MEX) files:
Would you like mex to locate installed compilers [y]/n?
(2)選擇y,MATLAB將自動搜尋計算機上已安裝的外部編譯器的型別、版本及所在路徑,並列出來讓使用者選擇:
Select a compiler:
[1] Borland C++Builder version 6.0 in C:\Program Files\Borland
[2] Digital Visual Fortran version 6.0 in C:\Program Files\Microsoft Visual Studio
[3] Lcc C version 2.4 in D:\MATLAB6P5P1\sys\lcc
[4] Microsoft Visual C/C++ version 6.0 in C:\Program Files\Microsoft Visual Studio
[0] None
Compiler:
(3)選擇其中一種(在這裡選擇了3),MATLAB讓使用者進行確認:
Please verify your choices:
Compiler: Lcc C 2.4
Location: D:\MATLAB6P5P1\sys\lcc
Are these correct?([y]/n):
(4)選擇y,結束MATLAB編譯器的配置。
3 一個簡單的MEX檔案例子
【例1】用m檔案建立一個1000×1000的Hilbert矩陣。
tic
m=1000;
n=1000;
a=zeros(m,n);
for i=1:1000
for j=1:1000
a(i,j)=1/(i+j);
end
end
toc
在matlab中新建一個Matlab_1.cpp 檔案並輸入以下程式:
#include "mex.h"
//計算過程
void hilb(double *y,int n)
{
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
*(y+j+i*n)=1/((double)i+(double)j+1);
}
//介面過程
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
double x,*y;
int n;
if (nrhs!=1)
mexErrMsgTxt("One inputs required.");
if (nlhs != 1)
mexErrMsgTxt("One output required.");
if (!mxIsDouble(prhs[0])||mxGetN(prhs[0])*mxGetM(prhs[0])!=1)
mexErrMsgTxt("Input must be scalars.");
x=mxGetScalar(prhs[0]);
plhs[0]=mxCreateDoubleMatrix(x,x,mxREAL);
n=mxGetM(plhs[0]);
y=mxGetPr(plhs[0]);
hilb(y,n);
}
該程式是一個C語言程式,它也實現了建立Hilbert矩陣的功能。在MATLAB命令視窗輸入以下命令:mex Matlab_1.cpp,即可編譯成功。進入該資料夾,會發現多了兩個檔案:Matlab_1.asv和Matlab_1.dll,其中Matlab_1.dll即是MEX檔案。執行下面程式:
tic
a=Matlab_1(1000);
toc
elapsed_time =
0.0470
由上面看出,同樣功能的MEX檔案比m檔案快得多。
4 MEX檔案的組成與引數
MEX檔案的原始碼一般由兩部分組成:
(1)計算過程。該過程包含了MEX檔案實現計算功能的程式碼,是標準的C語言子程式。
(2)入口過程。該過程提供計算過程與MATLAB之間的介面,以入口函式mxFunction實現。在該過程中,通常所做的工作是檢測輸入、輸出引數個數和型別的正確性,然後利用mx-函式得到MATLAB傳遞過來的變數(比如矩陣的維數、向量的地址等),傳遞給計算過程。
MEX檔案的計算過程和入口過程也可以合併在一起。但不管那種情況,都要包含#include "mex.h",以保證入口點和介面過程的正確宣告。注意,入口過程的名稱必須是mexFunction,並且包含四個引數,即:
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
其中,引數nlhs和nrhs表示MATLAB在呼叫該MEX檔案時等式左端和右端變數的個數,例如在MATLAB命令視窗中輸入以下命令:
[a,b,c]=Matlab_1(d,e,f,g)
則nlhs為3,nrhs為4。
MATLAB在呼叫MEX檔案時,輸入和輸出引數儲存在兩個mxArray*型別的指標陣列中,分別為prhs[]和plhs[]。prhs[0]表示第一個輸入引數,prhs[1]表示第二個輸入引數,…,以此類推。如上例中,d→prhs[0],e→prhs[1],f→prhs[2],f→prhs[3]。同時注意,這些引數的型別都是mxArray *。
介面過程要把引數傳遞給計算過程,還需要從prhs中讀出矩陣的資訊,這就要用到下面的mx-函式和mex-函式。
5 常用的mex-函式和mx-函式
在MATLAB6.5版本中,提供的mx-函式有106個,mex-函式有38個,下面我們僅介紹常用的函式。
5.1入口函式mexFunction
該函式是C MEX檔案的入口函式,它的格式是固定的:
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
說明:MATLAB函式的呼叫方式一般為:[a,b,c,…]=被呼叫函式名稱(d,e,f,…),nlhs儲存了等號左端輸出引數的個數,指標陣列plhs具體儲存了等號左端各引數的地址,注意在plhs各元素針向的mxArray記憶體未分配,需在介面過程中分配記憶體;prhs儲存了等號右端輸入引數的個數,指標陣列prhs具體儲存了等號右端各引數的地址,注意MATLAB在呼叫該MEX檔案時,各輸入引數已存在,所以在介面過程中不需要再為這些引數分配記憶體。
5.2出錯資訊釋出函式mexErrMsgTxt,mexWarnMsgTxt
兩函式的具體格式如下:
#include "mex.h"
void mexErrMsgTxt(const char *error_msg);
void mexWarnMsgTxt(const char *warning_msg);
其中error_msg包含了要顯示錯誤資訊,warning_msg包含要顯示的警告資訊。兩函式的區別在於mexErrMsgTxt顯示出錯資訊後即返回到MATLAB,而mexWarnMsgTxt顯示警告資訊後繼續執行。
5.3 mexCallMATLAB和mexString
兩函式具體格式如下:
#include "mex.h"
int mexCallMATLAB(int nlhs, mxArray *plhs[],
int nrhs, mxArray *prhs[], const char *command_name);
int mexString(const char *command);
mexCallMATLAB前四個引數的含義與mexFunction的引數相同,command_name可以MATLAB內建函式名、使用者自定義函式、M檔案或MEX檔名構成的字串,也可以MATLAB合法的運算子。
mexString用來操作MATLAB空間已存在的變數,它不返回任何引數。
mexCallMATLAB與mexString差異較大,請看下面的例子。
【例2】試用MEX檔案求5階完全圖鄰接矩陣 的特徵值及對應的特徵向量。
5階完全圖的鄰接矩陣為:(這裡找不到圖片了,抱歉。不過不會影響您對本文的理解。)
下面是求該矩陣的MEX檔案。
[Matlab_2.cpp]
#include "mex.h"
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
double x;
mxArray *y,*z,*w;
int n;
if (nrhs!=1)
mexErrMsgTxt("One inputs required.");
if (nlhs != 3)
mexErrMsgTxt("Three output required.");
if (!mxIsDouble(prhs[0])||mxGetN(prhs[0])*mxGetM(prhs[0])!=1)
mexErrMsgTxt("Input must be a scalar.");
x=mxGetScalar(prhs[0]);
plhs[0]=mxCreateDoubleMatrix(x,x,mxREAL);
plhs[1]=mxCreateDoubleMatrix(x,x,mxREAL);
plhs[2]=mxCreateDoubleMatrix(x,x,mxREAL);
n=mxGetM(plhs[0]);
y=plhs[0];
z=plhs[1];
w=plhs[2];
//利用mexCallMATLAB計算特徵值
mexCallMATLAB(1,&plhs[1],1,prhs,"ones");
mexCallMATLAB(1,&plhs[2],1,prhs,"eye");
mexCallMATLAB(1,&plhs[0],2,&plhs[1],"-");
mexCallMATLAB(2,&plhs[1],1,&plhs[0],"eig");
//演示mexString的功能
mexString("y=y*2");
mexString("a=a*2");
}
在MATLAB命令視窗輸入以下命令:
>> mex Matlab_2.cpp
>> clear
>> a=magic(5)
a =
17 24 1 8 15
23 5 7 14 16
4 6 13 20 22
10 12 19 21 3
11 18 25 2 9
>> [y,z,w]=Matlab_2(5)
??? Undefined function or variable 'y'.
a =
34 48 2 16 30
46 10 14 28 32
8 12 26 40 44
20 24 38 42 6
22 36 50 4 18
y =
0 1 1 1 1
1 0 1 1 1
1 1 0 1 1
1 1 1 0 1
1 1 1 1 0
z =
0.8333 -0.1667 -0.1667 0.2236 0.4472
-0.1667 0.8333 -0.1667 0.2236 0.4472
-0.1667 -0.1667 0.8333 0.2236 0.4472
-0.5000 -0.5000 -0.5000 0.2236 0.4472
0 0 0 -0.8944 0.4472
w =
-1 0 0 0 0
0 -1 0 0 0
0 0 -1 0 0
0 0 0 -1 0
0 0 0 0 4
由上面可以看出,K5的特徵值為–1和4,其中–1是四重根。MATLAB提供了mexGetVariable、mexPutVariable函式,以實現MEX空間與其它空間交換資料的任務,具體可以參看MATLAB幫助文件。
5.4建立二維雙精度矩陣函式mxCreateDoubleMatrix
其格式具體如下:
#include "matrix.h"
mxArray *mxCreateDoubleMatrix(int m, int n, mxComplexity ComplexFlag);
其中m代表行數,n代表列數,ComplexFlag可取值mxREAL 或mxCOMPLEX。如果建立的矩陣需要虛部,選擇mxCOMPLEX,否則選用mxREAL。
類似的函式有:
mxCreateCellArray |
建立n維元胞mxArray |
mxCreateCellMatrix |
建立二維元胞mxArray |
mxCreateCharArray |
建立n維字串mxArray |
mxCreateCharMatrixFromStrings |
建立二維字串mxArray |
mxCreateDoubleMatrix |
建立二維雙精度浮點mxArray |
mxCreateDoubleScalar |
建立指定值的二維精度浮點mxArray |
mxCreateLogicalArray |
建立n維邏輯mxArray,初值為false |
mxCreateLogicalMatrix |
建立二維邏輯mxArray,初值為false |
mxCreateLogicalScalar |
建立指定值的二維邏輯mxArray |
mxCreateNumericArray |
建立n維數值mxArray |
mxCreateNumericMatrix |
建立二維數值mxArray,初值為0 |
mxCreateScalarDouble |
建立指定值的雙精度mxArray |
MxCreateSparse |
建立二維稀疏mxArray |
mxCreateSparseLogicalMatrix |
建立二維稀疏邏輯mxArray |
MxCreateString |
建立指定字串的1 n的串mxArray |
mxCreateStructArray |
建立n維架構mxArray |
mxCreateStructMatrix |
建立二維架構mxArray |
5.5 獲取行維和列維函式mxGetM、mxGetN
其格式如下:
#include "matrix.h"
int mxGetM(const mxArray *array_ptr);
int mxGetN(const mxArray *array_ptr);
與之相關的還有:
mxSetM:設定矩陣的行維
mxSetN:設定矩陣的列維
5.6 獲取矩陣實部和虛部函式mxGetPr、mxGetPi
其格式如下:
#include "matrix.h"
double *mxGetPr(const mxArray *array_ptr);
double *mxGetPi(const mxArray *array_ptr);
與之相關的函式還有:
mxSetPr:設定矩陣的實部
mxSetPi:設定矩陣的虛部
【例3】實現字串的倒序輸出。
#include "mex.h"
void revord(char *input_buf,int buflen,char *output_buf)
{
int i;
//實現字串倒序
for(i=0;i<buflen-1;i++)
*(output_buf+i)=*(input_buf+buflen-i-2);
}
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
//定義輸入和輸出參量的指標
char *input_buf,*output_buf;
int buflen,status;
//檢查輸入引數個數
if(nrhs!=1)
mexErrMsgTxt("One input required.");
else if(nlhs>1)
mexErrMsgTxt("Too many output arguments.");
//檢查輸入引數是否是一個字串
if(mxIsChar(prhs[0])!=1)
mexErrMsgTxt("Input must be a string.");
//檢查輸入引數是否是一個行變數
if(mxGetM(prhs[0])!=1)
mexErrMsgTxt("Input must a row vector.");
//得到輸入字串的長度
buflen=(mxGetM(prhs[0])*mxGetN(prhs[0]))+1;
//為輸入和輸出字串分配記憶體
input_buf=mxCalloc(buflen,sizeof(char));
output_buf=mxCalloc(buflen,sizeof(char));
//將輸入參量的mxArray結構中的數值拷貝到C型別字串指標
status=mxGetString(prhs[0],input_buf,buflen);
if(status!=0)
mexWarnMsgTxt("Not enough space. String is truncated.");
//呼叫C程式
revord(input_buf,buflen,output_buf);
plhs[0]=mxCreateString(output_buf);
}
這個程式中需要注意的地方是mxCalloc函式,它代替了標準C程式中的calloc函式用於動態分配記憶體,而mxCalloc函式採用的是MATLAB的記憶體管理機制,並將所有申請的記憶體初始化為0,因此凡是C程式碼需要使用calloc函式的地方,對應的Mex檔案應該使用mxCalloc函式。同樣,凡是C程式碼需要使用realloc函式的地方,對應的Mex檔案應該使用mxRealloc函式。
在MATLAB命令視窗中對revord.cpp程式程式碼編譯連結:
>> mex revord.cpp
在MATLAB命令視窗中對C-MEX檔案revord.dll進行測試:
>> x='I am student.';
>> revord(x)
ans =
.tneduts ma I
終於寫完了,相信大家對mex檔案應該有點熟悉了,具體還要到實際應用中慢慢體會。