基于openAPI的C++/SYCL实现矩阵乘法并行加速

snow_and_sword

已于 2023-12-01 16:53:06 修改

阅读量1.3k

点赞数 35

分类专栏： oneAPI 文章标签： c++ 矩阵人工智能

于 2023-12-01 16:31:29 首次发布

本文链接：https://blog.csdn.net/snow_and_sword/article/details/134737583

版权

oneAPI 专栏收录该内容

3 篇文章 1 订阅

订阅专栏

基于openAPI的C++/SYCL实现矩阵乘法并行加速

一、任务需求
二、问题分析
三、oneAPI 介绍
四、oneAPI并行编程介绍
- 4.1 编程框架
- 4.2 编程流程
五、矩阵乘法代码实现
- 5.1 使用缓存和存储
- 5.2 使用统一共享内存（USM）
六、结论
七、收获与感悟
致谢

一、任务需求

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

项目简介: 该项目旨在使用 SYCL（基于 C++ 的高性能并行编程模型）来加速j矩阵乘法过程。通过利用 GPU 和其他硬件加速器的并行计算能力，大幅提高矩阵乘法的效率。
- 具体而言，在并行计算矩阵乘法时，可以利用线程块和线程的层次结构来优化计算。通过合理划分矩阵数据并利用共享内存来减少全局内存访问的次数，可以⼤幅提高计算效率。
- 还可以利用GPU上的多个计算单元并执行行矩阵乘法，进⼀步提高计算速度
技术栈及实现方案:
- SYCL：用于并行编程和加速计算
- Intel® oneAPI：提供工具和库来优化性能和简化开发过程
- Intel® AI Analytics Toolkit：包含加速AI和机器学习工作流的工具
- C++：项目的主要编程语言

二、问题分析

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

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

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

三、oneAPI 介绍

随着科学技术的迅猛发展，高性能计算在人工智能、药物研制、智慧医疗、计算化学等领域发挥着日益重要的作用。然而，随着后摩尔时代的来临，计算机系统结构进入了百花齐放、百家争鸣的繁荣时期，CPU、GPU、FPGA和AI芯片等相互补充。硬件的多样性带来了软件设计与开发的复杂性，高性能计算并行程序的计算效率以及在不同计算平台之间的可移植性日益变得重要。为了解决这一问题，Intel推出了oneAPI。

3.1 oneAPI诞生背景

随着科技的迅猛发展，高性能计算在人工智能、药物研制、智慧医疗、计算化学等领域的作用日益显著。然而，随着后摩尔时代的到来，计算机系统结构迎来了繁荣时期，CPU、GPU、FPGA和AI芯片等各异的硬件相互补充。这种硬件多样性导致了软件设计与开发的复杂性的增加，高性能计算并行程序的计算效率和在不同计算平台之间的可移植性变得日益关键。为了解决这一问题，Intel推出了oneAPI。

3.2 oneAPI 简介

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

使用Intel oneAPI编程具有以下好处：
1. 跨硬件平台：允许开发者针对不同类型的硬件（包括非Intel硬件）编写统一的代码，增加代码的可移植性和灵活性。
2. 性能优化：提供针对Intel硬件优化的库和工具，可以帮助开发者提高应用程序的性能。
3. 简化开发流程：通过统一的API和工具链，简化了对多种硬件的编程和优化过程，降低了学习曲线。
4. 支持现代编程语言：支持包括C++、Python在内的多种现代编程语言，适用于广泛的开发需求。
5. 促进创新：通过提供高级的编程模型和工具，使开发者能够更容易地探索和实现新的算法和应用。

简而言之，oneAPI 是一种通用的可以适配Intel的多种异构硬件的并行编程库，个人理解可以用于HPC场景以及深度学习的卷积优化、向量传播等场景，并降低适配不同硬件的编译适配成本。

OneAPI提供了云平台DevCloud可以免费注册后使用，即使笔记本不是Intel的也可以使用Intel OneAPI配套硬件的体验，还无需自行安装toolkit。

jupyter.oneapi.devcloud.intel.com/ 这是我在本次实验中使用的云平台链接，本次实验主要基于Jupyter Lab进行，使用简单，用Jupyter Notebook构建了可视化界面，即使不会使用命令行进行编译也可以直接复制官方教程中的脚本进行编译，免去了CPP的编译痛苦。

3.3 oneAPI整体架构

3.3.1 oneAPI 开放式规范

