OpenCL向量加法

http://blog.csdn.net/lwfcgz/article/details/44056607

简介

下面一个例子介绍了向量加法的OpenCL版，相当于学习C语言中的“Hello World”，本篇教程中的代码以及其余相关教程都可以通过OLCF github下载

vecAdd.c

[cpp]  view plain  copy

 

1. #include <stdio.h>  
2. #include <stdlib.h>  
3. #include <math.h>  
4. #include <CL/opencl.h>  
5.    
6. // OpenCL kernel. Each work item takes care of one element of c  
7. const char *kernelSource =                                       "\n" \  
8. "#pragma OPENCL EXTENSION cl_khr_fp64 : enable                    \n" \  
9. "__kernel void vecAdd(  __global double *a,                       \n" \  
10. "                       __global double *b,                       \n" \  
11. "                       __global double *c,                       \n" \  
12. "                       const unsigned int n)                    \n" \  
13. "{                                                               \n" \  
14. "    //Get our global thread ID                                  \n" \  
15. "    int id = get_global_id(0);                                  \n" \  
16. "                                                                \n" \  
17. "    //Make sure we do not go out of bounds                      \n" \  
18. "    if (id < n)                                                 \n" \  
19. "        c[id] = a[id] + b[id];                                  \n" \  
20. "}                                                               \n" \  
21.                                                                 "\n" ;  
22.    
23. int main( int argc, char* argv[] )  
24. {  
25.     // Length of vectors  
26.     unsigned int n = 100000;  
27.    
28.     // Host input vectors  
29.     double *h_a;  
30.     double *h_b;  
31.     // Host output vector  
32.     double *h_c;  
33.    
34.     // Device input buffers  
35.     cl_mem d_a;  
36.     cl_mem d_b;  
37.     // Device output buffer  
38.     cl_mem d_c;  
39.    
40.     cl_platform_id cpPlatform;        // OpenCL platform  
41.     cl_device_id device_id;           // device ID  
42.     cl_context context;               // context  
43.     cl_command_queue queue;           // command queue  
44.     cl_program program;               // program  
45.     cl_kernel kernel;                 // kernel  
46.    
47.     // Size, in bytes, of each vector  
48.     size_t bytes = n*sizeof(double);  
49.    
50.     // Allocate memory for each vector on host  
51.     h_a = (double*)malloc(bytes);  
52.     h_b = (double*)malloc(bytes);  
53.     h_c = (double*)malloc(bytes);  
54.    
55.     // Initialize vectors on host  
56.     int i;  
57.     for( i = 0; i < n; i++ )  
58.     {  
59.         h_a[i] = sinf(i)*sinf(i);  
60.         h_b[i] = cosf(i)*cosf(i);  
61.     }  
62.      size_t globalSize, localSize;  
63.     cl_int err;  
64.    
65.     // Number of work items in each local work group  
66.     localSize = 64;  
67.    
68.     // Number of total work items - localSize must be devisor  
69.     globalSize = ceil(n/(float)localSize)*localSize;  
70.    
71.     // Bind to platform  
72.     err = clGetPlatformIDs(1, &cpPlatform, NULL);  
73.    
74.     // Get ID for the device  
75.     err = clGetDeviceIDs(cpPlatform, CL_DEVICE_TYPE_GPU, 1, &device_id, NULL);  
76.    
77.     // Create a context    
78.     context = clCreateContext(0, 1, &device_id, NULL, NULL, &err);  
79.    
80.     // Create a command queue   
81.     queue = clCreateCommandQueue(context, device_id, 0, &err);  
82.    
83.     // Create the compute program from the source buffer  
84.     program = clCreateProgramWithSource(context, 1,  
85.                             (const char **) & kernelSource, NULL, &err);  
86.    
87.     // Build the program executable   
88.     clBuildProgram(program, 0, NULL, NULL, NULL, NULL);  
89.    
90.     // Create the compute kernel in the program we wish to run  
91.     kernel = clCreateKernel(program, "vecAdd", &err);  
92.    
93.     // Create the input and output arrays in device memory for our calculation  
94.     d_a = clCreateBuffer(context, CL_MEM_READ_ONLY, bytes, NULL, NULL);  
95.     d_b = clCreateBuffer(context, CL_MEM_READ_ONLY, bytes, NULL, NULL);  
96.     d_c = clCreateBuffer(context, CL_MEM_WRITE_ONLY, bytes, NULL, NULL);  
97.    
98.     // Write our data set into the input array in device memory  
99.     err = clEnqueueWriteBuffer(queue, d_a, CL_TRUE, 0,  
100.                                    bytes, h_a, 0, NULL, NULL);  
101.     err |= clEnqueueWriteBuffer(queue, d_b, CL_TRUE, 0,  
102.                                    bytes, h_b, 0, NULL, NULL);  
103.    
104.     // Set the arguments to our compute kernel  
105.     err  = clSetKernelArg(kernel, 0, sizeof(cl_mem), &d_a);  
106.     err |= clSetKernelArg(kernel, 1, sizeof(cl_mem), &d_b);  
107.     err |= clSetKernelArg(kernel, 2, sizeof(cl_mem), &d_c);  
108.     err |= clSetKernelArg(kernel, 3, sizeof(unsigned int), &n);  
109.    
110.     // Execute the kernel over the entire range of the data set    
111.     err = clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &globalSize, &localSize,  
112.                                                               0, NULL, NULL);  
113.    
114.     // Wait for the command queue to get serviced before reading back results  
115.     clFinish(queue);  
116.    
117.     // Read the results from the device  
118.     clEnqueueReadBuffer(queue, d_c, CL_TRUE, 0,  
119.                                 bytes, h_c, 0, NULL, NULL );  
120.    
121.     //Sum up vector c and print result divided by n, this should equal 1 within error  
122.     double sum = 0;  
123.     for(i=0; i<n; i++)  
124.          sum += h_c[i];  
125.     printf("final result: %f\n", sum/n);  
126.    
127.     // release OpenCL resources  
128.     clReleaseMemObject(d_a);  
129.     clReleaseMemObject(d_b);  
130.     clReleaseMemObject(d_c);  
131.     clReleaseProgram(program);  
132.     clReleaseKernel(kernel);  
133.     clReleaseCommandQueue(queue);  
134.     clReleaseContext(context);  
135.    
136.     //release host memory  
137.     free(h_a);  
138.     free(h_b);  
139.     free(h_c);  
140.    
141.     return 0;  
142. }  

