Matlab 教材：MEX 範例--迴圈

請讀者抓下這個檔案：

ex.m
這是一個 Matlab 的 function，需要兩個向量當輸入，最多會有兩個輸出。 假設 C、D 是兩個向量，而且 C 比 D 長很多，當你下指令
[A B]=ex(C,D)
A 是 C 的奇數項和 D 疊積的結果，B 是 C 的偶數項和 D 疊積的結果。 讀者不用在意疊積的定義，這不是本文的重點， 請直接看下面的說明就夠了。

這個程式的原始碼如下： 

function [A,B]=ex(C,D)
LN=length(C);
SN=length(D);
A=zeros(1,ceil((LN-2*SN+2)/2));
B=zeros(1,floor((LN-2*SN+2)/2));
for i=1:2:LN-2*SN+2
    for j=1:SN;
	A((i+1)/2)=A((i+1)/2)+C(i-2+2*j)*D(j);
    end;
end;
for i=2:2:LN-2*SN+2
    for j=1:SN;
	B(i/2)=B(i/2)+C(i-2+2*j)*D(j);
    end;
end;
因為這個程式用了大量的迴圈，所以當 C 或 D 的長度很長時﹝以百萬為單位﹞， 執行時間將會拖得很長。 為了改善效率，我們有必要利用 C 重寫這個程式。

這裡是同一個程式的 C 版本：

ex2.c
請讀者下載編譯後依序輸入以下的指令：
  1. C=rand(1,1000000);
  2. D=rand(1,35);
  3. tic;[A B]=ex(C,D);toc;
  4. tic;[A B]=ex2(C,D);toc;
在我的電腦上，ex2 的效率是 ex 的兩倍以上。 這並不是一個非常明顯的例子，以我的電腦而言，只是從四秒到兩秒的提升。 在我過去的個人經驗，曾經有從六個小時到三秒鐘的飛躍性提升的例子。 這說明了用把迴圈交由 C 來完成確實是比較有效率的。

現在，我們來看看 ex2.c 的原始碼： 

#include "mex.h"
  
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    int i,j,ln,sn;
    double *lv,*sv,*v1,*v2;
    if (nlhs>2)
	mexErrMsgTxt("Two return values");
    if (nrhs!=2)
	mexErrMsgTxt("Require 2 inputs");
    for (i=0;i<2;++i) {
	if (!mxIsDouble(prhs[i]))
	    mexErrMsgTxt("Only for double value.");
	if (mxIsComplex(prhs[i]))
	    mexErrMsgTxt("Not for complex value.");
	if (mxGetM(prhs[i])!=1)
	    mexErrMsgTxt("Only for row vector.");
    }
    ln=mxGetN(prhs[0]);
    sn=mxGetN(prhs[1]);
    if (ln<sn*2-1)
	mexErrMsgTxt("Empty result.");
    plhs[0]=mxCreateDoubleMatrix(1,(ln-2*sn+3)/2,mxREAL);
    plhs[1]=mxCreateDoubleMatrix(1,(ln-2*sn+2)/2,mxREAL);
    lv=mxGetPr(prhs[0]);
    sv=mxGetPr(prhs[1]);
    v1=mxGetPr(plhs[0]);
    v2=mxGetPr(plhs[1]);
    for(i=0;i<(ln-2*sn+3);i+=2){
	for(j=0;j<sn;++j)
	    *(v1+i/2)+=(*(lv+i+2*j))*(*(sv+j));
    }
    for(i=0;i<(ln-2*sn+2);i+=2){
	for(j=0;j<sn;++j)
	    *(v2+i/2)+=(*(lv+i+2*j+1))*(*(sv+j));
    }
}
這裡用到了大量 Matlab 提供的函式。 Matlab 提供給 C 的函式分成兩類：
mx 開頭
這一類型都是屬於矩陣屬性的操作，例如為矩陣配置記憶體。
mex 開頭
提供有關於 Matlab 的操作，如在 Matlab 上顯示錯誤訊息。
所以 mexFunction 和 Matlab 的操作相關─它提供了一個進入點。

接著來介紹 mexFunction。 mexFunction 提供了四個參數：nlhs、plhs、nrhs 及 prhs。 在 nrhs 和 prhs 中存放著傳給 subroutine 的參數及參數的資訊 (rhs 就是等號右邊的意思 right hand side)， 而 nlhs 和 plhs 則存放著要回傳給 matlab 的結果及關於結果的資訊 (所以 lhs 就是等號左邊)。

現在看一下這些參數代表的意義， nlhs 是呼叫 subroutine 時使用者希望得到的輸出數目， nrhs 則是傳進 subroutine 的參數數目。例如下指令

