基于Intel SYCL的矩阵乘法

1. 问题阐述

任务1：熟悉SYCL编程环境 

       登录SYCL编程环境，支持oneAPI的英特尔或英伟达图形图像处理器（简称GPU），或支持英伟达图像处理器的Codeplay插件(详见https://developer.codeplay.com/products/oneapi/nvidia)， 确保必需的编译器和函数库都已安装。 

任务2：矩阵乘法运算 

       按以下步骤写一个SYCL程序来实现包含两个矩阵的矩阵乘法运算： 

组建2个序列（向量）的浮点数，每个序列的规格是N（如N=1024*1024），构成矩阵输入值。用随机值初始化序列。使用缓存和存储来实现对设备（GPU）上矩阵内存的分配并运行。运行SYCL Kernel实现两个矩阵的并行运算，这里你需要运用SYCL nd_range概念来定义Kernel的运行范围。使用SYCL排队系统来运行设备上的Kernel。 Kernel运行结束后，使用存储将结果从设备（GPU）检索回主机。 

通过比较主机上的预设结果来验证你的程序是否正确。 

任务3：统一共享内存（简称“USM”） 

        统一共享内存在主机和设备（GPU）之间提供了一个统一的存储模式。修改你的程序以使用统一共享内存来实现内存分配和数据转换，从而替代缓存和存储。如果使用不同种类的USM将获得加分。重写代码以使用统一共享内存来实现内存分配和数据转换（如需要），照常运行Kernel来验证程序运行结果。 

任务4：程序运行结果分析（视运行结果而加分） 

       当使用缓存/存储和统一共享内存（简称“USM”）（或不同种类的USM）时，你可以通过测量或比较程序运行时间来分析或讨论程序运行结果可能存在的任何差异。

2. 软硬件及开发环境

直接使用英特尔oneAPI Developer Cloud 服务，直接利用Develop Cloud平台中的CPU和GPU硬件，启动JupyterLab递交代码到计算节点运行。平台已经提供了DPC++ Compiler的编译环境，也模拟了一些硬件。

3. 主要实现方案

使用缓存和存储进行矩阵乘法运算

准备阶段：设置矩阵的大小以及block_size的大小，在主机端定义浮点类型容器分别为矩阵A的输入，矩阵B的输入，通过GPU计算后的结果，通过CPU计算后的结果。为矩阵A，B用赋随机值，也为GPU和CPU计算结果分别准备容器。

进入外部计算加速器阶段：创建SYCL设备选择队列，这里先选择GPU。创建输入输出缓冲区二维的buffer，将主机上的存储数据放在buffer中，矩阵A和B的buffer设置为read_only模式，结果矩阵部分buffer设置为write_only。队列Q提交指令任务，这里要通过中间对象handler &cgh，以描述runtime对象的并行信息和计算顺序。然后为每个buffer创造进程管理器accessor，通过与handler的联系将buffer的数据传入设备。下面就是基于accessor进行的在设备（这里是GPU）中的并行计算，并行计算函数在

handler.parallel_for({gloabal,local},[=](id)){

//code

})中进行实际的并行计算，最后队列Q.wait()等待命令组完成。

核心代码如下：

//创建SYCL设备序列,选择设备GPU

    sycl::queue deviceQueue(sycl::gpu_selector_v);

    //创建输入和输出缓存区，# Create Buffer

    buffer<float, 2> a_buf(vector1.data(), range<2>(N, N)), b_buf(vector2.data(), range<2>(N, N)), c_buf(vector3.data(),

            range<2>(N, N));

    //提交一个SYCL命令组

    deviceQueue.submit([&](sycl::handler & cgh) {

        //定义访问器以访问输入和输出缓冲区

        accessor a{a_buf, cgh};

        accessor b{b_buf, cgh};

        accessor c{c_buf, cgh};



        range global{N, N};

        range local{B, B};

        //定义内核函数来执行矩阵乘法

        cgh.parallel_for(nd_range{global, local}, [ = ](nd_item<2> it) {

           int j = it.get_global_id(0);

           int i = it.get_global_id(1);

           for (int k = 0; k < N; k++) {

               c[j][i] += a[j][k] * b[k][i];

           }

        });

    });

    //等待命令组完成

    deviceQueue.wait();

数据回收阶段：使用host_accessor(buffer,read_only)将设备端的buffer中的数据返回到主机上预先开辟的结果矩阵空间中

host_accessor h_a(c_buf, read_only);

并行计算部分（Kernel函数）：这一部分有主要的是三个类，range（计算范围）、id（所属索引）、item（计算点位）由于本实验是矩阵计算，维度为2，这里使用nd_range。并行计算的范围有分global_size和work_group_size，类似于将一个整体的计算范围又划分成了更小的计算单元以加速计算。代码中id(0)表示行方向（多行可并行计算），id(1)表示列方向（多列可并行计算），根据矩阵乘法行乘列的方法再进行累和。核心代码如下：

cgh.parallel_for(nd_range{global, local}, [ = ](nd_item<2> it) {

           int j = it.get_global_id(0);

           int i = it.get_global_id(1);

           for (int k = 0; k < N; k++) {

               c[j][i] += a[j][k] * b[k][i];

           }

        });