oneAPI 这一开放式规范包括一种跨架构的编程语言 Data Parallel C++（DPC++）、一套用于API编程的函数库以及底层硬件接口（oneAPI Level Zero），如下图1所示。有了这些组件，英特尔和其它企业就能创建他们自己的 oneAPI 实现来支持他们自己的产品，或基于 oneAPI 进行新产品开发。
Alt

3.3.2 SYCL规范

SYCL第一次是在2014年引入，它是一种基于C++异构平行编程框架，用来加速高性能计算，机器学习，内嵌计算，以及在相当宽泛的处理器构架之上的计算量超大的桌面应用。这些处理器包括了CPU, GPU, FPGA, 和张量加速器。

2021年2月9号 , 科纳斯组织（Khronos® Group），作为一个由工业界主流公司组成的创建先进的互联标准的开放协会，宣布了SYCL 2020最终版规范的批准和发布。这个规范是单源C++并行编程的开放标准。作为多年来规范开发的一个主要的里程碑，SYCL 2020是在SYCL 1.2.1的功能的基础之上建立的，用以进一步改善可编程性，更小的代码尺寸，和高效的性能。基于C++17之上的SYCL 2020, 使得标准C++应用的加速更为容易，而且推动使之与ISO C++的路线图变得更为一致。SYCL 2020 将会进一步加速在多平台上的采用和部署，包括使用除了OpenCLTM之外的多样的加速API 后端。
SYCL 2020集成了超过40项新的特征，包括了为简化代码所做的更新，和更小的代码尺寸。一些主要增加的内容包括：

统一的共享存储（USM）使得带有指向器的代码，可以在不需要缓冲与存取器的情况下自然地工作
工作组和子工作组算法，在工作项目中增加了有效率的并行操作
类模板参数推导(CTAD)与模版减量指南，简化了类模板实例化
扩展的互操作性通过各种后端加速API实现高效加速
统一的共享存储（USM）使得带有指向器的代码，可以在不需要缓冲与存取器的情况下自然地工作;
并行的减量增加了一种内置的减量操作，来减少样板代码以达到具有内置的减量操作加速硬件上的最大性能。

3.3.3 DPC++介绍

oneAPI包含一种全新的跨架构编程语言 DPC++，DPC++基于 C++编写，由一组C++类、模板与库组成，同时兼容 Kronos 的 SYCL 规范，图2给出了DPC++与SYCL、C++关系。同时，intel DPC++兼容性工具可以实现将CUDA代码迁移到DPC++上，其中大约会有80%-90%的代码实现了自动迁移并提供内联注释，很大程度上帮助开发人员减轻代码移植的负担。
Alt DPC++是一种单一源代码语言，其中主机代码和异构加速器内核可以混合在同一源文件中。在主机上调用 DPC++程序，并将计算加载到加速器。程序员使用熟悉的C++和库结构，并添加诸如工作目标队列、数据管理缓冲区和并行性并行的函数，以指导计算和数据的哪些部分应该被加载。

3.3.4 oneAPI工具包

oneAPI 编程模式兼容性堪称达到了历史最强。目前在各个领域应用比较广泛的高性能计算开发工具如 Fortran，在 AI 领域的 Python，以及像 OpenMP 这样不同领域使用的语言都可以做到无缝对接，同时，oneAPI 也支持一些主流的 AI 工具包，包括 Hadoop、Spark、TensorFlow、PyTorch、PaddlePaddle、OpenVINO 等等，形成更适合人工智能时代的软件栈。oneAPI有六个工具包，几乎涵盖了高性能计算、物联网、渲染、人工智能、大数据分析这些领域。

Intel® oneAPI Base Toolkit：这个工具包是 oneAPI 其他产品的基础，包含了几个在 Parallel Studio中常用的软件以及 icc 编译器、MPI、DPCPP 等。这个工具包使开发人员都可以跨CPU、GPU和FPGA构建、测试和部署以性能为中心、以数据为中心的应用程序。
Intel® oneAPI HPC Toolkit ：这个工具包提供可扩展的快速C ++、Fortran、OpenMP和MPI应用程序。从某种程度上来说 Intel® oneAPI Base Toolkit 加 Intel® oneAPI HPC Toolkit 基本就包含Intel Parallel Studio XE的功能了。
Intel® oneAPI IoT Toolkit ：这个工具包主要用于建立可在网络边缘运行的高性能、高效、可靠的解决方案，属于物联网领域。
Intel® AI Analytics Toolkit ：这个工具包提供优化的深度学习框架和高性能Python库，加速端到端机器学习和数据科学库。这些组件是使用 oneAPI 库构建的，用于低级计算优化。这可以最大化从预处理到机器学习的性能。
Intel® oneAPI Rendering Toolkit：它主要用于创建高性能、高保真的可视化应用程序，适用于各种渲染领域。
Intel® Distribution of OpenVINO™ Toolkit：这个工具包用于从设备到云部署高性能推理应用程序。该工具包基于卷积神经网络（CNN），可将工作负载扩展到整个英特尔®硬件（包括加速器），并最大限度地提高性能。该工具包可以使深度学习推理从边缘到云，加速人工智能工作负载,包括计算机视觉、音频、演讲,语言,和推荐系统。支持异构执行在英特尔架构和AI加速器CPU、iGPU,英特尔Movidius视觉处理单元(VPU)、FPGA,和英特尔高斯 & 神经加速器(Intel® GNA)。

