基于因特尔OneAPI实现矩阵并行乘法运算-Toy模板网

这篇具有很好参考价值的文章主要介绍了基于因特尔OneAPI实现矩阵并行乘法运算。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

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

OneAPI介绍

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

问题描述

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

实现思路

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

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

分配内存：在主机端分配内存空间用于存储输⼊矩阵和输出矩阵，同时在GPU端分配内存空间用于存储相应的输入和输出数据。
数据传输：将输入矩阵数据从主机端内存传输到GPU端内存中。
核函数调用：在SYCL中，矩阵乘法的计算通常会在GPU上使用核函数来实现并行计算。核函数会分配线程块和线程来处理不同的数据块。
并行计算：在核函数中，每个线程负责计算输出矩阵的⼀个单独的元素。为了最大限度地利用 GPU的并行计算能力，通常会使用⼆维线程块和线程网格的方式来处理矩阵的乘法计算。
数据传输：计算完成后，将输出矩阵数据从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);