Matlab混合编程
Matlab混合编程
混合编程目的
在Matlab中采用混合编程目的主要包括
- 利用已有的函数库，避免重复工作
- 加速计算，特别是减少循环所用时间
- 利用GPU等进行异构编程
混合编程方法—mex函数
目前已有的方法包括两种：（1）将c/Fortran源程序改写为mex函数，然后编译为二进制文件进行调用；（2）将c/Fortran源程序编译为动态链接库然后进行调用。在这两种方法中，前者计算速度更快，更适合混合编程方法。
mex函数是一类特殊的c/Fortran函数接口，编译后的二进制文件可以被matlab直接调用。
1.1.mex函数声明
首先看下mex函数声明格式
void mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
其中包含四个参数，其中nlhs与nrhs分别是输出参数与输入参数个数（Num of arguments on Left/Right Hand Side），plhs与prhs则是输入参数指针数组，mxArray为 Matlab 中数据结构体，其声明包含在头文件mex.h内。
1.2.mex函数参数调用
传入mex函数的参数都保存在mxArray结构体数组内，在使用参数进行计算时，需首先将其转换为C语言对应类型。
常用几个获得参数函数包括：
double mxGetScalar(const mxArray *pm);
double *mxGetPr(const mxArray *pm);
其中pm为对应mxArray类型的输入参数指针，如prhs[0]等。
- mxGetScalar
获取标量参数； - mxGetPr
获取向量/矩阵参数；
在使用mxGetPr或mxGetScalar函数获得输入参数对应的double类型指针之后，我们便可以直接对参数的值进行相应操作。需要注意的是，在Matlab中数组是按照列优先进行排列的，即是说当我们采用double类型指针*p指向一个[m x n]大小的矩阵时，p(i,j)对应的元素应为
p[(i-1)+(j-1)*m]
另外，在C函数中一般需要将矩阵维度也作为参数传入，但是在mex函数中，结构体mxArray则包含了相应的矩阵维度信息，采用函数
size_t mxGetM(const mxArray *pm);
size_t mxGetN(const mxArray *pm);
可以分别获得矩阵行数与列数，减少参数个数。
Matlab中自定义类型mxArray实际是一个包含多个矩阵信息的结构体，函数mxGetM，mxGetN和mxGetPr等分别返回结构体包含的对应元素。
1.3.mex函数返回参数
mex返回参数个数一般是确定的，可以用如下语句进行判断
if (nlhs != 2)
mexErrMsgTxt("Wrong number of output arguments");
在对输出变量元素进行操作之前，需要首先为其申请内存
mxArray *mxCreateDoubleMatrix(mwSize m, mwSize n, mxComplexity ComplexFlag);
其中m和n为矩阵大小，ComplexFlag为元素类型，包括mxREAL与mxCOMPLEX，前者代表申请大小为[m x n]的实数矩阵，后者申请相同大小的复数矩阵。
在为输出参数申请内存完毕后，即可用mxGetPr函数来获取输出参数对应double类型指针，然后进行赋值操作。
示例
目标：使用混合编译的方法，在Matlab中调用 Polylib 函数库
1.采用mex函数
使用c添加mex接口函数
#include "mex.h"
void mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double *p_c, *p_d;
double *p_a, *p_b;
int c_rows, c_cols;
int d_rows, d_cols;
int numEl;
int n;
mxAssert(nlhs==2 && nrhs==2, "Error: number of variables");
c_rows = mxGetM(prhs[0]);// get rows of c
c_cols = mxGetN(prhs[0]);// get cols of c
d_rows = mxGetM(prhs[1]);// get rows of d
d_cols = mxGetN(prhs[1]);// get cols of d
mxAssert(c_rows==d_rows && c_cols==d_cols, "Error: cols and rows");
// create output buffer
plhs[0] = mxCreateDoubleMatrix(c_rows, c_cols, mxREAL);
plhs[1] = mxCreateDoubleMatrix(c_rows, c_cols, mxREAL);
// get buffer pointers
p_a = (double*)mxGetData(plhs[0]);
p_b = (double*)mxGetData(plhs[1]);
p_c = (double*)mxGetData(prhs[0]);
p_d = (double*)mxGetData(prhs[1]);
// compute a = c + d; b = c - d;
numEl = c_rows*c_cols;
for (n = 0; n < numEl; n++)
{
p_a[n] = p_c[n] + p_d[n];
p_b[n] = p_c[n] - p_d[n];
}
}
2.采用动态链接库
使用gcc编译为动态链接库，使用Matlab写接口对应的函数
CC=gcc-4.6
target:
${CC} -c polylib.c
${CC} -shared -fPIC -o libpolylib.dyld polylib.o
在Matlab中查看动态库包含指令
>> loadlibrary('libpolylib.dyld');
>> list = libfunctions('libpolylib','-full')
list =
'[doublePtr, doublePtr, doublePtr] Dgj(doublePtr, doublePtr, doublePtr, int32, double, double)'
'[doublePtr, doublePtr, doublePtr] Dglj(doublePtr, doublePtr, doublePtr, int32, double, double)'
'[doublePtr, doublePtr, doublePtr] Dgrjm(doublePtr, doublePtr, doublePtr, int32, double, double)'
'[doublePtr, doublePtr, doublePtr] Dgrjp(doublePtr, doublePtr, doublePtr, int32, double, double)'
'[doublePtr, doublePtr, doublePtr] Imgj(doublePtr, doublePtr, doublePtr, int32, int32, double, double)'
'[doublePtr, doublePtr, doublePtr] Imglj(doublePtr, doublePtr, doublePtr, int32, int32, double, double)'
'[doublePtr, doublePtr, doublePtr] Imgrjm(doublePtr, doublePtr, doublePtr, int32, int32, double, double)'
'[doublePtr, doublePtr, doublePtr] Imgrjp(doublePtr, doublePtr, doublePtr, int32, int32, double, double)'
'[double, doublePtr] hgj(int32, double, doublePtr, int32, double, double)'
'[double, doublePtr] hglj(int32, double, doublePtr, int32, double, double)'
'[double, doublePtr] hgrjm(int32, double, doublePtr, int32, double, double)'
'[double, doublePtr] hgrjp(int32, double, doublePtr, int32, double, double)'
'[doublePtr, doublePtr] jacobd(int32, doublePtr, doublePtr, int32, double, double)'
'[doublePtr, doublePtr, doublePtr] jacobfd(int32, doublePtr, doublePtr, doublePtr, int32, double, double)'
'[doublePtr, doublePtr] zwgj(doublePtr, doublePtr, int32, double, double)'
'[doublePtr, doublePtr] zwglj(doublePtr, doublePtr, int32, double, double)'
'[doublePtr, doublePtr] zwgrjm(doublePtr, doublePtr, int32, double, double)'
'[doublePtr, doublePtr] zwgrjp(doublePtr, doublePtr, int32, double, double)'
想要调用的是zwgj函数，输入参数类型为(doublePtr, doublePtr, int32, double, double)，注意Matlab中对应变量类型
| C Type | Equivalent Matlab Type |
|---|---|
| char,byte | int8 |
| unsigned char,byte | uint8 |
| short | int16 |
| unsigned short | uint16 |
| int | int32 |
| long(Windows®) | int32,long |
| long(Linux®) | int64,long |
| unsigned int | unit32 |
| unsigned long(Windows) | uint32,long |
| unsigned long (Linux) | uint64,long |
| float | single |
| double | double |
| char * | 1xn char array |
| *char[] | cell array of strings |
其中libfunctions函数显示的变量类型对应为
| C Pointer Type | Argument Data Type | Equivalent Matlab Type |
|---|---|---|
| double * | doublePtr | double |
| float * | singlePtr | single |
| integer pointer types (int *) | (u)int(size)Ptr | (u)int(size) |
| Matrix of signed bytes | int8Ptr | int8 |
| Null-terminated string passed by value | cstring | 1xn char array |
| Array of pointers to strings (or one **char) | stringPtrPtr | cell array of strings |
| enum | enumPtr | |
| type ** | Same as typePtr with an added Ptr (for example, double ** is doublePtrPtr) | lib.pointer object |
| void * | voidPtr | |
| void ** | voidPtrPtr | lib.pointer object |
| C-style structure | structure | MATLAB struct |
| mxArray * | MATLAB array | MATLAB array |
| mxArray ** | MATLAB arrayPtr | lib.pointer object |
对应的Matlab接口函数为
loadlibrary('libpolylib.dyld');
np = int32(5);
z = zeros(1,np); w = zeros(1,np);
[z, w] = calllib('libpolylib', 'zwgj', z, w,np,alpha,beta);
其中np转换为对应的int32类型变量