四、oneAPI并行编程介绍

4.1 编程框架

**平台模型：**oneAPI的平台模型基于SYCL*平台模型。它指定控制一个或多个设备的主机。主机是计算机，通常是基于CPU的系统，执行程序的主要部分，特别是应用范围和命令组范围。主机协调并控制在设备上执行的计算工作。设备是加速器，是包含计算资源的专门组件，可以快速执行操作的子集，通常比系统中的CPU效率更高。每个设备包含一个或多个计算单元，可以并行执行多个操作。每个计算单元包含一个或多个处理元素，充当单独的计算引擎。图3所示为平台模型的可视化描述。一个主机与一个或多个设备通信。每个设备可以包含一个或多个计算单元。每个计算单元可以包含一个或多个处理元素。
在这里插入图片描述

**执行模型：**执行模型基于SYCL*执行模型。它定义并指定代码（称为内核kernel）如何在设备上执行并与控制主机交互。主机执行模型通过命令组协调主机和设备之间的执行和数据管理。命令组（由内核调用、访问器accessor等命令组成）被提交到执行队列。访问器(accessor)形式上是内存模型的一部分，它还传达执行的顺序要求。使用执行模型的程序声明并实例化队列。可以使用程序可控制的有序或无序策略执行队列。有序执行是一项英特尔扩展。设备执行模型指定如何在加速器上完成计算。从小型一维数据到大型多维数据集的计算通过ND-range、工作组、子组（英特尔扩展）和工作项的层次结构中进行分配，这些都在工作提交到命令队列时指定。需注意的是，实际内核代码表示为一个工作项执行的工作。内核外的代码控制执行的并行度多大；工作的数量和分配由 ND-range和工作组的规格控制。下图4描述了ND-range、工作组、子组和工作项之间的关系。总工作量由ND-range的大小指定。工作的分组由工作组大小指定。本例显示了ND-range的大小X*Y*Z，工作组的大小 X’* Y’*Z’，以及子组的大小X’。因此，有X*Y*Z工作项。有(X*Y*Z)/ (X’*Y’*Z’)工作组和(X*Y*Z)/X’子组。
在这里插入图片描述

**内存模型：*oneAPI 的内存模型基于 SYCL 内存模型。它定义主机和设备如何与内存交互。它协调主机和设备之间内存的分配和管理。内存模型是一种抽象化，旨在泛化和适应不同主机和设备配置。在此模型中，内存驻留在主机或设备上，并由其所有，通过声明内存对象来指定如图5所示。内存对象有两种：缓冲器和图像。这些内存对象通过访问器在主机和设备之间进行交互，访问器传达期望的访问位置（如主机或设备）和特定的访问模式（如读或写)。
在这里插入图片描述

在oneAPI内存模型中的Buffer Model：缓冲区(Buffer)将数据封装在跨设备和主机的 SYCL 应用中。访问器(Accessor)是访问缓冲区数据的机制。设备(device)和主机(host)可以共享物理内存或具有不同的内存。当内存不同时，卸载计算需要在主机和设备之间复制数据。DPC++不需要您管理数据复制。通过创建缓冲区(buffer)和访问器(accessor)，DPC++能够确保数据可供主机和设备使用，而无需您介入。DPC++ 还允许您明确地显式控制数据移动，以实现最佳性能。需要注意的是，在这种内存模式下，若有多个内核程序使用相同的缓冲区，访问器需要根据依赖关系，以对内核执行进行排序以避免争用缓冲区而出现错误（通过主机访问器或缓冲区破坏实现。

**编程模型：*面向 oneAPI 的内核编程模式基于 SYCL 内核编程模型。它支持主机和设备之间的显式并行性。并行性是显式的，因为程序员决定在主机和设备上执行什么代码；它不是自动的。内核代码在加速器上执行。采用 oneAPI 编程模型的程序支持单源，这意味着主机代码和设备代码可以在同一个源文件中。但主机代码中所接受的源代码与设备代码在语言一致性和语言特性方面存在差异。SYCL 规范详细定义了主机代码和设备代码所需的语言特性。

