描述
编写一个基于oneAPI的C++/SYCL程序来执行矩阵乘法操作。需要考虑大尺寸矩阵的乘法操作以及不同线程之间的数据依赖关系。通常在实现矩阵乘法时，可以使用块矩阵乘法以及共享内存来提高计算效率。
环境配置
1、访问英特尔统一登录账户注册页面 https://idzcn.com/devcloud.htm 开始注册一个新的英特尔用户账号。点击“创建账户”开始创建一个全新的英特尔账号。 填写基本联络信息，电子邮件、姓氏、名字信息、密码等。完成后，点击“下一步：验证您的电子邮件”按钮递交表单。 发送验证码并验证完毕后，点击“创建账户”完成账号创建。
2、通过访问 https://devcloud.intel.com/oneapi/home/ 页面中，“Get Free Access” 再次登录账号后，继续相关信息的填写。填写完成“Business or Institution Name”中企业或学校机构名称，并在下一个拉菜单“What type of user are you?”中选择您的用户类型。点击“Next: Terms and Conditions”递交表单。
3、阅读 Intel® Developer Cloud 相关服务条款内容，点击“I Accept”进入下一步。选择感兴趣的英特尔产品及技术领域，订阅相关技术/产品更新或简报（可选项）。点击“Next: Submit”进入下一步。核验确认相关的信息，如确认无误，点击“Submit”进入下一步。系统显示“Thank you for applying for Intel® Developer Cloud access” 则表示账号注册已完成。
4、我们实验使用的硬件环境需要通过访问 https://devcloud.intel.com/oneapi/home/ 页面中的远程云服务来进入 DevCloud for oneAPI 专项登录页面。点击页面右上角的“Sign In”或“Enroll” 或下方的“Sign In” 或“Get Free Access”链接，输入并完成登录及后续服务激活步骤。
5、勾选“I accept these terms”(我同意这些条款)，点击“Submit”进递交表单后，会在网页右上角显示账户信息及服务有效期。可以通过点击网页左侧“Get Started”进入该选项的页面，在“Get Started”选项页面中， 可以通过点击页面最左下角 Connect with Jupyter* Lab 中的“Launch JupyterLab*”按钮直接启动 Jupypter 服务。
矩阵乘法
主函数
由于要对矩阵乘法不断优化，所以选择写一个主函数，在主函数中记录运行时长，随后只需要为不同的方法实现矩阵乘法的函数即可。
#include <CL/sycl.hpp>
#include <ctime>
#include <chrono>
using namespace sycl;
// 不同的实现方法
void mm_kernel(queue &q, std::vector<float> &matrix_a, std::vector<float> &matrix_b, std::vector<float> &matrix_c, size_t N, size_t M);
// 验证计算结果
bool almost_equal(float a, float b){
float tolerance = 1e-6;
float diff = fabs(a - b);
a = fabs(a);
b = fabs(b);
float bigger = (b > a) ? b : a;
if(diff <= bigger * tolerance) return true;
return false;
}
int main(int argc, char *argv[]) {
size_t N = 1024;
size_t M = 16;
int VERIFY = 0;
int PRINT_OUTPUT_MATRIX = 0;
// 定义矩阵，a、b为输入，c为输出，d用于检测
std::vector<float> matrix_a(N*N);
std::vector<float> matrix_b(N*N);
std::vector<float> matrix_c(N*N);
std::vector<float> matrix_d(N*N);
// 初始化
float v1 = 2.f;
float v2 = 3.f;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++){
matrix_a[i*N+j] = v1++;
matrix_b[i*N+j] = v2++;
matrix_c[i*N+j] = 0.f;
matrix_d[i*N+j] = 0.f;
}
queue q(property::queue::enable_profiling{});
std::cout << "Offload Device : " << q.get_device().get_info<info::device::name>() << "\n";
std::cout << "max_work_group_size : " << q.get_device().get_info<info::device::max_work_group_size>() << "\n";
// 记录起始时间
auto start = std::chrono::high_resolution_clock::now().time_since_epoch().count();
// 调用不同方式的实现方法
mm_kernel(q, matrix_a, matrix_b, matrix_c, N, M);
// 计算运行时间
auto duration = std::chrono::high_resolution_clock::now().time_since_epoch().count() - start;
std::cout << "Compute Duration : " << duration / 1e+9 << " seconds\n";
// 输出结果
if (PRINT_OUTPUT_MATRIX){
for (int i=0; i<N; i++){
for (int j=0; j<N; j++){
std::cout << matrix_c[i*N+j] << " ";
}
std::cout << "\n";
}
} else {
std::cout << " [0][0] = " << matrix_c[0] << "\n";
}
// 与普通的乘法进行比较
if (VERIFY){
int fail = 0;
for(int i=0; i<N; i++){
for (int j = 0; j < N; j++) {
for(int k=0; k<N; k++){
matrix_d[i*N+j] += matrix_a[i*N+k] * matrix_b[k*N+j];
}
if(!almost_equal(matrix_c[i*N+j], matrix_d[i*N+j])) fail = 1;
}
}
if(fail == 1){
std::cout << "FAIL\n";
} else {
std::cout << "PASS\n";
}
}
return 0;
}
使用KML库
#include <CL/sycl.hpp>
#include "oneapi/mkl/blas.hpp" // oneMKL DPC++ 接口
using namespace sycl;
void mm_kernel(queue &q, std::vector<float> &matrix_a, std::vector<float> &matrix_b, std::vector<float> &matrix_c, size_t N, size_t M) {
std::cout << "Configuration : MATRIX_SIZE= " << N << "x" << N << "\n";
// 创建矩阵缓冲
buffer a(matrix_a);
buffer b(matrix_b);
buffer c(matrix_c);
// 乘法引子
float alpha = 1.f, beta = 1.f;
// 参与乘法运算的矩阵的转置
oneapi::mkl::transpose transA = oneapi::mkl::transpose::nontrans;
oneapi::mkl::transpose transB = oneapi::mkl::transpose::nontrans;
// 调用函数
oneapi::mkl::blas::gemm(q, transA, transB, N, N, N, alpha, b, N, a, N, beta, c, N);
c.get_access<access::mode::read>();
}
Basic Parallel Kernel Implementation
void mm_kernel(queue &q, std::vector<float> &matrix_a, std::vector<float> &matrix_b, std::vector<float> &matrix_c, size_t N, size_t M) {
std::cout << "Configuration : MATRIX_SIZE= " << N << "x" << N << "\n";
// 创建缓冲
buffer a(matrix_a);
buffer b(matrix_b);
buffer c(matrix_c);
// 提交指令
auto e = q.submit([&](handler &h){
// 获取缓冲数据
auto A = a.get_access<access::mode::read>(h);
auto B = b.get_access<access::mode::read>(h);
auto C = c.get_access<access::mode::write>(h);
// 遍历c的所有坐标
h.parallel_for(range<2>{N,N}, [=](item<2> item){
const int i = item.get_id(0);
const int j = item.get_id(1);
for (int k = 0; k < N; k++) {
C[i*N+j] += A[i*N+k] * B[k*N+j];
}
});
});
c.get_access<access::mode::read>();
// 统计运行时间
auto kernel_duration = (e.get_profiling_info<info::event_profiling::command_end>() - e.get_profiling_info<info::event_profiling::command_start>());
std::cout << "Kernel Execution Time : " << kernel_duration / 1e+9 << " seconds\n";
}ND-range
只需要将上述并行代码的for循环进行以下修改即可：
// 定义ND-Range和work-group的大小
range<2> global_size(N,N);
range<2> work_group_size(M,M);
// 以group为单位进行矩阵乘法
h.parallel_for(nd_range<2>{global_size, work_group_size}, [=](nd_item<2> item){
const int i = item.get_global_id(0);
const int j = item.get_global_id(1);
for (int k = 0; k < N; k++) {
C[i*N+j] += A[i*N+k] * B[k*N+j];
}
};
Local Memory
在ND-Range中的for循环中需要多次写入全局内存，所以可以考虑使用一个临时变量记录A和B矩阵work-group相乘的结果，在for循环执行结束后再写入C矩阵，这样就可以减少对全局内存的读写，也可以提高运行速度。
// 定义ND-Range和work-group的大小
range<2> global_size(N,N);
range<2> work_group_size(M,M);
// 以group为单位进行矩阵乘法
h.parallel_for(nd_range<2>{global_size, work_group_size}, [=](nd_item<2> item){
const int i = item.get_global_id(0);
const int j = item.get_global_id(1);
float temp = 0.f;
for (int k = 0; k < N; k++) {
temp += A[i*N+k] * B[k*N+j];
}
C[i*N+j] = c;
};
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
