Fortran嵌套例程向OpenACC移植咨询：GPU并行化方案抉择

背景概述 #

我刚接触OpenACC，计划将一套采用MPI+OpenMP混合模式的大型Fortran代码移植到GPU，编译器使用nvfortran。代码核心是3D数组（索引(i1,i2,i3)）的时间演化，所有操作仅沿单一坐标进行，不会混合坐标，例如：

A(i1,i2,i3) = B(i1+1,i2,i3) + C(i1-1,i2,i3)
! OR
A(i1,i2,i3) = B(i1,i2+1,i3) + C(i1,i2-1,i3)

不会出现A(i1,i2,i3) = B(i1+1,i2,i3) + C(i1,i2-1,i3)这类跨坐标操作。现有优化逻辑为：先沿i1执行所有操作，转置数组为(i2,i1,i3)后沿i2操作，再转置为(i3,i1,i2)沿i3操作，最终转置回(i1,i2,i3)。

现有代码框架（简化版）：

integer :: m1,n1,m2,n2,m3,n3
integer :: m1c,n1c,m2c,n2c,m3c,n3c
integer :: Nmax, nchunks, ichunk
double precision, allocatable(:,:,:) :: A,B,C
double precision, allocatable(:,:,:) :: sliceA, sliceB, sliceC
!
allocate(A(m1:n1,m2:n2,m3:n3),B(m1:n1,m2:n2,m3:n3),C(m1:n1,m2:n2,m3:n3))
... ! A和B初始化
! 获取分块数量
call getChunkSize(Nmax,m1,n1,m2,n2,m3,n3,1,m2c,n2c,m3c,n3c,nchunks)
! OpenMP并行区域
!$OMP PARALLEL DEFAULT(SHARED), &
!$OMP          PRIVATE(sliceA,sliceB,sliceC, &
!$OMP                  m2c, n2c, m3c, n3c, &
!$OMP                  ichunk)
... ! 其他操作
!$OMP DO
do ichunk=1,nchunks
  ! 获取当前块的索引范围
  call getChunkSize(Nmax,m1,n1,m2,n2,m3,n3,ichunk,m2c,n2c,m3c,n3c)
  !
  allocate(sliceA(m1:n1,m2c:n2c,m3c:n3c),sliceB(m1:n1,m2c:n2c,m3c:n3c),sliceC(m1:n1,m2c:n2c,m3c:n3c))
  sliceA(:,:,:) = A(m1:n1,m2c:n2c,m3c:n3c)
  sliceB(:,:,:) = B(m1:n1,m2c:n2c,m3c:n3c)
  !
  call operations(sliceA,sliceB,sliceC)
  !
  C(m1:n1,m2c:n2c,m3c:n3c) = sliceC(:,:,:)
  deallocate(sliceA,sliceB,sliceC)
end do
!$OMP END DO
... ! 其他操作
!$OMP END PARALLEL
!
! 重复类似逻辑，分别处理(i2,i1,i3)和(i3,i1,i2)维度的切片

其中Nmax控制切片最大元素数以优化缓存，getChunkSize将大数组拆分为nchunks个块，每个块元素数≤Nmax。外层用OpenMP并行处理每个切片，调用operations计算后写回结果。

operations例程包含三重循环，仅i1索引有邻域访问（i1±1、i1±2），i2、i3索引无跨值操作：

do i3=m3c,n3c
  do i2=m2c,n2c
    do i1=m1,n1
      sliceC(i1,i2,i3) = ...
    end do
  end do
end do

两种并行化方案的优缺点分析 #

方案1：保留现有分块逻辑，映射到GPU gangs/workers/vectors #

实现要点

1. 外层块循环用!$ACC PARALLEL LOOP GANG替代OpenMP的!$OMP DO，可通过预处理指令区分OpenMP/OpenACC模式。
2. operations例程需声明为OpenACC设备例程：
   • 若例程内手动控制并行，用!$ACC ROUTINE PARALLEL；
   • 若让编译器自动并行内部循环，用!$ACC ROUTINE SEQUENTIAL并在循环上添加!$ACC LOOP WORKER/VECTOR。

Nmax取值建议

