使用因特尔OneAPI工具实现并行计算

因特尔OneApi介绍

对于大多数应用程序开发者来说，使用高级语言进行编程已经成为一种常见的实践。现有的高级语言编译器已经很好地将程序开发与底层计算机体系结构分离开来，使开发者能够专注于算法和应用程序的开发，而无需深入了解底层处理器指令的操作码。

随着应用程序复杂性的增加和对计算能力的要求，使用各种硬件加速器（如GPU、FPGA等）成为提高性能和满足实时需求的关键。然而，针对这些硬件加速器进行开发需要深入了解底层硬件体系结构以及特定的开发技巧和工具。这确实需要更专业的知识和经验，超出了普通软件工程师或算法工程师的范畴。

英特尔oneAPI是一个综合性的软件开发工具集，旨在支持跨多种硬件架构的高性能计算。它提供了一套统一的编程模型和工具，使开发人员能够轻松地利用不同类型的处理器和加速器来加速应用程序的执行。oneAPI的目标是实现代码的可移植性和可扩展性，使开发人员能够更高效地利用现代硬件。

OneAPI的目标是为开发者提供一致的编程体验，无论使用哪种硬件加速器，都能够以统一的方式进行开发。这使得开发者能够更高效地利用异构计算平台的优势，提高应用程序的性能和效率。英特尔oneAPI Developer Cloud 服务，可免安装额外环境，直接利用Developer Cloud平台中的CPU与GPU硬件完成相应的作业。

1.并行矩阵乘法

1.1 暴力计算

我们首先使用暴力计算直接进行矩阵相乘，代码如下：

 %%writefile src/matrix_btruteforce.cpp
 #include <CL/sycl.hpp>
 #include <iostream>
 #include <vector>
 #include <chrono>
 ​
 constexpr size_t N = 1024;
 ​
 // 矩阵乘法函数，计算矩阵A和矩阵B的乘积结果存储在矩阵C中
 void matrixMultiplication(const std::vector<float>& matrixA, const std::vector<float>& matrixB, std::vector<float>& matrixC) {
     for (int i = 0; i < N; ++i) {
         for (int j = 0; j < N; ++j) {
             float sum = 0.0f;
             for (int k = 0; k < N; ++k) {
                 sum += matrixA[i * N + k] * matrixB[k * N + j];
             }
             matrixC[i * N + j] = sum;
         }
     }
 }
 ​
 int main() {
     std::vector<float> matrixA(N * N, 2.0f);
     std::vector<float> matrixB(N * N, 3.0f);
     std::vector<float> matrixC(N * N, 0.0f);
 ​
     auto start = std::chrono::high_resolution_clock::now();
 ​
     matrixMultiplication(matrixA, matrixB, matrixC);
 ​
     auto stop = std::chrono::high_resolution_clock::now();
     auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(stop - start);
 ​
     std::cout << "暴力求解时间: " << duration.count() << " milliseconds" << std::endl;
 ​
     // 打印结果矩阵C的前2x2部分，也可以打印更多信息
     for (size_t i = 0; i < 2; ++i) {
         for (size_t j = 0; j < 2; ++j) {
             std::cout << matrixC[i * N + j] << " "; // 打印矩阵C中对应位置的元素
         }
         std::cout << std::endl;
     }
 ​
     return 0;
 }

• 执行时间：1781ms

1.2 并行矩阵计算

1.2.1 代码思路：

首先，代码中定义了一些常量：N表示矩阵的维度（这里是1024x1024），BLOCK_SIZE表示在局部内存中处理块的大小（这里是16x16）。

代码创建了三个std::vector<float>来代表矩阵A、B和C。矩阵A被初始化所有元素为2.0f，矩阵B被初始化所有元素为3.0f，矩阵C初始化为0。

在try块中，代码首先创建了一个SYCL队列，选择了一个GPU设备。

定义了全局和局部的工作项尺寸。全局尺寸（global_size）定义了整个工作负载的尺寸，而局部尺寸（local_size）定义了工作组内的工作项数量。

使用矩阵A、B和C的数据创建了三个SYCL缓冲区（bufferA，bufferB和bufferC），以便在设备上进行计算。

通过myQueue.submit提交了一个lambda表达式，它定义了如何执行矩阵乘法的计算。

在lambda中，首先获取了缓冲区的访问器，这些访问器允许在设备上读写缓冲区。

