基于因特尔OneAPI实现矩阵并行乘法运算

基于因特尔OneAPI实现矩阵并行乘法运算

OneAPI介绍

Intel oneAPI 是一个跨行业、开放、基于标准的统一的编程模型，旨在提供一个适用于各类计算架构的统一编程模型和应用程序接口。其核心思想是使开发者只需编写一次代码，便可在跨平台的异构系统上运行，支持的底层硬件架构包括 CPU、GPU、FPGA、神经网络处理器以及其他专为不同应用设计的硬件加速器等。这意味着，oneAPI不仅提高了开发效率，同时具备一定的性能可移植性。通过采用这一编程模型，开发者能够更灵活地利用不同类型的硬件，充分发挥各种计算资源的潜力，从而更好地适应不同应用场景的需求。

问题描述

编写⼀个基于oneAPI的C++/SYCL程序来执行矩阵乘法操作，需要考虑大尺寸矩阵的乘法操作以及不同线程之间的数据依赖关系。

实现思路

运行代码使用因特尔的云平台Get Started | Intel® DevCloud，进入oneapi/essential，然后进入02-SYCL-Program-Structure

编码方面利用基于SYCL的编程模型在GPU上实现矩阵乘法的计算，步骤如下：

1. 分配内存：在主机端分配内存空间用于存储输⼊矩阵和输出矩阵，同时在GPU端分配内存空间用于存储相应 的输入和输出数据。

2. 数据传输：将输入矩阵数据从主机端内存传输到GPU端内存中。

3. 核函数调用：在SYCL中，矩阵乘法的计算通常会在GPU上使用核函数来实现并行计算。核函数 会分配线程块和线程来处理不同的数据块。

4. 并行计算：在核函数中，每个线程负责计算输出矩阵的⼀个单独的元素。为了最大限度地利用 GPU的并行计算能力，通常会使用⼆维线程块和线程网格的方式来处理矩阵的乘法计算。

5. 数据传输：计算完成后，将输出矩阵数据从GPU端内存传输回主机端内存中，以便进⼀步处理或 分析。

在并行计算矩阵乘法时，可以利用线程块和线程的层次结构来优化计算。通过合理划分矩阵数据并利用共享内 存来减少全局内存访问的次数，可以⼤幅提高计算效率。此外，还可以利用GPU上的多个计算单元并执行行矩 阵乘法，进⼀步提高计算速度

代码详解

 %%writefile lab/vector_add.cpp
 //==============================================================
 // Copyright © Intel Corporation
 //
 // SPDX-License-Identifier: MIT
 // =============================================================
 #include <sycl/sycl.hpp>
 #include <iostream>
 #include <vector>
 #include <random>
  
 using namespace sycl;
 using namespace std;
  
 class CustomDeviceSelector {
  public:
   CustomDeviceSelector(std::string vendorName) : vendorName_(vendorName){};
   int operator()(const device &dev) {
     int device_rating = 0;
     if (dev.is_gpu() & (dev.get_info<info::device::name>().find(vendorName_) !=
                         std::string::npos))
       device_rating = 3;
     else if (dev.is_gpu())
       device_rating = 2;
     else if (dev.is_cpu())
       device_rating = 1;
     return device_rating;
   };
  
  private:
   std::string vendorName_;
 };
  
 int main() {
     const int N = 1024;
     
     //二维数组
     vector<vector<float>> matrix1(N, vector<float>(N));
     vector<vector<float>> matrix2(N, vector<float>(N));
     vector<vector<float>> matrix3(N, vector<float>(N));
     //初始化
     random_device rd;
     mt19937 rng(rd());
     uniform_int_distribution<int> dist(0, 2); 
  
     for (size_t i = 0; i < N; i++) {
         for (size_t j = 0; j < N; j++) {
             matrix1[i][j] = dist(rng);
             matrix2[i][j] = dist(rng);
         }
     }
     
     cout<<"\nmatrix1:\n";
      for (size_t i=0;i<N;i++)
 {
     for (size_t j=0;j<N;j++){
         cout<<matrix1[i][j]<<" ";
     }
     cout<<"\n";
 }
     cout<<"\nmatrix2:\n";
      for (size_t i=0;i<N;i++)
 {
     for (size_t j=0;j<N;j++){
         cout<<matrix2[i][j]<<" ";
     }
     cout<<"\n";
 }
     //Create buffers
     buffer<float, 2> Matrix1_buffer(reinterpret_cast<float*>(matrix1.data()), range<2>(N, N));
     buffer<float, 2> Matrix2_buffer(reinterpret_cast<float*>(matrix2.data()), range<2>(N, N));
     buffer<float, 2> Output_buffer(reinterpret_cast<float*>(matrix3.data()), range<2>(N, N));
  
  
     //Choose device
     std::string vendor_name = "Intel";
     CustomDeviceSelector selector(vendor_name);
     
  
     //Submit task to multiply matrices
     queue q(selector);
     q.submit([&](handler &h) {
       //# Create accessors for buffers
       accessor M1 (Matrix1_buffer,h,read_only);
       accessor M2 (Matrix2_buffer,h,read_only);
       accessor M3 (Output_buffer,h,write_only);
  
       h.parallel_for(nd_range<2>({N, N}, {64, 64}), [=](nd_item<2> item) {
           //# Multiplication
           size_t row = item.get_global_id(0);
           size_t col = item.get_global_id(1);
           for (size_t k = 0; k < N; ++k) {
               M3[row][col] += M1[row][k] * M2[k][col];
           }
         });
     });
  
  
     //Create accessor
     host_accessor h_a(Output_buffer, read_write);
  
     
     cout << "\nmatrix3:\n";
     for (size_t i = 0; i < N; i++) {
         for (size_t j = 0; j < N; j++) {
             cout << matrix3[i][j] << " ";
         }
         cout << "\n";
     }
     return 0;
 }

