matlab gpu程序效率,提高MATLAB程序的运行效率

在使用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中的函数在编译时就收纳到自己的程序里了。