国产加速器海光DCU&GPGPU深算处理器异构编程实战（下）

技术瘾君子1573

已于 2024-08-26 07:00:55 修改

阅读量686

点赞数 27

分类专栏： Linux并行计算&HPC高性能计算文章标签：开发语言 OpenMP Fortran DCU 异构编程

于 2024-08-26 06:30:00 首次发布

本文链接：https://blog.csdn.net/qq_27815483/article/details/141326281

版权

Linux并行计算&HPC高性能计算专栏收录该内容

54 篇文章 4 订阅

订阅专栏

四、使用DCU加速Fortran程序

4.1 概述

4.2 在Fortran程序中使用HIP编程接口

第一二章内容请看国产加速器海光DCU&GPGPU深算处理器异构编程实战（上）-CSDN博客

第三章内容请看国产加速器海光DCU&GPGPU深算处理器异构编程实战（中）-CSDN博客

四、使用DCU加速Fortran程序

Fortran语言自20世纪50年代诞生至今，一直被广泛应用于计算物理学、计算化学、流体力学、气候预测、有限元分析等计算密集的学科领域。同C/C++语言一样，Fortran语言也是科学/工程计算领域的主流编程语言，使用Fortran语言编写的应用程序运行在世界各地大大小小的计算集群上，为人类的各种研究与创新提供数据支撑。

本章将介绍如何让使用Fortran语言编写的程序也能利用DCU进行并行加速。

4.1 概述

在上一章，我们介绍了如何使用C/C++语言编写DCU程序，在3.1.7节中我们总结了编写DCU程序主要有两方面工作要做，其一是使用HIP编程接口实现DCU的设备管理、内存管理、DCU内存与主机内存间的数据传输等功能；另一方面是编写DCU核函数让DCU执行计算任务。我们之所以能够做到这些，是基于HIP提供了一套C/C++的编程接口，并且对C/C++语言做了扩展，使之能够编写DCU核函数。

对于Fortran编写的程序，虽然DTK中也提供了HIP的Fortran编程接口，但没有对Fortran语言提供扩展，不能使用Fortran语言编写核函数。因此如果我们想使用DCU加速Fortran编写的程序，如果是直接编写DCU核函数，只能使用C/C++来编写，再让Fortran程序调用核函数，换言之使用DCU加速Fortran程序实际上是一个Fortran/C互用性（interoperability）的问题；另一条路线则是使用OpenMP，用导语的方式指导编译器将并行计算部分卸载到DCU上。

本章主要介绍使用C/C++编写核函数并在Fortran程序中调用的方式，使用OpenMP的方式将在第5章进行介绍。

4.2 在Fortran程序中使用HIP编程接口

我们知道在C++程序中使用HIP编程接口时需要通过#include在代码中插入头文件hip_runtime.h，使用其它DTK中的库时还需要插入其他相应的头文件，如hiprand.h、hipblas.h等等（这些头文件都在DTK中的include目录下）。在Fortran程序中使用HIP编程接口和其他DTK库时也需要在程序通过USE语句使用相应的Fortran模块。

这些Fortran模块的源代码在23.04版本的DTK中的src/hipfor路径下，主要模块总结在表4-1中。

表4-1 hipfor中的主要模块

模块	主要内容
hipfort	HIP运行时编程接口
hipfort_types	dim3类型和一些枚举类型
hipfort_check	错误检查接口
hipfort_hipfft	hipFFT库的接口
hipfort_hipblas	hipBLAS库的接口
hipfort_hiprand	hipRAND库的接口
hipfort_hipsolver	hipSOLVER库的接口
hipfort_hipsparse	hipSPARSE库的接口

查看这些模块的源代码，你会发现这些模块包含着使用Fortran 2003标准中的内置模块iso_c_binding编写的大量Fortran接口，这些接口绑定到HIP的相关C语言接口，在Fortran中使用HIP接口就是以这种方式实现的。

我们结合Fortran程序来介绍，在3.5.1节中的C程序use_libs.cc中我们实现了使用DCU随机生成两个矩阵并计算矩阵乘。在程序中我们没有编写DCU核函数，只使用了HIP运行时编程接口和hipRAND、hipBLAS两个库中的函数，刚好适合作为本节例子。

我们仿照use_libs.cc的流程，用Fortran语言来实现如下(use_libs.f90)：

