本文详细介绍了一种高效的并行计算方法——数值归并运算，通过OpenCL实现标量和向量值的归并求和，并对比了二者的运算效率。文章提供了完整的OpenCL主机端和设备端代码，深入解析了归并算法的高效原因及OpenCL开发心得。
数值归并运算实质上就是求N项常数的总和的一种并行计算方法。
设数组的长度为N，若依次去运算数组中的每个元素(sum+=array[i])，显然时间复杂度就是O(N)；
对于数值归并运算还有一种更快捷的方法，可以通过将数组第i个元素与第（i+N/2）个元素相加，得到一个长度为N/2的数组；
得到的新数组中的每个元素，即为前一数组中元素两两相加所得（若原数组为奇数个，则保留最后一个元素至新数组）；
将这个过程重复(log2N)次，即可得到最终的计算结果，时间复杂度也为O(log2N)。
比较原始运算方法和归并求和运算方法，二者的计算次数均为(N-1)次。但是原始方法需要等前一工作项得到结果后才能进行下一次运算，浪费了其他工作项的运算能力。通过归并运算，可以使几乎全部（若总工作项<N/2的话则为全部）工作项投入并行化运算，从而减少了运算层次，提高程序效率。
OpenCL的归并计算在《OpenCL实战》这本书有一个C语言的实现，同时该作者比较了在归并运算中标量值和向量值运算效率上的区别，附上源代码。
OpenCL主机端完整代码：
#include "pch.h"
#include <iostream>
#define _CRT_SECURE_NO_WARNINGS
/* 执意使用老版本、非安全性的函数，可以使用 _CRT_SECURE_NO_WARNINGS 标记来忽略这些警告问题*/
#define PROGRAM_FILE "reduction.cl"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
/* 设置输出数组的大小，在预编译阶段会申请一部分内存，申请内存的大小若超过限制(一般为1M)则会报错*/
#define ARRAY_SIZE 1048576
#define NUM_KERNELS 2
/* 忽略No.4996警告，即OpenCL引起的版本冲突（许多指令由于安全性等已被禁用或已被更新为其他指令）*/
#pragma warning(disable : 4996)
#ifdef MAC
#include <OpenCL/cl.h>
#else
#include <CL/cl.h>
#endif
/* Find a GPU or CPU associated with the first available platform */
/* 发现可用平台下的GPU或CPU设备*/
cl_device_id create_device() {
cl_platform_id platform;
cl_device_id dev;
int err;
/* Identify a platform */
/* 识别平台*/
err = clGetPlatformIDs(1, &platform, NULL);
if (err < 0) {
perror("Couldn't identify a platform");
exit(1);
}
/* Access a device */
/* 使用一个设备*/
err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &dev, NULL);
if (err == CL_DEVICE_NOT_FOUND) {
err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_CPU, 1, &dev, NULL);
}
if (err < 0) {
perror("Couldn't access any devices");
exit(1);
}
return dev;
}
/* Create program from a file and compile it */
/* 从一个文件创建项目并进行编译*/
cl_program build_program(cl_context ctx, cl_device_id dev, const char* filename) {
cl_program program;
FILE *program_handle;
char *program_buffer, *program_log;
size_t program_size, log_size;
int err;
/* Read program file and place content into buffer */
program_handle = fopen(filename, "r");
if (program_handle == NULL) {
perror("Couldn't find the program file");
exit(1);
}
fseek(program_handle, 0, SEEK_END);
program_size = ftell(program_handle);
rewind(program_handle);
program_buffer = (char*)malloc(program_size + 1);
program_buffer[program_size] = '\0';
fread(program_buffer, sizeof(char), program_size, program_handle);
fclose(program_handle);
/* Create program from file */
program = clCreateProgramWithSource(ctx, 1,
(const char**)&program_buffer, &program_size, &err);
if (err < 0) {
perror("Couldn't create the program");
exit(1);
}
free(program_buffer);
/* Build program */
err = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
if (err < 0) {
/* Find size of log and print to std output */
clGetProgramBuildInfo(program, dev, CL_PROGRAM_BUILD_LOG,
0, NULL, &log_size);
program_log = (char*)malloc(log_size + 1);
program_log[log_size] = '\0';
clGetProgramBuildInfo(program, dev, CL_PROGRAM_BUILD_LOG,
log_size + 1, program_log, NULL);
printf("%s\n", program_log);
free(program_log);
exit(1);
}
return program;
}
using namespace std;
int main() {
/* OpenCL structures */
/* OpenCL 结构*/
cl_device_id device;
cl_context context;
cl_program program;
/* 数组形式存放kernel变量，NUM_KERNELS是内核数量*/
cl_kernel kernel[NUM_KERNELS];
cl_command_queue queue;
/* prof_event用于*/
cl_event prof_event;
cl_int i, j, err;
size_t local_size, global_size;
/* 若存在多个kernel函数，则命名放于这里*/
char kernel_names[NUM_KERNELS][20] =
{ "reduction_scalar", "reduction_vector" };
/* Data and buffers */
/* 数据和缓存*/
float data[ARRAY_SIZE];
float sum, actual_sum, *scalar_sum, *vector_sum;
/* 定义三个内存对象*/
cl_mem data_buffer, scalar_sum_buffer, vector_sum_buffer;
cl_int num_groups;
/* 定义计时变量*/
cl_ulong time_start, time_end, total_time;
/* Initialize data */
/* 初始化数据*/
for (i = 0; i < ARRAY_SIZE; i++) {
data[i] = 1.0f*i;
}
/* Create device and determine local size */
/* 创建设备并决定本地大小*/
device = create_device();
err = clGetDeviceInfo(device, CL_DEVICE_MAX_WORK_GROUP_SIZE,
sizeof(local_size), &local_size, NULL);
if (err < 0) {
perror("Couldn't obtain device information");
exit(1);
}
/* Allocate and initialize output arrays */
/* 申请内存并初始化输出数组*/
num_groups = ARRAY_SIZE / local_size;
scalar_sum = (float*)malloc(num_groups * sizeof(float));
vector_sum = (float*)malloc(num_groups / 4 * sizeof(float));
for (i = 0; i < num_groups; i++) {
scalar_sum[i] = 0.0f;
}
for (i = 0; i < num_groups / 4; i++) {
vector_sum[i] = 0.0f;
}
/* Create a context */
/* 创建一个上下文*/
context = clCreateContext(NULL, 1, &device, NULL, NULL, &err);
if (err < 0) {
perror("Couldn't create a context");
exit(1);
}
/* Build program */
/* 构建程序*/
program = build_program(context, device, PROGRAM_FILE);
/* Create data buffer */
/* 创建数据buffer缓存，OpenCL一共可以创建Buffer和Image两种内存对象类型，实际应用中具体用途有所区别*/
data_buffer = clCreateBuffer(context, CL_MEM_READ_ONLY |
CL_MEM_COPY_HOST_PTR, ARRAY_SIZE * sizeof(float), data, &err);
scalar_sum_buffer = clCreateBuffer(context, CL_MEM_READ_WRITE |
CL_MEM_COPY_HOST_PTR, num_groups * sizeof(float), scalar_sum, &err);
vector_sum_buffer = clCreateBuffer(context, CL_MEM_READ_WRITE |
CL_MEM_COPY_HOST_PTR, num_groups * sizeof(float), vector_sum, &err);
if (err < 0) {
perror("Couldn't create a buffer");
exit(1);
};
/* Create a command queue */
/* 创建一个命令队列*/
queue = clCreateCommandQueue(context, device,
CL_QUEUE_PROFILING_ENABLE, &err);
if (err < 0) {
perror("Couldn't create a command queue");
exit(1);
};
for (i = 0; i < NUM_KERNELS; i++) {
/* Create a kernel */
/* 创建内核*/
kernel[i] = clCreateKernel(program, kernel_names[i], &err);
if (err < 0) {
perror("Couldn't create a kernel");
exit(1);
};
/* Create kernel arguments */
/* 多个kernel程序进行内核参数设定，i=1为kernel系数*/
err = clSetKernelArg(kernel[i], 0, sizeof(cl_mem), &data_buffer);
if (i == 0) {
global_size = ARRAY_SIZE;
err |= clSetKernelArg(kernel[i], 1, local_size * sizeof(float), NULL);
err |= clSetKernelArg(kernel[i], 2, sizeof(cl_mem), &scalar_sum_buffer);
}
else {
global_size = ARRAY_SIZE / 4;
err |= clSetKernelArg(kernel[i], 1, local_size * 4 * sizeof(float), NULL);
err |= clSetKernelArg(kernel[i], 2, sizeof(cl_mem), &vector_sum_buffer);
}
if (err < 0) {
perror("Couldn't create a kernel argument");
exit(1);
}
/* Enqueue kernel */
/* 将命令入列，若依次执行多个kernel函数，则可以设置一个循环循环入列执行操作，clFinish是否进入循环待验证*/
err = clEnqueueNDRangeKernel(queue, kernel[i], 1, NULL, &global_size,
&local_size, 0, NULL, &prof_event);
if (err < 0) {
perror("Couldn't enqueue the kernel");
exit(1);
}
/* Finish processing the queue and get profiling information */
/* 完成命令队列的处理并得到概述信息*/
clFinish(queue);
/*通过clGetEventProfilingInfo得到事件发生的时间*/
clGetEventProfilingInfo(prof_event, CL_PROFILING_COMMAND_START,
sizeof(time_start), &time_start, NULL);
clGetEventProfilingInfo(prof_event, CL_PROFILING_COMMAND_END,
sizeof(time_end), &time_end, NULL);
/*计算每个过程执行时间*/
total_time = time_end - time_start;
/* Read the result */
/* 读取结果*/
if (i == 0) {
/*将标量传入kernel计算*/
err = clEnqueueReadBuffer(queue, scalar_sum_buffer, CL_TRUE, 0,
num_groups * sizeof(float), scalar_sum, 0, NULL, NULL);
if (err < 0) {
perror("Couldn't read the buffer");
exit(1);
}
sum = 0.0f;
for (j = 0; j < num_groups; j++) {
sum += scalar_sum[j];
}
}
else {
/*将向量传入kernel计算*/
err = clEnqueueReadBuffer(queue, vector_sum_buffer, CL_TRUE, 0,
num_groups / 4 * sizeof(float), vector_sum, 0, NULL, NULL);
if (err < 0) {
perror("Couldn't read the buffer");
exit(1);
}
sum = 0.0f;
for (j = 0; j < num_groups / 4; j++) {
sum += vector_sum[j];
}
}
/* Check result */
/* 校验运算结果*/
printf("%s: ", kernel_names[i]);
actual_sum = 1.0f * ARRAY_SIZE / 2 * (ARRAY_SIZE - 1);
if (fabs(sum - actual_sum) > 0.01*fabs(sum))
printf("Check failed.\n");
else
printf("Check passed.\n");
std::cout << "Total time = " << total_time << std::endl;
/* Deallocate event */
/* 释放事件——为什么要每次释放事件? 实时监测每次运算中的报错，prof_event作为参数输入不同cl函数可以输出不同信息*/
clReleaseEvent(prof_event);
}
/* Deallocate resources */
/* 释放资源*/
free(scalar_sum);
free(vector_sum);
for (i = 0; i < NUM_KERNELS; i++) {
clReleaseKernel(kernel[i]);
}
/* 在循环执行的最后释放所有的内存对象，这些内存对象可以循环使用不需要中途释放重新建立，否则太影响效率*/
clReleaseMemObject(scalar_sum_buffer);
clReleaseMemObject(vector_sum_buffer);
clReleaseMemObject(data_buffer);
clReleaseCommandQueue(queue);
clReleaseProgram(program);
clReleaseContext(context);
return 0;
}
设备端 reduction.cl文件完整代码：
__kernel void reduction_scalar(__global float* data,
__local float* partial_sums, __global float* output) {
int lid = get_local_id(0);
int group_size = get_local_size(0);
partial_sums[lid] = data[get_global_id(0)];
barrier(CLK_LOCAL_MEM_FENCE);
for(int i = group_size/2; i>0; i >>= 1) {
if(lid < i) {
partial_sums[lid] += partial_sums[lid + i];
}
barrier(CLK_LOCAL_MEM_FENCE);
}
if(lid == 0) {
output[get_group_id(0)] = partial_sums[0];
}
}
__kernel void reduction_vector(__global float4* data,
__local float4* partial_sums, __global float* output) {
int lid = get_local_id(0);
int group_size = get_local_size(0);
partial_sums[lid] = data[get_global_id(0)];
barrier(CLK_LOCAL_MEM_FENCE);
for(int i = group_size/2; i>0; i >>= 1) {
if(lid < i) {
partial_sums[lid] += partial_sums[lid + i];
}
barrier(CLK_LOCAL_MEM_FENCE);
}
if(lid == 0) {
output[get_group_id(0)] = dot(partial_sums[0], (float4)(1.0f));
}
}
该代码的运算速度对比：
归并算法十分高效的原因：
(1) 未使用取模运算符%。该运算符在CPU和GPU上的耗时差异很大，常用于确定数的奇偶性。一般情况下，在并行异构开发中不考虑使用该运算符；
(2) 与内存有关，局部内存物理层面上成组分布于各个内存条上。这些内存条中存在内存交织(Memory Interleave)的技术，因此工作项同一时间访问不同内存条不会占用额外的时间，其读写操作可并行完成。
使用OpenCL一段时间的小总结：
- CPU与GPU之间的传输消耗是一个很值得重视的问题，在使用OpenCL之前必须对整个项目进行评估，是否该工程的运算量达到了需要使用GPU的程度。
- 虽然OpenCL是一种跨平台的标准，但是在AMD显卡的设备条件下使用OpenCL进行开发还是应该首先搭建Linux+Eclipse(C++)+CodeXL的环境，使用Intel显卡的设备可以使用Windows+Visual Studio(latest version)+Nsight的环境，而对于Apple OS下的开发可以使用Xcode环境。搭建这些环境能极大提高使用OpenCL的效率，因为能实时反馈OpenCL在设备上的性能。
- 读取Kernel, 发现设备，建立平台等等步骤，不同项目之间可以直接copy+paste，因为初始化的小流程是亘古不变的。
- OpenCL重点关注如何对齐内存（如何设置全局内存、本地内存、私有内存和常量内存的问题），使用什么样的内存对象，使用什么样的运算函数简便kernel中的计算过程，如何将串行计算转化为并行计算的方法（例如将简单求和运算变为归并运算）以及内存对象释放的问题。
- 如何提高OpenCL的应用能力：需要多自己动手实验，尝试去写需要分配local memory的程序；同时对OpenCL中用到的数据类型要进行总结，在使用时才可以游刃有余；多去掌握可以进行通过并行化，节约计算时间的算法，多提问多思考；多参与社区中的讨论（很多社区已经不活跃了，但是例如GPU World中还是有很多久远的内容可以借鉴），关注博客中的评论中的水友回复+楼主回答，也可以关注其他GPU编程相关社区（Cuda/C++ AMP）以及GPU处理的OpenGL的相关内容；多和比自己强的人交流，设定一个要超越的目标不断努力。
- 使用Git hub上传自己的代码，通过代码管理巩固自己已有的知识体系。
扫一扫
3443
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
