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

最新推荐文章于 2023-10-27 09:48:22 发布
最新推荐文章于 2023-10-27 09:48:22 发布
阅读量1.1w 收藏 23
点赞数 8
分类专栏： OpenCL
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/cloud_desktop/article/details/19822025
版权


一、 案例简述


　　本案例讲述使用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))
	{
		...
	}










Leon Lu
关注
  • 8
    点赞
  • 踩
  • 23
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
专栏目录
opencl 多gpu矩阵数组运算
12-06
1：opencl 多gpu运算； 2：opencl 矩阵运算； 3：opencl 数组运算；
OpenCL例程3-矩阵相乘
HAIHUAYEN的专栏
03-31 1999
cpp文件代码 // OpenCl-3.cpp : Defines the entry point for the console application. // #include "stdafx.h" #include #include #include using namespace std; #define MATRIX_A_ROW 8 // 矩阵A行数 #defin
openCL-矩阵相乘
数学屌丝走在it的道路上
11-20 4583
这是AMD异构大赛发的书上的例子，我自己加了一些东西，实现了一下。不过还没有学会如何分析运行的时间，先贴上代码 这是第一次用amd APP KernelAnalyzer感觉还可以吧， simpleMultiply.cl // Enter your kernel in this window __kernel v
6 OpenCL matmpy矩阵乘法.docx
11-24
6 OpenCL matmpy矩阵乘法
AM5728/AM57XX openCL使用案例三(矩阵转置)
07-18
下面对4x4矩阵进行转置，调用openCL，4X4矩阵采用二维数组进行存储，在程序设计上，让转置过程分10次转置完成，就是一次转一行，因此OpenCL的工作维数是二维。 openCL程序分为两个部份，一部份是内核代码，负责具体算法。另一部份是主程序负责初始化OpenCL和准备数据。主程序加载内核代码，并按照既定方法进行运算。
C++实战之OpenCL矩阵相乘
天马行空，隐姓埋名
11-06 5325
简单概念理解在opencl中，有个索引空间NDRange的概念，NDRange是一个N维的索引空间，N可以是1，2，3。NDRange由一个长度为N的整数阵列来定义，他指定了索引空间各个维度的宽度，每个work-item的全局id和局部id，都是N维元组。有多个work-item构成的叫做work-group，作業組的 ID 跟作業項的全局 ID差不多。一個長度為 N 的陣列定義了每個維度上作業組的
OpenCL矩阵乘法
张小殊的博客
10-19 3672
OpenCL (Open Computing Language) 是由Khronos Group制定的开放式并行计算标准。它允许开发者利用异构计算平台上的多核CPU和GPU等处理器进行并行计算。OpenCL可以用于各种应用领域，包括科学计算、图形处理、机器学习等。
OpenCL之矩阵乘法实现
热门推荐
I AM BACK
05-01 1万+
本篇博客写的是矩阵乘法的opencl实现，算是一个入门级的程序，附上kernel和源码，很适合初学者看，一起学习opencl吧！
openCL矩阵乘法
最新发布
weixin_49562509的博客
10-27 77
【代码】openCL矩阵乘法。
C++实战之OpenCL矩阵相乘优化（二）
天马行空，隐姓埋名
11-21 9258
前言上一篇文章，分析了简单的矩阵相乘在opencl里面的优化kernel代码，每个work-item只负责计算结果矩阵的一个元素。下一步准备每次计算出结果矩阵的块元素，看看计算时间是如何。这个矩阵系列参考国外一个大神的教程： https://cnugteren.github.io/tutorial/pages/page4.html 有每个kernel的详解工程，还有github代码工程。具体分析
opencl gpu 调用实现矩阵运算和迪杰斯特拉算法
07-13
opencl gpu 调用实现矩阵运算和迪杰斯特拉算法
OpenCL学习stepbystep(一)
03-03
从今天开始学习OpenCL……安装AMDOpenCLAPP的过程我这里就不仔细说明了。现在，我们开始写一个简单的OpenCL程序，计算两个数组相加的和，放到另一个数组中去。程序用cpu和gpu分别计算，最后验证它们是否相等。OpenCL程序的流程大致如下：下面是sourcecode中的主要代码：注意：如果我们系统中安装不止一个opencl平台，比如我的os中，有intel和amd两家opencl平台，用上面这行代码，有可能会出错，因为它得到了intel的opencl平台，而intel的平台只支持cpu，而我们后面的操作都是基于gpu，这时我们可以用下面的代码，得到AMD的opencl平台上面这
OpenCL 初实践（1）矩阵相乘
u010741691的博客
08-11 900
对于设定global item size 对kernel函数 是很有影响的，不同的item size 会影响构造
PyOpenCL－02.向量与矩阵运算
Sun zi chao的博客
06-21 450
Python是目前较流行的一种科学计算语言。语法简洁,上手快，易于维护的优点。但其运算速度是真心的慢。那我们能否利用Python的简洁+OpenCL的运算能力呢？答案是可以的，那就是PyOpenCL。
关于《OPENCL异构并行计算》中矩阵乘法的分析
qq_20028731的博客
04-15 2498
《OPENCL异构并行计算》中讲了如何利用OPENCL进行矩阵乘法运算，并给出了使用局部存储器优化、使用向量加载指令以及一个工作项同时计算多个输出的例子，这里对其进行简单分析
OpenCL概述 续篇
jaccen的专栏
05-10 1042
摘要：Kernel就是在程序代码中的一个函数，这个函数能在OpenCL设备上执行。一个Kernel对象就是kernel函数以及其相关的输入参数。 Kernel对象： Kernel就是在程序代码中的一个函数，这个函数能在OpenCL设备上执行。一个Kernel对象就是kernel函数以及其相关的输入参数。 Kernel对象通过程序对象以及指定的函数名字创建。注意：函数必
OpenCL矩阵乘法的例子
2Know的博客
10-13 2695
前两篇博客中介绍过矩阵转置的两种方法，矩阵乘法可以先对矩阵做转置运算，然后再对应相乘； 矩阵大小是65*8192 先对8192*65的矩阵转置为65*8192； 然后由65个工作项，每个工作项负责一行数据的乘法； 完成65行数据的乘法； 转置是由8192个工作项完成的； 代码中给出了C语言实现矩阵转置的测试代码； 感兴趣的朋友可以测试一下C语言代码的运算时间和GPU运算时间的差别；
opencl学习（一）
random_repick的专栏
11-22 2461
细分了流程，并注释了相应函数 对get_global_id的理解：如果clEnqueueNDRangeKernel设置的是一维且global_item_size = 9，那么get_global_id(0)返回0到8，对应9个work-group。如果clEnqueueNDRangeKernel设置的是二维且global_item_size[2] = {9, 8}，那么get_global_id(0)和get_global_id(1)的结果是(0,j)(1,j)(2,j)(3,j)(4,j)(5,j)(6.
opencl矩阵乘法
10-25
OpenCL是一种跨平台的并行计算框架，可以用于GPU、CPU和FPGA等不同类型的处理器。矩阵乘法是一个常见的并行计算问题，可以使用OpenCL来加速计算。 在OpenCL中，矩阵乘法可以通过将矩阵分成多个小块，然后在不同的处理器上并行计算来实现。具体来说，可以将输入矩阵A和B分别划分成多个小块，然后在不同的处理器上并行计算每个小块的结果。最后将所有小块的结果合并起来得到最终的输出矩阵C。 在实现OpenCL矩阵乘法时，需要编写OpenCL内核函数来执行矩阵乘法计算。内核函数可以使用OpenCL提供的向量和矩阵数据类型来简化计算过程。同时，还需要编写主机代码来管理OpenCL设备和内存，并将数据传输到设备上进行计算。 下面是一个简单的OpenCL矩阵乘法示例代码，仅供参考： ``` __kernel void matrix_multiply(__global float* A, __global float* B, __global float* C, int M, int N, int K) { int i = get_global_id(0); int j = get_global_id(1); float sum = 0.0f; for (int k = 0; k < K; k++) { sum += A[i*K+k] * B[k*N+j]; } C[i*N+j] = sum; } int main() { // 初始化OpenCL设备和内存 cl_device_id device; cl_context context; cl_command_queue queue; cl_mem bufferA, bufferB, bufferC; // ... // 将数据传输到设备上 clEnqueueWriteBuffer(queue, bufferA, CL_TRUE, 0, sizeof(float)*M*K, A, 0, NULL, NULL); clEnqueueWriteBuffer(queue, bufferB, CL_TRUE, 0, sizeof(float)*K*N, B, 0, NULL, NULL); // 执行内核函数 size_t global_size[2] = {M, N}; clEnqueueNDRangeKernel(queue, kernel, 2, NULL, global_size, NULL, 0, NULL, NULL); // 将结果传输回主机 clEnqueueReadBuffer(queue, bufferC, CL_TRUE, 0, sizeof(float)*M*N, C, 0, NULL, NULL); // 清理OpenCL设备和内存 // ... return 0; } --相关问题--:

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值