PROGRAM MATRIX_MULTIPLICATION
    USE hipfort
    USE hipfort_check
    USE hipfort_hiprand
    USE hipfort_hipblas
    IMPLICIT NONE
    INTEGER :: argc, N
    INTEGER(C_LONG_LONG) :: SEED
    INTEGER(C_SIZE_T) :: num_elements
    CHARACTER(LEN=80) :: argv
    TYPE(C_PTR) :: gen, start, stop_, handle
    REAL(4), TARGET :: dt
    REAL(8) :: alpha, beta
    REAL(8), ALLOCATABLE :: A(:,:), B(:,:), C(:,:)
    REAL(8), POINTER :: dA(:,:), dB(:,:), dC(:,:)
    INTEGER,PARAMETER :: max_size = 5 ! 打印的最大维度
    INTEGER :: TIME ! 对于flang需要声明

    ! 获取命令行参数个数
    argc = COMMAND_ARGUMENT_COUNT()

    IF (argc /= 1) THEN
        CALL GET_COMMAND_ARGUMENT(0,argv)
        WRITE (*,*) "Usage: ", TRIM(argv), " N"
        STOP 1
    ELSE
        CALL GET_COMMAND_ARGUMENT(1,argv)
        READ (argv, *) N
        IF (N < 2) THEN
            WRITE (*,*) "N must larger than 2."
            STOP 1
        END IF
    END IF

    ! 初始化hipRAND随机数发生器
    CALL hipCheck(hiprandCreateGenerator &
                  (gen, HIPRAND_RNG_PSEUDO_DEFAULT))
    SEED = TIME()
    CALL hipCheck(hiprandSetPseudoRandomGeneratorSeed(gen, SEED))

    ! 计时器
    CALL hipCheck(hipEventCreate(start))
    CALL hipCheck(hipEventCreate(stop_))

    ! 分配矩阵内存空间
    ALLOCATE(A(N,N), B(N,N), C(N,N))
    CALL hipCheck(hipMalloc(dA, N, N))
    CALL hipCheck(hipMalloc(dB, N, N))
    CALL hipCheck(hipMalloc(dC, N, N))

    ! 随机生成矩阵A和B
    CALL hipCheck(hipEventRecord(start,C_NULL_PTR))
    num_elements = N * N
    CALL hipCheck(hiprandGenerateUniformDouble &
                  (gen, C_LOC(dA), num_elements))
    CALL hipCheck(hiprandGenerateUniformDouble &
                  (gen, C_LOC(dB), num_elements))
    CALL hipCheck(hipEventRecord(stop_,C_NULL_PTR))
    CALL hipCheck(hipEventSynchronize(stop_))
    CALL hipCheck(hipEventElapsedTime(C_LOC(dt),start,stop_))
    WRITE(*,'("Generating Matrices took ", F8.3, " ms")') dt

    ! 打印矩阵A,B
    CALL hipCheck(hipMemcpy(A, dA, N * N, hipMemcpyDeviceToHost))
    CALL hipCheck(hipMemcpy(B, dB, N * N, hipMemcpyDeviceToHost))
    WRITE(*,*) "Matrix A:"
    CALL print_matrix(N, A)
    WRITE(*,*)
    WRITE(*,*) "Matrix B:"
    CALL print_matrix(N, B)
    WRITE(*,*)

    ! 创建hipBLAS句柄
    CALL hipCheck(hipblasCreate(handle))

    ! 矩阵乘法，将结果存储在矩阵C中
    alpha = 1.0
    beta = 0.0
    ! 预热
    CALL hipCheck(hipblasDgemm(handle, &
                  HIPBLAS_OP_N, HIPBLAS_OP_N, 2, 2, 2, &
                  alpha, dA, N, dB, N, beta, dC, N))
    ! 正式计算
    CALL hipCheck(hipEventRecord(start,C_NULL_PTR));
    CALL hipCheck(hipblasDgemm(handle, &
                  HIPBLAS_OP_N, HIPBLAS_OP_N, N, N, N, &
                  alpha, dA, N, dB, N, beta, dC, N))
    CALL hipCheck(hipEventRecord(stop_,C_NULL_PTR))
    CALL hipCheck(hipEventSynchronize(stop_))
    CALL hipCheck(hipEventElapsedTime(C_LOC(dt), start, stop_))
    WRITE(*,'("Matrix multiplication (hipblasDgemm) ", &
          "took " F8.3, " ms")') dt

    ! 打印矩阵C
    CALL hipCheck(hipMemcpy(C, dC, N * N, hipMemcpyDeviceToHost))
    WRITE(*,*) "Matrix C:"
    CALL print_matrix(N,C)

    ! 清理
    CALL hipCheck(hipFree(dA))
    CALL hipCheck(hipFree(dB))
    CALL hipCheck(hipFree(dC))
    CALL hipCheck(hipEventDestroy(start))
    CALL hipCheck(hipEventDestroy(stop_))
    CALL hipCheck(hiprandDestroyGenerator(gen))
    CALL hipCheck(hipblasDestroy(handle))
    DEALLOCATE(A,B,C)

    CONTAINS

    SUBROUTINE print_matrix(N, M)
        IMPLICIT NONE
        REAL(8) :: M(N,N)
        INTEGER :: N, i, j
        INTEGER,PARAMETER :: max_size = 5

        DO i = 1, MIN(N, max_size)
            DO j = 1, MIN(N, max_size)
                WRITE(*,'(F13.6)',ADVANCE='NO') M(i,j)
            ENDDO
            IF (N > max_size) THEN
                WRITE(*,'(" ...")',ADVANCE='NO')
            ENDIF
            WRITE(*,*)
        ENDDO
        IF (N > max_size) THEN
            WRITE(*,'("."/ "." / ".")')
        ENDIF
    END SUBROUTINE
END PROGRAM

在上面的代码中，我们使用了hipfort、hipfort_check、hipfort_hiprand和hipfort_hipblas四个Fortran模块，此后就可以像在C语言中一样方便的使用HIP运行时API和hipRAND、hipBLAS两个库中的函数了。

4.3 hipfor的使用方法

我们先编译上节中编写的Fortran程序use_libs.f90，在程序中我们通过USE语句使用了hipfor中的4个模块，但我们需要首先编译出这些模块。此时我们有两种方式，其一是在代码的第一行添加INCLUDE语句来插入hipfor的源代码：

INCLUDE `hipfort.f`

因为我们的程序以及hipfor的代码都是Fortran的标准代码，因此不需要像C代码那样只能使用hipcc进行编译，可以使用GCC中的gfortran、Intel的ifort或DTK中包含的flang。

DTK中之所以没有提供编译好的hipfor模块，原因在于Fortran模块在不同编译器甚至不同编译器版本间不通用。

我们使用gfortran编译如下(需要加载DTK-23.04环境，如使用其他DTK版本，需要将$ROCM_PATH替换为23.04版本根目录的绝对路径)：

gfortran use_libs.f90 -o use_libs -I$ROCM_PATH/src/hipfor -L$ROCM_PATH/lib -lamdhip64 -lhiprand -lhipfft -lhipblas -lhipsolver -lhipsparse

这里的-I$ROCM_PATH/src/hipfor是让编译器能找到hipfort.f文件，此外hipfor模块中绑定了大量的HIP编程接口，因此还需要链接一系列DTK中的库。

我们的程序只包含一个文件，按这种方式使用hipfor就足够了，如果是包含多个源文件的项目，则至少需要保证编译的第一个使用hipfor的文件中插入hipfort.f，否则编译器会报找不到模块的错误。

因此，更推荐按第二种方式使用hipfor，即首先将hipfor构建好。将DTK-23.04中src路径下的hipfor目录整个拷贝到工作路径下，比如说use_lib.f90所在的路径。最好将hipfort.f重命名为hipfort.f90，这样告诉编译器hipfor中的代码是使用自由格式编写的。

用下面的命令将hipfor编译出静态库libhipfort.a以及一系列.mod文件（以使用gfortran为例，换成ifort或flang同样可以）：

gfortran hipfort.f90 -c -o libhipfort.a

再进一步编译我们的程序，并与libhipfort.a链接在一起，注意到这次的-I./hipfor是为了让编译器找到.mod后缀的模块文件：

gfortran use_libs.f90 ./hipfor/libhipfort.a -o use_libs -I./hipfor -lamdhip64 -lhiprand -lhipfft -lhipblas -lhipsolver -lhipsparse

我们运行程序，结果与C语言编写的版本是一样的，都用上了DCU：

./use_libs 4096
Generating Matrices took    5.505 ms
 Matrix A:
     0.008903     0.629933     0.709649     0.638452     0.644310 ...
     0.342506     0.346941     0.301005     0.960350     0.240107 ...
     0.249814     0.120858     0.016834     0.924211     0.464020 ...
     0.722379     0.389470     0.720806     0.880512     0.681257 ...
     0.501467     0.831721     0.861388     0.547661     0.694420 ...
.
.
.

 Matrix B:
     0.005949     0.260324     0.646069     0.458477     0.029385 ...
     0.675923     0.769810     0.819317     0.350884     0.333020 ...
     0.277287     0.791301     0.326005     0.550576     0.047234 ...
     0.380112     0.276836     0.119294     0.472814     0.783367 ...
     0.479333     0.519103     0.606373     0.648004     0.932445 ...
.
.
.

Matrix multiplication (hipblasDgemm) took   28.362 ms
 Matrix C:
  1011.616693  1028.564987  1023.063039  1033.902850  1031.037938 ...
  1023.053500  1039.319391  1027.577098  1039.415822  1041.772753 ...
  1025.272843  1034.920788  1024.315198  1046.460633  1021.596934 ...
  1007.171633  1026.780537  1018.048457  1020.943283  1021.284208 ...
  1018.007430  1034.625873  1025.411990  1030.159926  1039.294852 ...
.
.
.

最后做一些补充说明：