4.2 编程流程

DPC++程序设计大致可分为以下5个步骤：
(1) 申请Host内存
(2) 创建SYCL缓冲区并为其定义访问缓冲区内存的方法。
设备(device)和主机(host)可以共享物理内存或具有不同的内存。当内存不同时，卸载计算需要在主机和设备之间复制数据。而通过创建缓冲区(buffer)和访问器(accessor)的方式，DPC++就不需要您管理数据复制，其能够确保数据可供主机和设备使用，而无需您介入。DPC++ 还允许您明确地显式控制数据移动，以实现最佳性能。

//创建vector1向量的SYCL缓冲区；
buffer vector1_buffer(vector1,R);
定义了访问缓冲区内存的accessor；
accessor v1_accessor (vector1_buffer,h,read_only);

（3）创建队列以向Device（包括Host）提交工作（包括选择设备和排队）

q.submit([&](handler& h) {
    //COMMAND GROUP CODE
});

可以通过选择器(selector)选择 CPU、GPU、FPGA和其他设备。使用默认的 q，这意味着 DPC++ 运行时会使用默认选择器（default selector）选择功能最强大的设备。
（4）调用oneAPI的核函数（Kernel）在Device上完成指定的运算。
该内核将应用于索引空间中的每个点，内核封装在C++ lambda函数中。DPC++中内核的形式如下：

h.parallel_for(range<1>(1024), [=](id<1> i){
    A[i] =  B[i] + C[i];
});

在该循环中，每个迭代都是完全独立的，并且不分顺序。使用 parallel_for 函数表示并行内核。
（5）将SYCL缓冲区的数据读到Host端。

五、矩阵乘法代码实现

本实验使用的头文件如下：

#include <sycl/sycl.hpp>
#include <iostream>
#include <vector>
#include <random> //用于随机初始化矩阵
#include <chrono> //用于记录运行时间

本实验同时使用了JupyterLab中给出的一种选择最佳设备的方法，其类定义如下：

class CustomDeviceSelector {
 public:
  CustomDeviceSelector(std::string vendorName) : vendorName_(vendorName){};
  int operator()(const device &dev) {
    int device_rating = 0;
    //We are querying for the custom device specific to a Vendor and if it is a GPU device we
    //are giving the highest rating as 3 . The second preference is given to any GPU device and the third preference is given to
    //CPU device. 
    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_;
};

5.1 使用缓存和存储

首先组建2个序列（向量）的浮点数，每个序列的大小是 $N * N$ （ $N = 1024$ ），构成矩阵输入值。用随机值初始化序列。使用缓存和存储来实现对设备（GPU）上矩阵内存的分配并运行。

    auto start = std::chrono::high_resolution_clock::now();
    
    const int N = 1024;
    const int Size = N * N;
    
    //# Create two vectors for the input matrices and one for the output
    std::vector<float> A(Size),B(Size),C(Size,0);
    
    //# 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++) {
        A[i] = dist(rng);
        B[i] = dist(rng);
    }
    
    //# Create buffers
    buffer Matrix1_buffer(A);
    buffer Matrix2_buffer(B);
    buffer Output_buffer(C);

    //# Choose the best device
    // Pass in the name of the vendor for which the device you want to query
    std::string vendor_name = "Intel";
    // std::string vendor_name = "AMD";
    // std::string vendor_name = "Nvidia";
    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);

随后，运行SYCL Kernel实现两个矩阵的并行运算。此处我们使用SYCL nd_range概念来定义Kernel的运行范围，并使用SYCL排队系统来运行设备上的Kernel。

      //继承上文的代码部分
      h.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];
          }
        });
    });

在Kernel运行结束后，使用存储将结果从设备（GPU）检索回主机，并输出运行结果和运行时间。

    //# Create a host accessor to copy data from device to host
    host_accessor h_a(Output_buffer,read_only);

    //# Print Input Values
    std::cout<<"\nInput Matrix1:\n";    
    for (size_t i = 0; i < N; i++) {
        for (size_t j = 0; j < N; j++) {
            std::cout << A[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 << B[i*N+j] << " ";
        }
        std::cout << "\n";
    }
    //# 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 << C[i*N+j] << " ";
        }
        std::cout << "\n";
    }
    
    //# Record Endtime
    auto end = std::chrono::high_resolution_clock::now();
    
    auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
    std::cout << "\nThe run time is: " << duration.count() << " ms.\n";

5.2 使用统一共享内存（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);