3-5 BzGx} MEX 仵袙d

本節將介紹兩個範例，以便對二維矩陣進行處理，

第一個範例是 arrayx2.c，可將輸入二維矩陣（可能包含實部和虛部）乘以 2，其內容如下：

原始檔（03-應用程式介面/arrayx2.c）：（灰色區域按兩下即可拷貝）

/* 此函式為 MATLAB 的 MEX 檔案，其輸入為一個向量，輸出則為2乘以此向量。 */

#include "mex.h"	/* mex.h 包含 mxArray 結構的定義，以及其他相關資訊 */

/* prhs = pointer to the right-hand-side arguments，即指向輸入變數列的指標 */
/* prls = pointer to the  left-hand-side arguments，即指向輸出變數列的指標 */
#define IN  prhs[0]	/* 定義 IN  為此函式的第一個輸入變數，以便後續取用 */
#define OUT plhs[0]	/* 定義 OUT 為此函式的第一個輸出變數，以便後續取用 */

/* 此函式的功能為將向量 x 的每一個元素乘以2，再將結果送到 y 向量。 */
/* 此函式將會被 mexFunction 所呼叫。 */
void timestwo(double y[], double x[], int length) {
	int i;
	for (i=0; i<length; i++)
		y[i] = 2.0*x[i];
}

/* 此函式為和 MATLAB 溝通的主要函式 */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] ) {
	double *xr, *xi, *yr, *yi;
	int length = mxGetM(IN)*mxGetN(IN);
  
	/* 檢查輸出和輸入變數個數是否都是1，其中		  */
	/* nrhs = no. of right-hand-side arguments（輸入變數個數）*/
	/* nlhs = no. of  left-hand-side arguments（輸出變數個數）*/
	if(nrhs!=1)	/* 檢查輸入變數個數是否是1 */
		mexErrMsgTxt("One input required.");
	if(nlhs>1)	/* 檢查輸出變數個數是否是1 */
		mexErrMsgTxt("Too many output arguments");
  
	/* 檢查輸入變數是否合格 */
	if(!mxIsDouble(IN))		/* 檢查輸入變數是否為 double */
		mexErrMsgTxt("Input must be a double.");
  
	/* 配置記憶體給輸出變數 */
	OUT = mxCreateDoubleMatrix(mxGetM(IN), mxGetN(IN), mxCOMPLEX);
  
	/* 取得輸入和輸出變數的實部指標 */
	xr = mxGetPr(IN);
	yr = mxGetPr(OUT);
	/* 取得輸入和輸出變數的虛部指標 */
	if (mxIsComplex(IN)) {
		xi = mxGetPi(IN);
		yi = mxGetPi(OUT);
	}
  
	/* 執行實際的運算：將輸入向量的實部乘以2 */
	timestwo(yr, xr, length);
	/* 執行實際的運算：將輸入向量的虛部乘以2 */
	if (mxIsComplex(IN))
		timestwo(yi, xi, length);
}

Hint

由於所有的陣列在 MATLAB 內部都表示成一個長向量（例如，對二為矩陣來說，就是一個由各個行向量所串接而成的長向量），因此在上述範例中，我們可由 length = mxGetM(IN)*mxGetN(IN) 來取得此長向量的長度。

接著，我們可以編譯以上程式：

>> mex arrayx2.c

欲確認可執行檔的存在，可輸入如下：

>> which arrayx2 D:\matlabBook\MATLAB程式設計：進階篇\03-應用程式介面\arrayx2.mexw64

欲進行測試，可輸入如下：

>> arrayx2([1 2 3]) ans = 2 4 6 >> arrayx2([1+i 2+2i 3+3i; 4 5 6]) ans = 2.0000 + 2.0000i 4.0000 + 4.0000i 6.0000 + 6.0000i 8.0000 10.0000 12.0000 >> arrayx2('Test string') ??? Input must be a double. 本節的第二個 MEX 範例是 pairdist.c（檔名為「pairwise distance」的簡稱），可接受兩個輸入矩陣 A 和 B，矩陣的大小分別是 p×m 及 p×n，此兩個矩陣分別代表由 m 個直行向量及 n 個直行向量所形成的集合，而每一個向量的長度都是 p。此 MEX 檔案傳回一個距離矩陣 C，其中 C(i, j)為向量 A(:, i) 和向量 B(:, j) 的矩離。

Hint

