Intel® DevCloud for oneAPI SYCL编程项目实践

问题陈述

实验所用的硬件环境和软件环境

本次实验使用oneAPI中支持SYCL编程模型的C++编译器，使用英特尔oneAPI Developer Cloud服务，可以免安装额外环境，利用CPU作为主机（Host），同时利用GPU作为设备（Device）完成作业

通过队列q输出GPU设备名称：

Intel® DevCloud编程环境注册

注册登录Intel账号

进入如图位置，往下拉到最底

连接并进入jupyter lab环境

进入箭头所指目录

在此目录可以学习sycl和oneapi相关的官方文档

本次实验的环境使用oneAPI_Essentials/02_SYCL_Program_Structure/目录中SYCL_Program_Structure.ipynb的Lab Exercise:Vector Add实验，在这个位置可以直接写自己的代码并运行，无需额外安装操作

如何获取运行结果

在上述Vector Add的cell中填写代码

运行

此时代码会保存在另一个位置的cpp文件中：

再运行下面这个cell即可

在Dev C++平台上面运行矩阵乘法

代码如下：

#include<vector>
#include <chrono>
#include <iostream>
using namespace std;
//cpu上的计算 
double cpu_kernel(std::vector<float> &A, std::vector<float> &B, std::vector<float> &C, int M, int N, int K) {double duration = 0.0;std::chrono::high_resolution_clock::time_point s, e;s = std::chrono::high_resolution_clock::now();for(int i = 0; i < M; i++) {for(int j = 0; j < N; j++) {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;}}e = std::chrono::high_resolution_clock::now();duration = std::chrono::duration<float, std::milli>(e - s).count();return duration;
}
void compute(const int M,const int N,const int K,const int block_size,const int iterations) {//matrix A:(1024*1024)  matrix B:(1024*1024)  matrix C=A*B(1024*1024)cout << "C矩阵尺寸： c(" << M << "," <<  N << ") ="<< "A矩阵 a(" << M << "," << K << ") *"<< " B矩阵b(" << K << "," << N << ")\n";std::vector<float> A(N*K);std::vector<float> B(K*M);std::vector<float> C(N*M);for(int i=0; i < M * K; i++) {A[i] = rand() / double(RAND_MAX);}for(int i=0; i < K * N; i++) {B[i] = rand() / double(RAND_MAX);}for(int i=0; i < M * N; i++) {C[i] = 0.0f;}//CPU平均执行时间 double duration_cpu = 0.0f;int warmup = 2;for(int run = 0; run < iterations + warmup; run++) {float duration = cpu_kernel(A, B, C, M, N, K);if(run >= warmup) duration_cpu += duration;}duration_cpu = duration_cpu / iterations;//运行时间输出printf("\nCPU计算时间：%lf (ms); \n",duration_cpu);
}
int main() {compute(1024,1024,1024, 4, 5);return 0;
}

运行结果

使用缓存存储实现矩阵乘法

在buffer模型中，使用缓存存储将结果从设备（GPU）检索回主机有两个方法，第一个是使用host_accessor，第二个是利用buffer destruction机制

buffer_destruction机制：缓冲区创建发生在单独的函数作用域中。当执行超出此函数作用域时，调用缓冲区析构函数，从而放弃数据的所有权并将数据复制回主机内存

完整代码以及注释：

