OpenCL案例研究之一<矩阵乘法>

一、 案例简述

　　本案例讲述使用OpenCL计算矩阵乘法：C = A * B 。

　　设A、B、C分别是大小为N*P、P*M和N*M的矩阵，那么顺序实现的C代码可以如下所示：

 // C Functionvoid 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)){...}