使用统一共享内存（USM）

统一共享内存（USM）在主机和设备（GPU）之间提供了一个统一的存储模式。通过用以下代码修改程序以使用统一共享内存来实现内存分配和数据转换，从而替代缓存和存储。（主要内容不变，仅仅改变了存储方式）

     //# USM allocation using malloc_shared
     float *M1 = malloc_shared<float>(Size, q);
     float *M2 = malloc_shared<float>(Size, q);
     float *M3 = malloc_shared<float>(Size, q);
     
     //# Initialize matrices A and B with random values
     std::mt19937 rng(42);
     uniform_int_distribution<int> dist(0, 4); 
     //std::uniform_real_distribution<float> dist(0.0, 1.0);
     for (size_t i = 0; i < Size; i++) {
         M1[i] = dist(rng);
         M2[i] = dist(rng);
         M3[i]=0;
     }
     std::cout<<"\nInput Matrix1:\n";    
     for (size_t i = 0; i < N; i++) {
         for (size_t j = 0; j < N; j++) {
             std::cout << M1[i*N+j] << " ";
         }
         std::cout << "\n";
     }
     std::cout<<"\nInput Matrix2:\n";    
     for (size_t i = 0; i < N; i++) {
         for (size_t j = 0; j < N; j++) {
             std::cout << M2[i*N+j] << " ";
         }
         std::cout << "\n";
     }
 ​
     q.parallel_for(nd_range<2>({N, N}, {16, 16}), [=](nd_item<2> item) {
       //# Multiplication
       size_t row = item.get_global_id(0);
       size_t col = item.get_global_id(1);
       for (size_t k = 0; k < N; ++k) {
           M3[row*N+col] += M1[row*N+k] * M2[k*N+col];
       }
     });
 ​
     //# Print Output values 
     std::cout<<"\nOutput Values:\n";
     for (size_t i = 0; i < N; i++) {
         for (size_t j = 0; j < N; j++) {
             std::cout << M3[i*N+j] << " ";
         }
         std::cout << "\n";
     }
     free(M1, q);
     free(M2, q);
     free(M3, q);
 ​
 ​

总结

通过本次实验学会了如何运用SYCL来进行并行编程，学会了因特尔OneAPI的工具和库函数来进行矩阵乘法运算，了解了因特尔的平台的运用方法，很方便，直接在网站上就可以运行自己的并行运算代码，增深了对高性能计算和并行处理的理解。

感谢因特尔相关团队提供的技术支持和云平台的支持。