1. 引言
1.1 起因
这周去见应用数学的毕设老师Z，闲聊时Z说另一个同学代码里某个for循环要跑几小时，但这个for循环里操作并不复杂，主要就是转换和赋值。然后我就想，可不可以使用C/C++加速呢？
作为坚定的C/C++党，正好也想学学python中使用C/C++，C/C++中使用asm汇编这些技术，索性在github上建了一个新仓库Extend-Embed，专门学习怎样在一种语言编写的代码中使用另一种语言。
这篇文章记录怎样通过MEX方式，在MATLAB中使用C/C++，最主要的参考资料是：Writing MATLAB C/MEX Code
1.2 环境
操作系统：Windows 7
MATLAB：R2013a
Visual Studio：Community 2013
2. 准备
2.1 基本概念
简单粗暴地讲，在MATLAB中使用C/C++就是三步走的流程：
- 编写C/C++文件
- 用mex命令把C/C++文件编译成MEX(Matlab Executable)文件
- 像普通MATLAB函数一样直接使用
举个例子，我们写了一个hello.c，功能是输出一句”Hello, world!”，那么使用：
mex hello.c
后根据平台不同会得到hello.mex或者hello.mexw64，然后使用：
hello
就会输出一句”Hello， world！“
mex版本的Hello World
2.2 配置MEX
不得不说配置过程还是挺坑的。自己试验加上网上查资料试了无数遍，终于发现配置过程要经常反其道而行之，也就是说按正常思路配置基本都要GG。最终选了一堆no反而配好了环境时：
我跟你讲我就是这个表情
输入mex -setup，会询问是否让mex自己找编译器，
注意选择y是不行的：
找不到编译器
这时要选择n，得到各种支持的编译器名称：
各种支持的编译器
我安装了Visual Studio Community 2013，所以选择21。因为它给的默认路径不对，所以选择n后输入正确路径：
输入正确的vs路径
输入路径后选择y确认，输出信息里有Done就是配置成功了：
配置成功
3. 混合编程
在MATLAB安装目录下有一个extern\examples\mex文件夹，里面有很多例子；也可以参考Writing MATLAB C/MEX Code。不过必须知道的新概念也就几条，其它都可以边用边学，所以实际上手非常容易。
3.1 新概念
头文件方面，使用MEX文件时都需要include “mex.h”。
main函数方面，入口不再是main函数，而是mexFunction，声明格式如下：
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
其中mxArray是mex.h里自己定义的数据结构，直接看四个奇怪变量名的含义：
lhs = left hand side，rhs = right hand side
nlhs = 输出变量的数量
plhs = 指向输出变量的指针数组
nrhs = 输入变量的数量
prhs = 指向输入变量的指针数组（注意const）
用MEX方式实现了函数test，在MATLAB里面使用：
F = test(X, Y)
那么输入方面，nrhs为2，prhs[0]是X，prhs[1]是Y；输出方面，nlhs为1，plhs[0]是F。
需要注意的是两件事：第一，要自己分配输出的存储空间；第二，要自己检查输入是否合法。这两件事直接讲不好讲清楚，看后面的例子就很容易理解了。
此外，mex.h里提供了不少辅助函数，输出信息用mexPrintf、mexErrMsgTxt，得到矩阵尺寸用mxGetM和mxGetN，得到指向矩阵数据的指针用mxGetPr，etc
3.2 例子
写了一个discrete_cpp.cpp，功能是传入一个矩阵A返回一个矩阵B，B和A的形状相同，B中元素B(i, j)和A中元素A(i, j)的对应关系是：
- 若 A(i, j) 小于1/3，则 B(i, j) = 0
- 若 A(i, j) 大于等于1/3 且小于2/3，则 B(i, j) = 1
- 若 A(i, j) 大于2/3，则 B(i, j) = 2
类似于做了离散化。
#include "mex.h"
// 使用MEX必须包含的头文件
int discrete(double d);
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[]){
if (nrhs != 1)
mexErrMsgTxt("Wrong number of input arguments.\n");
// 检查输入变量数量是否正确，否则报错
if (nlhs > 1)
mexErrMsgTxt("Too many output argumnents.\n");
// 检查输出变量数量是否正确，否则报错
#define A_IN prhs[0]
#define B_OUT plhs[0]
int M = mxGetM(A_IN);
int N = mxGetN(A_IN);
// 得到输入矩阵A的行数和列数
B_OUT = mxCreateDoubleMatrix(M, N, mxREAL);
// 为输出矩阵B分配存储空间
double *A = mxGetPr(A_IN);
double *B = mxGetPr(B_OUT);
// 取得输入矩阵A和输出矩阵B的数据指针
for (int i = 0; i < M * N; ++i)
B[i] = discrete(A[i]);
// 调用discrete，根据A(i, j)计算B(i, j)
}
int discrete(double d) {
if (d < 1.0 / 3.0)
return 0;
else if (d < 2.0 / 3.0)
return 1;
return 2;
}
4. 测试
编写discrete的MATLAB版本discrete_matlab.m：
function [ B ] = discrete_matlab( A )
[M, N] = size(A);
B = zeros(M, N);
for i = 1:M
for j = 1:N
if A(i, j) < 1/3
B(i, j) = 0;
elseif A(i, j) < 2/3
B(i, j) = 1;
else
B(i, j) = 2;
end
end
end
end
编写测试脚本test_discrete：
MATRIX_M = 100;
MATRIX_N = 100;
MATRIX_NUMBER = 1000;
data = rand(MATRIX_M, MATRIX_N, MATRIX_NUMBER);
tic
for i = 1:MATRIX_NUMBER
discrete_cpp(data(:, :, i));
end
toc
tic
for i = 1:MATRIX_NUMBER
discrete_matlab(data(:, :, i));
end
toc
也就是生成MATRIX_NUMBER个MATRIX_M * MATRIX_N的随机矩阵，然后分别用mex版本和普通matlab函数做discrete，看各自消耗的时间。
测试结果……第一行是MEX版本的时间，第二行是MATLAB版本的时间，可以看出MEX方式在数据量大时还是挺有优势的。
M = 4，N = 4，NUMBER = 1000
M = 20，N = 20，NUMBER = 1000
M = 50，N = 50，NUMBER = 1000
M = 100，N = 100，NUMBER = 1000
所以这图还真是各种生动形象：
the spooky beast
5. 参考资料
- Writing MATLAB C-MEX Code：
http://www.docin.com/p-684614383.html - 解决”mex setup或mbuild setup找不到编译器”问题
http://blog.csdn.net/sanqima/article/details/49514889