#include <chrono>
#include<vector>
#include <iostream>
#include <CL/sycl.hpp>
using namespace std;
using namespace sycl;
//gpu上的计算 
double gpu_kernel(std::vector<float> &A, std::vector<float> &B, std::vector<float> &C, int M, int N, int K, int block_size, sycl::queue &q) {//确保行数和列数被块大小block_size整除auto grid_rows = (M + block_size - 1) / block_size * block_size;auto grid_cols = (N + block_size - 1) / block_size * block_size;//局部范围决定了每个工作组的大小auto local_ndrange  = range<2>(block_size, block_size);//全局范围定义了整个计算的范围 auto global_ndrange = range<2>(grid_rows, grid_cols);//缓冲区 buffer buf1(A);buffer buf2(B);buffer buf3(C);double duration = 0.0f;auto e = q.submit([&](sycl::handler &h) {//访问器 accessor A1(buf1,h);accessor B1(buf2,h);accessor C1(buf3,h);h.parallel_for<class k_name_t>(sycl::nd_range<2>(global_ndrange, local_ndrange), [=](sycl::nd_item<2> index) {//通过index获取行列号 int row = index.get_global_id(0);int col = index.get_global_id(1);float sum = 0.0f;for (int i = 0; i < K; i++) {sum += A1[row * K + i] * B1[i * N  + col];}C1[row * N + col] = sum;});});e.wait();duration += (e.get_profiling_info<info::event_profiling::command_end>() -e.get_profiling_info<info::event_profiling::command_start>()) /1000.0f/1000.0f;return duration;
}
//cpu上的计算 
double cpu_kernel(std::vector<float> &A, std::vector<float> &B, std::vector<float> &C, int M, int N, int K) {double duration = 0.0;std::chrono::high_resolution_clock::time_point s, e;s = std::chrono::high_resolution_clock::now();for(int i = 0; i < M; i++) {for(int j = 0; j < N; j++) {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;}}e = std::chrono::high_resolution_clock::now();duration = std::chrono::duration<float, std::milli>(e - s).count();return duration;
}
//验证cpu计算结果和gpu计算结果的偏差 
int verify(std::vector<float> &cpu_res, std::vector<float> &gpu_res, int length){int err = 0;for(int i = 0; i < length; i++) {if( fabs(cpu_res[i] - gpu_res[i]) > 1e-3) {err++;printf("\n%lf, %lf", cpu_res[i], gpu_res[i]);}}return err;
}
int compute(const int M,const int N,const int K,const int block_size,const int iterations,sycl::queue &q) {//matrix A:(1024*1024)  matrix B:(1024*1024)  matrix C=A*B(1024*1024)cout << "C矩阵尺寸： c(" << M << "," <<  N << ") ="<< "A矩阵 a(" << M << "," << K << ") *"<< " B矩阵b(" << K << "," << N << ")\n";std::vector<float> A(N*K);std::vector<float> B(K*M);std::vector<float> C(N*M); std::vector<float> C2(N*M);for(int i=0; i < M * K; i++) {A[i] = rand() / double(RAND_MAX);}for(int i=0; i < K * N; i++) {B[i] = rand() / double(RAND_MAX);}for(int i=0; i < M * N; i++) {C[i] = 0.0f;C2[i] = 0.0f;}//GPU平均执行时间 double duration_gpu = 0.0f;//CPU平均执行时间 double duration_cpu = 0.0f;// 代码“预热”次数 ，即总共计算20次，但只统计11--20次的时候的平均运行时间 int warmup = 10;for (int run = 0; run < iterations + warmup; run++) {float duration = gpu_kernel(A, B, C, M, N, K, block_size, q);if(run >= warmup) duration_gpu += duration;}//计算gpu平均运行时间 duration_gpu = duration_gpu / iterations;//计算cpu平均运行时间，cpu运行较慢，所以取iteration值为原来一半 warmup = 2;for(int run = 0; run < iterations/2 + warmup; run++) {float duration = cpu_kernel(A, B, C2, M, N, K);if(run >= warmup) duration_cpu += duration;}duration_cpu = duration_cpu / iterations/2;//错误结果验证int errCode = 0;errCode = verify(C2, C, M*N);printf("\n总共有%d处错误\n", errCode);//运行时间输出printf("\n1\nGPU计算时间：%lf (ms); \nCPU计算时间：%lf (ms); \n",duration_gpu, duration_cpu);return errCode;
}
int main() {//enable_profiling() 是一个SYCL属性（property），用于在SYCL队列上启用性能分析auto propList = cl::sycl::property_list {cl::sycl::property::queue::enable_profiling()};//队列q 使用 gpu设备选择器 queue q( cl::sycl::gpu_selector_v , propList);int errCode = compute(1024,1024,1024, 4, 10, q);return errCode;
}

代码解释

本代码除了main函数之外，还包含四个函数，分别是gpu_kernel、cpu_kernel、verify、compute

其中gpu_kernel负责矩阵在gpu上的乘法运算和计算运算时间的功能，cpu_kernel负责矩阵在cpu上的乘法运算和计算运算时间的功能，verify负责矩阵在gpu和cpu上运算后的结果的对比，compute负责矩阵随机初始化、前三个函数的调用、输出结果等等多个功能

