基于英特尔oneAPI的C++/SYCL编程

本文将介绍Intel oneAPI的功能与作用,以及我如何从零开始使用Intel DevCloud for oneAPI进行C++/SYCL编程,并对三道作业题进行了作答。

一、英特尔oneAPI介绍

Intel oneAPI是一套全面的开发工具集,旨在简化在多种架构(如CPU、GPU、FPGA等)上进行应用程序和解决方案的开发。这套工具集的核心理念是提供统一的编程模型,以支持多种硬件类型,从而使开发人员能够编写一次代码,然后在不同的设备和架构上运行,无需为每种设备单独编写和优化代码。

oneAPI包括几个关键组件:

1. DPC++(Data Parallel C++):这是一种基于C++的编程语言,扩展了C++以支持数据并行编程。它允许开发人员利用不同类型硬件的并行计算能力。

2. oneAPI库:这些高性能库覆盖了多个领域,如数学、机器学习、数据分析和图像处理等。它们经过优化,可在多种硬件上高效运行。

3. Intel Advisor、VTune Profiler和其他工具:这些工具帮助开发者优化性能,分析代码效率,并确保在不同平台上的可移植性。

4. 兼容性和可扩展性:oneAPI旨在与现有的软件和工具兼容,同时也支持未来的硬件创新。

通过使用oneAPI,开发人员可以针对Intel的各种硬件产品(如CPU、GPU和FPGA)开发应用程序,而无需深入了解每种硬件的内部机制。这不仅提高了开发效率,还有助于实现更广泛的硬件支持和更好的性能优化。

二、DevCloud使用

首先在浏览器中访问下列快捷网址 https://idzcn.com/devcloud.htm 或直接访问英特尔统一登录账户注册页面 https://www.intel.com/content/www/cn/zh/secure/forms/devcloud/enrollment.html 开始注册一个新的英特尔用户账号。

然后访问 Get Started | Intel® DevCloud

在页面底部点击 Lanch JupyterLab

进入该页面后即可开始进行编程

可以参考 oneAPI_Essentials/01_oneAPI_Intro/oneAPI_Intro.ipynb 中的教程进行,其中 Simple Exercise 部分指导了如何进行 Build and Run:

! chmod 755 q; chmod 755 run_simple.sh;if [ -x "$(command -v qsub)" ]; then ./q run_simple.sh; else ./run_simple.sh; fi

三、作业1:并行矩阵乘法

1. 问题描述

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

2. 思路分析

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

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

② 数据传输:将输入矩阵数据从主机端内存传输到GPU端内存中。

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

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

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

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

3. 代码实现

CPU计算函数如下,返回运行时间:

double cpu_cacu(float *cA, float *cB, float *cC, int M, int N, int K) {
    double duration = 0.0;
    std::chrono::high_resolution_clock::time_point s, e;
    s = std::chrono::high_resolution_clock::now();
    for(int i = 0; i < M; i++) {
        for(int j = 0; j < N; j++) {
            float sum = 0.0f;
            for(int k = 0; k < K; k++) {
                sum +=  cA[i * K + k] * cB[k * N  + j];
            }
            cC[i * N + j] = sum;
        }
    }
    e = std::chrono::high_resolution_clock::now();
    duration = std::chrono::duration<float, std::milli>(e - s).count();
    return(duration);
}

GPU计算函数如下,在GPU上使用核函数来实现并行计算,核函数会分配线程块和线程来处理不同的数据块,计算完成后将输出矩阵数据从GPU端内存传输回主机端内存中:

