Matlab混合编程

Matlab混合编程

混合编程目的

在Matlab中采用混合编程目的主要包括

1. 利用已有的函数库，避免重复工作
2. 加速计算，特别是减少循环所用时间
3. 利用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]等。

1. mxGetScalar
   获取标量参数；
2. 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类型变量