在compute函数当中，使用了两个参数iteration和warmup

前者iteration是取gpu/cpu的iteration次的运算平均时间值，warmup起到“热身”作用，比如当iteration=10，warmup=10时，gpu会运算20次，而最终结果取11--20次计算时间的平均值

运行结果

使用USM显式数据移动实现矩阵乘法

数据的显式移动是指使用malloc_device的USM实现，其中主机和设备之间的数据移动可以由开发人员使用memcpy显式地完成

下面展示使用这种思想的矩阵乘法实现代码:

#include <chrono>
#include <iostream>
#include <CL/sycl.hpp>
using namespace std;
using namespace sycl;
//验证cpu计算结果和gpu计算结果的偏差 
int verify(float *cpu_res, float *gpu_res, int length){int err = 0;for(int i = 0; i < length; i++) {if( fabs(cpu_res[i] - gpu_res[i]) > 1e-3) {err++;printf("\n%lf, %lf", cpu_res[i], gpu_res[i]);}}return(err);
}
int main() {//enable_profiling() 是一个SYCL属性（property），用于在SYCL队列上启用性能分析auto propList = cl::sycl::property_list {cl::sycl::property::queue::enable_profiling()};//队列q 使用 gpu设备选择器 queue q( cl::sycl::gpu_selector_v , propList);int block_size=4; int M=1024;int N=1024;int K=1024;//gpu计算//确保行数和列数被块大小block_size整除auto grid_rows = (M + block_size - 1) / block_size * block_size;auto grid_cols = (N + block_size - 1) / block_size * block_size;//局部范围决定了每个工作组的大小auto local_ndrange  = range<2>(block_size, block_size);//全局范围定义了整个计算的范围 auto global_ndrange = range<2>(grid_rows, grid_cols);//主机上的数据 float *A=static_cast<float *>(malloc(M*K*sizeof(float)));float *B=static_cast<float *>(malloc(K*N*sizeof(float)));float *C=static_cast<float *>(malloc(M*N*sizeof(float)));float *C2=static_cast<float *>(malloc(M*N*sizeof(float)));for(int i=0; i < M * K; i++) {A[i] = rand() / double(RAND_MAX);}for(int i=0; i < K * N; i++) {B[i] = rand() / double(RAND_MAX);}for(int i=0; i < M * N; i++) {C[i] = 0.0f;C2[i] = 0.0f;}//malloc_device:数据可以在设备和主机显示地迁移 float *A_device=malloc_device<float>(M*K,q);float *B_device=malloc_device<float>(K*N,q);float *C_device=malloc_device<float>(M*N,q);//把主机上的数据复制到设备上面q.memcpy(A_device,A,sizeof(float)*M*K).wait(); q.memcpy(B_device,B,sizeof(float)*K*N).wait(); q.memcpy(C_device,C,sizeof(float)*M*N).wait(); double duration=0.0;for(int i=0;i<20;i++){std::chrono::high_resolution_clock::time_point s, e;s = std::chrono::high_resolution_clock::now();auto e1 = q.submit([&](sycl::handler &h) {h.parallel_for<class k_name_t>(sycl::nd_range<2>(global_ndrange, local_ndrange), [=](sycl::nd_item<2> index) {//通过index获取行列号 int row = index.get_global_id(0);int col = index.get_global_id(1);float sum = 0.0f;for (int i = 0; i < K; i++) {sum += A_device[row * K + i] * B_device[i * N  + col];}C_device[row * N + col] = sum;});});e1.wait();//将处理后的、位于设备上的C矩阵数据复制到主机上的C中q.memcpy(C,C_device,sizeof(float)*M*N).wait();e = std::chrono::high_resolution_clock::now();if(i>=10)duration = duration+std::chrono::duration<float, std::milli>(e - s).count();}printf("\nGPU计算时间：%lf (ms);  \n",duration/10.0);duration = 0.0;for(int i=0;i<7;i++){std::chrono::high_resolution_clock::time_point s, e;s = std::chrono::high_resolution_clock::now();for(int i = 0; i < M; i++) {for(int j = 0; j < N; j++) {float sum = 0.0f;for(int k = 0; k < K; k++) {sum +=  A[i * K + k] * B[k * N  + j];}C2[i * N + j] = sum;}}e = std::chrono::high_resolution_clock::now();if(i>=2)duration = duration+std::chrono::duration<float, std::milli>(e - s).count();}printf("\nCPU计算时间：%lf (ms);  \n",duration/5.0);//错误结果验证int errCode = 0;errCode = verify(C2, C, M*N);printf("\n总共有%d处错误\n", errCode);free(A_device,q);free(B_device,q);free(C_device,q);free(A);free(B);free(C);return errCode;
}

