/ Coding-Sucks

Matlab 連結 C Library

決定寫 C library 之前,建議是先用 profile reporter 確認程式的效能瓶頸。針對造成瓶頸的函式再花時間去打造其 C library 比較有經濟效益。

供 Matlab 連結使用之 C DLL 程式碼可分為兩個部分:第一部分為介面程式,用來處理 Matlab 與 C 之間的變數型態轉換;另一部分為功能主體,用來做函數所需要的運算。該 DLL 必須遵循 mexFunction 架構,並適當轉換輸入及輸出變數型態。

mexFunction 架構

典型的 mexFunction 程式碼的架構為

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    /* parse input arguments */
    /* create output arguments */
    /* base function */
}

首先必須要引入標頭檔 mex.h,這裡面有 Matlab 矩陣型態的定義以及各個輔助函式的 prototype。接著宣告進入點函數 mexFunction(...),差不多就是 main(...) 的意思,這函數的 prototype 是不可以改變的:

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);

mexFunction 包含四個引數,分別為

nlhs:呼叫函式時,等號左邊的變數數目,即 plhs 之個數
plhs:呼叫函式時,等號左邊之變數群
nrhs:呼叫函式時,等號右邊的變數數目,即 prhs 之個數
prhs:呼叫函式時,等號右邊之變數群

由於 Matlab 的函數呼叫允許多個輸出參數,如:

[M, I] = larger(A, B);

呼叫時 nlhs 以及 plhs 即用來存放等號左邊的輸出參數 M 與 I,nrhs 與 prhs 即存放等號右邊的輸入參數 A 及 B。此時,nlhs 及 nrhs 均為 2 。但如果使用者以

M = larger(A, B);

的形式呼叫,這時候 nlhs 會為 1。DLL 裡面可以依照 nlhs 及 nrhs 的不同,判斷使用者呼叫的輸入及輸出型態,以便進行 function overloading。

Matlab 中所有的變數,無論是矩陣或純量型態,在 mexFunction 中都以 mxArray 型態表示,純量就視為一個「大小為1x1之矩陣」。plhs 及 prhs 之變數型態為 mxArray **,亦即一「指向 Matlab 矩陣之指標陣列」。以前述 [M, I]=larger(A, B); 的例子來看,plhs[0] 指向 M 而 plhs[1] 指向 I,prhs[0] 指向 A 而 prhs[1] 則指向 B。

Matlab 矩陣與 C 陣列轉換

mexFunction 中的 mxArray 只是一個指標,實際仍然是 Matlab 之資料型態,並無法直接被 C 程式碼使用。要用 Matlab 提供的一些輔助函式進行轉換,從 mxArray 中取得可供 C 使用的數值。詳細可以參考 Matlab 手冊中關於 MX-Functions 的部分。

要從 mxArray 裡面取得資料,我們至少必須知道兩件事:矩陣的大小和矩陣的內容。矩陣大小可由 mxGetM 和 mxGetN 取得,矩陣內容則由 mxGetPr 以及 mxGetPi 取得。

int mxGetM(const mxArray *array_ptr);
int mxGetN(const mxArray *array_ptr);

mxGetM 及 mxGetN 分別用來取得矩陣之第一維以及第二維尺寸。函數傳回一個 int 整數,即為對應之尺寸大小。若矩陣超過二維,則必須使用 mxGetNumberOfDimensions 及 mxGetDimensions 函數。

double *mxGetPr(const mxArray *array_ptr);
double *mxGetPi(const mxArray *array_ptr);

mxGetPr 及 mxGetPi 則分別用來取出矩陣的實數與虛數部分,並回傳一個 double * 指向所取出資料陣列的開頭。無論原始矩陣的維度(dimension)是多少,所取出之資料均為一個一維陣列,因此必須配合由 mxGetM 及 mxGetN 取得的矩陣大小並計算正確之索引位置。

配合上述輔助函式,可利用下列程式碼,把 [M, I] = larger(A, B); 中矩陣 A 的實數部分變成 C 的變數:

double *A;
int mA, nA;
A = mxGetPr(prhs[0]);
mA = mxGetM(prhs[0]);
nA = mxGetN(prhs[0]);

接著便可以使用 A[i] 來存取矩陣,如果是二維陣列就使用 A[i+j*mA] 之方式。

建立輸出矩陣

