英特尔oneAPI——统一编程模型和跨架构编程语言DPC++详细介绍与模板匹配算法实例

      本文对oneAPI进行了详细的介绍，包括其出现的背景、模型架构（开发式规范、SYCL、DPC++、oneAPI工具包）。然后对基于跨架构的DPC++编程语言进行详细介绍，包括其四大编程模型、编程流程、简单的矩阵乘法示例。最后，本文介绍了模板匹配算法的基本原理，并基于oneAPI进行了编程实现。

目录

一、oneAPI诞生背景

二、oneAPI是什么

三、oneAPI整体架构

1.oneAPI 开放式规范

2.SYCL规范

3.DPC++介绍

4.oneAPI工具包

四、oneAPI并行编程介绍

1.编程框架

2.编程流程

3.简单示例-矩阵加法

五、基于oneAPI的模板匹配算法

1.模板匹配介绍

2.模板匹配算法

3.基于oneAPI的实现

致谢 

---

一、oneAPI诞生背景

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

二、oneAPI是什么

       Intel oneAPI 是一个跨行业、开放、基于标准的统一的编程模型，旨在提供一个适用于各类计算架构的统一编程模型和应用程序接口。也就是说，应用程序的开发者只需要开发一次代码，就可以让代码在跨平台的异构系统上执行，底层的硬件架构可以是CPU、GPU、FPGA、神经网络处理器，或者其他针对不同应用的硬件加速器等等。由此可见，oneAPI既提高开发效率，又可以具有一定的性能可移植性。

三、oneAPI整体架构

1.oneAPI 开放式规范

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

 图1  oneAPI 开放式规范架构

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项新的特征，包括了为简化代码所做的更新，和更小的代码尺寸。一些主要增加的内容包括：

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

3.DPC++介绍

       oneAPI包含一种全新的跨架构编程语言 DPC++，DPC++基于 C++编写，由一组C++类、模板与库组成，同时兼容 Kronos 的 SYCL 规范，图2给出了DPC++与SYCL、C++关系。同时，intel DPC++兼容性工具可以实现将CUDA代码迁移到DPC++上，其中大约会有80%-90%的代码实现了自动迁移并提供内联注释，很大程度上帮助开发人员减轻代码移植的负担。

图2 DPC++与SYCL、C++关系

       DPC++是一种单一源代码语言，其中主机代码和异构加速器内核可以混合在同一源文件中。在主机上调用 DPC++程序，并将计算加载到加速器。程序员使用熟悉的C++和库结构，并添加诸如工作目标队列、数据管理缓冲区和并行性并行的函数，以指导计算和数据的哪些部分应该被加载。

4.oneAPI工具包

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

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

四、oneAPI并行编程介绍

1.编程框架

       oneAPI编程框架和OpenCL类似，包含平台模型、执行模型、内存模型、编程模型等四个模型，下面分别说明。

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

 图 3  oneAPI平台模型

       执行模型：执行模型基于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’子组。

 图4  ND-range、工作组、子组和工作项之间的关系图

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

图5   oneAPI内存模型 

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

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

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的核函数在Device上完成指定的运算。

       该内核将应用于索引空间中的每个点，内核封装在C++ lambda函数中。DPC++中内核的形式如下：

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

        在该循环中，每个迭代都是完全独立的，并且不分顺序。使用 parallel_for 函数表示并行内核。

（5）将SYCL缓冲区的数据读到Host端。

3.简单示例-矩阵加法

       下面给出一个oneAPI程序的例子vectorAdd_dpcpp.cpp，其功能为计算两个一维向量的相加。编译器使用dpcpp，具体编译命令为：dpcpp vectorAdd_dpcpp.cpp -o vectorAdd_dpcpp。