double gpu_cacu(float *A, float *B, float *C, int M, int N, int K, int BLOCK, sycl::queue &q)
{
  auto grid_rows = M / tileY;
  auto grid_cols = N / tileX;
  auto local_ndrange  = range<2>(BLOCK, BLOCK);
  auto global_ndrange = range<2>(grid_rows, grid_cols);
  double duration = 0.0f;
  auto e = q.submit([&](sycl::handler &h) {
      h.parallel_for<class k_name_t>(
          sycl::nd_range<2>(global_ndrange, local_ndrange), [=](sycl::nd_item<2> index) {
              int row = tileY * index.get_global_id(0);
              int col = tileX * index.get_global_id(1);
              float sum[tileY][tileX] = {0.0f};
              float subA[tileY] = {0.0f};
              float subB[tileX] = {0.0f};
              for (int k = 0; k < N; k++) {
                for(int m = 0; m < tileY; m++) {
                    subA[m] = A[(row + m) * N + k];
                }
                for(int p = 0; p < tileX; p++) {
                    subB[p] = B[k * N + p + col];
                }
                for (int m = 0; m < tileY; m++) {
                  for (int p = 0; p < tileX; p++) {
                    sum[m][p] += subA[m] * subB[p];
                  }
                }
              }
              for (int m = 0; m < tileY; m++) {
                for (int p = 0; p < tileX; p++) {
                  C[(row + m) * N + col + p] = sum[m][p];
                }
              }
          });
    });
    e.wait();
    duration += (e.get_profiling_info<info::event_profiling::command_end>() -
    e.get_profiling_info<info::event_profiling::command_start>()) /1000.0f/1000.0f;
    return(duration);
}

其他部分: 

int verify(float *cpu_res, float *gpu_res, int length){
    int err = 0;
    for(int i = 0; i < length; i++) {
       if( fabs(cpu_res[i] - gpu_res[i]) > 1e-3) {
          err++;
          printf("\n%lf, %lf", cpu_res[i], gpu_res[i]);
       }
    }
    return(err);

}

int cacu(const int M, const int N, const int K, const int block_size, const int iterations, sycl::queue &q) {
  auto A = malloc_shared<float>(M * K, q);
  auto B = malloc_shared<float>(K * N, q);
  auto C = malloc_shared<float>(M * N, q);
  auto C_host = malloc_host<float>(M * N, q);
  for(int i=0; i < M * K; i++) {
      A[i] = random_float();
  }
  for(int i=0; i < K * N; i++) {
      B[i] = random_float();
  }
  for(int i=0; i < M * N; i++) {
      C[i] = 0.0f;
      C_host[i] = 0.0f;
  }
  double flopsPerMatrixMul = 2.0 * static_cast<double>(M) * static_cast<double>(N) * static_cast<double>(K);
  double duration_gpu = 0.0f;
  double duration_cpu = 0.0f;
  int warmup = 10;
  for (int run = 0; run < iterations + warmup; run++) {
    float duration = gpu_cacu(A, B, C, M, N, K, block_size, q);
    if(run >= warmup) duration_gpu += duration;
  }
  duration_gpu = duration_gpu / iterations;
  warmup = 2;
  for(int run = 0; run < iterations/2 + warmup; run++) {
      float duration = cpu_cacu(A, B, C_host, M, N, K);
      if(run >= warmup) duration_cpu += duration;
  }
  duration_cpu = duration_cpu / iterations/2;
  int errCode = 0;
  errCode = verify(C_host, C, M*N);
  if(errCode > 0)
      printf("\nThere are %d errors\n", errCode);
  cout << "Matrix size: c(" << M << "," <<  N << ") ="
       << " a(" << M << "," << K << ") *"
       << " b(" << K << "," << N << ")\n";
  printf("GPU Computation Time = %lf (ms); \n"
          "CPU Computaiton Time = %lf (ms); \n",
          duration_gpu, duration_cpu);
  free(A, q);
  free(B, q);
  free(C, q);
  free(C_host, q);
  return(errCode);
}

int main() {
  auto propList = cl::sycl::property_list {cl::sycl::property::queue::enable_profiling()};
  queue my_gpu_queue( cl::sycl::gpu_selector_v , propList);
  int errCode = cacu(512, 512, 512, 4, 10, my_gpu_queue);
  return(errCode);
}