接着，定义了两个局部访问器：localA和localB。这些被用于在局部内存中存储矩阵A和B的块以加速计算。

parallel_for构造定义了一个命名为matrix_mul的类，用来并行执行矩阵乘法的内核。

在内核中，每个工作项计算自己负责的全局行和列的索引。然后，通过迭代k的方式，在块内进行局部内存的加载、计算和存储。

内核结束时，计算得到的sum被写入结果矩阵C的对应位置。

等待队列完成所有计算后，测量并输出执行时间。

代码最后打印出矩阵C的前两行和两列的元素，以验证计算结果。

1.2.2 并行矩阵计算代码

 %%writefile src/matrix.cpp
 #include <CL/sycl.hpp>
 #include <iostream>
 #include <vector>
 #include <chrono>
 ​
 constexpr size_t N = 1024;
 constexpr size_t BLOCK_SIZE = 16;
 ​
 int main() {
     std::vector<float> matrixA(N * N, 2.0f); //初始话矩阵 2.0
     std::vector<float> matrixB(N * N, 3.0f); // 初始化为 3.0
     std::vector<float> matrixC(N * N, 0.0f);
 ​
     try {
         sycl::queue myQueue(sycl::gpu_selector{});
         sycl::range<2> global_size(N, N);
         sycl::range<2> local_size(BLOCK_SIZE, BLOCK_SIZE);
         sycl::buffer<float, 2> bufferA(matrixA.data(), sycl::range<2>(N, N));
         sycl::buffer<float, 2> bufferB(matrixB.data(), sycl::range<2>(N, N));
         sycl::buffer<float, 2> bufferC(matrixC.data(), sycl::range<2>(N, N));
 ​
         auto start = std::chrono::high_resolution_clock::now();
         myQueue.submit([&](sycl::handler& cgh) {
             auto accessorA = bufferA.get_access<sycl::access::mode::read>(cgh);
             auto accessorB = bufferB.get_access<sycl::access::mode::read>(cgh);
             auto accessorC = bufferC.get_access<sycl::access::mode::write>(cgh);
             sycl::accessor<float, 2, sycl::access::mode::read_write, sycl::access::target::local> localA(sycl::range<2>(BLOCK_SIZE, BLOCK_SIZE), cgh);
             sycl::accessor<float, 2, sycl::access::mode::read_write, sycl::access::target::local> localB(sycl::range<2>(BLOCK_SIZE, BLOCK_SIZE), cgh);
 ​
             cgh.parallel_for<class matrix_mul>(sycl::nd_range<2>(global_size, local_size), [=](sycl::nd_item<2> item) {
                 const int globalRow = item.get_global_id(0);
                 const int globalCol = item.get_global_id(1);
                 const int localRow = item.get_local_id(0);
                 const int localCol = item.get_local_id(1);
                 float sum = 0.0f;
 ​
                 for (int k = 0; k < N; k += BLOCK_SIZE) {
                     localA[localRow][localCol] = accessorA[globalRow][k + localCol];
                     localB[localRow][localCol] = accessorB[k + localCol][globalCol]; // 使用 k + localCol 访问
                     item.barrier(sycl::access::fence_space::local_space);
 ​
                     for (int n = 0; n < BLOCK_SIZE; ++n) {
                         sum += localA[localRow][n] * localB[n][localCol];
                     }
                     item.barrier(sycl::access::fence_space::local_space);
                 }
 ​
                 accessorC[globalRow][globalCol] = sum;
             });
         });
 ​
         myQueue.wait_and_throw();
         auto stop = std::chrono::high_resolution_clock::now();
         auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(stop - start);
         std::cout << "运行时间: " << duration.count() << " milliseconds" << std::endl;
     } catch (sycl::exception const& e) {
         std::cerr << "出现异常 " << e.what() << std::endl;
         return 1;
     }
 ​
     // 打印前两行和前两列的结果
     std::cout << "输出结果:" << std::endl;
     for (size_t i = 0; i < 2; ++i) { // 打印矩阵的前两行
         for (size_t j = 0; j < 2; ++j) { // 打印矩阵的前两列
             std::cout << matrixC[i * N + j] << " ";
         }
         std::cout << std::endl;
     }
 ​
     return 0;
 }

1.2.3 运行结果

可以看到使用并行计算，所花费的时间是158ms，这里打印了得到的矩阵的前两行和前两列，结果和暴力计算是一样的。

