Matlab混合编程

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

1. 利用已有的函数库，避免重复工作
2. 加速计算，特别是减少循环所用时间
3. 利用GPU等进行异构编程

目前已有的方法包括两种：（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类型变量