代码分析

内核（kernel）：

kernel是OpenCL代码的核心部分，整个内核必须通过C字符串的形式读入，最简单的办法是像代码一样定义一个长长的字符串，在真实的项目代码中通常都会从单独的文件中读入内核。

[cpp]  view plain  copy

 

1. // OpenCL kernel. Each work item takes care of one element of c  
2. const char *kernelSource =                                      "\n" \  
3. "#pragma OPENCL EXTENSION cl_khr_fp64 : enable                    \n" \  
4. "__kernel void vecAdd(  __global double *a,                       \n" \  
5. "                       __global double *b,                       \n" \  
6. "                       __global double *c,                       \n" \  
7. "                       const unsigned int n)                    \n" \  
8. "{                                                               \n" \  
9. "    //Get our global thread ID                                  \n" \  
10. "    int id = get_global_id(0);                                  \n" \  
11. "                                                                \n" \  
12. "    //Make sure we do not go out of bounds                      \n" \  
13. "    if (id < n)                                                 \n" \  
14. "        c[id] = a[id] + b[id];                                  \n" \  
15. "}                                                               \n" \  
16.                                                                 "\n" ;  

下面是内核的函数声明：

__kernel void  vecAdd(  __global  double  *a, __global  double  *b,                         __global double  *c,  const  unsigned  int  n)

__kernel是一个定义OpenCL内核的关键字，__global则定义函数指针指向全局设备内存空间，否则可以使用一般的C语言函数声明语法。内核的返回值必须为空void

int  id = get_global_id(0);

通过get_global_id函数可以获得当前工作单元（work item）的全局id，参数为0表示获取X维上的ID。