运行结果如下:

可以看出GPU的速度是CPU的约4倍,体现了并行计算的速度优势。

四、作业2:并行排序算法

1. 问题描述

使用基于oneAPI的C++/SYCL实现⼀个高效的并行归并排序。需要考虑数据的分割和合并以及线程之间的协作。

2. 思路分析

归并排序是⼀种分治算法,其基本原理是将待排序的数组分成两部分,分别对这两部分进行排序,然后将已排序的子数组合并为⼀个有序数组。可考虑利用了异构并行计算的特点,将排序和合并操作分配给多个线程同时执行,以提高排序效率。具体实现过程如下:

① 将待排序的数组分割成多个较小的子数组,并将这些⼦数组分配给不同的线程块进行处理。

② 每个线程块内部的线程协作完成子数组的局部排序。

③ 通过多次迭代,不断合并相邻的有序⼦数组,直到整个数组有序。

在实际实现中,归并排序可使用共享内存来加速排序过程。具体来说,可以利用共享内存来存储临时数据,减少对全局内存的访问次数,从而提高排序的效率。另外,在合并操作中,需要考虑同步机制来保证多个线程之间的数据⼀致性。 需要注意的是,在实际应用中,要考虑到数组大小、线程块大小、数据访问模式等因素,来设计合适的算法和参数设置,以充分利用目标计算硬件GPU的并行计算能力,提高排序的效率和性能。

3. 代码实现

并行归并排序函数如下,将待排序的数组分割成多个较小的子数组,并将这些⼦数组分配给不同的线程块进行处理:

void merge(std::vector<int>& arr, int l, int m, int r) {
    int n1 = m - l + 1;
    int n2 = r - m;

    std::vector<int> L(n1), R(n2);

    std::copy(arr.begin() + l, arr.begin() + l + n1, L.begin());
    std::copy(arr.begin() + m + 1, arr.begin() + m + 1 + n2, R.begin());

    int i = 0, j = 0, k = l;
    while (i < n1 && j < n2) {
        if (L[i] <= R[j]) {
            arr[k] = L[i];
            i++;
        } else {
            arr[k] = R[j];
            j++;
        }
        k++;
    }

    while (i < n1) {
        arr[k] = L[i];
        i++;
        k++;
    }

    while (j < n2) {
        arr[k] = R[j];
        j++;
        k++;
    }
}

void parallelMergeSort(std::vector<int>& arr, int l, int r, queue& q) {
    if (l >= r) {
        return;
    }
    int m = l + (r - l) / 2;

    q.submit([&](handler& h) {
        h.single_task([=, &arr]() {
            parallelMergeSort(arr, l, m, q);
        });
    }).wait();

    q.submit([&](handler& h) {
        h.single_task([=, &arr]() {
            parallelMergeSort(arr, m + 1, r, q);
        });
    }).wait();

    merge(arr, l, m, r);
}

主函数: 

int main() {
    queue q(default_selector{}, property::queue::enable_profiling());

    std::vector<int> arr;
    std::ifstream in_file("problem-2.txt", std::ios::in);
    while (!in_file.eof()){
        int t;
        in_file>>t;
        arr.emplace_back(t);
    }
    in_file.close();

    parallelMergeSort(arr, 0, arr.size() - 1, q);

    for (int i = 0; i < arr.size(); i++) {
        std::cout << arr[i] << " ";
    }
    std::cout << std::endl;

    return 0;
}

五、作业3:图像卷积并行加速

1. 问题描述

使用基于oneAPI的C++/SYCL实现一个用于计算图像的卷积操作。输入为一个图像矩阵和一个卷积核矩阵,输出为卷积后的图像。

2. 思路分析

