利用oneapi进行矩阵相乘的并行优化

作业要求

并⾏矩阵乘法

描述
编写⼀个基于oneAPI的C++/SYCL程序来执行矩阵乘法操作。需要考虑大尺寸矩阵的乘法操作以及不同线程之
间的数据依赖关系。通常在实现矩阵乘法时，可以使用块矩阵乘法以及共享内存来提高计算效率。
分析
利用基于SYCL的编程模型在GPU上实现矩阵乘法的计算，步骤如下：

1. 分配内存：在主机端分配内存空间用于存储输⼊矩阵和输出矩阵，同时在GPU端分配内存空间用于存储相应的输入和输出数据。
2. 数据传输：将输入矩阵数据从主机端内存传输到GPU端内存中。
3. 核函数调用：在SYCL中，矩阵乘法的计算通常会在GPU上使用核函数来实现并行计算。核函数会分配线程块和线程来处理不同的数据块。
4. 并行计算：在核函数中，每个线程负责计算输出矩阵的⼀个单独的元素。为了最大限度地利用GPU的并行计算能力，通常会使用⼆维线程块和线程网格的方式来处理矩阵的乘法计算。
5. 数据传输：计算完成后，将输出矩阵数据从GPU端内存传输回主机端内存中，以便进⼀步处理或分析。在并行计算矩阵乘法时，可以利用线程块和线程的层次结构来优化计算。通过合理划分矩阵数据并利用共享内存来减少全局内存访问的次数，可以⼤幅提高计算效率。此外，还可以利用GPU上的多个计算单元并执行行矩阵乘法，进⼀步提高计算速度。

代码详情

首先导入oneapi的包

#include <sycl/sycl.hpp>
using namespace sycl;

设置矩阵的唯独大小，A是M*N的矩阵，B是N*Z的矩阵
W为工作组的大小

int M = 1024;
int N = 1024;
int Z = 1024;
int W = 16;

用一维数组存储二维矩阵 假设二维矩阵维度为M*N,则二维矩阵的元素A_2[x][y]存储在一维数组A_1[x*N+y]处

    std::vector<float> A(M\*N), B(N\*Z), C(M\*Z);
    std::cout << "A_MATRIX_SIZE= " << M << "x" << N "\n";
    std::cout << "B_MATRIX_SIZE= " << N << "x" << Z "\n";

通过主机向量创建缓冲区

缓冲区在SYCL应用程序中封装数据,以跨设备和主机

    buffer a(matrix_A);
    buffer b(matrix_B);
    buffer c(matrix_C);

一个队列提交命令组以在SYCL运行时执行

队列是一种机制，用于向设备提交工作。

    queue q;

    auto e = q.submit([&](handler &h){
        //......
    });

将三个数组读取到设备上

        accessor A(a, h, read_only);
        accessor B(b, h, read_only);
        accessor C(c, h, write_only);

ND-Range设置为结果矩阵的维度M*Z

        range<2> global_size(M,Z);
        range<2> work_group_size(W,W);

并行计算矩阵的主逻辑部分代码，每个C[x][y]并行计算

        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*Z+j] += A[i*N+k] * B[k*Z+j];
            }
        });

逐行打印结果矩阵C

    for(int i=0;i<M;i++){
        for(int j=0;j<Z;j++)
            std::cout<<C[i*Z+j]<<" ";
        std::cout<<std::endl;
    }

完整代码如下

#include <sycl/sycl.hpp>
using namespace sycl;

int M = 1024;
int N = 1024;
int Z = 1024;
int W = 16;

int main(){
    std::vector<float> A(M*N), B(N*Z), C(M*Z);
    std::cout << "A_MATRIX_SIZE= " << M << "x" << N "\n";
    std::cout << "B_MATRIX_SIZE= " << N << "x" << Z "\n";

    buffer a(matrix_A);
    buffer b(matrix_B);
    buffer c(matrix_C);

    queue q;

    auto e = q.submit([&](handler &h){
        accessor A(a, h, read_only);
        accessor B(b, h, read_only);
        accessor C(c, h, write_only);

        range<2> global_size(M,Z);
        range<2> work_group_size(W,W);

        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*Z+j] += A[i*N+k] * B[k*Z+j];
            }
        });
    });

    for(int i=0;i<M;i++){
        for(int j=0;j<Z;j++)
            std::cout<<C[i*Z+j]<<" ";
        std::cout<<std::endl;
    }
    return 0;
}

收获

学到了一些关于使用 SYCL 进行并行计算的基本概念和用法

1. 缓冲区和访问器：代码中使用了 buffer 来创建用于在设备和主机之间传递数据的缓冲区。通过 accessor，可以在内核中访问这些缓冲区的数据，分别用于只读（read_only）和写入（write_only）。

   buffer a(A);
   accessor A(a, h, read_only);
   

2. 队列和命令组：通过 queue 和 submit，代码将命令组提交到设备执行。这包含了并行计算的具体实现，其中包括了 ND-Range 的定义、工作组大小的设置以及并行计算的内核函数。

   queue q;
   auto e = q.submit([&](handler &h) {
       // 内核函数和设置
   });
   

3. 并行计算内核：内核函数使用 parallel_for 启动了一个并行计算，其中的工作项负责计算矩阵乘法的一部分。在这个内核函数中，使用了 nd_range 和 range 来定义全局 ND-Range 的大小和工作组大小。

   h.parallel_for(nd_range<2>{global_size, work_group_size}, [=](nd_item<2> item){
       // 并行计算的具体实现
   });
   

4. 结果输出：最后，代码在主函数中打印了结果矩阵 C，展示了并行计算的结果。

   for(int i = 0; i < M; i++){
       for(int j = 0; j < Z; j++)
           std::cout << C[i*Z+j] << " ";
       std::cout << std::endl;
   }
   

通过这次作业，学到了如何使用 SYCL 进行基本的并行计算，包括数据传递、内核编写、队列提交等方面的内容。 SYCL 提供了一种相对高层次的抽象，使得编写并行代码更加方便，并且代码具有可移植性，可以在不同的支持 SYCL 的设备上运行。