OpenACC中!$acc update self位置影响数组更新结果的原因咨询

我是OpenACC新手，在用它加速粒子代码时发现，主机端更新数组eta时，!$acc update self指令的位置会导致不同结果，以下是复现代码及模块：

主程序代码 #

program approximateFun
use funs
use paras_mod

    integer :: Nx
    real(dp) :: dx
    real(dp), dimension(:), allocatable :: x, eta
    real(dp), dimension(5) :: xp, fAtxp

    !$acc declare create(Nx)
    !$acc declare create(x)
    !$acc declare create(eta)
    !$acc declare create(dx)
    !$acc declare create(fAtxp)
    !$acc declare create(xp)

    Nx = 16
    !$acc update device(Nx)

    xp = (/3.9, 4.1, 4.5, 5.0, 5.6/)
    !$acc update device(xp)

    allocate(x(1 : Nx))
    allocate(eta(1 : Nx))
    eta = 0.0d0

    dx = 2 * pi / (Nx - 1)
    !$acc update device(dx)

    do i = 1, Nx
        x(i) = (i - 1.0d0) * dx
    end do
    !$acc update device(x)


    call calc_etaVec(x, Nx, eta)
    !$acc update self(eta)        ! gives the correct results

    !$acc parallel loop present(dx, xp, eta, fAtxp)
    do i = 1, 5
        call calcFunAtx(xp(i), dx, eta, fAtxp(i))
    end do

    !$acc update self (fAtxp)
    !!$acc update self(eta) !---> gives wrong result

    write(6, *) 'eta', eta
    do i = 1, 5
        write(6, *) 'xp, fAtxp', xp(i), fAtxp(i)
    end do

    deallocate(x)
    deallocate(eta)

end program approximateFun

依赖模块代码 #

funs模块 #

MODULE funs
    use paras_mod
    implicit none
    CONTAINS

    subroutine calc_etaVec(x, nx, eta)
        integer, intent(in) :: nx
        real(dp), dimension(:), intent(in) :: x
        real(dp), dimension(:), intent(out) :: eta

        integer :: i
        !$acc parallel loop present(x, eta)
        do i = 1, nx
            eta(i) = sin(x(i))
        end do
    end subroutine

    subroutine calcFunAtx(xp, dx, eta, fAtx)
        real(dp), intent(in) :: xp, dx
        real(dp), dimension(:), intent(in) :: eta
        real(dp), intent(out) :: fAtx
        integer :: idx

        !$acc routine seq

        idx = 1 + floor(xp / dx)
        fAtx = eta(idx)
    end subroutine calcFunAtx

END MODULE

paras_mod模块 #

module paras_mod
    implicit none
    save
    INTEGER, PARAMETER :: dp = selected_real_kind(14,300)
    REAL(dp), PARAMETER :: PI=4.0d0*atan(1.0d0)
end module paras_mod

两种情况的输出结果 #

情况1：!$acc update self(eta)紧跟calc_etaVec后（结果正确） #

0.000000000000000        0.4067366430758002
   0.7431448254773941        0.9510565162951535        0.9945218953682734 
   0.8660254037844388        0.5877852522924732        0.2079116908177597 
  -0.2079116908177591       -0.5877852522924730       -0.8660254037844384 
  -0.9945218953682732       -0.9510565162951536       -0.7431448254773946 
  -0.4067366430758009       -2.4492935982947064E-016

情况2：!$acc update self(eta)放在后续循环之后（结果错误） #

0.000000000000000        0.4067366430758002
    0.7431448254773941        0.9510565162951535        0.9945218953682734
    0.000000000000000         0.000000000000000         0.000000000000000
    0.000000000000000         0.000000000000000         0.000000000000000
    0.000000000000000         0.000000000000000         0.000000000000000
    0.000000000000000         0.000000000000000

原因分析 #

核心问题出在OpenACC的隐式拷贝规则和present子句的行为：

• 你用!$acc declare create(eta)仅在设备端分配了eta的存储空间，并未将主机端初始的全0值拷贝到设备（create只分配空间，不做数据初始化拷贝）。
• 调用calc_etaVec后，设备端的eta被正确计算为sin(x(i))的完整结果，但主机端的eta仍然是初始的全0值，此时主机和设备端的eta数据不一致。
• 当后续执行!$acc parallel loop present(dx, xp, eta, fAtxp)时，present(eta)会触发OpenACC runtime检查数据一致性：由于主机端eta从未同步过设备端的修改，runtime会认为主机端数据是“最新”的，自动将主机端的全0eta拷贝回设备端，覆盖掉之前计算的正确值。
• 而calcFunAtx只访问了eta的前5个元素（xp对应的idx都在前5位），这部分数据在被覆盖前已经完成读取计算，所以fAtxp结果正确，但设备端剩下的eta元素已被0覆盖，最后同步回主机时就只有前5个元素正确，其余为0。
• 反之，如果在calc_etaVec后立刻执行!$acc update self(eta)，会先把设备端完整的eta同步回主机，此时主机和设备端数据一致，后续present(eta)不会触发隐式拷贝，设备端eta保持正确，最终同步回主机的就是完整的正确结果。

优化建议 #

• 用!$acc enter data copyin(eta)替代!$acc declare create(eta)，初始化时就将主机端的eta拷贝到设备，避免后续隐式拷贝触发；
• 保持calc_etaVec后立刻同步eta到主机的操作，确保主机设备数据一致性；
• 在后续并行循环中使用deviceptr(eta)替代present(eta)，明确指定使用设备端数据，规避隐式拷贝逻辑。

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