图像卷积是一种常见的图像处理操作,用于应用各种滤波器和特征检测器。其原理可以简单地描述为在图像的每个像素上应用一个小的矩阵(通常称为卷积核或滤波器),并将卷积核中的元素与图像中对应位置的像素值相乘,然后将所有乘积的和作为结果。这个过程可以看作是对图像进行了平滑、锐化、边缘检测等操作。 假设有一个大小为 M×N 的输入图像 I 和一个大小为 m×n 的卷积核 K 。图像卷积操作可以用下面的数学公式来表示:

其中,S(i, j)是卷积操作的结果图像中位置 (i, j) 处的像素值。I(i + k, j + l) 是图像中位置 (i + k, j + l) 处的像素值,K(k, l) 是卷积核中位置 (k, l) 处的权重。卷积核通常是一个小的⼆维矩阵,用于捕捉图像中的特定特征。在异构计算编程中,可以使用并行计算来加速图像卷积操作。通过将图像分割成小块,然后在GPU上并行处理这些块,可以实现高效的图像卷积计算。通过合理的块大小和线程组织方式,可以最大限度地利用GPU的并行计算能力来加速图像处理过程。 基于GPU的图像卷积操作的原理基于并行处理和矩阵乘法的基本原理,通过将图像数据和卷积核数据分配给不同的线程块和线程,利用GPU的并行计算能力实现对图像的快速处理。

3. 代码实现

#include <CL/sycl.hpp>
#include <vector>
#include <iostream>

using namespace sycl;

void conv2d(const std::vector<float>& input, std::vector<float>& output, const std::vector<float>& kernel, int width, int height, int kernelSize) {
    int pad = kernelSize / 2;

    queue q(default_selector{});

    buffer<float, 1> inputBuffer(input.data(), range<1>(width * height));
    buffer<float, 1> outputBuffer(output.data(), range<1>(width * height));
    buffer<float, 1> kernelBuffer(kernel.data(), range<1>(kernelSize * kernelSize));

    q.submit([&](handler& h) {
        auto inputAcc = inputBuffer.get_access<access::mode::read>(h);
        auto outputAcc = outputBuffer.get_access<access::mode::write>(h);
        auto kernelAcc = kernelBuffer.get_access<access::mode::read>(h);

        h.parallel_for(range<2>(height, width), [=](id<2> idx) {
            int x = idx[1];
            int y = idx[0];
            float sum = 0.0f;

            for (int ky = -pad; ky <= pad; ky++) {
                for (int kx = -pad; kx <= pad; kx++) {
                    int ix = x + kx;
                    int iy = y + ky;

                    if (ix >= 0 && ix < width && iy >= 0 && iy < height) {
                        sum += inputAcc[iy * width + ix] * kernelAcc[(ky + pad) * kernelSize + (kx + pad)];
                    }
                }
            }

            outputAcc[y * width + x] = sum;
        });
    });

    q.wait();
}

int main() {
    int width = 5, height = 5, kernelSize = 3;
    std::vector<float> input(width * height, 1.0f);
    std::vector<float> output(width * height, 0.0f);
    std::vector<float> kernel = {1, 1, 1, 1, 1, 1, 1, 1, 1};

    conv2d(input, output, kernel, width, height, kernelSize);

    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            std::cout << output[y * width + x] << " ";
        }
        std::cout << std::endl;
    }

    return 0;
}

sample运行结果如下:

六、总结心得

学习了oneAPI的使用以及直接利用Developer Cloud平台中的CPU与GPU硬件,体验到了利用oneAPI进行C++/SYCL编程的强大之处。英特尔 oneAPI 基础工具包包括编译器、性能库、分析和 debug 工具以及一个兼容性工具,可以帮助开发者把在 CUDA 上编写的代码迁移到 DPC++,加速了在 XPU 架构(如 CPU、GPU、FPGA 和 AI)上高性能计算(HPC)、机器学习、嵌入式计算和计算密集型应用程序的异构并行编程。

  • 1
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值