代码解释：

这个代码使用USM显式数据移动方案实现矩阵乘法，首先分别创建主机上的内存再使用malloc_device创建设备上的内存，在主机上的矩阵内存被随机赋值之后，通过队列q的memcpy把主机上的数据复制到设备的内存中，让这些数据通过队列提交到加速器上处理之后再重新赋值给主机的内存中，最后通过verify进行结果验证

这个代码在CPU和GPU计算之前，也先进行了“代码预热”，在CPU跑够2轮、GPU跑够10轮之后，再让CPU跑5轮，取时间平均值、让GPU跑10轮，取时间平均值

运行结果：

使用USM隐式数据移动实现矩阵乘法

USM隐式数据移动是指使用malloc_host、malloc_shared进行内存的分配和管理，下面演示使用这种思想实现矩阵乘法的代码：

#include <chrono>
#include <iostream>
#include <CL/sycl.hpp>
using namespace std;
using namespace sycl;
//gpu上的计算 
double gpu_kernel(float *A, float *B, float *C, int M, int N, int K, int block_size, sycl::queue &q) {//确保行数和列数被块大小block_size整除auto grid_rows = (M + block_size - 1) / block_size * block_size;auto grid_cols = (N + block_size - 1) / block_size * block_size;//局部范围决定了每个工作组的大小auto local_ndrange  = range<2>(block_size, block_size);//全局范围定义了整个计算的范围 auto global_ndrange = range<2>(grid_rows, grid_cols);double duration = 0.0f;auto e = q.submit([&](sycl::handler &h) {h.parallel_for<class k_name_t>(sycl::nd_range<2>(global_ndrange, local_ndrange), [=](sycl::nd_item<2> index) {//通过index获取行列号 int row = index.get_global_id(0);int col = index.get_global_id(1);float sum = 0.0f;for (int i = 0; i < K; i++) {sum += A[row * K + i] * B[i * N  + col];}C[row * N + col] = sum;});});e.wait();duration += (e.get_profiling_info<info::event_profiling::command_end>() -e.get_profiling_info<info::event_profiling::command_start>()) /1000.0f/1000.0f;return(duration);
}
//cpu上的计算 
double cpu_kernel(float *A, float *B, float *C, int M, int N, int K) {double duration = 0.0;std::chrono::high_resolution_clock::time_point s, e;s = std::chrono::high_resolution_clock::now();for(int i = 0; i < M; i++) {for(int j = 0; j < N; j++) {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;}}e = std::chrono::high_resolution_clock::now();duration = std::chrono::duration<float, std::milli>(e - s).count();return(duration);
}
//验证cpu计算结果和gpu计算结果的偏差 
int verify(float *cpu_res, float *gpu_res, int length){int err = 0;for(int i = 0; i < length; i++) {if( fabs(cpu_res[i] - gpu_res[i]) > 1e-3) {err++;printf("\n%lf, %lf", cpu_res[i], gpu_res[i]);}}return(err);
}
int compute(const int M,const int N,const int K,const int block_size,const int iterations,sycl::queue &q) {//matrix A:(1024*1024)  matrix B:(1024*1024)  matrix C=A*B(1024*1024)cout << "C矩阵尺寸： c(" << M << "," <<  N << ") ="<< "A矩阵 a(" << M << "," << K << ") *"<< " B矩阵b(" << K << "," << N << ")\n";//malloc_shared:数据可以在设备和主机隐式地迁移 auto A = malloc_shared<float>(M * K, q);auto B = malloc_shared<float>(K * N, q);auto C = malloc_shared<float>(M * N, q);//malloc_host: 数据显式地分配在主机上 auto C2 = malloc_host<float>(M * N, q);for(int i=0; i < M * K; i++) {A[i] = rand() / double(RAND_MAX);}for(int i=0; i < K * N; i++) {B[i] = rand() / double(RAND_MAX);}for(int i=0; i < M * N; i++) {C[i] = 0.0f;C2[i] = 0.0f;}//GPU平均执行时间 double duration_gpu = 0.0f;//CPU平均执行时间 double duration_cpu = 0.0f;// 代码“预热”次数 ，即总共计算20次，但只统计11--20次的时候的平均运行时间 int warmup = 10;for (int run = 0; run < iterations + warmup; run++) {float duration = gpu_kernel(A, B, C, M, N, K, block_size, q);if(run >= warmup) duration_gpu += duration;}//计算gpu平均运行时间 duration_gpu = duration_gpu / iterations;//计算cpu平均运行时间，cpu运行较慢，所以取iteration值为原来一半 warmup = 2;for(int run = 0; run < iterations/2 + warmup; run++) {float duration = cpu_kernel(A, B, C2, M, N, K);if(run >= warmup) duration_cpu += duration;}duration_cpu = duration_cpu / iterations/2;//错误结果验证int errCode = 0;errCode = verify(C_host, C, M*N);printf("\n总共有%d处错误\n", errCode);//运行时间输出printf("\n\nGPU计算时间：%lf (ms); \nCPU计算时间：%lf (ms); \n",duration_gpu, duration_cpu);//释放内存free(A, q);free(B, q);free(C, q);free(C_host, q);return errCode;
}
int main() {//enable_profiling() 是一个SYCL属性（property），用于在SYCL队列上启用性能分析auto propList = cl::sycl::property_list {cl::sycl::property::queue::enable_profiling()};//队列q 使用 gpu设备选择器 queue q( cl::sycl::gpu_selector_v , propList);int errCode = compute(1024,1024,1024, 4, 10, q);return errCode;
}

