一、 案例简述
本案例讲述使用OpenCL计算矩阵乘法:C = A * B 。
设A、B、C分别是大小为N*P、P*M和N*M的矩阵,那么顺序实现的C代码可以如下所示:
// C Function
void mat_mul(
int Ndim, int Mdim, int Pdim,
float* A, float* B, float* C)
{
int i, j, k;
float tmp;
for (i = 0; i < Ndim; i++) {
for (j = 0; j < Mdim; j++) {
tmp = 0.0;
for (k = 0; k < Pdim; k++)
tmp += A[i*Pdim + k] * B[k*Mdim + j];
C[i*Mdim + j] = tmp;
}
}
}
二、 OpenCL实现矩阵乘法
1. 内核函数实现
// OpenCL Kernel Function
__kernel void HelloOpenCL(
const int Ndim,
const int Mdim,
const int Pdim,
__global const float* A,
__global const float* B,
__global float* C)
{
int i = get_global_id(0);
int j = get_global_id(1);
int k;
float tmp;
if ((i < Ndim) && (j < Mdim)) {
tmp = 0.0;
for (k = 0; k < Pdim; k++)
tmp += A[i*Pdim + k] * B[k*Mdim + j];
C[i*Mdim + j] = tmp;
}
}
为每个工作项分配一个要计算的乘法矩阵的元素。将针对i,j的外层循环删除,替换为函数调用,查找这两维中对应工作项的全局ID。要特别当心,必须保证得到的工作项ID在矩阵C的范围内。这三个矩阵都留在全局内存中。
2. 宿主机代码实现
下面是在《基于CUDA的OpenCL开发环境搭建与入门程序示例》中main.cpp宿主机代码为基础的补丁文件。测量运行时间的部分:首先,在clCreateCommandQueue()函数中设置CL_QUEUE_PROFILING_ENABLE标志;然后,在clEnqueueNDRangeKernel()函数中设置事件对象;最后,通过clGetEventProfilingInfo()函数获取命令入队时间和命令执行结束时间。注意:时间的单位是纳秒,在最后打印时转换为秒显示。
--- /root/Desktop/main.cpp
+++ /root/Desktop/main_new.cpp
@@ -143,8 +143,10 @@
}
// 4. Choose the first device
- commandQueue = clCreateCommandQueue(context,
- devices[0], 0, NULL);
+ commandQueue = clCreateCommandQueue(context,
+ devices[0],
+ CL_QUEUE_PROFILING_ENABLE,
+ NULL);
if (commandQueue == NULL) {
perror("Failed to create commandQueue for device 0.");
exit(1);
@@ -183,14 +185,33 @@
/******** 第四部分 创建内核和内存对象 ********/
- #define ARRAY_SIZE 10
+ const int Ndim = 3;
+ const int Mdim = 4;
+ const int Pdim = 5;
+
+ int szA = Ndim * Pdim;
+ int szB = Pdim * Mdim;
+ int szC = Ndim * Mdim;
cl_kernel kernel = 0;
cl_mem memObjects[3] = {0, 0, 0};
- float a[ARRAY_SIZE];
- float b[ARRAY_SIZE];
- float result[ARRAY_SIZE];
+ float *A;
+ float *B;
+ float *C;
+
+ A = (float *)malloc(szA * sizeof(float));
+ B = (float *)malloc(szB * sizeof(float));
+ C = (float *)malloc(szC * sizeof(float));
+
+ int i, j;
+
+ for (i = 0; i < szA; i++)
+ A[i] = i + 1;
+
+ for (i = 0; i < szB; i++)
+ B[i] = i + 1;
+
// 8. Create the kernel
kernel = clCreateKernel(program, "HelloOpenCL", NULL);
@@ -200,23 +221,18 @@
}
// 9. Create memory objects
- for (int i = 0; i < ARRAY_SIZE; i++) {
- a[i] = (float)i + 1;
- b[i] = (float)i + 1;
- }
-
memObjects[0] = clCreateBuffer(context, CL_MEM_READ_ONLY |
CL_MEM_COPY_HOST_PTR,
- sizeof(float) * ARRAY_SIZE,
- a, NULL);
+ sizeof(float) * szA,
+ A, NULL);
memObjects[1] = clCreateBuffer(context, CL_MEM_READ_ONLY |
CL_MEM_COPY_HOST_PTR,
- sizeof(float) * ARRAY_SIZE,
- b, NULL);
+ sizeof(float) * szB,
+ B, NULL);
memObjects[2] = clCreateBuffer(context, CL_MEM_READ_WRITE |
CL_MEM_COPY_HOST_PTR,
- sizeof(float) * ARRAY_SIZE,
- result, NULL);
+ sizeof(float) * szC,
+ C, NULL);
if (memObjects[0] == NULL || memObjects[1] == NULL ||
memObjects[2] == NULL) {
perror("Error in clCreateBuffer.\n");
@@ -225,48 +241,98 @@
/******** 第五部分 执行内核 ********/
- size_t globalWorkSize[1] = { ARRAY_SIZE };
- size_t localWorkSize[1] = { 1 };
// 10. Set the kernel arguments
- errNum = clSetKernelArg(kernel, 0, sizeof(cl_mem), &memObjects[0]);
- errNum |= clSetKernelArg(kernel, 1, sizeof(cl_mem), &memObjects[1]);
- errNum |= clSetKernelArg(kernel, 2, sizeof(cl_mem), &memObjects[2]);
+ errNum = clSetKernelArg(kernel, 0, sizeof(int), &Ndim);
+ errNum |= clSetKernelArg(kernel, 1, sizeof(int), &Mdim);
+ errNum |= clSetKernelArg(kernel, 2, sizeof(int), &Pdim);
+ errNum |= clSetKernelArg(kernel, 3, sizeof(cl_mem), &memObjects[0]);
+ errNum |= clSetKernelArg(kernel, 4, sizeof(cl_mem), &memObjects[1]);
+ errNum |= clSetKernelArg(kernel, 5, sizeof(cl_mem), &memObjects[2]);
if (errNum != CL_SUCCESS) {
perror("Error in clSetKernelArg.\n");
exit(1);
}
// 11. Queue the kernel up for execution across the array
- errNum = clEnqueueNDRangeKernel(commandQueue, kernel, 1, NULL,
- globalWorkSize, localWorkSize,
- 0, NULL, NULL);
+ size_t global[2];
+ cl_event prof_event;
+ cl_ulong ev_start_time = (cl_ulong)0;
+ cl_ulong ev_end_time = (cl_ulong)0;
+ double rum_time;
+
+ global[0] = (size_t)Ndim;
+ global[1] = (size_t)Mdim;
+
+ errNum = clEnqueueNDRangeKernel(commandQueue, kernel, 2, NULL,
+ global, NULL, 0, NULL, &prof_event);
if (errNum != CL_SUCCESS) {
perror("Error in clEnqueueNDRangeKernel.\n");
exit(1);
}
+ clFinish(commandQueue);
+ errNum = clWaitForEvents(1, &prof_event);
+ if (errNum != CL_SUCCESS) {
+ perror("Error in clWaitForEvents.\n");
+ exit(1);
+ }
+
+ errNum = clGetEventProfilingInfo(prof_event,
+ CL_PROFILING_COMMAND_QUEUED,
+ sizeof(cl_ulong),
+ &ev_start_time,
+ NULL);
+
+ errNum |= clGetEventProfilingInfo(prof_event,
+ CL_PROFILING_COMMAND_END,
+ sizeof(cl_ulong),
+ &ev_end_time,
+ NULL);
+
+ if (errNum != CL_SUCCESS) {
+ perror("Error in clGetEventProfilingInfo.\n");
+ while(1);
+ exit(1);
+ }
+
// 12. Read the output buffer back to the Host
errNum = clEnqueueReadBuffer(commandQueue, memObjects[2],
CL_TRUE, 0,
- ARRAY_SIZE * sizeof(float), result,
+ sizeof(float) * szC, C,
0, NULL, NULL);
if (errNum != CL_SUCCESS) {
perror("Error in clEnqueueReadBuffer.\n");
exit(1);
}
+ rum_time = (double)(ev_end_time - ev_start_time);
+
/******** 第六部分 测试结果 ********/
- printf("\nTest: a * b = c\n\n");
-
- printf("Input numbers:\n");
- for (int i = 0; i < ARRAY_SIZE; i++)
- printf("a[%d] = %f, b[%d] = %f\n", i, a[i], i, b[i]);
-
- printf("\nOutput numbers:\n");
- for (int i = 0; i < ARRAY_SIZE; i++)
- printf("a[%d] * b[%d] = %f\n", i, i, result[i]);
+
+ printf("\nArray A:\n");
+ for (i = 0; i < Ndim; i++) {
+ for (j = 0; j < Pdim; j++)
+ printf("%.3f\t", A[i*Pdim + j]);
+ printf("\n");
+ }
+
+ printf("\nArray B:\n");
+ for (i = 0; i < Pdim; i++) {
+ for (j = 0; j < Mdim; j++)
+ printf("%.3f\t", B[i*Mdim + j]);
+ printf("\n");
+ }
+
+ printf("\nArray C:\n");
+ for (i = 0; i < Ndim; i++) {
+ for (j = 0; j < Mdim; j++)
+ printf("%.3f\t", C[i*Mdim + j]);
+ printf("\n");
+ }
+
+ printf("\n\nRunning Time: %f s\n", rum_time*1.0e-9);
while(1);
3. 运行结果
(1). N = 3,M = 4,P = 5。
(2). N = 1000,M = 1000,P = 1000。
三、 代码优化 <工作项分组和减少数据移动>
矩阵乘法的核心是一个乘加计算。大多数处理器中的ALU都有足够的带宽,可以保证这个计算可以接近峰值性能运行,不过只有隐藏数据移动时的开销时才能做到这一点。因此,优化矩阵运算的根本就是尽量减少数据移动。上面的代码中的矩阵乘法内核,三个矩阵都保留在全局内存中。这意味着每次乘法都要通过内存层次结构反复地传递行和列(全局内存到私有内存)。
1. 第一次优化
在此优化版本中,每个工作项计算矩阵中的一行。NDRange从一个2D range(分区)集(匹配矩阵C的维度)变为一个1D range(分区)集(匹配矩阵C的行数)。如下内核代码中,每个工作项管理C中一整行的更新,不过在完成这个更新之前,要把矩阵A中相关的行从全局内存复制到私有内存。注:优化代码以(N = 1000,M = 1000,P = 1000)为例。
// OpenCL Kernel Function
__kernel void HelloOpenCL(
const int Ndim,
const int Mdim,
const int Pdim,
__global const float* A,
__global const float* B,
__global float* C)
{
int i = get_global_id(0);
//int j = get_global_id(1);
int j, k;
float tmp;
float Awrk[1000];
if (i < Ndim) {
for (k = 0; k < Pdim; k++)
Awrk[k] = A[i*Pdim + k];
for (j = 0; j < Mdim; j++) {
tmp = 0.0;
for (k = 0; k < Pdim; k++)
tmp += Awrk[k] * B[k*Mdim + j];
C[i*Mdim + j] = tmp;
}
}
}
在宿主机代码中,只需将clEnqueueNDRangeKernel()函数中的维度数由"2"改为"1"。
errNum = clEnqueueNDRangeKernel(commandQueue, kernel, 1, NULL,
global, NULL, 0, NULL, &prof_event);
最终运行结果: 0.295128 秒,是上面代码中0.692887秒的42.6% !
本次优化通过增加内核工作量,减少数据移动来提高运行速度。其实如果是在GPU中而不是CPU中,单独使用上述方式增加内核工作量,运行时间反而变为原来几倍。
2. 第二次优化
首先获取平台上的设备数及其最大计算单元数,使用以下代码。
//OpenCL设备信息
cl_uint numDevices;
cl_device_id deviceIds[1];
size_t maxComputeUnits;
errNum = clGetDeviceIDs (platformIds[0],
CL_DEVICE_TYPE_GPU,
0, NULL, &numDevices);
if ((errNum != CL_SUCCESS) || (numDevices < 1)) {
perror("Error in clGetDeviceIDs or no GPU deivce.");
exit(1);
}
errNum = clGetDeviceIDs (platformIds[0],
CL_DEVICE_TYPE_GPU,
1, &deviceIds[0], NULL);
if ((errNum != CL_SUCCESS) || (numDevices < 1)) {
perror("Error in clGetDeviceIDs.");
exit(1);
}
errNum = clGetDeviceInfo(deviceIds[0],
CL_DEVICE_MAX_COMPUTE_UNITS,
sizeof(cl_uint),
&maxComputeUnits,
NULL);
printf("numDevices = %d, maxComputeUnits = %d\n", numDevices, maxComputeUnits);
运行结果: 平台中由一个设备,最大计算单元数为7。在此,我们划分为4个工作组,每个工作组大小设置为250。具体代码如下:
size_t global[1];
size_t local[1];
global[0] = (size_t)Ndim;
local[0] = (size_t)250;
errNum = clEnqueueNDRangeKernel(commandQueue, kernel, 1, NULL,
global, local, 0, NULL, &prof_event);
最终运行时间为: 0.208867 秒,约为原有基础上的2/3 !
4. 附录:关于二维数据的工作组分组
在宿主机代码中,根据clGetDeviceInfo()函数查找CL_DEVICE_MAX_WORK_GROUP_SIZE标志所对应工作组内最多元素个数,在我的电脑上,为1024。工作组的大小通常是64的倍数,最好不超过256。localx,localy也有要求,根据我的试验,必须是4的倍数。
size_t local[2];
local[0] = 8;
local[1] = 128;
errNum = clEnqueueNDRangeKernel(commandQueue, kernel, 2, NULL, global, local, 0, NULL, &prof_event);
由于globalx与globaly必须恰好整除localx与localy,因此通过下述写法确保在图像宽与高不能整除对应localx和localy值时,能够满足条件。
global[0] = (nWidth + local[0] - 1) / local[0];
global[1] = (nHeight + local[1] - 1) / local[1];
global[0] *= local[0];
global[1] *= local[1];
为此,在相应内核代码中,必须判断所计算出的globalx与globaly的值是否还在图像范围之内。如下:
if((x < nWidth) && (y < nHeight))
{
...
}