#include <CL/sycl.hpp>
using namespace sycl;
static const size_t numElements = 50000;
void work(queue &q) {
  std::cout << "Device : "
            << q.get_device().get_info<info::device::name>()
            << std::endl;
  float vector1[numElements] , vector2[numElements] , vector3[numElements];
  auto R = range(numElements);
   for (int i = 0; i < numElements; ++i) {
        vector1[i] = rand()/(float)RAND_MAX;
        vector2[i] = rand()/(float)RAND_MAX;
    }
//2.创建vector1、vector2、vector3向量的SYCL缓冲区；
  buffer vector1_buffer(vector1,R);
  buffer vector2_buffer(vector2,R);
  buffer vector3_buffer(vector3,R);

//3.向Device提交工作（定义了访问缓冲区内存的accessor；）
  q.submit([&](handler &h) {
    accessor v1_accessor (vector1_buffer,h,read_only);
    accessor v2_accessor (vector2_buffer,h,read_only);
    accessor v3_accessor (vector3_buffer,h);
   //4. 调用oneAPI的核函数在Device上完成指定的运算；
        h.parallel_for (range<1>(numElements), [=](id<1> index) {
  //核函数部分，若单独写一个函数，直接使用函数名（参数表）调用即可
      if (index < numElements)
        v3_accessor [index] = v1_accessor [index] + v2_accessor [index];
    });
  }).wait(); //排队等待
   // 5. 将SYCL缓冲区的数据读到Host端，检查误差 
    host_accessor h_c(vector3_buffer,read_only);
    for (int i = 0; i < numElements; ++i) {
     if (fabs(vector1[0] + vector2[0] - vector3[0] ) > 1e-8 ) {
            fprintf(stderr, "Result verification failed at element %d!\n", i);
            exit(EXIT_FAILURE);
        }
    }
}
int main() {  
  try {
    queue q;
    work(q);
  } catch (exception e) {
    std::cerr << "Exception: " << e.what() << std::endl;
    std::terminate();
  } catch (...) {
    std::cerr << "Unknown exception" << std::endl;
    std::terminate();
  }
}

五、基于oneAPI的模板匹配算法

1.模板匹配介绍

       模板匹配是图像处理中最基本、最常用的匹配方法。模板匹配是一项在一幅图像中寻找与模板图像最相似的区域的技术，该项技术可用于物体的定位、识别。由于模板匹配计算量庞大，多在PC机或工控机中实现，存在成本高、体积大、功耗高等缺点，限制了模板匹配应用场景。同时，在现有嵌入式平台进行模板匹配时，计算时间较长，难以满足对系统进行实时性响应的要求。在传统的模板匹配算法基础之上，结合并行计算方面的有关知识设计出一种并行的模板匹配算法，能在很大程度上减少模板匹配算法的执行时间。

2.模板匹配算法

       模板匹配的过程，简单的来讲就是通过模板图像与待匹配图像之间相似度的比较，然后在图像中找到模板图像所在位置的过程。其具体执行过程大致可以描述如下:首先按照像素来比较模板图像与待搜索图像之间的相似度，接着找到其中最大的相似度量区域作为我们要找的匹配位置。模板匹配算法的具体过程：通过把图像块在待搜索的图像上进行滑动的方法，对待搜索的图像块和模板图像进行一步一步的匹配。为了便于理解，可以把算法简单描述成下面这几个步骤：这里我们假设待搜索的图像是一张200×200的图像，而模板图像则用另外一张10×10的图像，那么模板匹配算法的具体步骤可以描述如下：

1. 从待搜索图像的左上角点(0，0)开始，切割出来一块(0，0)至(10，10)的临时图像；
2. 用当前临时图像与模板图像进行比较，比较结果记为c；
3. 对比结果c，记作结果图像(0，0)处的对比值；
4. 向右滑动一个像素，切割待搜索图像从(0，1)至 (10，11)的临时图像，再与模板图像进行比较，并把比较结果记录到结果中；
5. 重复（1）～（4）步直到临时图像切割到待搜索图像的右下角。

过程如图6所示：

 

图 6  模板匹配算法原理图

       设S(x,y)是大小为mxn的匹配图像，T(x,y)是MxN的模板图像，目前常规的模板匹配度量值计算方法有以下几种：

平均绝对差算法（MAD算法）：平均绝对差D(i,j)越小，表明越相似，故只需找到最小的D(i,j)即可确定能匹配的子图位置

