GPU&VS2012&CUDA&matlab&Arrayfire杂记（二）——matlab

Matlab&GPU&CUDA并行加速学习心得

1所需软件：

l  NVIDIA GPU一块或多块，允许不同型号的GPU共存（本机配置了华硕GTX960）；

l  NVIDIA开发驱动程序（安装板卡的时候会自动提示安装，没有提示安装用附带光盘）；

l  CUDA toolkit安装，这里选择的是cuda toolkit7.5安装，网址如下，里面包含了GPU计算软件开发工具包，并行Nsight调试器，默认安装就好了；https://developer.nvidia.com/cuda-toolkit；

l  CUDA SDK安装，好像NVIDIA已经不更新了，下载的V2.21版本，主要包含了很实用例子和函数，网址如下：http://www.nvidia.com.tw/object/cuda_get_tw_old.html

l  安装GPU-Z，主要用来查看GPU的运行状态；

l  VS2010\VS2012和Matlab2014a的下载安装这里就不多做介绍了；

2安装步骤

l  装CUDA之前电脑上已经有matlab2014a了，我首先安装了VS2010；

l  更新GPU驱动到最新版本；

l  安装CUDAtoolkit7.5，一路默认下去；

l  安装SDK，也是默认安装；

至此，完成了基本安装，下面就是一些环境变量配置了；

3配置环境变量

3.1基本变量

安装完成Toolkit 和SDK 后，已自动配置好系统环境变量。保险起见，手动

配置环境变量。在系统环境变量中新建如下项：

CUDA_SDK_PATH = C:\ProgramData\NVIDIACorporation\CUDA

Samples\v5.0\common

CUDA_PATH = C:\Program Files\NVIDIA GPUComputing

Toolkit\CUDA\v5.0

CUDA_LIB_PATH = %CUDA_PATH%\lib\Win32

CUDA_BIN_PATH = %CUDA_PATH%\bin

CUDA_SDK_LIB_PATH =%CUDA_SDK_PATH%\common\lib\Win32

CUDA_SDK_BIN_PATH =%CUDA_SDK_PATH%\bin\Win32

在系统环境变量 Path 后添加如下内容：

;%CUDA_LIB_PATH%;%CUDA_BIN_PATH%;%CUDA_SDK_LIB_PATH%;%CUD

A_SDK_BIN_PATH%;

3.2NVCC编译器变量设置（困扰了一段时间）

nvcc –ptx ××.cu配置（用来将.cu文件转换成ptx文件使得matlab可以在GPU的内核上并行执行程序），新建用户变量PATH，路径为C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v7.5\bin;D:\ProgramFiles (x86)\Microsoft Visual Studio 10.0\VC\bin，然后在dos下执行nvcc –ptx××.cu就ok了，生成的.ptx文件会保存在用户名下面，找不到可以在计算机中搜素一下就知道了。

3.3MEX-setup编译器变量设置（困扰了一段时间）

最后是重新装了VS2012搭配matlab2014解决了问题。

4CUDA C/C++关键字及函数高亮显示

4.1.cu文件中C\C++关键字高亮

这个设置是让VisualStudio2010 在编辑.cu 文件时，把.cu 文件里的C/C++语法高亮。设置方法：在Visual Studio 2010 的菜单依次选“ Tools|Options|Text Editor|File Extension（工具|选项|文本编辑器|文件扩展名） ” ， 在该窗口中将 “Editor （编辑器） ” 下拉框选择 “Microsoft Visual C++”，在“Extension（扩展名）”文本框中输入cu点击“Add（添加）”按钮，重复工作把cuh 添加为Visual C++类型，添加完成后点击“OK（确定）”按钮，如图 1所示。

图 1

重启VisualStudio 2010 后，.cu 文件C++关键字就高亮了。然而此时CUDA的关键字还是黑色的，下一步把CUDA 关键自高亮显示。

4.2CUDA关键字高亮设置

为了让CUDA 的关键字，如__device__、dim3 等的文字高亮，需按如下步骤设置：将C:\ProgramFiles\NVIDIACorporation\SDK\doc\syntax_highlighting\visual_studio_8目录下的usertype.dat（这个在下载的SDK文件中）文件复制到D:\ProgramFiles\Microsoft Visual Studio10.0\Common7\IDE\ 目录下 （ 对x64 位Win7 系统为X:\ProgramFiles(X86)\Microsoft Visual Studio 10.0\Common7\IDE\）。重启VisualStudio 2010 后打开.cu 文件，CUDA 的关键字应该变成蓝色了。

