【Intel校企合作课程】基于oneAPI的C++/SYCL（DPC++）的并行矩阵乘法操作

最新推荐文章于 2024-06-12 20:44:27 发布

Forever·流年

最新推荐文章于 2024-06-12 20:44:27 发布

阅读量1k

点赞数 20

文章标签： oneapi c++

本文链接：https://blog.csdn.net/qq_53020182/article/details/134699634

版权

1：作业简介：并行矩阵计算

1.1 问题陈述：

实现一个矩阵乘法的操作，需要考虑大尺寸矩阵的乘法操作以及不同线程之间的数据依赖关系。

1.2 预期解决方案

利用Intel图形处理器对矩阵进行乘法操作，采用DPC++并行线程模型，利用线程块和线程的层次结构来优化计算。

1.3 团队名称

何雨川

1.4 结果展示

我们使用了一个C++的矩阵乘法函数来和DPC++的矩阵乘法实现函数进行对比实验，如果输出“Pass”字样，代表结果一致，即DPC++算法实现是正确的。经过我们的实验，程序输出了正确的结果。

2：算法具体实现

2.1 分配内存和数据传输

我们使用了DPC++中的设备队列queue，队列策略采用了in order属性防止冲突。分配的矩阵有4个，其中两个矩阵是作乘法的两个矩阵；另外两个矩阵是结果矩阵，分别用于储存在DPC++设备上执行乘法之后的结果和在普通主机上执行乘法之后的结果，4个矩阵均为1000 × 1000大小。

然后，为4个矩阵在主机上分配内存，然后为a和b矩阵（两个执行乘法的矩阵）使用随机数进行赋值：

然后为了实现高效计算，我们需要利用DPC++的设备进行处理。于是再在DPC++管理的设备上分配了三块内存，通过malloc_device函数实现，然后通过memcpy函数将主机上的a和b矩阵数据传输到设备上，为了保证数据传输安全，使用了wait函数进行等待：

2.2 并行结构

由于结果矩阵的大小为m × k，使用无符号整型数据结构保存线程网格的列大小和行大小，大小如图所示：

其中BLOCK_SIZE是一个宏，大小设置为16。然后我们采用sycl命名空间里的range数据结构保存线程网格和线程块的规格，range定义并行调度中单个工作组的迭代域，或调度的总体维度。在这里将range的维度设置为3维：

然后，通过调用设备队列的parallel_for函数对核函数MartixM进行循环调用，采用ND-Range并行模型：

2.3 核函数算法

核函数MartixM的形参分别为：3个设备上分配内存的矩阵指针、矩阵的大小（m、n、k，均为1000）以及一个循环的句柄item_ct1。核函数的算法如下：

可以看出核函数的每一个线程最终计算结果矩阵里的一个值。线程通过get_group确定工作组id，通过get_local_range函数确定子组大小，二者相乘确定偏移量大小；通过get_local_id确定线程在子组中的id，与前面的偏移量大小相加确定全局id。

2.4 对比实验

核函数计算完毕后，结果矩阵的数据将从设备上拷贝到主机。这之后为了确定算法是否正确，需要和传统的C++矩阵乘法实现代码进行比较。将a、b矩阵和cpu_c矩阵（储存主机上矩阵相乘代码的执行结果）以及三个矩阵的大小传入cpu_matrix_mult函数中，计算得到结果后将cpu_c矩阵和c矩阵（存储从设备上传回的结果矩阵数据）进行比较，如果二者对应元素误差较大，则将标识位ok置为0：

如果ok不为0则表示算法正确，输出“Pass”；否则输出“error”。最后通过sycl空间里的free函数释放设备上分配的内存。

三、使用的oneAPI技术总结

使用了Intel oneAPI DPC++/C++ Compiler技术。

四、心得体会

使用DPC++可以很方便地将任务通过线程部署到GPU或者其他异构设备上，它是一个十分强大的技术。

#include <sycl/sycl.hpp>
#include <dpct.hpp>
#include <iostream>
#include<math.h>

#define BLOCK_SIZE 16

void  MartixM(int* a, int* b, int* c, int m, int n, int k,
	const sycl::nd_item<3>& item_ct1) {
	int row = item_ct1.get_group(1) * item_ct1.get_local_range(1) +
		item_ct1.get_local_id(1);
	int col = item_ct1.get_group(2) * item_ct1.get_local_range(2) +
		item_ct1.get_local_id(2);
	int sum = 0;
	if (col < k && row < m)
	{
		for (int i = 0; i < n; i++)
		{
			sum += a[row * n + i] * b[i * k + col];
		}
		c[row * k + col] = sum;
	}

}

void cpu_matrix_mult(int* h_a, int* h_b, int* h_result, int m, int n, int k) {
	for (int i = 0; i < m; ++i)
	{
		for (int j = 0; j < k; ++j)
		{
			int tmp = 0;
			for (int h = 0; h < n; ++h)
			{
				tmp += h_a[i * n + h] * h_b[h * k + j];
			}
			h_result[i * k + j] = tmp;
		}
	}
}

int main() {
	sycl::device dev_ct1;
	sycl::queue q_ct1(dev_ct1,
		sycl::property_list{ sycl::property::queue::in_order() });
	int m = 1000,n = 1000, k = 1000;
	int ok = 1;

	int* a = (int*)malloc(sizeof(int) * m * n);
	int* b = (int*)malloc(sizeof(int) * n * k);
	int* c = (int*)malloc(sizeof(int) * m * k);
	int* cpu_c = (int*)malloc(sizeof(int) * m * k);

	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			a[j + i * n] = rand() % 1024;
		}
	}

	for (int i = 0; i < n; i++) {
		for (int j = 0; j < k; j++) {
			b[j + i * k] = rand() % 1024;
		}
	}

	int* d_a, * d_b, * d_c;
	d_a = (int*)sycl::malloc_device(sizeof(int) * m * n, q_ct1);
	d_b = (int*)sycl::malloc_device(sizeof(int) * n * k, q_ct1);
	d_c = (int*)sycl::malloc_device(sizeof(int) * m * k, q_ct1);

	q_ct1.memcpy(d_a, a, sizeof(int) * m * n).wait();
	q_ct1.memcpy(d_b, b, sizeof(int) * n * k).wait();

	unsigned int grid_rows = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;
	unsigned int grid_cols = (k + BLOCK_SIZE - 1) / BLOCK_SIZE;
	sycl::range<3> dimGrid(1, grid_rows, grid_cols);
	sycl::range<3> dimBlock(1, BLOCK_SIZE, BLOCK_SIZE);


	q_ct1.parallel_for(sycl::nd_range<3>(dimGrid * dimBlock, dimBlock),
		[=](sycl::nd_item<3> item_ct1) {
			MartixM(d_a, d_b, d_c, m, n, k, item_ct1);
		});
	q_ct1.memcpy(c, d_c, sizeof(int) * m * k).wait();

	cpu_matrix_mult(a, b, cpu_c, m, n, k);
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < k; j++) {
			if (fabs(cpu_c[i * k + j] - c[i * k + j]) > (1.0e-10)) {
				ok = 0;
			}
		}
	}

	if (ok) {
		std::cout << "Pass" << std::endl;
	}
	else {
		std::cout << "error" << std::endl;
	}

	sycl::free(d_a, q_ct1);
	sycl::free(d_b, q_ct1);
	sycl::free(d_c, q_ct1);

	return 0;
}