由于hipfor中的接口都是Fortran中的FUNCTION函数，返回值是HIP错误类型，如hipSuccess，因此在Fortran中直接调用是不符合语法的：
```
 hipMalloc(dA, N, N) !不符合Fortran语言的语法
```
为此应该像我们的代码中那样，使用hipfort_check模块（在文件hipfort_check.f）中的SUBROUTINE子例程hipCheck，它会在接口的返回值不是hipSuccess时报出错误码并退出程序。或者至少定义一个整型变量来接收接口的返回值：
```
 ierr = hipMalloc(dA, N, N)
```
对于hipfor中的接口，DTK中暂时没有提供像C语言版本的接口那样的详细说明文档，因此在使用接口时可能需要在hipfor中搜索相关的源代码，比如说在use_libs.f90中，需要创建hipBLAS的句柄handle时，我们搜索hipfor的源码，找到hipblasCreate的接口：

```fortran interface hipblasCreate
```
 function hipblasCreate_orig(handle) bind(c, name="hipblasCreate")
   use iso_c_binding
   use hipfort_hipblas_enums
   implicit none
   integer(kind(HIPBLAS_STATUS_SUCCESS)) :: hipblasCreate_orig
   type(c_ptr) :: handle
 end function
```

  end interface
```

从而了解到需要将handle的类型声明为`iso_c_bingding`中的C指针类型`C_PTR`。

hipfor中的接口普遍是泛型接口（generic interface），以此来实现接口风格更贴近于Fortran语言，比如在我们的代码中hipMalloc接口的参数可以是REAL(8)类型的二维数组指针：
```
 CALL hipCheck(hipMalloc(dA, N, N))
```
而不必使用以下C语言风格的接口（虽然效果是一样的）：
```
 ! 通过C_LOC取指向dA首地址的C指针，或者直接将dA定义为C_PTR类型
 ! 8_c_size_t 意思是kind为c_size_t的整数8
 ! 因为双精度浮点数是8个字节
 CALL hipCheck(hipMalloc(C_LOC(dA), N * N * 8_c_size_t))
```
对于hipfor中每个接口可以接收何种参数，目前也需要查看源代码，并灵活使用。

4.4 在Fortran程序中调用核函数

下面介绍在Fortran程序中如何调用核函数。回顾前面的内容，我们之所以能通过hipfor在Fortran中调用HIP编程接口，是因为hipfor的模块中包含了使用iso_c_binding内置模块编写的一系列Fortran接口，绑定到HIP的C语言编程接口和库函数上。因此如果我们想在Fortran程序中调用核函数，建议也采用这种方式。

在3.5.4节的程序more_elements.cc中，我们编写了性能不错的核函数。在本节我们就编写一个Fortran版本的，并调用之前编写的核函数。

我们先将more_elements.cc中的设备端函数以及核函数写到一个文件kernels.cc中，由于宏定义既要在核函数中使用，又要在主程序中使用，因此也写到一个独立的头文件config.h，并分别插入kernels.cc以及Fortran源文件中。config.h中直接采用最优配置：

#ifndef CONFIG_H
#define CONFIG_H

#define Nsub  128
#define Nsub_ 16
#define BLOCKDIM_X 16
#define BLOCKDIM_Y 16
#define Nx Nsub / BLOCKDIM_X
#define Ny Nsub / BLOCKDIM_Y

#endif

此外由于Fortran中无法使用三个尖括号<<<>>>的方式调用核函数，因此核函数的调用也要写在kernels.cc中：

extern "C"
{
    void matrix_multiply_kernel(dim3* grid, dim3* block, int N,
                                int Nblk, double* A, double* B,
                                double* C)
    {
        matrix_multiply<<<*grid,*block>>>(N, Nblk, A, B, C);
    }
}

因为more_elements.cc中原有的核函数命名为matrix_multiply，为了避免重复，将启动核函数的函数命名为matrix_multiply_kernel，extern "C"是告诉编译器禁掉C++的特性，并按C语言的规则来命名，这样保证能和我们的Fortran程序正确链接。我们将kernels.cc编译成目标文件，因为文件中是HIP核函数，因此需要使用hipcc进行编译：

hipcc -c kernels.cpp

这样就生成了包含核函数调用的kernels.o文件。现在编写Fortran程序如下（more_elements.f90）：

