使用oneAPI和C++/SYCL进行高效矩阵乘法

简介：

随着高性能计算（HPC）和机器学习应用的不断增长，利用GPU的强大计算能力变得越来越重要。oneAPI提供了一种使用Data Parallel C++（DPC++）编写跨平台代码的方法，其核心是SYCL标准，旨在简化并行编程。本博客将介绍如何使用oneAPI和SYCL程序编写一个基本的矩阵乘法程序，展示如何利用GPU进行高效计算。

OneAPI简介：

oneAPI是一个跨平台编程模型，旨在编写高性能的多平台并行计算。它由英特尔公司开发，在并行计算领域应用广泛。oneAPI提供了一个统一的编程模型，兼容多种编程语言，提供了各种方便的工具和库。通过使用oneAPI，开发者能够在不同的硬件架构上编写高性能的应用程序。

SYCL简介：

SYCL是一种基于C++的高级编程模型，用于编写异构计算的代码。它是由Khronos Group开发的开放标准，旨在提供一种统一的方式来编写代码，该代码可以在各种硬件上执行，包括CPU、GPU、DSP和FPGA。SYCL是一种单一源代码模型允许在同一源文件中混合传统的C++代码和专门的并行计算代码，使得开发和调试过程更加简单，与现有的C++库和框架兼容性很好，为需要跨多种硬件平台进行开发的项目提供了一种统一和高效的解决方案。

代码实现：

我们的目标是编写一个程序，它可以高效地计算两个1024x1024的浮点矩阵相乘。这个程序分为几个主要部分：

环境设置和初始化： 

引入必要的头文件，如 <sycl/sycl.hpp>，并定义了一个基于SYCL的自定义设备选择器。

代码：

#include <sycl/sycl.hpp>
#include <iostream>
#include <vector>
#include <random>
#include <chrono>

class CustomDeviceSelector : public sycl::device_selector {
 public:
  CustomDeviceSelector(const std::string& vendorName) : vendorName_(vendorName) {}
  int operator()(const sycl::device &dev) const override {
    if (dev.is_gpu() && dev.get_info<sycl::info::device::vendor>().find(vendorName_) != std::string::npos)
      return 3; // 优先选择指定供应商的GPU
    if (dev.is_gpu())
      return 2; // 其次选择其他GPU
    return dev.is_cpu() ? 1 : -1; // 最后选择CPU，不选择其他设备
  }
 private:
  std::string vendorName_;
};

数据准备和设备选择：

在main函数中，我们先创建三个大小为1024x1024的浮点向量（A, B, C），分别用于存储两个输入矩阵和一个输出矩阵，并对两个输入矩阵进行初始化。使用自定义的设备选择器，优先选择指定供应商（例如“Intel”）的GPU设备，配置队列。

代码：

auto start = std::chrono::high_resolution_clock::now();

const int N = 1024;
std::vector<float> A(N * N), B(N * N), C(N * N);

// 初始化矩阵A和B
std::mt19937 rng(42);
std::uniform_real_distribution<float> dist(0.0, 1.0);
for (int i = 0; i < N * N; ++i) {
    A[i] = dist(rng);
    B[i] = dist(rng);
}

// 选择设备并执行矩阵乘法
CustomDeviceSelector selector("Intel"); // 更改为所需的制造商
sycl::queue q(selector);

核心代码：

以下是实现矩阵乘法的核心部分，利用SYCL缓冲区（sycl::buffer）自动管理内存，包括在GPU端的内存分配和数据传输，定义了一个核函数，使用二维网格和块来并行计算矩阵乘法。

{
    sycl::buffer<float> bufA(A.data(), sycl::range<1>(N * N));
    sycl::buffer<float> bufB(B.data(), sycl::range<1>(N * N));
    sycl::buffer<float> bufC(C.data(), sycl::range<1>(N * N));

    q.submit([&](sycl::handler &h) {
        auto accA = bufA.get_access<sycl::access::mode::read>(h);
        auto accB = bufB.get_access<sycl::access::mode::read>(h);
        auto accC = bufC.get_access<sycl::access::mode::write>(h);

        h.parallel_for(sycl::nd_range<2>(sycl::range<2>(N, N), sycl::range<2>(16, 16)), [=](sycl::nd_item<2> item) {
            int row = item.get_global_id(0);
            int col = item.get_global_id(1);
            float sum = 0.0f;
            for (int k = 0; k < N; ++k)
                sum += accA[row * N + k] * accB[k * N + col];
            accC[row * N + col] = sum;
        });
    });
} 

结果获取和打印：

打印程序的执行时间和输出矩阵，测量执行时间，通常用于性能分析。

auto end = std::chrono::high_resolution_clock::now();
auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
std::cout << "运行时间: " << duration.count() << " ms" << std::endl;
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";    
}

结论：

在本例中，我们使用了16x16的二维线程块来进行矩阵乘法，这个简单的例子展示了如何使用oneAPI和SYCL进行GPU加速的矩阵乘法。通过将复杂的内存管理和设备选择抽象化，SYCL使得跨平台的高性能并行编程变得更加容易。