1.3 结果对比

可以使用并行矩阵计算比直接使用暴力计算快了11倍，同时两者计算得到的结果都是一样的。

2.并行排序

2.1 编程思路

以下代码是一个并行归并排序算法的实现，使用了SYCL库，这是一个能够实现异构计算的C++标准库。代码包括了串行和并行排序的过程，并能在主机和GPU上运行。

mergeSequences函数：该函数用于合并两个已排序的子序列。它首先创建一个临时的SYCL缓冲区（tempBuf），用于存储合并的结果。然后，它将序列合并的工作提交给SYCL队列（deviceQueue），该队列负责在设备上异步执行合并操作。最后，它将合并后的数据复制回原始数据缓冲区。

parallelMergeSort函数：这是主要的排序函数，它使用分而治之的方法将数组分成更小的部分进行排序。首先检查是否到达递归的基础条件（start < end）。当递归深度maxDepth降到0以下时，使用串行排序；否则，继续递归划分数组，并在每个部分上调用自己，最后合并这些已排序的部分。

main函数：这是程序的入口点。它首先创建了一个大小为2的20次幂的数组，并用随机数填充。然后，它创建了一个SYCL队列和一个数据缓冲区，用于在设备上运行排序算法。调用parallelMergeSort函数进行排序，并验证排序是否正确。

2.2 并行归并排序代码

 %%writefile src/parallelsort.cpp
 #include <CL/sycl.hpp>
 #include <vector>
 #include <iostream>
 #include <algorithm>
 // 不再使用整个命名空间
 using namespace sycl;
 ​
 // 串行归并排序，使用标准库中的排序算法
 void serialMergeSort(std::vector<int>& data) {
     std::sort(data.begin(), data.end());
 }
 ​
 // 合并两个有序子序列
 void mergeSequences(queue &deviceQueue, buffer<int, 1> &bufferData, int start, int middle, int end) {
     std::vector<int> tempBuffer(end - start + 1);
     buffer<int, 1> tempBuf(tempBuffer.data(), tempBuffer.size());
     deviceQueue.submit([&](handler &cgh) {
         auto dataAcc = bufferData.get_access<access::mode::read_write>(cgh);
         auto tempAcc = tempBuf.get_access<access::mode::write>(cgh);
         cgh.parallel_for(range<1>(end - start + 1), [=](id<1> idx) {
             int leftIndex = start;
             int rightIndex = middle + 1;
             int tempIndex = idx[0];
             if (leftIndex <= middle && (rightIndex > end || dataAcc[leftIndex] < dataAcc[rightIndex])) {
                 tempAcc[tempIndex] = dataAcc[leftIndex++];
             } else {
                 tempAcc[tempIndex] = dataAcc[rightIndex++];
             }
         });
     });
     deviceQueue.wait();
     deviceQueue.submit([&](handler &cgh) {
         auto dataAcc = bufferData.get_access<access::mode::write>(cgh);
         auto tempAcc = tempBuf.get_access<access::mode::read>(cgh);
         cgh.parallel_for(range<1>(end - start + 1), [=](id<1> idx) {
             dataAcc[start + idx[0]] = tempAcc[idx[0]];
         });
     });
     deviceQueue.wait();
 }
 ​
 // 并行归并排序
 void parallelMergeSort(queue &deviceQueue, buffer<int, 1> &bufferData, int start, int end, int maxDepth) {
     if (start < end) {
         if (maxDepth <= 0) {
             auto hostData = bufferData.get_access<access::mode::read_write>();
             std::vector<int> localData(hostData.get_pointer() + start, hostData.get_pointer() + end + 1);
             serialMergeSort(localData);
             for (int i = start; i <= end; ++i) {
                 hostData[i] = localData[i - start];
             }
         } else {
             int middle = (start + end) / 2;
             parallelMergeSort(deviceQueue, bufferData, start, middle, maxDepth - 1);
             parallelMergeSort(deviceQueue, bufferData, middle + 1, end, maxDepth - 1);
             mergeSequences(deviceQueue, bufferData, start, middle, end);
         }
     }
 }
 ​
 int main() {
     // 设置数组的大小
     const int numElements = 1 << 20;
     std::vector<int> inputData(numElements);
     // 生成随机数据
     std::generate(inputData.begin(), inputData.end(), [] { return std::rand() % 100; });
     queue deviceQueue;
     buffer<int, 1> dataBuf(inputData.data(), range<1>(numElements));
     parallelMergeSort(deviceQueue, dataBuf, 0, numElements - 1, 4);
     // 获取并输出排序后的数据
     auto sortedBuffer = dataBuf.get_access<access::mode::read>();
     for (int i = 0; i < numElements; ++i) {
         if(i > 5000 && i < 5100) {
             std::cout << sortedBuffer[i] << " ";
         }
     }
     std::cout << "\n";
     // 验证数组是否已正确排序
     bool isSorted = std::is_sorted(sortedBuffer.get_pointer(), sortedBuffer.get_pointer() + numElements);
     std::cout << "Array is sorted: " << (isSorted ? "Yes" : "No") << std::endl;
     
     return 0;
 }