if (id < n)      c[id] = a[id] + b[id];

工作组（work group）的个数必定是整数，由于工作组的大小不一定是需要的线程数的整数倍，因此通常使用的线程数比需要的线程数要多，在程序设计时可以将无用的线程简单丢弃掉。

内存（Memory）

[cpp]  view plain  copy

 

1. // Host input vectors  
2. double *h_a;  
3. double *h_b;  
4. // Host output vector  
5. double *h_c;  
6.     
7. // Device input buffers  
8. cl_mem d_a;  
9. cl_mem d_b;  
10. // Device output buffer  
11. cl_mem d_c;  

主机CPU和GPU有不同的内存空间，因此需要分别定义，上面的代码中前半部分定义主机（host）CPU的内存指针，后半部分定义设备（device）内存的handle，分别用h_和d_前缀来区分。

线程映射（Thread Mapping）

[cpp]  view plain  copy

 

1. // Number of work items in each local work group  
2. localSize = 64;  
3.     
4. // Number of total work items - localSize must be devisor  
5. globalSize = ceil(n/(float)localSize)*localSize;  

为了将我们要解决的问题映射到底层硬件结构，必须定义局部尺寸（local size）和全局尺寸（global size）。局部尺寸定义了每个工作组中的工作单元数，在NVIDIA GPU上等价于每个线程块（thread block）中的线程数。全局尺寸定义了工作单元的总数目。局部尺寸必须是全局尺寸的倍数。

OpenCL前期准备（setup）

[cpp]  view plain  copy

 

1. // Bind to platform  
2. err = clGetPlatformIDs(1, &cpPlatform, NULL);  

每个硬件厂商都会绑定一个不同的平台（platform），在这里clGetPlatformIDs会将cpPlatform设置成包含系统可用平台的变量。举个例子，如果一个系统包含AMD CPU以及NVIDIA GPU，并且安装了恰当的OpenCL驱动，那么两个OpenCL平台会被返回。

[cpp]  view plain  copy

 

1. // Get ID for the device  
2. err = clGetDeviceIDs(cpPlatform, CL_DEVICE_TYPE_GPU, 1, &device_id, NULL);  

可以询问每一个平台都包含哪些设备，在这里我们通过使用CL_DEVICE_TYPE_GPU来查询GPU设备。

[cpp]  view plain  copy

 

1. // Create a context  
2. context = clCreateContext(0, 1, &device_id, NULL, NULL, &err);  

在使用OpenCL设备之前，必须先准备一个上下文（context），上下文对象用来管理命令队列（command queue）、内存（memory）、内核操作（Kernel activity），一个上下文对象可一般含多个设备。

[cpp]  view plain  copy

 

1. // Create a command queue  
2. queue = clCreateCommandQueue(context, device_id, 0, &err);  

命令队列（command queue）用来流式地将命令从主机送到指定的设备，可以把数据传输和内核操作命令放到命令队列上，当条件适宜的时候命令就会被执行。

编译内核（Compile Kernel）

[cpp]  view plain  copy

 

1. program = clCreateProgramWithSource(context, 1,  
2.                         (const char **) & kernelSource, NULL, &err);  
3.     
4. // Build the program executable  
5. clBuildProgram(program, 0, NULL, NULL, NULL, NULL);  
6.     
7. // Create the compute kernel in the program we wish to run  
8. kernel = clCreateKernel(program, "vecAdd", &err);  

为了保证OpenCL代码可以移植到许多不同的设备上，运行kernel的默认方式是JIT（Just-in-time, 实时编译）。首先创建一个program对象（包含一系列内核代码），然后再创建一系列的内核。

准备数据（prepare data）

[cpp]  view plain  copy

 

