Matlab混合编程

时间:2021-11-13 17:24:44

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[])

其中包含四个参数,其中nlhsnrhs分别是输出参数与输入参数个数(Num of arguments on Left/Right Hand Side),plhsprhs则是输入参数指针数组,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

    获取向量/矩阵参数;

在使用mxGetPrmxGetScalar函数获得输入参数对应的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实际是一个包含多个矩阵信息的结构体,函数mxGetMmxGetNmxGetPr等分别返回结构体包含的对应元素。


1.3.mex函数返回参数

mex返回参数个数一般是确定的,可以用如下语句进行判断

if (nlhs != 2)
mexErrMsgTxt("Wrong number of output arguments");

在对输出变量元素进行操作之前,需要首先为其申请内存

mxArray *mxCreateDoubleMatrix(mwSize m, mwSize n, mxComplexity ComplexFlag);

其中m和n为矩阵大小,ComplexFlag为元素类型,包括mxREALmxCOMPLEX,前者代表申请大小为[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类型变量