5创建工程

创建CUDA 工程在VisualStudio 2010 菜单选择“file|new|project（文件|新建|工程）”，在打开的新建项目窗口的“已安装的模板”一栏中选“NVIDIA|CUDA”，类型选择为“CUDA 7.5 Runtime”。在“名称”中输入工程名后，点击确定。可对系统提供的kernel.cu 示例进行编译运行。更多运行实例可找C:\Program Files\NVIDIA Corporation\Installer2\CUDASamples_7.5或者打开SDK；

6MATLAB&VS2010&CUDA的并行加速

在按照上述方法进行了正确的环境配置以后可以进行实例操作并查看运行时间。下面主要说明三种GPU的加速方法及实例http://cn.mathworks.com/help/distcomp/examples/illustrating-three-approaches-to-gpu-computing-the-mandelbrot-set.html

1.        使用现有的算法执行以GPU格式输入的数据（gpuArrary）；

2.        使用arrayfun函数独立对每个元素执行算法，arrayfun函数中执行的多个向量大小要一样；

3.        使用matlab/CUDA接口运行一些现有的CUDA/C++代码，将.cu文件生成ptx文件，然后通过matlab进行内核调用，方法较为简单，但是想要写出高效的内核程序必须对底层特别了解；

4.        GPU与C语言混合编译为MEX（动态链接库）进行加速执行。MEX-file可以通过C来调用NVIDIA写好的一些高级库，主要是写mex-function，包括fft等算法，熟悉c的话用这个方法能实现最快的计算。

5.        使用matlab并行工具箱，打开并行池（start parallel pool），调用了cpu的多线程进行计算，如果使用了for语句，改为 parfor进行并行计算

运行实例

6.1方法一（Using gpuArray）

% Setup

t = tic();

x = linspace( xlim(1), xlim(2),gridSize );

y = linspace( ylim(1), ylim(2),gridSize );

[xGrid,yGrid] = meshgrid( x, y );

z0 = xGrid + 1i*yGrid;

count = ones( size(z0) );

 

% Calculate

z = z0;

for n = 0:maxIterations

   z = z.*z + z0;

   inside = abs( z )<=2;

   count = count + inside;

end

count = log( count );

 

% Show

cpuTime = toc( t );

fig = gcf;

fig.Position = [200 200 600 600];

imagesc( x, y, count );

axis image

colormap( [jet();flipud( jet() );00 0] );

title( sprintf( '%1.2fsecs (without GPU)', cpuTime ) );

运行结果如图 2所示：

图 2CPU测试

6.2方法二（Element-wise Operation）

建立一个函数文件

function count = pctdemo_processMandelbrotElement(x0,y0,maxIterations)

z0 = complex(x0,y0);

z = z0;

count = 1;

while (count <= maxIterations) && (abs(z) <= 2)

    count = count + 1;

    z = z*z + z0;

end

count = log(count);

M文件

% Setup

t = tic();

x = gpuArray.linspace( xlim(1), xlim(2), gridSize );

y = gpuArray.linspace( ylim(1), ylim(2), gridSize );

[xGrid,yGrid] = meshgrid( x, y );

% Calculate

count = arrayfun( @pctdemo_processMandelbrotElement, ...

                  xGrid, yGrid, maxIterations );

% Show

count = gather( count ); % Fetch the data back from the GPU

gpuArrayfunTime = toc( t );

imagesc( x, y, count )

axis image

title( sprintf( '%1.3fsecs (GPU arrayfun) = %1.1fx faster', ...

gpuArrayfunTime, cpuTime/gpuArrayfunTime ) );

运行结果如图 3所示：

图 3GPU测试

6.3方法三（Working with CUDA）

CUDA\C++code

__device__