對於大量資料的處理，一般的慣例都是將各個向量以直行向量的方式堆成一個矩陣，再進行各種處理。

pairdist.c 的程式碼如下：

原始檔（03-應用程式介面/pairdist.c）：（灰色區域按兩下即可拷貝）

#include <math.h>
#include "mex.h"	/* mex.h 包含 mxArray 結構的定義，以及其他相關資訊 */

/* prhs = pointer to the right-hand-side arguments，即指向輸入變數列的指標 */
/* prls = pointer to the  left-hand-side arguments，即指向輸出變數列的指標 */
#define	MAT1 prhs[0]	/* 定義 MAT1 為此函數的第一個輸入變數，以便後續取用 */
#define	MAT2 prhs[1]	/* 定義 MAT1 為此函數的第二個輸入變數，以便後續取用 */
#define	OUT  plhs[0]	/* 定義 OUT  為此函數的第一個輸出變數，以便後續取用 */

/* 此函數為和 MATLAB 溝通的主要函數 */
void mexFunction(
	int nlhs,	mxArray *plhs[],
	int nrhs, const mxArray *prhs[]) {
	double	*out, *mat1, *mat2, squareSum;
	int m, p, n, p2, i, j, k;

	/* 檢查輸出變數個數是否為2、輸入變數個數是否為1，其中	  */
	/* nrhs = no. of right-hand-side arguments（輸入變數個數）*/
	/* nlhs = no. of  left-hand-side arguments（輸出變數個數）*/
	if (nrhs!=2)
		mexErrMsgTxt("PAIRDIST requires two input arguments.");

	/* 檢查維度是否符合要求 */
	p  = mxGetM(MAT1);
	m  = mxGetN(MAT1);
	p2  = mxGetM(MAT2);
	n = mxGetN(MAT2);
	if (p != p2)
		mexErrMsgTxt("Matrix sizes mismatch!");

	/* 檢查輸入變數資料型態是否符合要求 */
	if (!mxIsNumeric(MAT1) || mxIsSparse(MAT1)  || !mxIsDouble(MAT1))
		mexErrMsgTxt("Input 1 is not a full numerical array!");
	if (!mxIsNumeric(MAT2) || mxIsSparse(MAT2)  || !mxIsDouble(MAT2))
		mexErrMsgTxt("Input 2 is not a full numerical array!");

	/* 配置記憶體給輸出變數 */
	OUT = mxCreateDoubleMatrix(m, n, mxREAL);

	/* 取得輸入和輸出變數的實部指標 */
	out = mxGetPr(OUT);
	mat1 = mxGetPr(MAT1);
	mat2 = mxGetPr(MAT2);

	/* 實際計算部份 */
	/* MAT1(i+1, j+1) 的值是 mat1[j*m+i] */
	/* MAT2(i+1, j+1) 的值是 mat2[j*n+i] */
	/*  OUT(i+1, j+1) 的值是  out[j*m+i] */

	for (i=0; i<m; i++)
		for (j=0; j<n; j++) {
			squareSum = 0;
			for (k=0; k<p; k++)
				squareSum += pow(mat1[i*p+k]-mat2[j*p+k], 2.0);
			out[j*m+i] = sqrt(squareSum);
		}
}

Hint

在上述範例中，請特別注意矩陣的索引，若使用 0-based indexing，則 MATLAB 矩陣 OUT(i+1, j+1) 所對應的元素值是 out[j*m+i]。
因為距離矩陣是對稱的，因此若考慮計算效能，對角線以下的元素可以直接拷貝對角線以上的元素，而不必重算一遍。

接著，我們可以對以上程式進行編譯：

>> mex pairdist.c 編譯完此程式後，可測試如下： >> which pairdist D:\matlabBook\MATLAB程式設計：進階篇\03-應用程式介面\pairdist.mexw64 >> pairdist([2 3], [2 3 4]) ans = 0 1 2 1 0 1 >> A = rand(2, 3); >> B = rand(2, 4); >> C = pairdist(A, B) C = 0.5104 0.9806 0.4544 0.5649 0.3562 0.2347 0.5911 0.2628 0.2101 0.5891 0.2675 0.1639 >> pairdist('string1', 'string2') ??? Input 1 is not a full numerical array! 本節所介紹的範例，均是對二維矩陣進行處理，其它更複雜的範例，可見 MATLAB API 的相關手冊。

MATLAB程式設計：進階篇