代码解释：

和第一个代码一样，这个代码除了main函数之外，还包含四个函数，分别是gpu_kernel、cpu_kernel、verify、compute

不同的是，这个代码使用molloc_shared和malloc_host实现了统一共享内存

运行结果：

几种运行结果的对比分析

通过上面几种不同的运算，可以发现devcloud平台具有明显的高性能优势，除此之外，在devcloud的三种方式当中，可以总结出USM的优势如下：

SYCL*缓冲区功能很强大，但是，在c++程序中用缓冲区替换所有指针和数组可能会给程序员带来负担，因此在这种情况下可以考虑使用统一共享内存USM

1.当把c++代码移植到sycl时，想要尽可能更改少的代码

2.当需要控制数据移动时，使用显式USM分配

项目遇到的挑战，解决措施与心得体会

SYCL作为一个异构编程框架，对于从未了解过的初学者来说，了解SYCL的概念、编程模型方法是一项挑战，为了解决这个问题，我认真学习了老师提供的音频ppt以及oneapi官方文档的前四章，主要包括三个方面：对oneapi的介绍、sycl编程框架和sycl的统一共享内存

为了加深对oneAPI的理解，针对上述的学习内容，自己写出了四个总结性csdn博客，链接如下：

oneAPI学习笔记一：SYCL简介

oneAPI学习笔记二：jupyter官方文档（oneAPI_Intro）学习笔记

oneAPI学习笔记三：jupyter官方文档（SYCL Program Structure）学习笔记

oneAPI学习笔记四：jupyter官方文档（Unified Shared Memory）学习笔记

尽管如此，对于sycl的学习仍处于入门阶段，依然面临着诸多挑战，包括但不限于处理不同硬件兼容性的挑战、如何在不同设备上优化代码的挑战和针对一些特定于异构编程的错误处理的挑战等

总而言之，本次oneAPI实践让我接触到了Intel公司前沿的技术，受益匪浅。进行SYCL实验是具有挑战性但也是富有成就感的过程，在这个过程当中，提升了我应对处理一个陌生问题的能力，同时锻炼了我阅读学习英文官方技术文档的能力。除此之外，接触了之前没有听说过的词汇比如高性能计算（HPC）、异构编程等等。随着时间和实践的累积，希望能够更加熟练地应用SYCL，并从中获得更多的收获和经验