• 每个切片元素数需足够大，能填满GPU的流多处理器（SM），建议单切片包含至少几千到几万个元素；
• nchunks最好是GPU可用gang数的整数倍，避免负载不均衡；若无法满足，需确保剩余gang的工作量足够。

优缺点

• 优点：复用现有分块逻辑，代码改动量小；分块能控制GPU内存占用，避免大数组导致显存不足。
• 缺点：分块带来的主机-设备数据拷贝开销会随分块数量增加而上升；需手动协调gang/workers/vectors的映射，调试成本较高。

方案2：取消分块，仅并行化内部三重循环 #

实现要点

1. 设置Nmax=(n1-m1+1)*(n2-m2+1)*(n3-m3+1)，让nchunks=1，直接对整个数组调用operations。
2. 在operations的三重循环上添加!$ACC PARALLEL LOOP，或用!$ACC ROUTINE PARALLEL声明例程。
3. 用预处理指令隔离OpenMP/OpenACC代码：

   #ifdef _OPENACC
     ! OpenACC逻辑：直接处理全数组
     call operations(A,B,C)
   #else
     ! 原OpenMP分块逻辑
     !$OMP PARALLEL ...
     ! ...
     !$OMP END PARALLEL
   #endif
   

优缺点

• 优点：减少数据拷贝次数，仅需一次全数组的主机-设备传输；编译器可更高效地优化循环并行化，避免手动分块的额外开销。
• 缺点：若数组过大，会占用过多显存，可能导致显存溢出或分页（大幅降低性能）；完全抛弃现有分块逻辑，代码改动量较大。

更优方案建议 #

优先尝试方案2的自适应改进版：

1. 如果GPU显存能容纳整个数组，直接采用方案2，最大化并行效率；
2. 如果显存不足，根据GPU显存大小动态计算Nmax，让每个切片刚好能高效利用SM且不超出显存限制，实现自适应分块。
3. 结合MPI+OpenACC混合模式：每个MPI进程负责一部分数组，将该部分数据传到GPU处理，避免跨进程数据传输，进一步提升性能。

operations例程的OpenACC语法示例 #

以下示例包含临时数组的设备端分配/释放，以及循环并行化：

!$ACC ROUTINE PARALLEL
subroutine operations(sliceA, sliceB, sliceC)
  implicit none
  double precision, intent(in) :: sliceA(:,:,:), sliceB(:,:,:)
  double precision, intent(out) :: sliceC(:,:,:)
  double precision, allocatable :: temp(:,:,:)
  integer :: i1,i2,i3, m1,n1,m2,n2,m3,n3

  m1 = lbound(sliceA,1); n1 = ubound(sliceA,1)
  m2 = lbound(sliceA,2); n2 = ubound(sliceA,2)
  m3 = lbound(sliceA,3); n3 = ubound(sliceA,3)

  ! 在GPU设备上分配临时数组
  !$ACC ENTER DATA CREATE(temp(m1:n1,m2:n2,m3:n3))

  ! 合并i2、i3循环映射到gangs/workers，i1循环映射到vectors
  !$ACC PARALLEL LOOP COLLAPSE(2) DEFAULT(PRESENT)
  do i3=m3,n3
    do i2=m2,n2
      !$ACC LOOP VECTOR
      do i1=m1,n1
        ! 示例计算：仅i1邻域操作
        temp(i1,i2,i3) = sliceA(i1,i2,i3) + sliceB(i1+1,i2,i3) + sliceB(i1-1,i2,i3)
        sliceC(i1,i2,i3) = temp(i1,i2,i3) * 0.5d0
      end do
    end do
  end do

  ! 释放GPU设备上的临时数组
  !$ACC EXIT DATA DELETE(temp)
end subroutine operations

• !$ACC ROUTINE PARALLEL：声明该例程可在GPU上并行执行；
• !$ACC ENTER DATA CREATE(temp)：直接在设备上分配临时数组，避免主机-设备拷贝；
• COLLAPSE(2)：合并i2、i3循环，提升并行粒度；
• DEFAULT(PRESENT)：告知编译器所有数组已在设备上存在，无需额外拷贝。

内容的提问来源于stack exchange，提问作者Fabio Riva