#include "config.h"
PROGRAM MATRIX_MULTIPLICATION
    USE hipfort
    USE hipfort_check
    USE hipfort_hiprand
    USE hipfort_hipblas
    USE hipfort_types
    IMPLICIT NONE
    INTEGER :: argc, N, Nblk
    INTEGER(C_LONG_LONG) :: SEED
    INTEGER(C_SIZE_T) :: num_elements
    CHARACTER(LEN=80) :: argv
    TYPE(C_PTR) :: gen, start, stop_, handle
    REAL(4), TARGET :: dt
    REAL(8) :: alpha, beta
    REAL(8), ALLOCATABLE :: C_kernel(:,:), C_hipblas(:,:)
    REAL(8), POINTER :: dA(:,:), dB(:,:), dC(:,:)
    INTEGER,PARAMETER :: max_size = 5 ! 打印的最大维度
    INTEGER :: TIME ! 对于flang需要声明
    TYPE(dim3) :: numBlocks, blockSize

    INTERFACE
    SUBROUTINE matrix_multiply_kernel(grid, block, &
                                      N, Nblk, A, B, C) bind(c)
    USE hipfort_types
    IMPLICIT NONE
    TYPE(dim3) :: grid, block
    INTEGER, VALUE :: N, Nblk
    REAL(8) :: A(N,N), B(N,N), C(N,N)
    END SUBROUTINE
    END INTERFACE

    ! 获取命令行参数个数
    argc = COMMAND_ARGUMENT_COUNT()

    IF (argc /= 1) THEN
        CALL GET_COMMAND_ARGUMENT(0,argv)
        WRITE (*,*) "Usage: ", TRIM(argv), " N"
        STOP 1
    ELSE
        CALL GET_COMMAND_ARGUMENT(1,argv)
        READ (argv, *) N
        IF (N < 2) THEN
            WRITE (*,*) "N must larger than 2."
            STOP 1
        END IF
    END IF

    ! 初始化hipRAND随机数发生器
    CALL hipCheck(hiprandCreateGenerator &
                  (gen, HIPRAND_RNG_PSEUDO_DEFAULT))
    SEED = TIME()
    CALL hipCheck(hiprandSetPseudoRandomGeneratorSeed(gen, SEED))

    ! 计时器
    CALL hipCheck(hipEventCreate(start))
    CALL hipCheck(hipEventCreate(stop_))

    ! 分配矩阵内存空间
    ALLOCATE(C_kernel(N,N), C_hipblas(N,N))
    CALL hipCheck(hipMalloc(dA, N, N))
    CALL hipCheck(hipMalloc(dB, N, N))
    CALL hipCheck(hipMalloc(dC, N, N))

    ! 随机生成矩阵A和B
    CALL hipCheck(hipEventRecord(start,C_NULL_PTR))
    num_elements = N * N
    CALL hipCheck(hiprandGenerateUniformDouble &
                  (gen, C_LOC(dA), num_elements))
    CALL hipCheck(hiprandGenerateUniformDouble &
                  (gen, C_LOC(dB), num_elements))
    CALL hipCheck(hipEventRecord(stop_,C_NULL_PTR))
    CALL hipCheck(hipEventSynchronize(stop_))
    CALL hipCheck(hipEventElapsedTime(C_LOC(dt),start,stop_))
    WRITE(*,'("Generating Matrices took ", F8.3, " ms")') dt

    ! 创建hipBLAS句柄
    CALL hipCheck(hipblasCreate(handle))

    ! 矩阵乘法，将结果存储在矩阵C中
    alpha = 1.0
    beta = 0.0
    ! 预热
    CALL hipCheck(hipblasDgemm(handle, &
                  HIPBLAS_OP_N, HIPBLAS_OP_N, 2, 2, 2, &
                  alpha, dA, N, dB, N, beta, dC, N))
    ! 正式计算
    CALL hipCheck(hipEventRecord(start,C_NULL_PTR));
    CALL hipCheck(hipblasDgemm(handle, &
                  HIPBLAS_OP_N, HIPBLAS_OP_N, N, N, N, &
                  alpha, dA, N, dB, N, beta, dC, N))
    CALL hipCheck(hipEventRecord(stop_,C_NULL_PTR))
    CALL hipCheck(hipEventSynchronize(stop_))
    CALL hipCheck(hipEventElapsedTime(C_LOC(dt), start, stop_))
    WRITE(*,'("Matrix multiplication (hipblasDgemm) ", &
          "took " F8.3, " ms")') dt

    ! 将hipBLAS计算的结果拷贝回主机
    CALL hipCheck(hipMemcpy(C_hipblas, dC, SIZE(dC), &
                            hipMemcpyDeviceToHost))

    ! 将dC清零
    CALL hipCheck(hipMemset(C_loc(dC), 0, SIZE(dC) * 8_c_size_t))

    ! 使用核函数计算
    blockSize = dim3(BLOCKDIM_X, BLOCKDIM_Y, 1)
    Nblk = (N + Nsub - 1) / Nsub
    numBlocks = dim3(Nblk, Nblk, 1)
    CALL hipCheck(hipEventRecord(start,C_NULL_PTR));
    CALL matrix_multiply_kernel(numBlocks, blockSize, &
                                N, Nblk, dA, dB, dC)
    CALL hipCheck(hipEventRecord(stop_,C_NULL_PTR))
    CALL hipCheck(hipEventSynchronize(stop_))
    CALL hipCheck(hipEventElapsedTime(C_LOC(dt), start, stop_))
    WRITE(*,'("Matrix multiplication (kernel) ", &
          "took " F8.3, " ms")') dt

    ! 将核函数计算的结果拷贝回主机并与hipBLAS的结果做比较
    CALL hipCheck(hipMemcpy(C_kernel, dC, SIZE(dC), &
                            hipMemcpyDeviceToHost))
    CALL compare_matrices(N, C_kernel, C_hipblas)

    ! 清理
    CALL hipCheck(hipFree(dA))
    CALL hipCheck(hipFree(dB))
    CALL hipCheck(hipFree(dC))
    CALL hipCheck(hipEventDestroy(start))
    CALL hipCheck(hipEventDestroy(stop_))
    CALL hipCheck(hiprandDestroyGenerator(gen))
    CALL hipCheck(hipblasDestroy(handle))
    DEALLOCATE(C_kernel,C_hipblas)

    CONTAINS

    SUBROUTINE compare_matrices(N, C1, C2)
        IMPLICIT NONE
        REAL(8) :: C1(N,N), C2(N,N)
        INTEGER :: N, i, j
        REAL(8) :: max_error, error

        max_error = 0.0
        DO j = 1, N
            DO i = 1, N
                error = ABS(C1(i,j) - C2(i,j))
                IF (error > max_error) max_error = error
            ENDDO
        ENDDO
        WRITE(*,'("Max error: ", ES11.4)') max_error
    END SUBROUTINE
END PROGRAM

程序的流程与之前的C语言版本完全一样，就是使用hipBLAS和我们自己编写的核函数分别计算矩阵dA和dB的乘积，并将结果拷贝到主机端进行比较。这里的重点就是INTERFACE语句块中matrix_multiply_kernel的接口，其中的bind(c)意思是让链接器在链接时去寻找按C语言规则命名的matrix_multiply_kernel函数；其次因为用到了dim3类型，所以需要USE模块hipfort_types；最后注意Fortran接口的参数类型要与kernels.cc中函数实现的参数类型一一对应，因为Fortran语言在函数传参时默认传地址，因此Fortran中的TYPE(dim3)类型的参数对应C中的dim3*、REAL(8)对应C中的double*，而像N、Nblk这些直接传值的参数，则需要为Fortran接口中的参数类型添加VALUE属性。

现在编译more_elements.f90，因为使用了宏定义，需要添加编译选项-cpp，并与kernels.o、静态库libhipfort.a（假设仍在上节的目录中，并已编译完成）以及DTK中的动态库链接在一起：

gfortran -cpp more_elements.f90 kernels.o ./hipfor/libhipfort.a -o more_elements -I./hipfor -L$ROCM_PATH/lib -lamdhip64 -lhiprand -lhipfft -lhipblas -lhipsolver -lhipsparse

运行如下：

./more_elements 4096
Generating Matrices took    5.561 ms
Matrix multiplication (hipblasDgemm) took   28.494 ms
Matrix multiplication (kernel) took   23.517 ms
Max error:  0.0000E+00

这样就实现了在Fortran程序中调用核函数。

4.5 小节

本节中介绍了如何在Fortran程序中调用HIP编程接口以及核函数，让Fortran程序也能利用DCU进行加速。主要的技术路线就是利用Fortran和C的互用性，以及Fortran2003标准中的内置模块iso_c_binding，无论是使用HIP的Fortran接口库hipfor调用HIP编程接口，还是调用我们自己编写的核函数，本质上都是按照这一技术路线来实现的。因此，掌握Fortran和C的互用性是比较重要的，本书限于篇幅，没有详细展开，读者如果需要可以多查阅相关的资料。

对于核函数的参数是类的情况，建议使用iso_c_binding建立C语言中类与Fortran语言中派生类型（derived type）的对应关系，避免C语言和Fortran语言使用不同编译器时类中的成员在内存中的布局不一致。

如需使用hipfor中尚未包含的HIP编程接口，需要在Fortran中自行编写接口，hipfor中的写法以及4.4节中的写法可以作为参考。

五、OpenMP编程简介

OpenMP(Open Multi-Processing) 是一种用于共享内存的多处理器并行系统程序设计的编程模型。共享内存是指多个CPU访问同一个共享存储器。在基于共享内存编程模型的OpenMP中,可以通过导语#pragma omp (Fortran语言导语为!$omp)调用库函数并设置相关环境参数后，将单个线程运行的任务（图3-37）变为多个线程执行的并行任务(图3-38)，这些线程间的通信通过对共享内存的读写来完成的。OpenMP目前支持C语言、C++语言和Fortran语言，而当选择忽略这些导语，或者所用编译器不支持OpenMP时，程序自动变为串行语句，仍可以正常的编译和运行。