绝对误差和算法（SAD算法）：绝对误差和D(i,j)越小，表明越相似，故只需找到最小的D(i,j)即可确定能匹配的子图位置

平均误差平方和算法（MSD算法）：计算子图与模板图的L2距离和的平均值

误差平方和算法（SSD算法）：计算子图与模板图的L2距离和

归一化积相关算法（NCC算法）：通过归一化的相关性度量公式来计算二者之间的匹配程度

3.基于oneAPI的实现

       我们的例子使用了第四种方法即误差平方和算法，计算子图与模板图的L2距离和。具体代码如下所示：

相关代码也可以直接访问git，这里贴上链接：

https://github.com/zly5/Parallel-Computing-Labhttps://github.com/zly5/Parallel-Computing-Lab

%%writefile lab/gpu_sample.cpp
#include <chrono>
#include <cmath>
#include <iostream>
#include "CL/sycl.hpp"
//#include "device_selector.hpp"
// dpc_common.hpp can be found in the dev-utilities include folder.
// e.g., $ONEAPI_ROOT/dev-utilities/<version>/include/dpc_common.hpp
#include "dpc_common.hpp"

// stb/*.h files can be found in the dev-utilities include folder.
// e.g., $ONEAPI_ROOT/dev-utilities/<version>/include/stb/*.h
#define STB_IMAGE_IMPLEMENTATION
#include "stb/stb_image.h"
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb/stb_image_write.h"

using namespace std;
using namespace sycl;

static void ReportTime(const string &msg, event e) {
  cl_ulong time_start =
      e.get_profiling_info<info::event_profiling::command_start>();

  cl_ulong time_end =
      e.get_profiling_info<info::event_profiling::command_end>();

  double elapsed = (time_end - time_start) / 1e6;
  cout << msg << elapsed << " milliseconds\n";
}


// SYCL does not need any special mark-up for functions which are called from
// SYCL kernel and defined in the same compilation unit. SYCL compiler must be
// able to find the full call graph automatically.
// always_inline as calls are expensive on Gen GPU.
// Notes:
// - coeffs can be declared outside of the function, but still must be constant
// - SYCL compiler will automatically deduce the address space for the two
//   pointers; sycl::multi_ptr specialization for particular address space
//   can used for more control
__attribute__((always_inline)) static void ApplyFilter(uint8_t *I,
                                                       uint8_t *T,
                                                       float *result,
                                                       int i,
                                                       int j,
                                                       int Iw,
                                                       int Ih, 
                                                       int Tw,
                                                       int Th)  {
    
  if (i >= Ih - Th + 1 || j >= Iw - Tw + 1) {
      return;
  }
  float sum = 0.0;
  for (int k = 0; k < Th; k++) {
    for (int s = 0; s < Tw; s++) {
        float diff = I[(i + k) * Iw + j + s] - T[k * Tw + s];
        sum += diff * diff;
    }
  }
  result[i * Iw + j] = sum;  
}