[A B C]=func(D,E,F,G)
nlhs 會是 3，nrhs 會是 4。 nlhs 並不是由 subroutine 決定的， 而是使用者呼叫 subroutine 時由 Matlab 根據使用者的呼叫決定的。 這使得 subroutine 可以決定是否可以省略某些計算。 例入我們把上面的指令改成
[A B]=func(D,E,F,G)
這時 nlhs 為 2，func 就可以選擇不算出 C。 同樣的情況也適用於 nrhs， subroutine 可以隨著參數數目的不同而有不同的計算方式。

plhs 和 prhs 都是指向一個序列的指標， 這個序列放的都是指向 mxArray 結構的指標。 在 Matlab 中只有一種物件型別，就是 array。 所以 Matlab 在呼叫 C 寫成的 subroutine 時， 傳進去的要是 array，得到的回傳也必須要是 array。 在 C 裡面，Matlab 傳進來的 array 全部都是 mxArray 的結構， 而傳回 Matlab 的 array 也同樣都是 mxArray 的結構。 所以 plhs 和 prhs 實際上就是指向一堆要傳回去的以及要拿進來的 array。 plhs 指向要回傳給 Matlab 的 array， prhs 則是指向 Matlab 傳給 subroutine 的 array。 以

[A B C]=func(D,E,F,G)
為例，prhs[0] 指向矩陣 D，prhs[1] 指向 E。 而 plhs[0] 指到的矩陣在 subroutine 結束後會自動變成矩陣 A。 所以只要為 plhs 配置好記憶體， subroutine 結束後 Matlab 會自動取得 plhs 指定的矩陣，而不用特別去回傳 plhs。 接著介紹一些用來操作 prhs 和 plhs 的函式。

在應用上，mxArray 內的資料型態多半會是 double，但是也會有例外。 如果希望知道一個 mxArray 是存放 double 的變數， 我們可以下指令 mxIsDouble(prhs[0]) (prhs[0] 是 mxArray* 型態)， 如果 prhs[0] 指到的 mxArray 代表的是 double 的 array，則這個函式傳回 true， 否則傳回 false。

在 mxArray 內有兩個指標，一個是 pr，一個是 pi，分別指向兩個維度相同的矩陣。 pr 指向的是 array 的實數部分，pi 指向 array 的虛數部分。 當 array 是複數時，pr 和 pi 指向的兩個矩陣相加代表整個完整的 array。 當 array 是實數時，pi 會指向 null，這時 pr 就足夠代表整個 array。 如果想要知道一個 array 是否為複數，可以用 mxIsComplex。

當你想要得到 plhs[0] 指向的 mxArray 的 pr 指標，可以下指令

pr=mxGetPr(plhs[0]);
如果想得到 pi 指標，可以下指令
pi=mxGetPi(plhs[0]);
pr 和 pi 多半是 double* 型態， 藉由這兩個指標，我們可以對 plhs[0] 指向的 array 做修改。 如果這個 array 的維度是 m 乘 n， 則第 i 列第 j 行會在 pr+m*(j-1)+i-1 以及 pi+m*(j-1)+i-1。 有了這些資訊，剩下的就是對指標的操作了。

如果要知道一個 plhs[0] 指到的 array 的列數，可以下指令

mxGetM(plhs[0]);
會傳回一個整數，表示列數。 而行數則可以用
mxGetN(plhs[0]);
來得到。

在 ex2.c 中，我不但用到了上面介紹的函式，同時也用了 mexErrMsgTxt。 這是一個檢查用的函式，它會強制中斷 subroutine， 回到 matlab 並在 matlab 顯示錯誤訊息 (使用者自己寫入)。 配合上面介紹的函式，我對我的程式做了以下的限制：

請讀者自行對照程式找出這些程式碼。

現在還剩下一件事，就是如何建立一個新的 array。 通常我們會直接將建立好的 array (在 C 中會是 mxArray 結構) 交給 plhs。 這代表我們需要做記憶體的配置，這個函式也由 matlab 提供：

plhs[0]=mxCreateDoubleMatrix(1,(ln-2*sn+3)/2,mxREAL);
mxCreateDoubleMatrix 會建立一個 double 型態的 array， array 的列數由第一個參數指定，行數由第二個參數指定。 第三個參數如果是 mxREAL， 則建立一個實數 array，mxCOMPLEX 則為複數 array。 mxCreateDoubleMatrix 會將指向這個 mxArray 的指標回傳， 到此一個 array 算是完全建立了。 這個新建立的 array 內的元素都是 0，所以我們不用再初始化這個 array。

李易霖 (2004/08/04) --- [Prev] [Next] [Up]