unsigned int doIterations( double const realPart0,

                           double const imagPart0,

                           unsigned int const maxIters ) {

   // Initialize: z = z0

   double realPart = realPart0;

   double imagPart = imagPart0;

   unsigned int count = 0;

   // Loop until escape

   while ( ( count <= maxIters )

          && ((realPart*realPart + imagPart*imagPart) <= 4.0) ) {

      ++count;

      // Update: z = z*z + z0;

      double const oldRealPart = realPart;

      realPart = realPart*realPart - imagPart*imagPart + realPart0;

      imagPart = 2.0*oldRealPart*imagPart + imagPart0;

   }

   return count;

}

Matlabcode

% Load the kernel

cudaFilename = 'pctdemo_processMandelbrotElement.cu';

ptxFilename = ['pctdemo_processMandelbrotElement.',parallel.gpu.ptxext];

kernel = parallel.gpu.CUDAKernel( ptxFilename, cudaFilename );

% Setup

t = tic();

x = gpuArray.linspace( xlim(1), xlim(2), gridSize );

y = gpuArray.linspace( ylim(1), ylim(2), gridSize );

[xGrid,yGrid] = meshgrid( x, y );

% Make sure we have sufficient blocks to cover all of the locations

numElements = numel( xGrid );

kernel.ThreadBlockSize = [kernel.MaxThreadsPerBlock,1,1];

kernel.GridSize = [ceil(numElements/kernel.MaxThreadsPerBlock),1];

% Call the kernel

count = zeros( size(xGrid), 'gpuArray' );

count = feval( kernel, count, xGrid, yGrid, maxIterations, numElements );

% Show

count = gather( count ); % Fetch the data back from the GPU

gpuCUDAKernelTime = toc( t );

imagesc( x, y, count )

axis image

title( sprintf( '%1.3fsecs (GPU CUDAKernel) = %1.1fx faster', ...

gpuCUDAKernelTime, cpuTime/gpuCUDAKernelTime ) );

运行结果如图 4所示：

图 4CUDA测试

7DOA估计测试

7.1CPU测试

代码：data =hilbert(data);

 

Gtheata=(theta);

Gdata=(data);

R1=(Gdata'*Gdata);

GR1=R1;

GE= (E);

%Gdata= (data);

for  i =1:length(Gtheata)

    a_s =exp(-jay*2*3.14*[0:7]*15000*0.05/1500*sin(Gtheata(i)*3.14/180));

    beam_single(i) =1/((a_s) *inv(GR1)*(a_s'));

end

运行结果如图 5所示：

图 5DOA（CPU）估计

7.2GPU测试（gpuArray）

代码：data =hilbert(data);

Gtheata= gpuArray(theta);

Gdata=gpuArray(data);

R1 =gpuArray(Gdata'*Gdata);

GR1=R1;

GE=gpuArray(E);

%Gdata= (data);

for  i =1:length(Gtheata)

    a_s =exp(-jay*2*3.14*[0:7]*15000*0.05/1500*sin(Gtheata(i)*3.14/180));

    beam_single(i) =1/((a_s) *(GE/GR1)*(a_s'));

end

 G =gather(beam_single);

运行结果如图 6所示：

图 6DOA（GPU）测试

7.3并行池测试（4核，parfor）

代码;

parfor  i =1:length(Gtheata)

    a_s =exp(-jay*2*3.14*[0:7]*15000*0.05/1500*sin(Gtheata(i)*3.14/180));

    beam_single(i) =1/((a_s) *inv(GR1)*(a_s'));

end

 

图 7DOA（并行池）测试

8结论

综上所述，matlab本身在计算简单的波束形成具有强大的优化功能，在这种基础上去调整并行或者GPU的运算反而不行，经过profile report分析出在GPU的运算中主要是

exp(-jay*2*3.14*[0:7]*15000*0.05/1500*sin(Gtheata(i)*3.14/180));

   beam_single(i) =1/((a_s) *(GE/GR1) *(a_s'));

这两句话占用了98%时间，使时间大大延长，分析有以下原因：

1.        没有给beam-single和data在GPU中预分配内存；

2.        这种简单的gpuArray用法较适合数值的迭代运算，对这种非迭代效果不好；

3.        改进方法：在数据的实时成像方法中，选用第三种方法（CUDA）较为可靠，在VS2010中编写GPU的内核程序，然后生成ptx文件或者mex成动态链接库在matlab中进行调用加速。