基于OpenMP的特点，使得它具有优化改动小、编程简单、移植性好、可扩展性高等优点。随着计算机工业的高速发展，GPU加速卡也相继问世，比起传统的CPU端并行加速，使用GPU计算卡可以使程序的计算速度提升几倍甚至几十倍。

2013年7月，OpenMP 4.0发布，可以支持向量化和异构计算，OpenMP目前已更新到OpenMP 5.1版本。使用OpenMP进行异构卡编程其流程如图3-39所示，LLVM13中的clang/flang编译器已支持OpenMP的绝大部分功能特性，且都可以在CPU-DCU异构平台上使用。篇幅所限，本节中主要介绍OpenMP的DCU异构用法，其他用法请参考资相关资料。本节中的代码示例使用C语言，OpenMP Fortran语言的用法类似，具体请参考相关资料。

图 3-42 OpenMP编程模型

使用OpenMP在异构计算过程中具体包括计算任务转移、设备和环境参数的设置、数据的传输和映射、建立目标设备的数据环境、共享、异步等多个步骤，下面将为大家一一介绍。

5.1 OpenMP调用DCU

来看下面一段数组相乘的程序，程序将v1和v2两个数组初始化后，相乘，并赋值给数组p：

    //FileName:normal_test.c
    void vec_mult(int N)
    {
        int i;
        float p[N], v1[N], v2[N];
        init(v1, v2, N);

        for (i=0; i<N; i++)
        {
            p[i] = v1[i] * v2[i];
        }
    }

由于循环中的数组元素的计算彼此独立，适合使用OpenMP改写成并行计算，改写后的程序如下所示：

//target_test1.c
void vec_mult(int N)
{
    int i;
    float p[N], v1[N], v2[N];
    init(v1, v2, N);

    #pragma omp target map(to:v1[0:N],v2[0:N]) map(from:p[0:N])
    {
        #pragma omp parallel for
        for (i=0; i<N; i++)
        {
            p[i] = v1[i] * v2[i];
        }
    }
}

