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
类型变量