在使用MATLAB的过程中，我对MATLAB的运行效率感到很头疼，就尝试了一些办法去提高之。现在把它们在这里作个总结，留作备忘和分享，之后有了新的想法也会补充进来。
- 使用矩阵运算替换循环语句
- CPU并行计算
- GPU并行计算
- 与C++协作
使用矩阵运算替换循环语句
这应该是老生常谈了；由于MATLAB处理矩阵挺方便的，所以一般也没人故意把矩阵运算拆成分量写。
有时候，for循环转成矩阵运算需要稍动一下脑筋。比如下面的例子中，plot_x是一个长度为size2的列向量，我们需要把n个plot_x连接起来，组成一个长长的列向量。
plot_x2 = zeros(size2 * n, 1);
for i = 1:n
plot_x2(((i - 1) * n + 1):(i * n), 1) = plot_x;
end
plot_x = plot_x2;
可以用这样的写法来替代：
plot_x = reshape(ones(n, 1) * (plot_x)', [size2 * n, 1]);
不知道这样效率能提高多少；就算不提高太多，看起来也更加令人心情愉悦。
不只是加减乘除，调用函数的时候也最好一次把所有的数据用矩阵喂给函数，而不是用for循环分开许多次调用。
例如，下面的例子中，plot_y是一个n行size2列的矩阵，plot_x是一个长度为size2的列向量，指示plot_y每一列画图时对应的横坐标。
for i = 1:n
plot(plot_x', plot_y(i, :));
end
可以修改为这样：
plot_x = ones(n, 1) * plot_x';
plot(plot_x, plot_y);
如果没记错的话，当时n大约取100，修改后执行时间从1分钟缩短到1秒。
有时候，从for循环到矩阵运算的转变需要用到arrayfun函数。需要先写一个用于处理单个元素（而不是矩阵）的函数，然后用arrayfun调用它。
下面的例子中，plot_x和plot_y的意义与上例相同，broken2是长度为size2的列向量，指示之前的程序在计算plot_y的每一列（也就是plot_x的每一行）时有没有出错（出错的点就不需要画出来了）。
plot_x = ones(n, 1) * plot_x';
for i = 1:size2
if ~broken(i)
plot(plot_x(:, i), plot_y(:, i));
end
end
可以改成这样：
plot_x = ones(n, 1) * (plot_x + arrayfun(@isbroken, broken))';
plot(plot_x, plot_y);
其中，isbroken是自己写的一个函数，用来处理broken的每一个元素，其定义为：
function a = isbroken(broken)
if broken
a=NaN;
else
a=0;
end
end
注意这里isbroken的形参broken所指并不是原来的broken，而是它的一个个分量。
把自己写的函数传给arrayfun时，必须在前面加上@，否则会报错。在很多其它场合是可加可不加的，比如使用ode45积分时，将微分方程传给ode45可以不加。
对于数据处理非常复杂的场合（比如积分），只要对矩阵中每一个元素的处理都是独立（互不影响）的，使用arrayfun改写，亲测都可以获得相当高的效率提升。
传给arrayfun的各个矩阵，要求各个维度的大小全部相同，arrayfun会把同一行、列位置上的元素传给同一个自定义函数；或者，一部分矩阵的各个维度大小相同，另一部分矩阵只有一个元素。
CPU并行计算
这个特别简单。
在命令窗口执行：
parpool(n)
这样，就在你的电脑上开了n线程的线程池，等会儿就可以使用parfor用这n个线程并行计算。我的CPU是四核八线程，我一般取n为4或者6。
parfor的用法和for差不多，只需要注意不同次循环之间的数据不能互相影响。举个例子：
parfor i = 1:s
c(i,:,:)=pcrpoint(sgm, r, b(i), n, ts, tu, t0, rand(3,1), apr, tol);
end
这样4线程或者6线程调用自己写的pcrpoint去计算了。
parpool和parfor还有一些其它功能，我没用过，这里不提。
GPU并行计算
MATLAB是原生支持GPU运算的，并且不需要编程者对GPU的结构、CUDA编程模型非常了解。使用gpuDevice命令可以看到自己GPU的一些参数。
MATLAB对GPU的支持有以下的限制：
- 必须是NVIDIA的显卡（支持CUDA），AMD的不行。
- 要求Compute Capability为2.0或更高。这个一般都能满足。
- 需要安装CUDA Toolkits，这个到NVIDIA官网下载，按照默认的选项安装就好了。
- 需要编程者有足够的耐心，这个不是在开玩笑。改写一些复杂的函数时，会相当相当麻烦。另外，GPU的整个环境特别容易出问题，更新一下显卡驱动都可能让整个环境崩溃，出问题后要有耐心地卸载掉重装CUDA Toolkits。
要使用GPU计算，需要先把数据放到显存上。如果已经有了一个矩阵a，可以这样把数据复制到显存上：
dev_a = gpuDevice(a);
或者直接在显存上创建矩阵：
dev_a = zeros(5, 5, 'gpuArray');
ones，rands等函数也可以类似地使用。
接下来的计算就是自然而然地发生在GPU上了。比如，可以计算：
dev_c = dev_b + dev_a;
dev_c = sin(dev_a);
加减乘除和常见的内置函数都可以用。具体哪些可用而哪些不可用，这里有详细的列表。
如果需要将显存上的数据复制到CPU内存中，可以使用gather函数。
除了使用内置的函数以外，还可以自己写一个处理单个元素的函数，然后使用arrayfun调用。要命的是，在自己写的那个函数中，不能有任何矩阵出现，必须全部是对单个元素的操作。
这是我写的一个例子：
[x(:, 1), y(:, 1), z(:, 1)] = arrayfun(@ode45_lrz_onlyResult_ew, t0_span, x0, y0, z0, rt_span, at_span, sgm_span, r_span, b_span);
其中，t0_span等参数已经被提前置于显存上。ode45_lrz_onlyResult_ew是我写的一个函数，对每个矩阵中的单个元素进行积分。它是这样定义的：
function [yout_x, yout_y, yout_z] = ode45_lrz_onlyResult_ew(tfinal, y0x, y0y, y0z, RelTol, AbsTol, sgm, r, b)
...
end
Windows不允许GPU有两秒以上的时间无响应；所以，如果你交给GPU的计算任务太重，MATLAB会报错。在注册表中修改以下键值为0就可以取消这个限制：
HKEY_LOCAL_MACHINE\SYSTEM\CurrentControlSet\Control\GraphicsDrivers\TdrLevel
如果没有这个键值，就新建一个。这样做的副作用是，GPU全力计算的时候电脑可能会卡住，等算完就好过来了。
到这一步，效率已经非常高了。我测试的结果是，同样的任务放到GPU上计算，MATLAB的效率可以达到C++的四分之一，并且相较于用C++写CUDA，学习成本低。
与C++协作
MATLAB支持与C++协作。我尝试过的有两种方式：调用C++写好的exe，或者调用C++写好的lib。前者兼容性好一些，后者传参方便一些。
举一个调用exe的例子。MATLAB端把param写入到'm2c.bin'里，然后启动exe并等待结束，再从'c2m.bin'里读返回的yout；
fclose('all');
fid = fopen('m2c.bin', 'wb');
fwrite(fid, param', 'double');
fclose(fid);
system('ode_solver_launcher.exe');
fid = fopen('c2m.bin', 'rb');
yout = fread(fid, [3, size2 * n], 'double');
C++端从"m2c.bin"里读取param，计算，然后将计算结果yout写入到"c2m.bin"里：
FILE *fin = fopen("m2c.bin", "rb"), *fout = fopen("c2m.bin", "wb");
fread(param, sizeof(double), 3 * size, fin);
...
fwrite(yout, sizeof(double), 3 * n * size, fout);
fflush(fout);
fcloseall();
这里使用二进制文件来传输数据。如果数据量小，也可以用文本文件来传输，MATLAB用fprintf、fscanf，C++用printf、scanf即可。不赘述。
再举一个调用dll的例子。
C++端写dll可以参考这里。需要注意的有几点：
- 按照教程的方法，最新版本vs2017创建项目时不应该选择“Win32控制台应用程序”而是“Windows Desktop Wizard”（Windows桌面应用程序向导），这样才会弹出下一步的框。
- dll给MATLAB的接口必须是C（而不是C++）格式的。也就是说，在这个接口处，不能使用C++有而C没有的东西（比如
bool），并且在源文件中引用申明导出函数的头文件时，需要告诉编译器按C的方式编译：_EXTERN_C #include "lrz_pcrgraph_gpu_c_dllept.h" _END_EXTERN_C_EXTERN_C和_END_EXTERN_C其实就是extern "C" {和}。 - 参数传递时不要使用高维数组（尽管MATLAB原则上支持高维数组）。我在尝试中发现高维数组会有问题。
MATLAB端则需要将param转为libpointer对象，并创建适当的libpointer对象yout_x容纳返回值，使用头文件'lrz_pcrgraph_dllept.h'加载'lrz_pcrgraph.dll'，调用C++函数'pcrgraph'计算，卸载dll，将yout转化为矩阵，再进行下一步的处理。
yout_x = libpointer('doublePtr', zeros(s * n, 1));
param = libpointer('doublePtr', param);
if ~libisloaded('lrz_pcrgraph')
loadlibrary('lrz_pcrgraph.dll', 'lrz_pcrgraph_dllept.h');
end
calllib('lrz_pcrgraph', 'pcrgraph', yout_x, param);
unloadlibrary('lrz_pcrgraph');
yout_x = yout_x.value;
在上例中，C++中的double*对应MATLAB的doublePtr。具体的对应在这里可以看到。需要指出的是，实际测试时发现char*到int8Ptr或charPtr的对应有问题，可以使用unsigned char*到uint8Ptr来代替。
即使是调用已经编译好的dll，也需要装有C语言编译器，否则MATLAB会报错。
如果需要把程序放到别人电脑上计算，无论哪种方法，都需要注意：C++写的部分可能需要依赖别的dll才能运行。可以使用dumpbin工具检查exe或者dll依赖哪些其它的dll，然后把它们和刚刚写的程序放到同一个目录下；或者在编译前将Runtime Library编译选项改为/MT（默认为/MD)，就可以把用到的别的dll中的函数在编译时就收纳到自己的程序里了。
网友评论