性能测试部分：这里主要通过计算不同设备的计算时间进行比较。使用C++，chrono类std::chrono::high_resolution_clock::time_point s, e;

    //开始计时，在队列提交SYCL命令组之前（.submmit）

    s = std::chrono::high_resolution_clock::now();

//结束计时，在队列等待命令组完成后（.wait）

    e = std::chrono::high_resolution_clock::now();

最后通过e-s，及单位转换为毫秒输出GPU的计算时间

性能对比部分：对于上述的从“进入外部加速器阶段“开始到最后运行过程全部相同，唯一的改动在于，在选择设备时选择了CPU。最后将GPU和CPU两个设备的计算时间进行了对比。

结果检验部分：对于准备阶段的两个结果容器进行逐个数据检验，差距在0.001以上的视为错误

使用统一共享存储（USM）

准备阶段：在主机端malloc方法开辟内存空间（输入矩阵A、B和GPU结果矩阵C，CPU结果矩阵C_C），对于两个输入矩阵的计算空间进行随机值赋值。

代码如下：

//开辟主机内存空间

    float *A = static_cast<float *>(malloc(N * N * sizeof(float)));

    float *B = static_cast<float *>(malloc(N * N * sizeof(float)));

    float *C = static_cast<float *>(malloc(N * N * sizeof(float)));

    float *C_C = static_cast<float *>(malloc(N * N * sizeof(float)));

    进入外部设备加速器（这里是GPU），创建设备队列（设备选择器选择GPU），同样使用malloc方法开辟设备内存(a_g,b_g,c_g)将主机上的数据复制到设备队列上面

    代码如下:   

//开辟设备内存

    float *a_g = malloc_device<float>(N * N, q);

    float *b_g = malloc_device<float>(N * N, q);

    float *c_g = malloc_device<float>(N * N, q);

    //把主机上的数据复制到设备上面

    q.memcpy(a_g, A, sizeof(float)*N * N).wait();

    q.memcpy(b_g, B, sizeof(float)*N * N).wait();

    q.memcpy(c_g, C, sizeof(float)*N * N).wait();

再运行内核函数，内核函数内容与缓存时的方法及其类似，也是定位到行和列再进行累和。内核函数代码如下：

//定义内核函数来执行矩阵乘法

    q.parallel_for(nd_range{global, local}, [ = ](nd_item<2> it) {

        int j = it.get_global_id(0);

        int i = it.get_global_id(1);



        float sum = 0.0f;

        for (int x = 0; x < N; x++) {

           sum += a_g[j * N + i] * b_g[i * N  + i];

        }

        c_g[j * N + i] = sum;



    }).wait();

数据返回阶段：最后将计算结果返回到主机上

q.memcpy(C, c_g, sizeof(float)*N * N).wait();

性能测试阶段与上述相同也采用chrono计时的方式,通过计算不同设备的计算时间进行比较。使用C++，chrono类,代码如下：

std::chrono::high_resolution_clock::time_point s, e;

    //开始计时，在队列提交SYCL命令组之前（.submmit）

    s = std::chrono::high_resolution_clock::now();

//结束计时，在队列等待命令组完成后（.wait）

    e = std::chrono::high_resolution_clock::now();

最后通过e-s，及单位转换为毫秒输出GPU的计算时间

性能对比部分：对于上述的从“进入外部加速器阶段“开始到最后运行过程全部相同，唯一的改动在于，在选择设备时选择了CPU。最后将GPU和CPU两个设备的计算时间进行了对比。

结果检验部分：对于准备阶段的两个结果空间进行逐个数据检验，差距在0.001以上的视为错误。

4. 系统运行结果及结果分析

将代码复制到DevCloud平台，启动JupyterLab递交代码到计算节点运行

对于缓存方式运行结果如图：

设备选择GPU，计算时间为314.428ms

设备选择CPU，计算时间为869.353ms

对准备阶段中储存两个结果的vector容器进行结果检验，两次计算结果相同

对于统一共享存储运行结果如图：

设备选择GPU，计算时间为76.2384ms

设备选择CPU，计算时间为483.643ms

对准备阶段中两个内存块进行结果检验，两次计算结果相同

对于两种不同的存储模式，本次实验结果显示统一共享内存的方式计算结果更快。个人分析，本次实验计算任务需要多个处理器协同工作，并且需要访问共享的内存数据，所以从结果上统一共享内存可能会更快。

5. 过程中的困难与问题

在缓存存储模式中，在设备端计算结束后一定要用host_accessor返回将设备缓冲区数据返回到主机端，这样才能有效输出。

在对缓冲区进行定义时，要定义其类型，read_only还是write_only，从实验过程中的测试结果来看这一过程是必需的。

在统一共享内存的模式中，在对于设备端用malloc方法进行空间开辟后，该空间无法进行正常赋值（会导致程序中断），赋值操作应提前到开辟主机内存部分，后再将该赋值数据复制到设备端。

6. 参考与收获

https://www.bilibili.com/video/BV1ku4y1s7za

团队通过学习蔺杰教授有关intel oneAPI的讲解，了解了oneAPI的目标，并行计算的发展现状，intel目前的硬件加速器产品以及sycl的基本编程包括基本流程，实现模式，程序架构，一些基础类。

本实验还有继续改进提高之处，也正在研究使用其他硬件设备进行实验获得更多的对比数据。同时还在探寻USM的其他实现方法。