常见异构程序设计语言

目录

一、OpenMP

二、MPI

三、CUDA/HIP

四、OpenACC

五、Athread

六、OpenCL

七、oneAPI

---

        20世纪80年代，异构计算技术就已经诞生了。异构就是CPU、DSP、GPU、ASIC、协处理器、FPGA等各种计算单元、使用不同的类型指令集、不同的体系架构的计算单元，组成一个混合的系统，执行计算的特殊方式，这种计算方式叫做“异构计算”。不同异构计算单元有各自的体系架构，如何利用异构计算单元的算力资源，需要通过对应体系架构的编程来实现。针对目前主流异构计算体系，相应的编程语言有OpenMP、MPI、CUDA/HIP、OpenACC、OpenCL等。下面对代表性异构编程语言进行简单介绍。

一、OpenMP

        OpenMP是共享内存型计算架构的多线程并行编程语言，主要由编译制导语句、运行时库函数和环境变量等组成，可提供一个可移植、可扩展的共享内存应用编程模型。OpenMP支持Fortran和C/C++。

共享内存访问模型

        OpenMP共享内存并行编程模型中，所有处理器都可以用访问共享内存并使用多个线程进行计算。OpenMP使用Fork-Join模型来实现并行，主要机制是使用一个线程作为主线程启动程序，主线程进入并行区域，派生多个从线程，与主线程组成一个工作组，共同执行并行代码。离开并行区域后，从线程将被置于睡眠状态，直到程序抵达下一个并行区域。

执行机制

        OpenMP使用编译制导语句提供一套对并行算法更高层次的抽象实现，包括指定并行区域、指定循环区域、指定临界区、设置同步等各种指令。OpenMP以#pragma omp作为编译制导语句的标识符，由编译器自动对程序实现并行化，降低并行编程的难度。但是OpenMP不适合处理复杂同步和互斥操作的场景。

二、MPI

        目前成规模的计算组织方式以集群形式为主，采用多节点的分布式存储结构，这种结构可扩展计算节点，节点之间并行操作时需要通过通信来完成，主流的方式是消息传递。在这种形式的计算集群中，并行编程主要采用消息传递接口（Message Passing Interface，MPI）。MPI是消息传递库的标准规范，提供了大量消息传递例程，MPI进程可以在分布式环境中实现跨节点数据通信。MPI最大的优点就是可移植性和易用，基于MPI实现的并行程序组件，允许应用程序、软件库和其他工具在不同机器之间方便快速的移植。

执行机制

        MPI-1协议中提供了多种点对点通信和集体通信的例程。点对点通信包括阻塞和非阻塞两种形式，根据数据管理及发送方、接收方之间的同步方式，可分为标准、就绪、缓冲、同步四种通信模式。集体通信则是通信域中多个进程的操作，包括广播、散发、收集、组收集等聚合通信操作。MPI-2协议中增加了三个MPI编程模型的新特性：其一是并行I/O操作接口，可以实现多个进程同时对文件进行操作；其二是远程内存操作，有效兼容共享内存模型；其三是动态进程管理，运行MPI程序启动后可以继续创建和取消新的进程。MPI-3协议新增一些特性，包括非阻塞式集合通信、近邻集合通信，进一步扩展了共享内存机制。最新的MPI-4协议新增了RAM/单边通信、持久化集合通信等。

三、CUDA/HIP

        CUDA/HIP是当前主流的异构编程模型，其基于C++，可实现通过编程接口访问和操作GPU，完成计算密集型的任务处理。CUDA/HIP不要求开发人员显式地进行线程管理，大幅简化编程模型。

        CUDA/HIP对于线程的管理模式为“网格（grid）——线程块（block）——线程（Thread）”的层次结构。如下图所示，只有在同一个线程块中的线程才可以通过共享内存和线程同步进行协作。线程是独立的执行单元，拥有自己的程序计数器、变量寄存器和处理器状态等。线程块是一组并行任务的集合，会被分配给GPU的流式多处理器（Streaming Multiprocessing，SM）执行，之后线程块会被进一步分成多个线程束（NVIDIA GPU的一个线程束一般有32个线程，AMD GPU一般有64个线程），一般以线程束为单位进行调度和执行。

线程层次结构

        CUDA/HIP编程模型中分为主机端（host）和设备端（device），两者分别维护各自的存储单元，分别称为主机内存和设备内存。CUDA/HIP程序至少需要执行下面三个步骤：

①数据从主机内存（Host memory）复制到设备内存（Device memory）；

②主机端调用设备端核函数（kernel），加载并CUDA/HIP程序；

③结果数据从设备内存复制到主机内存。

        下面是一个简单的CUDA程序：

#include <cuda_runtime.h>

//核函数
__global__ void kernel(DATATYPE* devData, DATATYPE* devRes)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    devRes[i] = devData[i] * 2.0;
}

int main(int argc, char* argv[])
{
    //分配设备内存
    DATATYPE* devData, *devRes;
    cudaMalloc(&devData, sizeof(DATATYPE) * SIZE);
    cudaMalloc($devRes,  sizeof(DATATYPE) * SIZE);
    //数据复制到设备内存
    cudaMemcpy(devData, hostData, sizeof(DATATYPE) * SIZE, cudaMemcpyHostToDevice);
    //核函数
    kernel<<<DIM, BLOCK>>>(devData, devRes);
    //结果复制到主机内存
    cudaMemcpy(hostRes, devRes, sizeof(DATATYPE) * SIZE, cudaMemcpyDeviceToHost);

    return 0;
}

        Thrust是一个提供CUDA/HIP C++模板函数的库，与STL相似，封装了主机端和设备端的多种基本类型和复杂容器的数据结构，如vector；同时还通过排序、规约、遍历等高性能并行算法。