2.3 实验结果：

可以看到程序对给出的样例数组已经排好序

3. 图像卷积并行加速

3.1 编程思路

为了简洁，假设图像和卷积核是正方形的，并且卷积核的尺寸较小。 定义一个convolve2D函数，该函数接受输入图像、卷积核、输出图像的向量，以及图像和卷积核的尺寸。

通过定义缓冲区将输入图像、卷积核和输出图像传递给设备。

使用q.submit将任务提交到SYCL队列。定义一个parallel_for内核，其范围为输出图像的尺寸，减去卷积核的大小加

3.2 实验代码

 %%writefile src/conv.cpp
 #include <CL/sycl.hpp>
 #include <vector>
 #include <iostream>
 ​
 using namespace sycl;
 ​
 // 函数：执行卷积操作
 void convolve2D(const std::vector<float>& inputImage, const std::vector<float>& kernel,
                 std::vector<float>& outputImage, int imageSize, int kernelSize, queue& q) {
     // 创建SYCL缓冲区
     buffer<float, 2> bufferInput(inputImage.data(), range<2>(imageSize, imageSize));
     buffer<float, 2> bufferKernel(kernel.data(), range<2>(kernelSize, kernelSize));
     buffer<float, 2> bufferOutput(outputImage.data(), range<2>(imageSize, imageSize));
 ​
     // 提交任务到队列
     q.submit([&](handler& cgh) {
         // 获取访问器
         auto accInput = bufferInput.get_access<access::mode::read>(cgh);
         auto accKernel = bufferKernel.get_access<access::mode::read>(cgh);
         auto accOutput = bufferOutput.get_access<access::mode::write>(cgh);
 ​
         // 定义卷积操作
         cgh.parallel_for(range<2>(imageSize - kernelSize + 1, imageSize - kernelSize + 1), [=](id<2> idx) {
             int x = idx[0];
             int y = idx[1];
             float sum = 0.0f;
 ​
             // 遍历卷积核
             for (int i = 0; i < kernelSize; ++i) {
                 for (int j = 0; j < kernelSize; ++j) {
                     sum += accInput[x + i][y + j] * accKernel[i][j]; // 应用卷积核
                 }
             }
 ​
             accOutput[x + kernelSize/2][y + kernelSize/2] = sum; // 存储卷积结果
         });
     }).wait(); // 等待队列中的任务完成
 }
 ​
 int main() {
     // 图像和卷积核的大小
     const int imageSize = 1024; // 假设为1024x1024
     const int kernelSize = 3;   // 假设为3x3
 ​
     // 初始化输入图像、卷积核和输出图像
     std::vector<float> inputImage(imageSize * imageSize, 1.0f); // 示例：所有值初始化为1
     std::vector<float> kernel = {
         0.0625f, 0.125f, 0.0625f,
         0.125f, 0.25f, 0.125f,
         0.0625f, 0.125f, 0.0625f
     }; // 示例：高斯模糊核
     std::vector<float> outputImage(imageSize * imageSize, 0.0f);
 ​
     // 创建SYCL队列
     queue q;
 ​
     // 执行卷积操作
     convolve2D(inputImage, kernel, outputImage, imageSize, kernelSize, q);
 ​
     // 输出结果（这里只输出前10个元素进行验证）
     for (int i = 0; i < 10; ++i) {
         std::cout << "Output[" << i << "] = " << outputImage[i] << std::endl;
     
     }
 ​
     return 0;
 }

3.3 实验结果

可以看到输出的验证结果是正确的、

4. 实验代码附件

实验代码附件见SYCL.ipynb，里面有三道实验的代码和运行结果