函式既然是給 Matlab 用的,輸出當然也必須是 Matlab 矩陣。輸出矩陣的大小及內容是由 C DLL 決定,所以輸出必須由 DLL 負責建立。利用 mxCreateDoubleMatrix 我們可以建立 double 型態的輸出矩陣。

mxArray *mxCreateDoubleMatrix(int m, int n, mxComplexity ComplexFlag);

這函數可創造一個任意大小、實數或複數型態的二維 Matlab 矩陣。使用時給定矩陣的大小(m-by-n)以及矩陣的型態(mxREAL或mxCOMPLEX擇一)。函式即會創造一個新的 Matlab 矩陣,配置足夠的空間,並回傳一個 mxArray 指標指向所創造的新矩陣。若要創造多維矩陣或非實數/虛數型態之矩陣,可參考 mxCreateNumericArray 以及 mxCreateNumericMatrix。

此處用 mxCreateDoubleMatrix 建立輸出矩陣 M 及 I,並利用前述 mxGetPr 取得輸出矩陣對應之 C 陣列:

double *M, *I;
plhs[0] = mxCreateDoubleMatrix(mA, nA, mxREAL);
M = mxGetPr(plhs[0]);
plhs[1] = mxCreateDoubleMatrix(mA, nA, mxREAL);
I = mxGetPr(plhs[1]);

上面程式碼將 M 與 I 設定為與 A 相同大小之實數矩陣,並取出其資料指標 M 及 I 以備寫入輸出資料用。做好輸入及輸出參數的設置,接下來就依照一般 C 語法撰寫函式的主體功能部分,寫完後存檔為 larger.c:

#include "mex.h"
 
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *A, *B, *M, *I;
    int m, n, s, i;

    /* parse input arguments */ 
    A = mxGetPr(prhs[0]);
    B = mxGetPr(prhs[1]);
    m = mxGetM(prhs[0]);
    n = mxGetN(prhs[0]);

    /* create output arguments */ 
    plhs[0] = mxCreateDoubleMatrix(m, n, mxREAL);
    M = mxGetPr(plhs[0]);
    if(nlhs>1) {
      plhs[1] = mxCreateDoubleMatrix(m, n, mxREAL);
      I = mxGetPr(plhs[1]);
    }

    /* base function */
    s = m*n;
    if(nlhs==1) {
      for(i=0; i<s; i++) {
        if(A[i]>=B[i]) M[i] = A[i];
        else M[i] = B[i];
      }
    } else {
      for(i=0; i<s; i++) {
        if(A[i]>=B[i]) {
          M[i] = A[i];
          I[i] = 0;
        } else {
          M[i] = B[i];
          I[i] = 1;
        }
      }
    }
}

設定編譯器

在 Matlab 下編譯 C DLL 的指令為 mex。使用 mex 之前必須設定編譯使用的 C 編譯器。在 Matlab 下輸入 mex -setup 指令進行設定:

>> mex -setup
Select a compiler:
[1] Lcc C version 2.4 in C:\\PROGRAM FILES\\MATLAB704\\syslcc 

[0] None

Compiler:

mex 會偵測電腦中已安裝適合的C編譯器,Lcc 是 Matlab 內建的,可以用不用怕。

編譯

選好編譯器後就可以開始編譯了,假設剛剛寫好的檔案存檔為 larger.c,將 Matlab 工作目錄切換至剛剛存檔 larger.c 的目錄,執行:

mex larger.c

執行完畢後就可以在該目錄下看見編譯好的函式庫 larger.dll,如果失敗就要依照錯誤訊息檢查 DLL 程式碼以及 mex 設定。

測試

現在可以來測試看看編譯好的 DLL:

>> A=[1 2 3; 4 5 6; 7 8 9];
>> B=[9 8 7; 6 5 4; 3 2 1];
>> [M,I]=larger(A, B);
>> M

M =

     9     8     7
     6     5     6
     7     8     9

>> I

I =

     1     1     1
     1     0     0
     0     0     0

>>

也可以測試一下 M=larger(A,B); 是否能正常工作。接著試試看

>> larger(A,B)

------------------------------------------------------------------------
       Segmentation violation detected at Sat Dec 10 09:52:22 2005
------------------------------------------------------------------------

(以下刪除)

出現 segmentation fault,這是因為程式碼並沒有考慮「未指定任何輸出」,即 nlhs 為 0 時的情形。一般函式在沒有指定輸出的情況下,應該要傳回一個輸出結果。把程式碼中,

if(nlhs==1) {

這一行,修正為

if(nlhs<=1)

就可以解決這個問題。