1. // Create the input and output arrays in device memory for our calculation  
2. d_a = clCreateBuffer(context, CL_MEM_READ_ONLY, bytes, NULL, NULL);  
3. d_b = clCreateBuffer(context, CL_MEM_READ_ONLY, bytes, NULL, NULL);  
4. d_c = clCreateBuffer(context, CL_MEM_WRITE_ONLY, bytes, NULL, NULL);  
5.     
6. // Write our data set into the input array in device memory  
7. err = clEnqueueWriteBuffer(queue, d_a, CL_TRUE, 0,  
8.                                bytes, h_a, 0, NULL, NULL);  
9. err |= clEnqueueWriteBuffer(queue, d_b, CL_TRUE, 0,  
10.                                bytes, h_b, 0, NULL, NULL);  
11.     
12. // Set the arguments to our compute kernel  
13. err  = clSetKernelArg(kernel, 0, sizeof(cl_mem), &d_a);  
14. err |= clSetKernelArg(kernel, 1, sizeof(cl_mem), &d_b);  
15. err |= clSetKernelArg(kernel, 2, sizeof(cl_mem), &d_c);  
16. err |= clSetKernelArg(kernel, 3, sizeof(unsigned int), &n);  

在启动内核之前，我们必须创建主机和设备之间的缓存（buffer），并将主机数据（host data）和这些新创建的设备缓存想绑定，最后再设定内核参数。

启动内核（Launch Kernel）

[cpp]  view plain  copy

 

1. // Execute the kernel over the entire range of the data set  
2. err = clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &globalSize, &localSize,  
3.                                                           0, NULL, NULL);  

将结果拷贝回主机（Copy results to host）

[cpp]  view plain  copy

 

1. // Wait for the command queue to get serviced before reading back results  
2. clFinish(queue);  
3.     
4. // Read the results from the device  
5. clEnqueueReadBuffer(queue, d_c, CL_TRUE, 0,  
6.                             bytes, h_c, 0, NULL, NULL );  

我们可以阻塞程序直到命令队列变为空，然后把结果拷贝回主机。

编译（Compile）

[plain]  view plain  copy

 

1. $ module load cudatoolkit  
2. $ cc -lOpenCL vecAdd.c -o vecAdd.out  

运行（Running）

[plain]  view plain  copy

 

1. $ aprun ./vecAdd.out  
2. final result: 1.000000  

VecAdd.cc

C++绑定在OpenCL的开发中非常常用，它比标准C接口更为流畅，下面是一个使用这些绑定的例子。

[cpp]  view plain  copy

 