int main(int argc, char **argv) {
  // loading the src image
  int src_img_width, src_img_height, src_channels;
  // 使用灰度图像
  // 加载图片 源图片
  uint8_t *src_image = stbi_load("./tmp_src_img.jpg", &src_img_width, &src_img_height, &src_channels, 1);
  if (src_image == NULL) {
    cout << "Error in loading the image\n";
    exit(1);
  }
  cout << "Loaded src image with a width of " << src_img_width << ", a height of "
       << src_img_height << " and " << src_channels << " channels\n";
  
  // loading the template image
  int template_img_width, template_img_height, template_channels;
  // 加载图片 模板图片
  uint8_t *template_image = stbi_load("./tmp_template_img.jpg", &template_img_width, &template_img_height, &template_channels, 1);
  if (template_image == NULL) {
    cout << "Error in loading the image\n";
    exit(1);
  }
  cout << "Loaded template image with a width of " << template_img_width << ", a height of "
       << template_img_height << " and " << template_channels << " channels\n";

  if (src_img_width < template_img_width || src_img_height < template_img_height) {
    cout << "Error: The template is larger than the picture\n";
    exit(1);
  }
  
  
  // 分配的结果内存
  size_t num_counts = src_img_height * src_img_width;
  size_t src_size = src_img_height * src_img_width;
  size_t template_size = template_img_width * template_img_height;
  // 分配输出图像的内存
  // allocating memory for output images
  float *result = new float[num_counts];
  // 初始化
  // memset(image_ref, 0, num_counts * sizeof(float));

  // Create a device selector which rates available devices in the preferred
  // order for the runtime to select the highest rated device
  // Note: This is only to illustrate the usage of a custom device selector.
  // default_selector can be used if no customization is required.
  // 选择合适的设备
  
    //device_selector sel;

  // Using these events to time command group execution
  event e1, e2;
  
  // Wrap main SYCL API calls into a try/catch to diagnose potential errors
  try {
    // Create a command queue using the device selector and request profiling
    // 选择最适合的设备
    auto prop_list = property_list{property::queue::enable_profiling()};
    queue q(default_selector{}, dpc_common::exception_handler, prop_list);

    // See what device was actually selected for this queue.
    cout << "Running on " << q.get_device().get_info<info::device::name>()
         << "\n";
    // 源图像buffer
    buffer src_image_buf(src_image, range(src_size));
    // 模板图像buffer
    // This is the output buffer device writes to
    buffer template_image_buf(template_image, range(template_size));
    // 结果的buffer
    buffer result_buf(result, range(num_counts));
    cout << "Submitting lambda kernel...\n";

    // Submit a command group for execution. Returns immediately, not waiting
    // for command group completion.
    // 得到输出的比较之后的结果
    e1 = q.submit([&](auto &h) {
      accessor src_image_acc(src_image_buf, h, read_only);
      accessor template_image_acc(template_image_buf, h, read_only);
      accessor result_acc(result_buf, h, write_only);
      // 使用二维线程数
      h.parallel_for(range<2>{(size_t)src_img_height, (size_t)src_img_width}, [=](id<2> index) {
		// 内核程序执行
        ApplyFilter(src_image_acc.get_pointer(), template_image_acc.get_pointer(), result_acc.get_pointer(), index[0], index[1], src_img_width, src_img_height, template_img_width, template_img_height);
      });
    });
    q.wait_and_throw();  
  } catch (sycl::exception e) {
    cout << "SYCL exception caught: " << e.what() << "\n";
    return 1;
  }

  // report execution times:
  ReportTime("Lambda kernel time: ", e1);
  // cout << result[0] << " " << result[1];
    
  // 得到匹配位置的最小值
  int x,y;
  float minresult = result[0];
  for (int i = 0; i < src_img_height - template_img_height + 1; i++) {
    for (int j = 0; j < src_img_width - template_img_width + 1; j++) {
        if (minresult > result[i * src_img_width + j]) {
            y = i;
            x = j;
            minresult = result[i * src_img_width + j];
        }
    }
  }
  
  int x1 = x;
  int x2 = x + template_img_width - 1;
  int y1 = y;
  int y2 = y + template_img_height - 1;
    
  cout << x1 << "  " << x2 << "  " << y1 << "  " << y2 << "  ";
  
  // 对图片进行保存
    
  // 先标记两条横线
  for (int i = x1; i <= x2; i++) {
      src_image[y1 * src_img_width + i] = 0;
      src_image[y2 * src_img_width + i] = 0;
  }
  for (int i = y1 + 1; i < y2; i++) {
      src_image[i * src_img_width + x1] = 0;
      src_image[i * src_img_width + x2] = 0;
  }
    
  stbi_write_png("sepia_ref.png", src_img_width, src_img_height, src_channels, src_image, src_img_width * src_channels);
  return 0;
}

致谢 

       最后，感谢教育部-英特尔产学合作专业综合改革项目提供的DevCloud平台支持，感谢英特尔亚太研发有限公司技术团队提供的技术支持。

       有关高性能计算课程及相关资料请参阅以下链接：

       https://faculty.xidian.edu.cn/hmzhu/zh_CN/article/336134/content/1891.htm#article

       https://faculty.xidian.edu.cn/hmzhu/zh_CN/article/336134/content/1221.htm#article