在上面的例子中，数组v1，v2，p最初驻留在主机内存中，当遇到OpenMP制导语句的target(#pragma omp target)结构时，会执行如下操作，首先将计算任务交给设备，即DCU加速卡，并由关键字map把相应的数据从主机存储映射到设备存储中，然后在设备中，使用#pragma omp parallel for语句用来分配for循环的迭代，等到计算完成后，数组p的数据由map子句控制从设备内存返回到主机内存中，完成一次异构并行计算。

与HIP的差异

我们把normal_test.c中通过OpenMP改写后的代码用HIP C来实现，以清晰展示编译器帮助我们做的工作。具体代码如下：

//hip_code.c
__global__ void add(float *d_p, float *d_v1, float *d_v2)
{
  int tid = threadIdx.x + blockIdx.x * blockDim.x;
  if (tid < N)
  {
        d_v2[tid] = d_p[tid] + d_v1[tid];
  }
}

void vec_mult(int N)
{
  int i;

  //申请数据空间
  float *p = (float *) malloc(N * sizeof(float));
  float *v1 = (float *) malloc(N * sizeof(float));
  float *v2 = (float *) malloc(N * sizeof(float));
  float *d_p = NULL;
  float *d_v1 = NULL;
  float *d_v2 = NULL;
  hipMalloc((void **) &d_p, N * sizeof(float));
  hipMalloc((void **) &d_v1, N * sizeof(float));
  hipMalloc((void **) &d_v2, N * sizeof(float));

  //数据初始化
  init(v1, v2, p, N);
  hipMemcpy(d_p, p, sizeof(float) * N, hipMemcpyHostToDevice);
  hipMemcpy(d_v1, v1, sizeof(float) * N, hipMemcpyHostToDevice);
  hipMemcpy(d_v2, v2, sizeof(float) * N, hipMemcpyHostToDevice);
  dim3 blocksize(256, 1);
  dim3 gridsize(N / 256 + 1, 1);

  // 进行数组相加
  add<<<gridsize, blocksize >>> (d_p, d_v1, d_v2);

  // 数据从设备端传回主机端
  hipMemcpy(v2, d_v2, sizeof(float) * N, hipMemcpyDeviceToHost);

  //释放申请空间
  free(p);
  free(v1);
  free(v2);
  hipFree(d_p);
  hipFree(d_v1);
  hipFree(d_v2);
}

通过对比改造成并行计算的OpenMP代码和HIP代码，可以很明显的看出OpenMP编程模型的便捷性，它只需要加几行导语，并不需要对代码进行非常大的改动，就可以达到一个不错的加速效果，反观HIP就略显繁琐。

在OpenMP代码的实际运行过程中，map子句会提醒系统，将子句中的参数在设备端创建一个临时变量并分配相应的空间，然后将该变量的值传递给临时变量，在计算完成中后，系统会根据map子句的关键字来选择是否要将设备端数据映射到主机端，并自动释放这些临时变量。

在OpenMP中同样设定blocksize和gridsize等变量，可以通过使用omp_set_num_threads()等函数或设置OMP_NUM_THREADS参数来设置可调用线程数。但是与HIP不同的是，OpenMP中只可以设置线程数的大小，并不能像HIP中设置block和grid等资源的维度。

总结一下，与HIP相比，使用OpenMP编程模型来进行代码的并行计算只需要加几行导语，操作简单，改动小是其优点。但是HIP上更加丰富的操作和语法，会使我们对程序的控制更加的合理和精细。这两种编程模型各有优点，在改动时还需要看具体的实现场景来做出选择。

下面详细介绍target结构，target结构由指令#pragma omp target和一个target计算区域组成，当程序执行到该语句，相关计算将由主机转移到指定设备中，等到设备端完成计算后，再将相应数据从设备端内存映射到主机内存中。

target结构的语法如下所示：

C/C++:
#pragma omp target [子句[[,]子句]...]
          结构体

fortran:
!omp target [clause[[,]clause]...]
          结构体
!omp end target

支持的子句:

C/C++和Fortran通用：

private(list), firstprivate(list)
in_reduction (reduction-identifier : list)
map ([[map-type-modifier[ ,] [map-type-modifier[,] ...]map-type: ] locator-list)
is_device_ptr (list)
defaultmap (implicit-behavior[:variable-category])
nowait
depend([depend-modifier, ] dependence-type : locator-list)
allocate ([allocator : ] list)
uses_allocators (allocator[ (allocator-traits-array)] [,allocator[ (allocator-traits-array) ] ...] )

C/C++:

if ([ target : ] scalar-expression)
device([device-modifier:] integer-expression)

Fortran:

if ([ target : ] scalar-logical-expression)
device ([device-modifier:] scalar-integer-expression)

map子句

默认情况下，在target结构中所涉及到的变量都会传入和传出，map子句显式的指定变量的映射形式，其参数有to、from和tofrom。to表示在target结构计算开始之前将相关参数从主机端内存映射到设备端内存，from表示在target结构计算完成之后将相关参数从设备端内存映射回主机端内存，tofrom则意味着相关参数在target结构开始计算前从主机端内存映射到设备端内存，当target结构完成后，数据再由设备端返回到主机端。

if子句

若使用了if子句，当if中的表达式为真时，target结构会在指定或默认设备上运行。当if表达式为假时，则target结构会在主机端执行。

private子句、firstprivate子句和defaultmap子句

这三个子句都将对在target 区域中出现的标量变量进行操作。private子句将一个或多个变量声明为线程的私有变量，每个线程都有它自己的变量私有副本，其他线程无法访问（类似于HIP编程模型中kernel里的局部变量，存在线程寄存器里）。需要注意的是，private子句变量是未初始化的变量，且在target结构结束之后，其值也不会传输到主机端。与此相对，firstprivate子句会将其标注的变量做一个初始化，然后传入设备端，且初始化的值是主机端该变量的值。注意，当没有任何子句时，默认在target结构中出现的标量变量都是firstprivate变量。

defaultmap子句的参数（implicit-behavior）有alloc, to, from, tofrom, firstprivate, none和 default，作用参数（variable-category）scalar, aggregate, pointer。

nowait子句

每个target结构指令都会带有一个隐式的barrier，使用nowait可以除掉barrier，使得程序在运行至并行循环的结束处不会进行线程同步，先执行完循环的线程会继续执行后续代码，从而加快运行速度。

depend子句

depend子句的依赖类型有in、out和inout

in:所产生的任务是同级任务的从属任务，它将引用至少一个出现前面的out 型或者inout 型同级任务变量列表中的变量。

out或inout :所产生的任务是同级任务的从属任务，它将引用至少一个出现在前面的in 型、out 型或者inout 型同级任务中的变量。

device子句

device子句用于指定具体哪一个DCU设备来进行计算，可选设备范围是[0:N-1]，其中N是设备个数。

此外，类似于与程序段target_test1.c中的#pragma omp parallel for指令的结构还有for、section、task等结构，下表列出了相关指令和描述。表中的一些指令可相互嵌套，如之前提到的#pragma omp parallel for等，具体请参考OpenMP相关资料。

指令	描述
parallel	放在一个代码段之前，表示这段代码将分配给多个线程进行并行执行
for	放在for 循环之前，将循环分配给多个线程并行执行，但必须保证每次循环之间没有相关性
sections	放在被并行执行的代码段之前
single	用在一段只被单个线程执行的代码段之前，表示后面的代码段将仅被一个线程执行
task	生成一个任务结构
teams	表不紧随其后的循环会分配给多个线程组

5.2 OpenMP数据传输

本小节讲一下OpenMP在异构计算中的几种数据传输方式，与hip编程模型类似，OpenMP编程模型在进行异构计算的过程中会频繁的进行主机端和设备端之间的数据拷贝。默认情况下，数据会自动的在主机端和设备端之间传输，但是某些数据可能不需要初始化或计算完成之后不需要将其传回主机，而这些非必要的数据传输就会影响程序速度。因此，除了一些默认的数据传输模式，OpenMP还提供了一些数据拷贝的操作接口方便更加精细化控制拷贝过程。

首先需要说明的是，在OpenMP异构计算过程中，如果没有任何显式说明的情况下，数组的默认映射类型是tofrom，指针的映射类型是具有tofrom 类型的零长度的指针，标量变量的默认类型为firstprivate ，这样在target语句结束时，标量变量的值不能回传给主机。

5.2.1 map子句

对于数组元素来讲，map子句是最直接的显式设置数组映射的方式，其语法格式如下：

map[映射类型：变量列表]

map子句的映射类型最主要的有to、 from和tofrom，变量列表可以是数组，也可以是结构元素。映射含义如图3-41所示，to表示在target等结构计算开始之前将变量从主机端内存映射到设备端内存，from表示在target等结构计算完成之后将相关参数从设备端内存映射回主机端内存，tofrom则意味着相关参数在target等结构开始计算前从主机端内存映射到设备端内存，当target结构完成后，数据再由设备端返回到主机端。

注：map子句适用范围有target结构，target data结构，target enter data结构，target exit data结构。

图3-43 map子句控制数据传输

5.2.2 结构化数据传输

使用OpenMP进行代码并行化的过程中，在设备端的数据传输和存储可分为结构化和非结构化两种，target data结构就是结构化的代表。在指令#pragma omp target data中通过map等子句所提到的变量，在整个target data区域内均持续存在，即一个target data 区域内的多个target结构共享target data 的设备数据环境。

我们重新编辑target_test1.c，加入数组q, 且数组q的值取决于数组p的最后计算结果，具体如下：

//target_test2.c
void vec_mult(int N)
{
  int i;
  float p[N], v1[N], v2[N];
  float q[N];
  init(v1, v2, N);

  #pragma omp target map(to:v1[0:N],v2[0:N]) map(from:p[0:N])
  {
    #pragma omp parallel for
    for (i=0; i<N; i++)
    {
      p[i] = v1[i] * v2[i];
    }
  }

  #pragma omp target map(to:p[0:N]) map(from:q[0:N])
  {
    #pragma omp parallel for
    for (i=0; i<N; i++)
    {
      q[i] =p[i] * i;
    }
  }
}

程序target_test2.c中，由于出现了两个target结构，因此数据的传入和传出会执行两次，会比较浪费时间，若使用target data结构则会省去一次数据的传输，具体如下：

void vec_mult(int N)
{
  int i;
  float p[N], v1[N], v2[N];
  float q[N];
  init(v1, v2,N);

  #pragma omp target data map(to:v1[0:N],v2[0:N]) map(from:p[0:N] , q[0:N])
  {
    #pragma omp target
    {
      #pragma omp parallel for
      for (i=0; i<N; i++)
      {
        p[i] = v1[i] * v2[i];
      }
    }

    #pragma omp target
    {
      #pragma omp parallel for
      for (i=0; i<N; i++)
      {
        q[i] =p[i] * i;
      }
    }
  }//end target map
}

程序中的两个target结构拥有各自的数据环境，同时他们都共享target data的存储空间，这样可以节省设备端和主机端的数据传输开销。需要注意的是，target data结构只是建立了一个设备数据环境，并不执行程序，且在设备端执行target区域的过程中，主机端无法访问设备端的数据。

target data结构的语法如下所示：

C/C++:
#pragma omp target data[子句[[,]子句]...]
          结构体

fortran:
!omp target data[clause[[,]clause]...]
          结构体
!omp end target data

支持的子句:

C/C++和Fortran通用：

map ([[map-type-modifier[ ,] [map-type-modifier[,] ...]map-type: ] locator-list)
use_device_ptr (ptr-list), use_device_addr (list)

C/C++:

if ([ target : ] scalar-expression)
device([device-modifier:] integer-expression)

Fortran:

if ([ target : ] scalar-logical-expression)
device ([device-modifier:] scalar-integer-expression)

5.2.3 非结构化数据传输

target enter data 和target exit data

target enter data和target exit data结构都是非结构化数据结构，它可以在主机代码内的任意位置生成和删除设备上的数据。

下面C++程序演示了在一个类中的构造函数和析构函数里，创建、删除设备端数据。构造函数使用target enter data创建一个数组变量v，并在map子句中使用alloc修饰符来避免将值复制到设备端。析构函数删除数据（target exit data），并在map子句中使用delete修饰符来避免将数据复制回主机端。需要说明的是独立的target enter data发生在创建主机向量之后，target exit data发生在删除主机端数据之前。

//target_ unstructured_data1.cpp
class Matrix
{

  Matrix(int n) {
    len = n;
    v = new double[len];
    #pragma omp target enter data map(alloc:v[0:len])
  }

  ~Matrix() {
    // NOTE: delete map type should be used, since the corresponding 
    // host data will cease to exist after the deconstructor is called.

    #pragma omp target exit data map(delete:v[0:len])
    delete[] v;
  }

  private:
  double* v;
  int len;
};

下面一段代码则演示了分配和释放Matrix结构的数据成员。init_target函数分配结构内存并使用target enter data指令将其映射到设备端。free_target函数删除设备端的映射参数并释放主机端相应内存。

//target_ unstructured_data2.c
typedef struct {
  double *A;
  int N;
} Matrix;

void init_target (Matrix *mat, int n)
{
  mat->A = (double *)malloc(n*sizeof(double));
  mat->N = n;
  #pragma omp target enter data map(alloc:mat->A[:n])
}

void free_target (Matrix *mat)
{
  #pragma omp target exit data map(delete:mat->A[:mat->N])
  mat->N = 0;
  free(mat->A);
  mat->A = NULL;
}

target enter data 和target exit data的语法如下：

C/C++:
#pragma omp target enter data [子句[[,]子句]...]
#pragma omp target exit data [子句[[,]子句]...]

fortran:
!omp target enter data [子句[[,]子句]...]
!omp target exit data [子句[[,]子句]...]

支持的子句:

C/C++和Fortran通用：

map ([[map-type-modifier[ ,] [map-type-modifier[,] ...]map-type: ] locator-list)
nowait
depend([depend-modifier, ] dependence-type : locator-list)

C/C++:

if ([ target : ] scalar-expression)
device([device-modifier:] integer-expression)

Fortran:

if ([ target : ] scalar-logical-expression)
device ([device-modifier:] scalar-integer-expression)

target update结构

target update指令是一个独立指令，它必须在target data结构内部。根据实际需要，显式的更新target data结构区域内的变量，主要的相关参数有to和from两个子句指定。

首先我们看来看一段代码，这段代码向我们展示了如何使用target update指令在target data结构中更新数据：

//target_ update.c
void vec_mult(int N)
{
  int i;
  float p[N], v1[N], v2[N];
  float out1,out2;

  init(v1, v2, N);
  #pragma omp target data map(to:v1,v2) map(from:p)
  {
    #pragma omp target 
      #pragma omp parallel for
      for (i=0; i<N; i++)
      {
        p[i] = v1[i];
        v2[i] = v2[i]+i;
      }

    init_again(v1,N);
    #pragma omp target update to(v1)        
    #pragma omp target update from(v2)

    #pragma omp target
      #pragma omp parallel for
      for (i=0; i<N; i++)
      {
        p[i] = v1[i] + p[i];
      }
  }
}

target_ update.c代码中在第一个target任务执行完之后，init_again函数将数组v1重新初始化，并将初始化后的值使用target update指令从主机端传输到设备端；同时，将在第一个target结构中计算后的数组变量v2使用target update指令从设备端更新到主机端。

target update指令用于在target data区域内显式的更新设备端或主机端的数据，更新的方向由关键字to和from来控制。更新的变量需要在主机端和设备端均有一个对应变量，to关键字表明将主机端变量传递给对应的设备端变量，同样的，from关键字表明将设备端变量更新到对应的主机端变量。

target update的语法如下：

C/C++:
#pragma omp target update [子句[[,]子句]...]

fortran:
!omp target update [子句[[,]子句]...]

支持的子句:

C/C++和Fortran通用：

map ([[map-type-modifier[ ,] [map-type-modifier[,] ...]map-type: ] locator-list)
nowait
depend([depend-modifier, ] dependence-type : locator-list)

C/C++:

if ([ target : ] scalar-expression)
device([device-modifier:] integer-expression)

Fortran:

if ([ target : ] scalar-logical-expression)
device ([device-modifier:] scalar-integer-expression)

motion-clause:

to ([mapper(mapper-identifier) : ] locator-list)
from ([mapper(mapper-identifier) : ] locator-list)

另外还有以下两点需要注意：

target update结构中必须存在一个to子句或from子句，且一个变量只能存在于这两个子句之一，相关变量不能同在两个子句中同时存在，子句to 或from 中的变量必须具有映射类型；
其次，对于出现在to或from子句中的变量，在主机端和设备端都必须存在一个对应的变量，否则将不进行数据更；

declare target结构

使用declare target指令声明，可以使target结构的函数和全局变量在默认设备上运行。以下代码展示了使用declare target指令声明设备端函数declare_test，声明后的declare_test函数可以直接在设备端调用。

//declare_target1.c
#pragma omp declare target
    extern void declare_test(int N);
#pragma omp end declare target
void fib_test(int n)
{
    #pragma omp target 
    {
        declare_test (n);
    }
}

下面程序演示了使用declare target指令将全局变量映射到设备端数据环境。

//declare_target2.c
#pragma omp declare target
float p[N], v1[N], v2[N];
#pragma omp end declare target

void vec_mult()
{
   int i;
   init(v1, v2, N);
   #pragma omp target update to(v1, v2)

   #pragma omp target
       #pragma omp parallel for
       for (i=0; i<N; i++)
           p[i] = v1[i] * v2[i];

   #pragma omp target update from(p)
   output(p, N);
}

程序中将全局变量v1,v2和p使用declare target指令映射到设备端，在函数vec_mult中可直接使用，并使用target update指令来更新数据。

declare target的语法如下：

C/C++:
#pragma omp declare target
      //结构体
#pragma omp end declare target
-or-
#pragma omp declare target(extended-list)
-or-
#pragma omp declare target[子句[[,]子句]...]

fortran:
!omp declare target (extended-list)
-or-
!omp declare target[子句[[,]子句]...]

支持的子句:

C/C++和Fortran通用：

to (extended-list), link (list)
device_type(host | nohost | any)

extended-list: A comma-separated list of named variables, procedure names, and named common blocks.

defaultmap子句

defaultmap子句可以显式的指定target结构中标量的数据传输，默认情况下，target结构中标量变量的类型firstprivate 变量，即标量会在传入设备端时被赋值，其值与该标量在主机端对应元素相等，但是等到计算结束后，这些在target结构区域内的标量值并不会被传递到主机端的同名变量。

defaultmap子句的语法如下：

defaultmap[显式操作：可变类别]

其中显式操作的关键字有to、from、tofrom、firstprivate、none、default、alloc、 present。可变类别为scalar、aggregate、pointer。此时，如果希望target结构中的标量在传入时被初始化，并在并行计算结束后将计算结果传回主机，我们就可以这样写：

#pragma omp target defaultmap(tofrom:scalar)
{
          //结构体
}

总结一下，map子句可以控制结构体内变量的映射方向，独立指令target update可以在target data区域内显式的指定某一变量更新方向，declare target结构声明的全局变量可以在默认的目标设备上运行，而defaultmap子句可以控制结构体内标量，指针的映射方向，下表对比了OpenMP和HIP在异构计算中的数据传输方式

OpenMP	HIP
map(to:dst) 或 target data to(dst)	hipError_t hipMalloc(void* dst, size_t size); hipError_t hipMemcpy (void dst, const void src, size_t sizeBytes, hipMemcpyHostToDevice); hipError_t hipFree(void dst);
map(from:dst) 或 target data from(dst)	hipError_t hipMalloc(void* dst, size_t size); hipError_t hipMemcpy (void dst, const void src, size_t sizeBytes, hipMemcpyDeviceToHost); hipError_t hipFree(void dst);
map(tofrom:dst) 或 target data tofrom(dst)	hipError_t hipMalloc(void* dst, size_t size); hipError_t hipMemcpy (void dst, const void src, size_t sizeBytes, hipMemcpyHostToDevice); hipError_t hipMemcpy (void dst, const void src, size_t sizeBytes, hipMemcpyDeviceToHost); hipError_t hipFree(void dst);
target update to(dst) 或 target enter data to(dst)	hipError_t hipMemcpy (void dst, const void src, size_t sizeBytes, hipMemcpyHostToDevice);
target update from(dst) 或 target exit data from(dst)	hipError_t hipMemcpy (void dst, const void src, size_t sizeBytes, hipMemcpyDeviceToHost);
#pragma omp declare target float dst[size]; #pragma omp end declare target	hipError_t hipMalloc(void* dst, size_t size); hipError_t hipMemcpy (void dst, const void *src, size_t sizeBytes, hipMemcpyHostToDevice);

5.3 执行配置

teams

teams指令可以创建一个由多个线程组构成的线程组群，线程组之间不可通信且无同步机制，每个线程组的主线程都会执行team区域，其他线程若无其他指令则处于空闲状态，teams指令类似于HIP中的资源调用的block设置。

teams指令必须包含在target结构内，在teams 区域内能以组合结构形式嵌套parallel区域。例如， distribute 、distribute simd 、distribute for 、distribute for simd 、parallel等。target_teams.c程序演示了target指令和teams指令的用法，程序创建了包含两个线程组的线程组群，且每个线程组群各完成一半的计算。

//target_teams.c
float teams_test(float B[], float C[], int N)
{
   float sum0 = 0.0;
   float sum1 = 0.0;
   #pragma omp target map(to: B[:N], C[:N]) map(tofrom: sum0, sum1) 
   #pragma omp teams num_teams(2)
   {
      int i;
      if (omp_get_num_teams() != 2)
         abort();
      if (omp_get_team_num() == 0)
      {
     #pragma omp parallel for reduction(+:sum0)
     for (i=0; i<N/2; i++)
        sum0 += B[i] * C[i];
      }
      else if (omp_get_team_num() == 1)
      {
     #pragma omp parallel for reduction(+:sum1)
     for (i=N/2; i<N; i++)
        sum1 += B[i] * C[i];
      }
   }
   return sum0 + sum1;
}

程序中omp_get_num_teams()返回线程组群定义的线程组数，omp_get_team_num()返回当前线程组号。每个线程组的主线程运行各自#pragma omp parallel for指令区域内的计算主体。

teams指令的语法如下：

C/C++
#pragma omp teams [子句[[,]子句]...]
      //结构体

fortran:
!omp teams [子句[[,]子句]...]
        !结构体
!omp end teams

支持的子句:

C/C++和Fortran通用：

private (list), firstprivate (list), shared (list)
reduction ([default ,] reduction-identifier : list)
allocate ([allocator : ]list)

C/C++:

num_teams (integer-expression)
thread_limit (integer-expression)
default (shared | none)

Fortran:

num_teams (scalar-integer-expression)
thread_limit (scalar-integer-expression)
default (shared | firstprivate | private | none)

在程序target_teams.c中我们不难发现，使用teams指令时，只有主线程在工作，其余线程处于空闲状态，这样浪费了不少资源，此时我们引入distribute指令来优化代码。

distribute

distribute指令的功能实工作共享，需要注意的是此指令必须包含在一个teams区域内。distribute指令将其紧随的for循环计算以隐式任务形式静态的分配给teams指令指定的个线程组的主线程执行。下面我们将以一段代码来讲解一下distribute指令的用法和运行方式。

//target_teams_distribute.c
void distribute_test()
{
    #pragma omp target 
        #pragma omp teams num_teams(2) num_threads(2)
        #pragma omp distribute
        {
            for(int i=0;i<6;i++)
            {
                #pragma omp parallel for
                for(int j=0;i<4;i++)
                    p[i,j] = v1[i,j];
            }
        }
}

可以看出，teams指令创立了一个含有两个线程组的线程组群，每个线程组含有两个线程。distribute指令将循环迭代的计算任务分配给teams区域内全部线程组的主线程，当遇到内循环parallel for指令时，线程组内的其余线程被激活如图3-42。

图 3-44 teams、distribute指令执行过程

5.4 常用函数

在使用OpenMP进行异构编程过程中，我们经常需要对运行环境进行查询和设定，下表给出了常用的一些内置函数:

表 3-x 内置函数

函数名	描述
omp_set_num_threads (int num_threads)	设置后续并行区域中并行执行的线程数
omp_get_num_threads(void)	确定当前并行区域内活动线程数量
omp_get_max_threads(void)	返回当前的并行区域内可用的最大线程数
omp_get_thread_num(void)	返回当前的线程号。线程号的值在0 (主线程)到线程总数减1 之间
omp_in_parallel(void)	确定正在运行的代码区域处于并行状态还是串行状态。如果处于并行状态，则返回值必真;否则为假
omp_get_thread_limit(void)	返回可用线程的最大数量
omp_get_team_size(int level)	对于指定的当前嵌套层，返回线程组的线程数量
omp_get_place_num_procs(void)	返回指定位置的执行环境的可用处理器的数
omp_set_default_device(int device_num)	通过默认目标设备号来控制默认的目标设备
omp_get_default_device(void)	返回默认目标设备的设备号
omp_get_num_teams(void)	返回当前区域的线程组的数量，如果在线程组区域之外调用则返问值为1
omp_get_team_num(void)	返回当前调用线程的线程组号
omp_is_initial_device(void)	如果当前任务在主机设备上执行，则返回true;否则false
omp_get_initial_device	返回代表主机设备的设备编号

5.5 总结

图 3-45 OpenMP异构计算的组成

注：目前OpenMP编程模型还不支持多流。

技术瘾君子1573

关注

27
点赞
踩
19

收藏

觉得还不错? 一键收藏
打赏
0
评论
国产加速器海光DCU&GPGPU深算处理器异构编程实战（下）

Fortran语言自20世纪50年代诞生至今，一直被广泛应用于计算物理学、计算化学、流体力学、气候预测、有限元分析等计算密集的学科领域。同C/C++语言一样，Fortran语言也是科学/工程计算领域的主流编程语言，使用Fortran语言编写的应用程序运行在世界各地大大小小的计算集群上，为人类的各种研究与创新提供数据支撑。
复制链接

扫一扫