1. #define __CL_ENABLE_EXCEPTIONS  
2.    
3. #include "cl.hpp"  
4. #include <cstdio>  
5. #include <cstdlib>  
6. #include <iostream>  
7. #include <math.h>  
8.    
9. // OpenCL kernel. Each work item takes care of one element of c  
10. const char *kernelSource =                                      "\n" \  
11. "#pragma OPENCL EXTENSION cl_khr_fp64 : enable                    \n" \  
12. "__kernel void vecAdd(  __global double *a,                       \n" \  
13. "                       __global double *b,                       \n" \  
14. "                       __global double *c,                       \n" \  
15. "                       const unsigned int n)                    \n" \  
16. "{                                                               \n" \  
17. "    //Get our global thread ID                                  \n" \  
18. "    int id = get_global_id(0);                                  \n" \  
19. "                                                                \n" \  
20. "    //Make sure we do not go out of bounds                      \n" \  
21. "    if (id < n)                                                 \n" \  
22. "        c[id] = a[id] + b[id];                                  \n" \  
23. "}                                                               \n" \  
24.                                                                 "\n" ;  
25.    
26.    
27. int main(int argc, char *argv[])  
28. {  
29.    
30.     // Length of vectors  
31.     unsigned int n = 1000;  
32.    
33.     // Host input vectors  
34.     double *h_a;  
35.     double *h_b;  
36.     // Host output vector  
37.     double *h_c;  
38.    
39.     // Device input buffers  
40.     cl::Buffer d_a;  
41.     cl::Buffer d_b;  
42.     // Device output buffer  
43.     cl::Buffer d_c;  
44.    
45.     // Size, in bytes, of each vector  
46.     size_t bytes = n*sizeof(double);  
47.    
48.     // Allocate memory for each vector on host  
49.     h_a = new double[n];  
50.     h_b = new double[n];  
51.     h_c = new double[n];  
52.    
53.     // Initialize vectors on host  
54.     for(int i = 0; i < n; i++ )  
55.     {  
56.         h_a[i] = sinf(i)*sinf(i);  
57.         h_b[i] = cosf(i)*cosf(i);  
58.     }  
59.    
60.     cl_int err = CL_SUCCESS;  
61.     try {  
62.     // Query platforms  
63.         std::vector<cl::Platform> platforms;  
64.         cl::Platform::get(&platforms);  
65.         if (platforms.size() == 0) {  
66.             std::cout << "Platform size 0\n";  
67.             return -1;  
68.          }  
69.    
70.         // Get list of devices on default platform and create context  
71.         cl_context_properties properties[] =  
72.            { CL_CONTEXT_PLATFORM, (cl_context_properties)(platforms[0])(), 0};  
73.         cl::Context context(CL_DEVICE_TYPE_GPU, properties);  
74.         std::vector<cl::Device> devices = context.getInfo<CL_CONTEXT_DEVICES>();  
75.    
76.         // Create command queue for first device  
77.         cl::CommandQueue queue(context, devices[0], 0, &err);  
78.    
79.         // Create device memory buffers  
80.         d_a = cl::Buffer(context, CL_MEM_READ_ONLY, bytes);  
81.         d_b = cl::Buffer(context, CL_MEM_READ_ONLY, bytes);  
82.         d_c = cl::Buffer(context, CL_MEM_WRITE_ONLY, bytes);  
83.    
84.         // Bind memory buffers  
85.         queue.enqueueWriteBuffer(d_a, CL_TRUE, 0, bytes, h_a);  
86.         queue.enqueueWriteBuffer(d_b, CL_TRUE, 0, bytes, h_b);  
87.    
88.         //Build kernel from source string  
89.         cl::Program::Sources source(1,  
90.             std::make_pair(kernelSource,strlen(kernelSource)));  
91.         cl::Program program_ = cl::Program(context, source);  
92.         program_.build(devices);  
93.    
94.         // Create kernel object  
95.         cl::Kernel kernel(program_, "vecAdd", &err);  
96.    
97.         // Bind kernel arguments to kernel  
98.         kernel.setArg(0, d_a);  
99.         kernel.setArg(1, d_b);  
100.         kernel.setArg(2, d_c);  
101.         kernel.setArg(3, n);  
102.    
103.         // Number of work items in each local work group  
104.         cl::NDRange localSize(64);  
105.         // Number of total work items - localSize must be devisor  
106.         cl::NDRange globalSize((int)(ceil(n/(float)64)*64));  
107.    
108.         // Enqueue kernel  
109.         cl::Event event;  
110.         queue.enqueueNDRangeKernel(  
111.             kernel,  
112.             cl::NullRange,  
113.             globalSize,  
114.             localSize,  
115.             NULL,  
116.             &event);  
117.    
118.         // Block until kernel completion  
119.         event.wait();  
120.      // Read back d_c  
121.         queue.enqueueReadBuffer(d_c, CL_TRUE, 0, bytes, h_c);  
122.         }  
123.     catch (cl::Error err) {  
124.          std::cerr  
125.             << "ERROR: "<<err.what()<<"("<<err.err()<<")"<<std::endl;  
126.     }  
127.    
128.     // Sum up vector c and print result divided by n, this should equal 1 within error  
129.     double sum = 0;  
130.     for(int i=0; i<n; i++)  
131.         sum += h_c[i];  
132.     std::cout<<"final result: "<<sum/n<<std::endl;  
133.    
134.     // Release host memory  
135.     delete(h_a);  
136.     delete(h_b);  
137.     delete(h_c);  
138.    
139.     return 0;  
140. }  

编译（Compile）

需要先下载cl.hpp

[plain]  view plain  copy

 

1. $ module load cudatoolkit  
2. $ CC vecAdd.cc -lOpenCL -o vecAdd.out  

运行（Running）

[plain]  view plain  copy

 

1. $ aprun ./vecAdd.out  
2. final result: 1.000000