四、OpenACC

        OpenACC和OpenMP类似，也是一组编译器指令，通过指定C/C++、FORTRAN程序中的循环代码区域将计算任务从主机端加载到设备端执行，从而实现跨操作系统、主机端CPU和设备端的可移植性。

执行机制​​​​​

        OpenACC通过编译制导的方式将指定的程序代码加载到加速核心上执行，其余代码仍在主机端执行。例如在神威超算上使用OpenACC加速的异构并行程序可以有效降低应用从核阵列对程序加速的编程难度。

        OpenACC编译制导语法格式：

#pragma acc 编译指示1 [子句···] 编译指示2 [子句···] ···

        每条编译制导语句都需要以“#pragma acc”开头，如下：

float* Vector_Addition(float *a, float *b, float *c, int m)
{
    #pragma acc kernels loop copyin(a[:m], b[0:m]) copyout(c[0:m])
    for(int i = 0; i < m; i++)
    {
        c[i] = a[i] + b[i];
    }
}

五、Athread

        Athread是针对神威众核处理器设计的主从加速编程模型，相比OpenACC编程难度大，但可以灵活对核组内的从核进行控制和调度，提供更小的并行粒度，充分发挥从核阵列的加速性能。Athread与OpenACC类似，使用Athread库进行加速同样是将计算任务加载到从核阵列上去执行。主核加速线程库提供了控制线程组初始化、创建、分配任务和终止环境等供主核程序使用的操作接口；从核加速线程库提供了从核线程标识、核组内同步和DMA读写等供从核程序使用的操作接口。

执行流程

六、OpenCL

        OpenCL（Open Computing Language）是可移植到不同异构平台（例如CPU/GPU/DSP/FPGA）的并行语言标准，可以对超算、云服务器、移动设备和嵌入式平台的加速设备进行跨平台的并行编程。OpenCL与OpenACC不同，需要显式地控制设备，并且全权负责并行化。

        OpenCL的程序一般由两部分组成：主机代码和设备代码。主机代码在主机端执行，将内核代码从主机端提交到OpenCL设备端执行。示例核函数代码如下：

const char *kernelSource = "\n" \
"#program OPENCL EXTENSION cl_khr_fp64:enable  \n" \
"__kernel void vecAdd(__global double *a, \n" \
"                      __global double *b, \n" \
"                      __global double *c, \n" \
"                      const unsigned int n) \n" \
"{                                            \n" \
"    int id = get_global_id(0);              \n" \
"    if(id < n)                               \n" \
"        c[id] = a[id] + b[id];               \n" \
"}                                            \n" \
                                              "\n" ;

        除核心的核函数代码，在执行时，需要首先根据运行时环境匹配OpenCL平台，然后选择一个或多个设备来创建上下文、分配内存、创建设备的命令队列，最后进行数据传输和任务计算。对于给定的上下文，应用程序可以：

        ①创建一个或多个命令队列；

        ②创建即将在指定设备上运行的程序；

        ③创建内核函数；

        ④在主机端或设备端分配内存；

        ⑤将数据从主机端复制到设备端；

        ⑥将核函数提交到命令队列并执行；

        ⑦将结果数据从设备端复制到主机端。

七、oneAPI

        oneAPI是Intel开发的并行编程平台和工具集，其包含了Intel编译器、高性能计算库（MKL）、性能分析和调试工具、异构编程语言。oneAPI面向各种处理器，提供了跨架构的编程语言——Data Parallel C++(DPC++)及编程环境。

        Intel针对高性能计算、深度学习、计算机视觉等，开发了多种配套的工具包，如oneDPL、oneMKL、oneTBB、oneDAL、oneDNN等。下面是一个使用oneMKL接口实现的oneAPI程序，通过调用数学库中的blas::axpy()函数实现快速浮点向量间的乘加操作，并支持在不同平台上的加速计算。

#include <CL/sycl.hpp>
#include <exception>
#include "oneapi/mkl.hpp"
#include "oneapi/mkl/blas.hpp"

int main(int argc, char *argv[])
{
    double alpha = 2.0;
    int len = 1024;
    std::vector<double> x(len);
    std::vector<double y(len*3);

    //省略x、y的初始化

    //初始化设备端信息
    sycl::device my_device;
    my_device = sycl::device(sycl::gpu_selector());

    //异常处理
    auto my_exception_handler = [](sycl::exception_list exceptions)
    {
        for(std::exception_ptr const& e : exceptions)
        {
            try
            {
                std::rethrow_exception(e);
            }
            catch (sycl::exception const& e)
            {
                std::cout<<"Caught asynchronous SYCL exception:\n"<<e.what()<<std::endl;
            }
            catch (std::exception const& e)
            {
                std::cout<<"Caught asynchronous stl exception:\n"<<e.what()<<std::endl;
            }
        }
    };

    //创建设备端执行队列
    sycl::queue my_queue(my_device, my_exception_handler);
    //创建主机端和设备端之间的数据缓存
    sycl::buffer<double, 1> x_buffer(x.data(), x.size());
    sycl::buffer<double, 1> y_buffer(y.data(), y.size());

    //调用oneMKL中的blas库接口，执行y=alpha*x + y
    oneapi::mkl::blas::axpy(my_queue, len, alpha, x_buffer, 1, y_buffer, 3);
    my_queue.wait_and_throw();

    //打印结果
    auto y_accessor = y_buffer.template get_access<sycl::access::mode::read>();
    std::cout<<"y"<<"=["<<y_accessor[0]<<"]"<<std::endl;
    std::cout<<"["<<y_accessor[3]<<"]"<<std::endl;

    return 0;
}