前言:
众所周知对于循环这种运算,matlab的运算速度是不快的
想起来个冷笑话,黑客帝国里的主角呢奥之所以能快过子弹,因为Matrix是matlab写成的hhhh
因此将循环部分用C语言写,然后编译成matlab可以调用的mex文件,无疑可以加快循环运算速度,从而加快整个代码速度
之前因为合并数据需要不会数据库,强行在matlab里利用类似excel的vlookup功能,无奈太慢,只能想办法最后知道了这种方法,最近偶然在桥哥的知识星球里又提到了,也顺便复习和记录一下。
ps:今天在闲鱼上买了小车车,高兴!
1.环境配置
我的是matlab2019a,win10-64bit,c语言编译器我选的是TDM-GCC(gcc/g++),安装起来很简单,推荐matlab2015及以后的选择这种,详细安装过程可以见文末链接1中的For Matlab 2015部分
2.如何编写可以被编译成matlab可执行文件的c语言代码
我是在matlab中文论坛(地址见reference第2条)里找到如何编写的,不过他写的有点乱,我来消化后写一下自己的理解
mex接口函数跟一般的C语言文件主要是两点区别,
第一是必须引入头文件mex.h
#include "mex.h"
第二是多了一个叫mexFunction的函数
void mexFunction(int nlhs,mxArray *plhs[], int nrhs,const mxArray *prhs[])
举个栗子:
#include "mex.h" // 使用MEX文件必须包含的头文件
// 执行具体工作的C函数
double add(double x, double y)
{
return x + y;
}
// MEX文件接口函数
void mexFunction(int nlhs,mxArray *plhs[], int nrhs,const mxArray *prhs[])
{
double *z;
double x, y;
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
z = mxGetPr(plhs[0]);
x = *(mxGetPr(prhs[0]));
y = *(mxGetPr(prhs[1]));
*z = add(x, y);
}
也可以写成:
#include "mex.h" // 使用MEX文件必须包含的头文件
// MEX文件接口函数
void mexFunction(int nlhs,mxArray *plhs[], int nrhs,const mxArray *prhs[])
{
double *z;
double x, y;
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
z = mxGetPr(plhs[0]);
x = *(mxGetPr(prhs[0]));
y = *(mxGetPr(prhs[1]));
*z=x+y;
}
也就是说执行具体功能可以通过编写一个独立的c语言函数在接口函数里调用,也可以直接写在mexFunction接口函数里
我个人推荐将具体功能独立出来,这样一来:
c语言函数负责执行具体功能
mexFunction接口函数负责数据的输入输出(将matlab数据输入c语言环境,运算后再输出回matlab)形式组织
再调用c语言函数
分工明确。
于是关键问题就是如何编写mexFunction接口函数了
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
下面解释一下mexFunction各个参数的含义(这里看到是void,即无返回值,因为传值是通过指针数组plhs[]传出的)
参数 | 意义 | 英文全称 | 类型 |
---|---|---|---|
nlhs | 左边输出参数的数目 | number of left-hand side | 整型(int) |
plhs | 指向输出参数的指针 | pointer of left-hand side | 指针数组 |
nrhs | 右边输入参数的数目 | number of right-hand side | 整型(int) |
prhs | 指向输入参数的指针 | pointer of right-hand side | 指针数组 |
第一个参数nlhs是输出参数的数目,第二个参数plhs是指向输出参数的指针
第三个参数nrhs是输入参数的数目,第四个参数prhs是指向输入参数的指针
其中plhs和prhs类型都是指向mxArray类型数据的指针,prhs这里还加了个const,这里的const跟之前c语言中是一样的,代表不改变,因为prhs是指向输入参数的指针,mxArray这个数据类型是在头文件mex.h中定义的,in fact ,在matlab中大多数数据都是以这种类型存在。
还是拿之前的栗子来解释:
#include "mex.h" // 使用MEX文件必须包含的头文件
// MEX文件接口函数
void mexFunction(int nlhs,mxArray *plhs[], int nrhs,const mxArray *prhs[])
{
double *z;
double x, y;
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
z = mxGetPr(plhs[0]);
x = *(mxGetPr(prhs[0]));
y = *(mxGetPr(prhs[1]));
*z=x+y;
}
在这个栗子中,输入参数是x,y,输入参数是z(虽然没有返回)
plhs和prhs其实都是指针数组
prhs[0]代表指向第一个输入数据的指针,mxGetPr代表获得这个指针,后面再加个*就是c语言里面再正常不过的按指针取内容了
所以就这样x,y通过下面这两句传入到了mexFunction中
x = *(mxGetPr(prhs[0]));
y = *(mxGetPr(prhs[1]));
mxGetPr函数的功能是从指向mxArray类型数据的指针prhs[0]、prhs[1]中获得了指向double类型的指针
另一个函数是mxGetScalar,Scalar即标量,功能是把通过prhs传递进来的mxArray类型数据的指针所指向的数据(标量)赋给c程序里的变量。前面的mxGetPr是传递矢量的,否则矩阵或者向量就没有办法传进来
这里由于传入的每个参数是单个值,于是也可以用mxGetScalar改写:
#include "mex.h" // 使用MEX文件必须包含的头文件
// MEX文件接口函数
void mexFunction(int nlhs,mxArray *plhs[], int nrhs,const mxArray *prhs[])
{
double *z;
double x, y;
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
z = mxGetPr(plhs[0]);
x = mxGetScalar(prhs[0]));
y = mxGetScalar(prhs[1]));
*z=x+y;
}
(注意:输出参数因为需要通过指针传出,因此这里z必须用mxGetPr这个函数来从mxArray指针获取double指针)
但是这样还是有个问题:如果输入的不是单个的数据,是向量或者矩阵,即使我们通过mxGetPr获得了矩阵或者向量的指针,但是我们假如不知道矩阵的shape,还是不好取值或者计算,所以这里又有两个函数mxGetM和mxGetN可以通过传入参数指针获得传入参数(矩阵)的行和列数。
需要注意的是:
在matlab中矩阵的第一行和第一列的序号是从1开始的,而在c语言中是从0开始的,而且在c语言中的一维数组存储矩阵的数据时,是一列列从上到下,从左到右来存储的,也即
matlab中一个m*n的矩阵M(i,j)对应于c语言中一维数组N中的N[j*m+i]
举个栗子:
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
double *data;
int M,N;
int i,j;
data=mxGetPr(prhs[0]); //获得指向矩阵的指针
M=mxGetM(prhs[0]); //获得矩阵的行数
N=mxGetN(prhs[0]); //获得矩阵的列数
for(i=0;i<M;i++)
{
for(j=0;j<N;j++)
mexPrintf("%4.3f ",data[j*M+i]);
mexPrintf("\n");
}
}
假如是一个形如[1,2,3;4,5,6]的矩阵,则执行上述后会先后打印出1,4,2,5,3,6
以上讲的都是输入参数,由于输入数据在函数调用前就已经在matlab里申请过内存,由于mex函数与matlab共用一个地址空间,因此通过prhs传递指针即可传入输入参数。但输出参数却需要在mex函数里申请内存,才能将指针放在plhs中传递出去。由于返回指针类似必须是mxArray,所以matlab专门提供了一个函数:mxCreateDoubleMatrix来实现内存申请,函数原型为:
mxArray * mxCreateDoubleMatrix(int m,int n,mxComplexity ComplexFlag)
m,n分别是待申请矩阵的行数和列数,需要注意的是为矩阵申请内存后得到的是mxArray类型的指针,就可以放在plhs[]中传递出去了。但是对这个矩阵的处理(包括赋值)却需要在mex函数或者c语言函数中完成,这就需要通过前面的mxGetPr或者mxGetScalar。使用mxGetPr获得指向这个矩阵的double类型指针后,就可以对这个矩阵进行各种操作和运算了。
还是拿这个栗子讲,这里由于输出是一个数,于是m,n都是1,通过mxGetPr进一步获得指向double数据类型指针z
double *z;
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
z = mxGetPr(plhs[0]);
x = *(mxGetPr(prhs[0]));
y = *(mxGetPr(prhs[1]));
*z=x+y;
当然,matlab里使用到的并不只是double类型这一种矩阵,还有字符串类型,结构类型矩阵等,并也提供了对应的处理函数。
references
(本文中的栗子来源于两文,部分知识讲述摘于文2)